能用的单纯形法python代码

 1 import numpy as np
 2 
 3 def pivot():
 4     l = list(d[0][:-2])
 5     jnum = l.index(max(l)) #转入编号
 6     m = []
 7     for i in range(bn):
 8         if d[i][jnum] == 0:
 9             m.append(0.)
10         else:
11             m.append(d[i][-1]/d[i][jnum])
12     inum = m.index(min([x for x in m[1:] if x!=0]))  #转出下标
13     s[inum-1] = jnum
14     r = d[inum][jnum]
15     d[inum] /= r
16     for i in [x for x in range(bn) if x !=inum]:
17         r = d[i][jnum]
18         d[i] -= r * d[inum]
19         
20 def solve():
21     flag = True
22     while flag:
23         if max(list(d[0][:-1])) <= 0: #直至所有系数小于等于0
24             flag = False
25         else:
26             pivot()
27             
28 def printSol():
29     for i in range(cn - 1):
30         if i in s:
31             print("x"+str(i)+"=%.2f" % d[s.index(i)+1][-1])
32         else:
33             print("x"+str(i)+"=0.00")
34     print("objective is %.2f"%(-d[0][-1]))
35 
36 d = np.loadtxt("data.txt", dtype=np.float)
37 (bn,cn) = d.shape
38 s = list(range(cn-bn,cn-1)) #基变量列表
39 solve()
40 printSol()
41 
42 # txt文件里面的格式
43 
44 1 14 6 0 0 0 0 0
45 
46 1 1 1 1 0 0 0 4
47 
48 1 0 0 0 1 0 0 2
49 
50 0 0 1 0 0 1 0 3
51 
52 0 3 1 0 0 0 1 6
53 
54 # 求解模型
55 max z = x0+14*x1+6*x2
56 
57 s.t.
58 
59 x0 + x1 + x2 <= 4
60 
61 x0 <= 2
62 
63 x2 <= 3
64 
65 3*x1 + x2 <= 6

  1 Max z =  18*x1 + 12.5*x2
  2 
  3 s.t. x1 + x2 +x3 = 20
  4 
  5 x1 + x4 = 12
  6 
  7 x2 + x5 = 16
  8 
  9 x1,x2,x3,x4,x5 >= 0  
 10 
 11 运行下面脚本(单纯形法的Python实现)
 12 
 13 import numpy as np
 14 
 15 # 实体类　Table
 16 # 控制类　Simplex
 17 
 18 class Table:
 19 
 20     def __init__(self):
 21 
 22         pass
 23 
 24     def set_para(self,A,b,c,base,z0):
 25 
 26         """
 27 
 28         输入LP必须已经化为标准形式
 29 
 30         """
 31 
 32         self.A=A
 33 
 34         self.b=b
 35 
 36         self.c=c
 37 
 38         self.z0=z0
 39 
 40         self.base=base
 41 
 42         self.m,self.n=self.A.shape
 43 
 44     def build(self):
 45 
 46         self.table=np.zeros((self.m+1,self.n+1))
 47 
 48         self.table[:-1,:1]=self.b.T
 49 
 50         self.table[-1 ,0]=self.z0
 51 
 52         self.table[:-1,1:]=self.A
 53 
 54         self.table[-1, 1:]=c
 55 
 56         self.baseVar=base
 57 
 58     def is_best(self):
 59 
 60         for sigma_index in range(self.n):
 61 
 62             if sigma_index not in self.baseVar:
 63 
 64                 sigma=self.table[-1,1+sigma_index]
 65 
 66                 if sigma>0:
 67 
 68                     return False
 69 
 70         return True
 71 
 72     def is_no_solution(self):
 73 
 74         for sigma_index in range(self.n):
 75 
 76             if sigma_index not in self.baseVar:
 77 
 78                 sigma=self.table[-1,1+sigma_index]
 79 
 80                 if sigma>0:
 81 
 82                     no_solution_flag=True
 83 
 84                     for a in self.table[:-1,1+sigma_index]:
 85 
 86                         if a>0:
 87 
 88                             no_solution_flag=False
 89 
 90                     if no_solution_flag==True:
 91 
 92                         return True
 93 
 94         return False
 95 
 96     def get_inVar(self):
 97 
 98         max_sigma=0
 99 
100         inVar=None
101 
102         for sigma_index in range(self.n):
103 
104             if sigma_index not in self.baseVar:
105 
106                 sigma=self.table[-1,1+sigma_index]
107 
108                 if sigma>max_sigma:
109 
110                     max_sigma=sigma
111 
112                     inVar=sigma_index
113 
114         return inVar
115 
116 
117     def get_outVar(self,inVar):
118 
119         rates=[]
120 
121         for nobaseVar in range(self.m):
122 
123             a=self.table[nobaseVar,inVar+1]
124 
125             b=self.table[nobaseVar,     0 ]
126 
127             if a>0:
128 
129                 rate=b/a
130 
131                 rates.append((rate,nobaseVar))
132 
133         return min(rates)[1]    
134 
135     def in_out(self,inVar,outVar):
136 
137         a=self.table[outVar,inVar+1]
138 
139         self.table[outVar,:]/=a
140 
141         for i in range(self.m+1):
142 
143             if i != outVar:
144 
145                 self.table[i,:]-=self.table[outVar,:]*self.table[i,inVar+1]
146 
147         self.baseVar[outVar]=inVar
148 
149  
150     def show(self):
151 
152         print ('基变量/取值：',self.baseVar,end='/')
153 
154         print (self.table[:-1,0])
155 
156         print ("单纯形表")
157 
158         for i in range(self.m+1):
159 
160             for j in range(self.n+1):
161 
162                 print ('%6.2f'%self.table[i,j],end=' ')
163 
164             print ()
165 
166         print ()
167 
168 class Simplex:
169 
170     def __init__(self):
171 
172         self.table=Table()
173 
174         # 0 正常，尚未得到最优解，继续迭代
175 
176         # 1 无解，无界解
177 
178         # 2 达到最优解
179 
180         self.status=0
181 
182         self.inVar=None
183 
184         self.outVar=None
185 
186     def set_para(self,A,b,c,base,z0=0):
187 
188         self.table.set_para(A,b,c,base,z0)
189 
190     def output_result(self):
191 
192         self._main()
193 
194         if self.status==1:
195 
196             print("此问题无界")
197 
198         elif self.status==2:
199 
200             print("此问题有一个最优解")
201 
202         elif self.status==3:
203 
204             print("此问题有无穷多个最优解")
205 
206     def _main(self):
207 
208         self._build_table()
209 
210         while 1:
211 
212             self.table.show()
213 
214             if self._is_best() or self._is_no_solution():
215 
216                 return
217  
218             self._get_inVar()
219 
220             self._get_outVar()
221 
222             self._in_out()
223 
224     def _build_table(self):
225 
226         self.table.build()
227 
228     def _is_best(self):
229 
230         if self.table.is_best():
231 
232             self.status=2
233 
234             return True
235 
236         return False
237 
238     def _is_no_solution(self):
239 
240         if self.table.is_no_solution():
241 
242             self.status=1
243 
244             return True
245 
246         return False
247 
248 
249     def _get_inVar(self):
250 
251         self.inVar=self.table.get_inVar()
252 
253 
254     def _get_outVar(self):
255 
256         self.outVar=self.table.get_outVar(self.inVar)
257  
258 
259     def _in_out(self):
260 
261         self.table.in_out(self.inVar,self.outVar)
262 
263  
264 
265 if __name__=="__main__":
266 
267     s=Simplex()
268 
269     A=np.matrix([[1,1,1,0,0],
270 
271                  [ 1, 0,0,1,0],
272 
273                  [ 0, 1,0,0,1]])
274 
275     b=np.matrix([[20,12,16]])
276 
277     c=np.matrix([[18,12.5,0,0,0]])
278 
279     base=[2,3,4]
280 
281     s.set_para(A,b,c,base)
282 
283 s.output_result()