相關理論知識參考 單純形理論知識

算法可以在給定一個包含線性規劃問題的標準形式的描述下，求解該線性規劃問題。
例如某一個 pro.txt 文件內容如下:

6
3
3 -1 1 -2 0 0
2 1 0 1 1 0
-1 3 0 -3 0 1
-3 4 12
-7 7 -2 -1 -6 0

執行算法之后得到結果：

x_1 = 0.000000,x_2 = 0.000000,x_3 = 0.000000,x_4 = 1.500000,x_5 = 2.500000,x_6 = 16.500000
objective value is : -16.500000
1-th line constraint value is : -3.000000
2-th line constraint value is : 4.000000
3-th line constraint value is : 12.000000

代碼如下：

#encoding=utf-8 __author__ = 'ysg' import numpy as np #python 矩陣操作libclass Simplex():def __init__(self):self._A = "" # 系數矩陣self._b = "" #self._c = '' #約束self._B = '' #基變量的下標集合self.row = 0 #約束個數def solve(self, filename):#讀取文件內容，文件結構前兩行分別為 變量數 和 約束條件個數#接下來是系數矩陣#然后是b數組#然后是約束條件c#假設線性規劃形式是標準形式（都是等式）A = []b = []c = []with open(filename,'r') as f:self.var = int(f.readline())self.row = int(f.readline())for i in range(self.row):x =map(int, f.readline().strip().split(' '))A.append(x)b=(map(int, list(f.readline().split(' '))))c=(map(int, list(f.readline().split(' '))))self._A = np.array(A, dtype=float)self._b = np.array(b, dtype=float)self._c = np.array(c, dtype=float)# self._A = np.array([[3,-1,1,-2,0,0],[2,1,0,1,1,0],[-1,3,0,-3,0,1]],dtype=float)# self._b = np.array([-3,4,12],dtype=float)# self._c = np.array([-7, 7, -2, -1, -6, 0],dtype=float)self._B = list()self.row = len(self._b)self.var = len(self._c)(x,obj) = self.Simplex(self._A,self._b,self._c)self.pprint(x,obj,A)def pprint(self,x,obj,A):px = ['x_%d = %f'%(i+1,x[i]) for i in range(len(x))]print ','.join(px)print ('objective value is : %f'% obj)print '------------------------------'for i in range(len(A)):print '%d-th line constraint value is : %f' % (i+1, x.dot(A[i]))#添加人工變量得到一個初始解def InitializeSimplex(self,A,b):b_min, min_pos = (np.min(b), np.argmin(b)) # 得到最小bi#將bi全部轉化成正數if (b_min < 0):for i in range(self.row):if i != min_pos:A[i] = A[i] - A[min_pos]b[i] = b[i] - b[min_pos]A[min_pos] = A[min_pos]*-1b[min_pos] = b[min_pos]*-1#添加松弛變量slacks = np.eye(self.row)A = np.concatenate((A,slacks),axis=1)c = np.concatenate((np.zeros(self.var),np.ones(self.row)),axis=1)# 松弛變量全部加入基,初始解為bnew_B = [i + self.var for i in range(self.row)]#輔助方程的目標函數值obj = np.sum(b)c = c[new_B].reshape(1,-1).dot(A) - cc = c[0]#entering basise= np.argmax(c)while c[e] > 0:theta = []for i in range(len(b)):if A[i][e] > 0:theta.append(b[i]/A[i][e])else:theta.append(float("inf"))l = np.argmin(np.array(theta))if theta[l] == float('inf'):print 'unbounded'return False(new_B, A, b, c, obj) = self._PIVOT(new_B, A, b, c, obj, l , e)e = np.argmax(c)#如果此時人工變量仍在基中，用原變量去替換之for mb in new_B:if mb >= self.var:row = mb-self.vari = 0while A[row][i] == 0 and i < self.var:i+=1(new_B, A, b, c, obj) = self._PIVOT(new_B,A,b,c,obj,new_B.index(mb),i)return (new_B, A[:,0:self.var], b)#算法入口def Simplex(self,A,b,c):B = ''(B, A ,b) = self.InitializeSimplex(A,b)#函數目標值obj = np.dot(c[B],b)c = np.dot(c[B].reshape(1,-1), A) - cc = c[0]# entering basise = np.argmax(c)# 找到最大的檢驗數，如果大于0，則目標函數可以優化while c[e] > 0:theta = []for i in range(len(b)):if A[i][e] > 0:theta.append(b[i] / A[i][e])else:theta.append(float("inf"))l = np.argmin(np.array(theta))if theta[l] == float('inf'):print 'unbounded'return False(B, A, b, c, obj) = self._PIVOT(B, A, b, c, obj, l, e)e = np.argmax(c)x = self._CalculateX(B,A,b,c)return (x,obj)#得到完整解def _CalculateX(self,B,A,b,c):x = np.zeros(self.var,dtype=float)x[B] = breturn x# 基變換def _PIVOT(self,B,A,b,c,z,l,e):# main element is a_le# l represents leaving basis# e represents entering basismain_elem = A[l][e]#scaling the l-th lineA[l] = A[l]/main_elemb[l] = b[l]/main_elem#change e-th column to unit arrayfor i in range(self.row):if i != l:b[i] = b[i] - A[i][e] * b[l]A[i] = A[i] - A[i][e] * A[l]#update objective valuez -= b[l]*c[e]c = c - c[e] * A[l]# change the basisB[l] = ereturn (B, A, b, c, z)s = Simplex() s.solve('pro.txt')

注釋差不多寫的比較清楚，可以參考上面的理論知識鏈接。

大佬請繞路。

有錯歡迎指出。

新人創作打卡挑戰賽發博客就能抽獎！定制產品紅包拿不停！

總結

以上是生活随笔為你收集整理的Simplex 单纯形算法的python实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。