上期談到用Python實現求解運輸問題 (Transportation Problem, TP) 的表上作業法的第一步——利用Vogel法尋找初始基可行解：

【Python實現】運輸問題的表上作業法（一）：利用伏格爾 (Vogel) 法尋找初始基可行解

這期來講講找到初始基可行解之后怎樣判斷當前解是否是最優解。如果當前解已達到最優，那么無需再進行操作；如果當前解非最優，那么還要對當前解進行調整以達到最優。調整解的操作放到下一期再講，本期先談談怎樣判斷最優性。

文章目錄

• 位勢法（對偶變量法）
• 位勢法應用實例
• 位勢法的Python語句

通常判斷TP問題解的最優性有兩種方法：閉回路法和位勢法。考慮到兩種方法代碼的復雜程度，我們重點講講更容易實現的位勢法。位勢法（Potentials）又稱對偶變量法，是利用線性規劃的對偶原理計算解的檢驗數，從而通過檢驗數判斷最優性的方法。關于對偶原理，請參考運籌學相關書籍，這里不進行闡述，而只對位勢法原理進行簡單介紹。

位勢法（對偶變量法）

位勢法應用實例

下面通過一個運輸問題的例題進一步說明位勢法原理。
Example：
已知該TP問題的運輸（成本）表如下（將中間3x4的成本表記為表1）：

通過最小元素法找到一組初始基可行解：

把初始基可行解中非零元素對應位置的成本系數剝離出來，并添加上位勢（圖中u₁-u₃、v₁-v₄），記為表2：

構建位勢方程組（需要注意的是，TP問題的可行解中非零變量個數為m+n-1，其中m、n分別是成本矩陣的行數和列數（本例中m=3,n=4），所以位勢方程組應該含有m+n-1個方程；但是位勢的個數為m+n，即方程組有m+n個未知量，多于方程的個數，可知方程組是解不出來的，因此我們隨便令某個位勢=0，即可求解這個方程組，并且此操作對最后檢驗數無影響）：

求出位勢后，將u_i+v_j填入成本表中對應解變量為0的位置，記為表3：

表1-表3（σ_ij=c_ij-(u_i+v_j)，對應位置元素相減），得到檢驗數表σ：

判斷最優性：由于本例是極小化問題，如果檢驗數表中存在σ_ij<0，說明當前解未達到最優，需要對解進行進一步的調整；否則，當前解達到最優（本例中當前初始基可行解并非最優解，因為σ₂₄=-1<0）。如果采用的檢驗數σ_ij=(u_i+v_j)-c_ij，則全部反向判斷即可。

位勢法的Python語句

同樣使用上面的例題，但在尋找初始基可行解時運用Vogel法，編制Python程序進行求解。

首先將運輸表寫入EXCEL表格并保存，文件名為TP_PPT_Sample1.xlsx：

然后執行以下代碼：

#Vogel法尋找初始基可行解+位勢法判斷解的最優性 import numpy as np import copy import pandas as pddef TP_split_matrix(mat):c=mat[:-1,:-1]a=mat[:-1,-1]b=mat[-1,:-1]return (c,a,b)def TP_vogel(var): #Vogel法代碼，變量var可以是以numpy.ndarray保存的運輸表，或以tuple或list保存的(成本矩陣,供給向量,需求向量)import numpytypevar1=type(var)==numpy.ndarraytypevar2=type(var)==tupletypevar3=type(var)==listif typevar1==False and typevar2==False and typevar3==False:print('>>>非法變量<<<')(cost,x)=(None,None)else:if typevar1==True:[c,a,b]=TP_split_matrix(var)elif typevar2==True or typevar3==True:[c,a,b]=varcost=copy.deepcopy(c)x=np.zeros(c.shape)M=pow(10,9)for factor in c.reshape(1,-1)[0]:while int(factor)!=M:if np.all(c==M):breakelse:print('c:\n',c)#獲取行/列最小值數組row_mini1=[]row_mini2=[]for row in range(c.shape[0]):Row=list(c[row,:])row_min=min(Row)row_mini1.append(row_min)Row.remove(row_min)row_2nd_min=min(Row)row_mini2.append(row_2nd_min)#print(row_mini1,'\n',row_mini2)r_pun=[row_mini2[i]-row_mini1[i] for i in range(len(row_mini1))]print('行罰數：',r_pun)#計算列罰數col_mini1=[]col_mini2=[]for col in range(c.shape[1]):Col=list(c[:,col])col_min=min(Col)col_mini1.append(col_min)Col.remove(col_min)col_2nd_min=min(Col)col_mini2.append(col_2nd_min)c_pun=[col_mini2[i]-col_mini1[i] for i in range(len(col_mini1))]print('列罰數：',c_pun)pun=copy.deepcopy(r_pun)pun.extend(c_pun)print('罰數向量：',pun)max_pun=max(pun)max_pun_index=pun.index(max(pun))max_pun_num=max_pun_index+1print('最大罰數：',max_pun,'元素序號：',max_pun_num)if max_pun_num<=len(r_pun):row_num=max_pun_numprint('對第',row_num,'行進行操作：')row_index=row_num-1catch_row=c[row_index,:]print(catch_row)min_cost_colindex=int(np.argwhere(catch_row==min(catch_row)))print('最小成本所在列索引：',min_cost_colindex)if a[row_index]<=b[min_cost_colindex]:x[row_index,min_cost_colindex]=a[row_index]c1=copy.deepcopy(c)c1[row_index,:]=[M]*c1.shape[1]b[min_cost_colindex]-=a[row_index]a[row_index]-=a[row_index]else:x[row_index,min_cost_colindex]=b[min_cost_colindex]c1=copy.deepcopy(c)c1[:,min_cost_colindex]=[M]*c1.shape[0]a[row_index]-=b[min_cost_colindex]b[min_cost_colindex]-=b[min_cost_colindex]else:col_num=max_pun_num-len(r_pun)col_index=col_num-1print('對第',col_num,'列進行操作：')catch_col=c[:,col_index]print(catch_col)#尋找最大罰數所在行/列的最小成本系數min_cost_rowindex=int(np.argwhere(catch_col==min(catch_col)))print('最小成本所在行索引：',min_cost_rowindex)#計算將該位置應填入x矩陣的數值（a,b中較小值）if a[min_cost_rowindex]<=b[col_index]:x[min_cost_rowindex,col_index]=a[min_cost_rowindex]c1=copy.deepcopy(c)c1[min_cost_rowindex,:]=[M]*c1.shape[1]b[col_index]-=a[min_cost_rowindex]a[min_cost_rowindex]-=a[min_cost_rowindex]else:x[min_cost_rowindex,col_index]=b[col_index]#填入后刪除已滿足/耗盡資源系數的行/列，得到剩余的成本矩陣，并改寫資源系數c1=copy.deepcopy(c)c1[:,col_index]=[M]*c1.shape[0]a[min_cost_rowindex]-=b[col_index]b[col_index]-=b[col_index]c=c1print('本次迭代后的x矩陣：\n',x)print('a:',a)print('b:',b)print('c:\n',c)if np.all(c==M):print('【迭代完成】')print('-'*60)else:print('【迭代未完成】')print('-'*60)total_cost=np.sum(np.multiply(x,cost))if np.all(a==0):if np.all(b==0):print('>>>供求平衡<<<')else:print('>>>供不應求，需求方有余量<<<')elif np.all(b==0):print('>>>供大于求，供給方有余量<<<')else:print('>>>無法找到初始基可行解<<<')print('>>>初始基本可行解x*：\n',x)print('>>>當前總成本：',total_cost)[m,n]=x.shapevarnum=np.array(np.nonzero(x)).shape[1]if varnum!=m+n-1:print('【注意：問題含有退化解】')return (cost,x)def create_c_nonzeros(c,x):import numpy as npimport copynonzeros=copy.deepcopy(x)for i in range(x.shape[0]):for j in range(x.shape[1]):if x[i,j]!=0:nonzeros[i,j]=1#print(nonzeros)c_nonzeros=np.multiply(c,nonzeros)return c_nonzerosdef if_dsquare(a,b):print('a:',a.shape,'\n','b:',b.shape)correct='>>>位勢方程組可解<<<'error='>>>位勢方程組不可解，請檢查基變量個數是否等于(m+n-1)及位勢未知量個數是否等于(m+n)<<<'if len(a.shape)==2:if len(b.shape)==2:if a.shape[0]==a.shape[1] and a.shape==b.shape:print(correct)if_dsquare=Trueelse:print(error)if_dsquare=Falseelif len(b.shape)==1 and b.shape[0]!=0:if a.shape[0]==a.shape[1] and a.shape[0]==b.shape[0]:print(correct)if_dsquare=Trueelse:print(error)if_dsquare=Falseelse:print(error)if_dsquare=Falseelif len(a.shape)==1:if len(b.shape)==2:if b.shape[0]==b.shape[1] and a.shape[0]==b.shape[0]:print('>>>位勢方程組系數矩陣與方程組值向量位置錯誤<<<')if_dsquare='True but not solvable'else:print(error)if_dsquare=Falseelif len(b.shape)==1:print(error)if_dsquare=Falseelse:print(error)if_dsquare=Falseelse:print(error)if_dsquare=Falsereturn if_dsquaredef TP_potential(cost,x):[m,n]=x.shapevarnum=np.array(np.nonzero(x)).shape[1]if varnum!=m+n-1:sigma=Noneprint('【問題含有退化解，暫時無法判斷最優性】')else:#print(c_nonzeros.shape)c_nonzeros=create_c_nonzeros(c,x)cc_nonzeros=np.array(np.nonzero(c_nonzeros))A=[]B=[]length=c_nonzeros.shape[0]+c_nonzeros.shape[1]zeros=np.zeros((1,length))[0]for i in range(cc_nonzeros.shape[1]):zeros=np.zeros((1,length))[0]zeros[cc_nonzeros[0,i]]=1zeros[cc_nonzeros[1,i]+c_nonzeros.shape[0]]=1A.append(zeros)B.append(c_nonzeros[cc_nonzeros[0,i],cc_nonzeros[1,i]])l=[1]for j in range(length-1):l.append(0) #補充一個x1=0的方程以滿足求解條件A.append(l)B.append(0)#print(A)#print(B)A=np.array(A)B=np.array(B)judge=if_dsquare(A,B)if judge==True:sol=np.linalg.solve(A,B) #求解條件：A的行數（方程個數)=A的列數(變量個數）=B的個數（方程結果個數）才能解#print(sol)var=[] #創建位勢名稱數組for i in range(c_nonzeros.shape[0]):var.append('u'+str(i+1))for j in range(c_nonzeros.shape[1]):var.append('v'+str(j+1))#print(var)solset=dict(zip(var,sol))print('>>>當前位勢:\n',solset)u=[]v=[][m,n]=c_nonzeros.shapezero_places=np.transpose(np.argwhere(c_nonzeros==0))for i in range(m):u.append(sol[i])for j in range(n):v.append(sol[j+m])for k in range(zero_places.shape[1]):c_nonzeros[zero_places[0,k],zero_places[1,k]]=u[zero_places[0,k]]+v[zero_places[1,k]]#print(c_nonzeros)sigma=cost-c_nonzerosprint('>>>檢驗表σ：\n',sigma)if np.all(sigma>=0):print('>>>TP已達到最優解<<<')else:print('>>>TP未達到最優解<<<')else:sigma=Noneprint('>>>位勢方程組不可解<<<')return sigmapath=r'C:\Users\spurs\Desktop\MCM_ICM\Data files\TP_PPT_Sample1.xlsx' mat=pd.read_excel(path,header=None).values #c=np.array([[3,11,3,10],[1,9,2,8],[7,4,10,5]]) #a=np.array([7,4,9]) #b=np.array([3,6,5,6]) [c,x]=TP_vogel(mat) #[c,x]=TP_vogel([c,a,b]) sigma=TP_potential(c,x)

運行結果：

c:[[ 3. 11. 3. 10.][ 1. 9. 2. 8.][ 7. 4. 10. 5.]] 行罰數： [0.0, 1.0, 1.0] 列罰數： [2.0, 5.0, 1.0, 3.0] 罰數向量： [0.0, 1.0, 1.0, 2.0, 5.0, 1.0, 3.0] 最大罰數： 5.0 元素序號： 5 對第 2 列進行操作： [11. 9. 4.] 最小成本所在行索引： 2 本次迭代后的x矩陣：[[0. 0. 0. 0.][0. 0. 0. 0.][0. 6. 0. 0.]] a: [7. 4. 3.] b: [3. 0. 5. 6.] c:[[3.e+00 1.e+09 3.e+00 1.e+01][1.e+00 1.e+09 2.e+00 8.e+00][7.e+00 1.e+09 1.e+01 5.e+00]] 【迭代未完成】 ------------------------------------------------------------ c:[[3.e+00 1.e+09 3.e+00 1.e+01][1.e+00 1.e+09 2.e+00 8.e+00][7.e+00 1.e+09 1.e+01 5.e+00]] 行罰數： [0.0, 1.0, 2.0] 列罰數： [2.0, 0.0, 1.0, 3.0] 罰數向量： [0.0, 1.0, 2.0, 2.0, 0.0, 1.0, 3.0] 最大罰數： 3.0 元素序號： 7 對第 4 列進行操作： [10. 8. 5.] 最小成本所在行索引： 2 本次迭代后的x矩陣：[[0. 0. 0. 0.][0. 0. 0. 0.][0. 6. 0. 3.]] a: [7. 4. 0.] b: [3. 0. 5. 3.] c:[[3.e+00 1.e+09 3.e+00 1.e+01][1.e+00 1.e+09 2.e+00 8.e+00][1.e+09 1.e+09 1.e+09 1.e+09]] 【迭代未完成】 ------------------------------------------------------------ c:[[3.e+00 1.e+09 3.e+00 1.e+01][1.e+00 1.e+09 2.e+00 8.e+00][1.e+09 1.e+09 1.e+09 1.e+09]] 行罰數： [0.0, 1.0, 0.0] 列罰數： [2.0, 0.0, 1.0, 2.0] 罰數向量： [0.0, 1.0, 0.0, 2.0, 0.0, 1.0, 2.0] 最大罰數： 2.0 元素序號： 4 對第 1 列進行操作： [3.e+00 1.e+00 1.e+09] 最小成本所在行索引： 1 本次迭代后的x矩陣：[[0. 0. 0. 0.][3. 0. 0. 0.][0. 6. 0. 3.]] a: [7. 1. 0.] b: [0. 0. 5. 3.] c:[[1.e+09 1.e+09 3.e+00 1.e+01][1.e+09 1.e+09 2.e+00 8.e+00][1.e+09 1.e+09 1.e+09 1.e+09]] 【迭代未完成】 ------------------------------------------------------------ c:[[1.e+09 1.e+09 3.e+00 1.e+01][1.e+09 1.e+09 2.e+00 8.e+00][1.e+09 1.e+09 1.e+09 1.e+09]] 行罰數： [7.0, 6.0, 0.0] 列罰數： [0.0, 0.0, 1.0, 2.0] 罰數向量： [7.0, 6.0, 0.0, 0.0, 0.0, 1.0, 2.0] 最大罰數： 7.0 元素序號： 1 對第 1 行進行操作： [1.e+09 1.e+09 3.e+00 1.e+01] 最小成本所在列索引： 2 本次迭代后的x矩陣：[[0. 0. 5. 0.][3. 0. 0. 0.][0. 6. 0. 3.]] a: [2. 1. 0.] b: [0. 0. 0. 3.] c:[[1.e+09 1.e+09 1.e+09 1.e+01][1.e+09 1.e+09 1.e+09 8.e+00][1.e+09 1.e+09 1.e+09 1.e+09]] 【迭代未完成】 ------------------------------------------------------------ c:[[1.e+09 1.e+09 1.e+09 1.e+01][1.e+09 1.e+09 1.e+09 8.e+00][1.e+09 1.e+09 1.e+09 1.e+09]] 行罰數： [999999990.0, 999999992.0, 0.0] 列罰數： [0.0, 0.0, 0.0, 2.0] 罰數向量： [999999990.0, 999999992.0, 0.0, 0.0, 0.0, 0.0, 2.0] 最大罰數： 999999992.0 元素序號： 2 對第 2 行進行操作： [1.e+09 1.e+09 1.e+09 8.e+00] 最小成本所在列索引： 3 本次迭代后的x矩陣：[[0. 0. 5. 0.][3. 0. 0. 1.][0. 6. 0. 3.]] a: [2. 0. 0.] b: [0. 0. 0. 2.] c:[[1.e+09 1.e+09 1.e+09 1.e+01][1.e+09 1.e+09 1.e+09 1.e+09][1.e+09 1.e+09 1.e+09 1.e+09]] 【迭代未完成】 ------------------------------------------------------------ c:[[1.e+09 1.e+09 1.e+09 1.e+01][1.e+09 1.e+09 1.e+09 1.e+09][1.e+09 1.e+09 1.e+09 1.e+09]] 行罰數： [999999990.0, 0.0, 0.0] 列罰數： [0.0, 0.0, 0.0, 999999990.0] 罰數向量： [999999990.0, 0.0, 0.0, 0.0, 0.0, 0.0, 999999990.0] 最大罰數： 999999990.0 元素序號： 1 對第 1 行進行操作： [1.e+09 1.e+09 1.e+09 1.e+01] 最小成本所在列索引： 3 本次迭代后的x矩陣：[[0. 0. 5. 2.][3. 0. 0. 1.][0. 6. 0. 3.]] a: [0. 0. 0.] b: [0. 0. 0. 0.] c:[[1.e+09 1.e+09 1.e+09 1.e+09][1.e+09 1.e+09 1.e+09 1.e+09][1.e+09 1.e+09 1.e+09 1.e+09]] 【迭代完成】 ------------------------------------------------------------ >>>供求平衡<<< >>>初始基本可行解x*：[[0. 0. 5. 2.][3. 0. 0. 1.][0. 6. 0. 3.]] >>>當前總成本： 85.0 >>>位勢方程組可解<<< >>>當前位勢:{'u1': 0.0, 'u2': -2.0, 'u3': -5.0, 'v1': 3.0, 'v2': 9.0, 'v3': 3.0, 'v4': 10.0} >>>檢驗數表σ：[[ 0. 2. 0. 0.][ 0. 2. 1. 0.][ 9. 0. 12. 0.]] >>>TP已達到最優解<<<

可見，這次我們用Vogel法找到的初始基可行解已經是最優解（檢驗數表σ中元素均大于0），再次證明Vogel法比最小元素法求解效率更高。

對于存在退化現象的問題，由于筆者能力原因，暫時沒有找到最優解決方案，因此程序會提示暫時無法判斷最優性：

Example:

Result:

>>>供求平衡<<< >>>初始基本可行解x*：[[1. 0. 0. 0. 0.][0. 0. 1. 0. 0.][0. 0. 0. 1. 0.][0. 0. 0. 0. 1.][0. 1. 0. 0. 0.]] >>>當前總成本： 0.8999999999999999 【注意：問題含有退化解】 【問題含有退化解，暫時無法判斷最優性】

對于此類問題的解決方案筆者會在將來進行補充。

函數使用方法：
TP_potential (c,x): c是成本矩陣，x是當前可行解矩陣。
本函數返回的是當前解的檢驗數矩陣。

感謝閱讀，歡迎指正。

總結

以上是生活随笔為你收集整理的【Python实现】运输问题的表上作业法（二）：利用位势法判断当前解的最优性的全部內容，希望文章能夠幫你解決所遇到的問題。

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