对偶单纯型法解线性规划问题

本程序实现的是单纯型法解线性规划问题，在检验向量e<=0的情况下，迭代得到最优解。

首先，要求你在data.txt文件中输入数据，输入方式，例：

 

 在data.txt中输入：

3
1 1 1
2 3
 3 1 1 >=  1
-1 4 1 >=  2
 

便可出现执行结果!

实现程序如下:

1. #include<iostream>  
2. #include<fstream>  
3. #include<string>  
4. #include<cstdlib>  
5. #define INF -100000  
6. #define MAX 100000  
7. using namespace std;  
8. int N1, N2;  
9. //约束条件系数   
10. double **a;  
11. double *b; // 向量b  
12. double *e;//检验数向量  
13. double RHS = 0;//目标函数值    
14. int *x; // 存放基本可行解   
15. int m;  // 添加的辅助变量的个数   
16. int JudgeSym(string symbol)  
17. {  
18.     if(symbol==">=")  
19.       return -1;  
20.     if(symbol=="<=")  
21.       return 1;  
22.     if(symbol=="=")  
23.       return 0;  
24. }  
25. // 找出初值解    
26. void Begin()  
27. {  
28.      //进行列搜索， 得到初始解   
29.     int index = 0;  
30.     for(int j=1; j<N2+m; ++j){  
31.       int cnt = 0;  
32.      for(int i=1; i<N1; ++i)  
33.       if(a[i][j]){  
34.          ++cnt;  
35.          index = i;  
36.        }  
37.       if(cnt==1&&x[index]==0){  
38.          if(a[index][j] != 1){  
39.              for(int k=1; k<N2; ++k)  
40.                if(k!=j)  
41.                  a[index][k] /= a[index][j];  
42.              b[index] /= a[index][j];  
43.              a[index][j] = 1;             
44.          }  
45.         x[index] = j;   
46.       }  
47.     }  
48.     cout<<"得到对偶问题的初始基本可行解为:"<<endl;  
49.     for(int i=1; i<N1; ++i)  
50.         cout<<"x("<<x[i]<<")"<<" = "<<b[i]/a[i][x[i]]<<endl;  
51. }  
52. void GetFirstForm()  
53. {  
54.      double tmp = 0;  
55.      for(int j=1; j<N1; ++j)  
56.          if(e[x[j]]!=0){  
57.             tmp = e[x[j]];  
58.          for(int i=1; i<N2; ++i){  
59.              /*May Have  Error */ 
60.              e[i] -=   tmp*a[j][i];  
61.             }  
62.          RHS -= tmp*b[j];          
63.      }  
64. }  
65. int CheckRvec()  
66. {  
67.     double min = MAX;  
68.     int index = 0;  
69.     for(int i=1; i<N1; ++i)  
70.       if(b[i]<0&&b[i]<min){  
71.        min =b[i];  
72.        index = i;  
73.       }  
74.     if(index!=0)  
75.       return index;  
76.     cout<<"迭代得到一个最优解:RHS = "<<RHS<<endl;  
77.     for(int i=1; i<N1; ++i)  
78.       cout<<"x("<<x[i]<<")"<<" = "<<b[i]/a[i][x[i]]<<endl;  
79.     return -1;  
80. }  
81. int MinEj(int Xiout)  
82. {  
83.     double min = MAX;  
84.     int index = -1;  
85.     for(int j=1; j<N2+m; ++j)  
86.       if(a[Xiout][j]<0&&e[j]/a[Xiout][j]<min){  
87.         index = j;  
88.         min = e[j]/a[Xiout][j];  
89.       }  
90.     if(index == -1)  
91.      cout<<"a(r,j)>=0;无最优解!"<<endl;  
92.     return index;  
93. }  
94. void InputandInit()  
95. {  
96.     void Standardlize(int,int);  
97.     ifstream infile("data.txt");   
98.     infile>>N2;   
99.     N2 += 1;  
100.     //存储没有添加辅助变量的   
101.     double *tpe = new double[N2];  
102.     for(int i=1; i<N2; ++i){  
103.      infile>>tpe[i];  
104.      tpe[i] = -tpe[i];  
105.      }/*修改，根据BeginOne()的返回值，判断是否进行e的再次录入*/ 
106.     infile>>N1>>N2;  
107.     m = N1; /*最多添加的辅助变量的个数*/ 
108.     N1 += 1;  
109.     N2 += 1;  
110.     //正式进行目标函数的录入   
111.     e = new double[N2+m];  
112.     memset(e,0,sizeof(double)*(N2+m));  
113.     memcpy(e,tpe,sizeof(double)*N2);  
114.     delete tpe;/**/ 
115.     //为 系数矩阵申请内存   
116.     a = new double*[N1];  
117.     for(int i=1; i<N1; ++i){  
118.       a[i] = new double[N2+m]; /*a : N1/N2+N1 */ 
119.       memset(a[i],0,sizeof(double)*(N2+m));  
120.     }  
121.     //为 右端向量申请内存   
122.     b = new double[N1];  
123.     memset(b,0,sizeof(double)*N1);  
124.       
125.     m = 0; /*申请内存结束后,重新置m为零,使m记录实际添加的人工变量的个数*/ 
126.        
127.     //输入a系数矩阵   
128.     for(int i=1; i<N1; ++i){  
129.      for(int j=1; j<N2; ++j){  
130.        infile>>a[i][j];  
131.      }  
132.      string Symbol;  
133.      infile>>Symbol;  
134.      infile>>b[i];  
135.      int numSym = JudgeSym(Symbol);  
136.      Standardlize(numSym,i);  
137.      }  
138. }  
139. void Standardlize(int numSym, int line)  
140. {  
141.      int row = N1 -1;  
142.      int col = N2 -1;  
143.      static int cnt = 1;  
144.      switch(numSym){  
145.       case -1:  
146.         for(int j=1; j<N2; ++j)  
147.           a[line][j] = -a[line][j];  
148.         b[line] = -b[line];  
149.         a[line][col+cnt] = 1;/*添加辅助变量*/ 
150.         ++cnt;  
151.         ++m;  /*辅助变量加1*/ 
152.         break;  
153.       case 1:  
154.         a[line][col+cnt] = 1;  
155.          ++cnt;  
156.          ++m; /*辅助变量加1*/ 
157.         break;  
158.       default: return;  
159.      }  
160.      return;  
161. }  
162. void print()  
163. {  
164.      static int cnt = 0;  
165.      cout<<left<<showpos;  
166.      cout.width(5);   
167.      cout.precision(3);  
168.      cout<<cnt<<"-Z  ";  
169.     for(int i=1; i<N2+m; ++i)  
170.       cout<<e[i]<<" ";  
171.       cout<<RHS<<endl;  
172.     for(int i=1; i<N1; ++i){  
173.       cout<<"x("<<x[i]<<") ";  
174.       for(int j=1; j<N2+m; ++j)  
175.         cout<<a[i][j]<<" ";  
176.       cout<<b[i]<<endl;  
177.     }  
178.     cout<<"***********************"<<endl;  
179.     ++cnt;  
180. }   
181. int main()  
182. {  
183.     InputandInit();  
184.     /* Until then; data input finished */ 
185.     x = new int[N1];  
186.     memset(x,0,sizeof(int)*N1);  
187.     Begin();  
188.     GetFirstForm();  
189.     //打印第一张单纯行表  
190.     cout<<"The First Form Of Simplex Method:"<<endl;  
191.     print();  
192.     /*循环得出最优解*/ 
193.     int Xiout = 0, ej = 0;           /*找到最小的e(j)/a(Xiout)(j)*/ 
194.     while((Xiout=CheckRvec())!=-1 &&(ej = MinEj(Xiout))!=-1){  
195.         x[Xiout] = ej;  
196.         for(int i=1; i<N2+m; ++i)  
197.           if(ej!=i){  
198.            a[Xiout][i] = a[Xiout][i]/a[Xiout][ej];  
199.            }  
200.         /*Error will be found*/ 
201.         b[Xiout]/=a[Xiout][ej];  
202.         a[Xiout][ej] = 1;  
203.         //化检验数向量  
204.         double tmp = e[ej];  
205.         for(int i=1; i<N2+m; ++i)  
206.            e[i] -= tmp*a[Xiout][i];  
207.         RHS -= tmp*b[Xiout];  
208.         /*化右端向量*/ 
209.         for(int i=1; i<N1; ++i){  
210.           if(i!=Xiout){  
211.            tmp = a[i][ej];  
212.           for(int j=1; j<N2+m; ++j)  
213.             a[i][j] -= tmp*a[Xiout][j];  
214.           b[i] -= tmp*b[Xiout];  
215.           }  
216.         }  
217.         print();  
218.     }  
219.     system("pause");  
220.     return 0;  
221. }  

 

转载于:https://blog.51cto.com/hustluy/412999