线性规划 单纯形算法
根据单联通区域求极值的想法得来的算法思想,故而叫做simplex algorithm。
暂且将写下的程序贴出,做个备份。
具体算法分析待看完证明后贴出。
此算法实现完全是对《算法导论》上算法框架的实现。
注:要使用g++编译器,vc要将变量定义做一些更改。
根据单联通区域求极值的想法得来的算法思想,故而叫做simplex algorithm。
暂且将写下的程序贴出,做个备份。
具体算法分析待看完证明后贴出。
此算法实现完全是对《算法导论》上算法框架的实现。
注:要使用g++编译器,vc要将变量定义做一些更改。
- #include <iostream>
- #include <set>
- #include <vector>
- #include <ctime>
- #include <climits>
- #include <cmath>
- using namespace std;
- int const BUF = 100;
- set<int> N,B;
- set<int>::iterator p,q;
- double bi[BUF];
- double ci[BUF];
- double v;
- double A[BUF][BUF];
- double res[BUF];
- int Num;
- //question:
- // request: max v+ci*x (ci=0,if and only if x[i] not in N)
- // constraint: y = bi+A*x (x[i] is in B and y[i] is in B)
- // solution: res and v
- // output the form of the question's equalities.
- // you can invoke this function after "pivot" function to see the changes.
- void output()
- {
- set<int>::iterator pt;
- printf("z = %.2lf ",v);
- p = N.begin();
- while(p!=N.end()){
- pt = p;
- pt++;
- if(pt!=N.end()){
- if(ci[*p]>0)
- printf(" + %.2lf*x(%d) ",ci[*p],*p);
- else
- printf(" - %.2lf*x(%d) ",fabs(ci[*p]),*p);
- }else {
- if(ci[*p]>0)
- printf(" + %.2lf*x(%d)/n",ci[*p],*p);
- else
- printf(" - %.2lf*x(%d)/n",fabs(ci[*p]),*p);
- }
- p++;
- }
- for(p=B.begin();p!=B.end();p++){
- printf("x(%d) = %.2lf ",*p,bi[*p]);
- q = N.begin();
- while(q!=N.end()){
- pt = q;
- pt++;
- if(pt!=N.end()){
- if(A[*p][*q]>0)
- printf(" - %.2lf*x(%d)",A[*p][*q],*q);
- else
- printf(" + %.2lf*x(%d)",fabs(A[*p][*q]),*q);
- }else {
- if(A[*p][*q]>0)
- printf(" - %.2lf*x(%d)/n",A[*p][*q],*q);
- else
- printf(" + %.2lf*x(%d)/n",fabs(A[*p][*q]),*q);
- }
- q++;
- }
- }
- printf("/n");
- }
- // input format:
- // first line one number Num representing
- // the num of all variables including the elements of N and B
- // next Num lines representing A[i][j] which is the coefficient of x[i]
- // next two lines each includes Num numbers --- the ci and bi respectively
- // for example
- /*
- input:@file = "t.in"
- 6
- 0 0 0 0 0 0
- 0 0 0 0 0 0
- 0 0 0 0 0 0
- 1 1 3 0 0 0
- 2 2 5 0 0 0
- 4 1 2 0 0 0
- 3 1 2 0 0 0
- 0 0 0 30 24 36
- 0
- representing the standard format:
- z = 3*x(0) + x(1) + 2*x(2)
- x(3) = 30 - 1*x(0) - x(1) - 3*x(2)
- x(4) = 24 - 2*x(0) - 2*x(1) - 5*x(2)
- x(5) = 36 - 4*x(0) - *x(1) - 2*x(2)
- */
- void init()
- {
- freopen("t.in","r",stdin);
- freopen("t.out","w",stdout);
- scanf("%d",&Num);
- for(int i=0;i<Num;i++)
- for(int j=0;j<Num;j++)
- scanf("%lf",&A[i][j]);
- for(int i=0;i<Num;i++)
- scanf("%lf",&ci[i]);
- for(int i=0;i<Num;i++)
- scanf("%lf",&bi[i]);
- scanf("%lf",&v);
- for(int i=0;i<Num;i++)
- ci[i]!=0 ? N.insert(i):B.insert(i);
- output();
- srand((unsigned)clock());
- }
- // relex each x[i]'s coefficient the main part of this algorithm.
- //x[e] is the new basic variable. x[l] is changed into nonbasic variable.
- void pivot(int l,int e)
- {
- bi[e] = bi[l]/A[l][e];
- for(p=N.begin();p!=N.end();p++)
- A[e][*p] = A[l][*p]/A[l][e];
- A[e][l] = 1/A[l][e];
- for(p=B.begin();p!=B.end();p++){
- if(*p==l)continue;
- bi[*p] -= A[*p][e]*bi[e];
- for(q=N.begin();q!=N.end();q++)
- if(*q!=e)
- A[*p][*q] -= A[*p][e]*A[e][*q];
- A[*p][l] = -1*A[*p][e]*A[e][l];
- }
- v += ci[e]*bi[e];
- for(p=N.begin();p!=N.end();p++)
- if(*p!=e)
- ci[*p] -= ci[e]*A[e][*p];
- ci[l] = -1*ci[e]*A[e][l];
- N.insert(l),N.erase(e);
- B.insert(e),B.erase(l);
- }
- // basing a little random selection on selecting
- // which x[i] will be changed into basic variables
- void simplex()
- {
- double temp = INT_MAX;
- vector<int> v;
- while(1){
- temp = INT_MAX;
- v.clear();
- for(p=N.begin();p!=N.end();p++)
- if(ci[*p]>0)
- v.push_back(*p);
- if(v.size()==0){
- temp = 0;
- break;
- }
- int t = rand()%v.size();
- int tn;
- for(q=B.begin();q!=B.end();q++)
- if(A[*q][t]>0){
- double tk = bi[*q]/A[*q][t];
- if(temp>tk)temp=tk,tn=*q;
- }
- printf("l = %d e = %d/n",tn,t);
- if(temp==INT_MAX)
- break;
- else
- pivot(tn,t);
- output();
- }
- if(temp==INT_MAX){
- printf("no answer/n");
- //which means unbounded,i.e. the question has no answer
- exit(0);
- }
- }
- void result()
- {
- printf("/n");
- for(p=B.begin();p!=B.end();p++)
- res[*p] = bi[*p];
- printf("ans seq: ");
- for(int i=0;i<Num;i++)
- printf("%.2lf ",res[i]);
- printf("/nmax: %.2lf/n",v);
- }
- int main()
- {
- init();
- simplex();
- result();
- return 0;
- }