数学建模——运输问题（Python实现）

目录

1、概述

（1）运输问题

（2）基本思想 

（3）表上作业法求解运输问题步骤 

2、知识点细讲 

（1）运输问题及其数学模型

​（2）表上作业法求解运输问题——思想 

（3）表上作业法求解运输问题——寻找基可行解 

（4）表上作业法求解运输问题——解的最优检验 

（5）表上作业法求解运输问题——解的改进

（6）产销不平衡的运输问题 

​（7）有转运的运输问题 

3、案例分析1——沃格尔法（寻找初始基可行解）

（1）案例分析

（2）Python实现 

（3）结果 

 4、案例分析2——位势法（解的最优检验）

（1）案例分析

​ （2）Python实现 

（3）结果

 5、致谢

---

1、概述

（1）运输问题

（2）基本思想 

（3）表上作业法求解运输问题步骤 

2、知识点细讲 

（1）运输问题及其数学模型

 

（2）表上作业法求解运输问题——思想 

（3）表上作业法求解运输问题——寻找基可行解 

（4）表上作业法求解运输问题——解的最优检验 

 

 

 

 

（5）表上作业法求解运输问题——解的改进

 

（6）产销不平衡的运输问题 

（7）有转运的运输问题 

     3、案例分析1——沃格尔法（寻找初始基可行解）

（1）案例分析

 

（2）Python实现 

#运输问题求解：使用Vogel逼近法寻找初始基本可行解
import numpy as np
import copy
import pandas as pd

def main():
    mat=pd.read_csv('表上作业法求解运输问题.csv',header=None).values
    #mat = pd.read_excel('表上作业法求解运输问题.xlsx', header=None).values
    #c=np.array([[4,12,4,11],[2,10,3,9],[8,5,11,6]]) #成本矩阵
    #a=np.array([16,10,22])                          #供给量
    #b=np.array([8,14,12,14])                        #需求量
    [c,x]=TP_vogel(mat)
    #[c,x]=TP_vogel([c,a,b])

def 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 numpy
    typevar1=type(var)==numpy.ndarray
    typevar2=type(var)==tuple
    typevar3=type(var)==list
    if 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]=var
        cost=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):
                    break
                else:
                    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+1
                    print('最大罚数：',max_pun,'元素序号：',max_pun_num)
                    if max_pun_num<=len(r_pun):
                        row_num=max_pun_num
                        print('对第',row_num,'行进行操作：')
                        row_index=row_num-1
                        catch_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-1
                        print('对第',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=c1
                    print('本次迭代后的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.shape
        varnum=np.array(np.nonzero(x)).shape[1]
        if varnum!=m+n-1:
            print('【注意：问题含有退化解】')
    return (cost,x)

if __name__ =='__main__':
    main()

[注意]（1）copy模块、（2）mat的操作、（3）在使用pandas导入文件数据的时候，运行read_csv正常，运行的read_excel 出现错误，直接按照提示安装openpyxl。

（3）结果 

c:
 [[ 4. 12.  4. 11.]
 [ 2. 10.  3.  9.]
 [ 8.  5. 11.  6.]]
行罚数： [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 列进行操作：
[12. 10.  5.]
最小成本所在行索引： 2
本次迭代后的x矩阵：
 [[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0. 14.  0.  0.]]
a: [16. 10.  8.]
b: [ 8.  0. 12. 14.]
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [8.0e+00 1.0e+09 1.1e+01 6.0e+00]]
【迭代未完成】
------------------------------------------------------------
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [8.0e+00 1.0e+09 1.1e+01 6.0e+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 列进行操作：
[11.  9.  6.]
最小成本所在行索引： 2
本次迭代后的x矩阵：
 [[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0. 14.  0.  8.]]
a: [16. 10.  0.]
b: [ 8.  0. 12.  6.]
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+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 列进行操作：
[4.e+00 2.e+00 1.e+09]
最小成本所在行索引： 1
本次迭代后的x矩阵：
 [[ 0.  0.  0.  0.]
 [ 8.  0.  0.  0.]
 [ 0. 14.  0.  8.]]
a: [16.  2.  0.]
b: [ 0.  0. 12.  6.]
c:
 [[1.0e+09 1.0e+09 4.0e+00 1.1e+01]
 [1.0e+09 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[1.0e+09 1.0e+09 4.0e+00 1.1e+01]
 [1.0e+09 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+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.0e+09 1.0e+09 4.0e+00 1.1e+01]
最小成本所在列索引： 2
本次迭代后的x矩阵：
 [[ 0.  0. 12.  0.]
 [ 8.  0.  0.  0.]
 [ 0. 14.  0.  8.]]
a: [4. 2. 0.]
b: [0. 0. 0. 6.]
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
行罚数： [999999989.0, 999999991.0, 0.0]
列罚数： [0.0, 0.0, 0.0, 2.0]
罚数向量： [999999989.0, 999999991.0, 0.0, 0.0, 0.0, 0.0, 2.0]
最大罚数： 999999991.0 元素序号： 2
对第 2 行进行操作：
[1.e+09 1.e+09 1.e+09 9.e+00]
最小成本所在列索引： 3
本次迭代后的x矩阵：
 [[ 0.  0. 12.  0.]
 [ 8.  0.  0.  2.]
 [ 0. 14.  0.  8.]]
a: [4. 0. 0.]
b: [0. 0. 0. 4.]
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
行罚数： [999999989.0, 0.0, 0.0]
列罚数： [0.0, 0.0, 0.0, 999999989.0]
罚数向量： [999999989.0, 0.0, 0.0, 0.0, 0.0, 0.0, 999999989.0]
最大罚数： 999999989.0 元素序号： 1
对第 1 行进行操作：
[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
最小成本所在列索引： 3
本次迭代后的x矩阵：
 [[ 0.  0. 12.  4.]
 [ 8.  0.  0.  2.]
 [ 0. 14.  0.  8.]]
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. 12.  4.]
 [ 8.  0.  0.  2.]
 [ 0. 14.  0.  8.]]
>>>当前总成本： 244.0

Process finished with exit code 0

 4、案例分析2——位势法（解的最优检验）

（1）案例分析

 （2）Python实现 

#Vogel法寻找初始基可行解+位势法判断解的最优性
import numpy as np
import copy
import pandas as pd

def 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 numpy
    typevar1=type(var)==numpy.ndarray
    typevar2=type(var)==tuple
    typevar3=type(var)==list
    if 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]=var
        cost=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):
                    break
                else:
                    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+1
                    print('最大罚数：',max_pun,'元素序号：',max_pun_num)
                    if max_pun_num<=len(r_pun):
                        row_num=max_pun_num
                        print('对第',row_num,'行进行操作：')
                        row_index=row_num-1
                        catch_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-1
                        print('对第',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=c1
                    print('本次迭代后的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.shape
        varnum=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 np
    import copy
    nonzeros=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_nonzeros

def 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=True
            else:
                print(error)
                if_dsquare=False
        elif 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=True
            else:
                print(error)
                if_dsquare=False
        else:
            print(error)
            if_dsquare=False
    elif 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=False
        elif len(b.shape)==1:
                print(error)
                if_dsquare=False
        else:
            print(error)
            if_dsquare=False
    else:
        print(error)
        if_dsquare=False
    return if_dsquare

def TP_potential(cost,x):
    [m,n]=x.shape
    varnum=np.array(np.nonzero(x)).shape[1]
    if varnum!=m+n-1:
        sigma=None
        print('【问题含有退化解，暂时无法判断最优性】')
    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]]=1
            zeros[cc_nonzeros[1,i]+c_nonzeros.shape[0]]=1
            A.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.shape
            zero_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_nonzeros
            print('>>>检验表σ：\n',sigma)
            if np.all(sigma>=0):
                print('>>>TP已达到最优解<<<')
            else:
                print('>>>TP未达到最优解<<<')
        else:
            sigma=None
            print('>>>位势方程组不可解<<<')
    return sigma

mat = pd.read_excel('表上作业法求解运输问题.xlsx', header=None).values
# mat = pd.read_csv('表上作业法求解运输问题.csv', header=None).values
# c=np.array([[4,12,4,11],[2,10,3,9],[8,5,11,6]])
# a=np.array([16,10,22])
# b=np.array([8,14,12,14])
[c, x] = TP_vogel(mat)
# [c,x]=TP_vogel([c,a,b])
sigma= TP_potential(c, x)

（3）结果

c:
 [[ 4. 12.  4. 11.]
 [ 2. 10.  3.  9.]
 [ 8.  5. 11.  6.]]
行罚数： [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 列进行操作：
[12. 10.  5.]
最小成本所在行索引： 2
本次迭代后的x矩阵：
 [[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0. 14.  0.  0.]]
a: [16. 10.  8.]
b: [ 8.  0. 12. 14.]
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [8.0e+00 1.0e+09 1.1e+01 6.0e+00]]
【迭代未完成】
------------------------------------------------------------
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [8.0e+00 1.0e+09 1.1e+01 6.0e+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 列进行操作：
[11.  9.  6.]
最小成本所在行索引： 2
本次迭代后的x矩阵：
 [[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0. 14.  0.  8.]]
a: [16. 10.  0.]
b: [ 8.  0. 12.  6.]
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[4.0e+00 1.0e+09 4.0e+00 1.1e+01]
 [2.0e+00 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+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 列进行操作：
[4.e+00 2.e+00 1.e+09]
最小成本所在行索引： 1
本次迭代后的x矩阵：
 [[ 0.  0.  0.  0.]
 [ 8.  0.  0.  0.]
 [ 0. 14.  0.  8.]]
a: [16.  2.  0.]
b: [ 0.  0. 12.  6.]
c:
 [[1.0e+09 1.0e+09 4.0e+00 1.1e+01]
 [1.0e+09 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[1.0e+09 1.0e+09 4.0e+00 1.1e+01]
 [1.0e+09 1.0e+09 3.0e+00 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+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.0e+09 1.0e+09 4.0e+00 1.1e+01]
最小成本所在列索引： 2
本次迭代后的x矩阵：
 [[ 0.  0. 12.  0.]
 [ 8.  0.  0.  0.]
 [ 0. 14.  0.  8.]]
a: [4. 2. 0.]
b: [0. 0. 0. 6.]
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 9.0e+00]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
行罚数： [999999989.0, 999999991.0, 0.0]
列罚数： [0.0, 0.0, 0.0, 2.0]
罚数向量： [999999989.0, 999999991.0, 0.0, 0.0, 0.0, 0.0, 2.0]
最大罚数： 999999991.0 元素序号： 2
对第 2 行进行操作：
[1.e+09 1.e+09 1.e+09 9.e+00]
最小成本所在列索引： 3
本次迭代后的x矩阵：
 [[ 0.  0. 12.  0.]
 [ 8.  0.  0.  2.]
 [ 0. 14.  0.  8.]]
a: [4. 0. 0.]
b: [0. 0. 0. 4.]
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
【迭代未完成】
------------------------------------------------------------
c:
 [[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]
 [1.0e+09 1.0e+09 1.0e+09 1.0e+09]]
行罚数： [999999989.0, 0.0, 0.0]
列罚数： [0.0, 0.0, 0.0, 999999989.0]
罚数向量： [999999989.0, 0.0, 0.0, 0.0, 0.0, 0.0, 999999989.0]
最大罚数： 999999989.0 元素序号： 1
对第 1 行进行操作：
[1.0e+09 1.0e+09 1.0e+09 1.1e+01]
最小成本所在列索引： 3
本次迭代后的x矩阵：
 [[ 0.  0. 12.  4.]
 [ 8.  0.  0.  2.]
 [ 0. 14.  0.  8.]]
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. 12.  4.]
 [ 8.  0.  0.  2.]
 [ 0. 14.  0.  8.]]
>>>当前总成本： 244.0
a: (7, 7) 
 b: (7,)
>>>位势方程组可解<<<
>>>当前位势:
 {'u1': 0.0, 'u2': -2.0, 'u3': -5.0, 'v1': 4.0, 'v2': 10.0, 'v3': 4.0, 'v4': 11.0}
>>>检验表σ：
 [[ 0.  2.  0.  0.]
 [ 0.  2.  1.  0.]
 [ 9.  0. 12.  0.]]
>>>TP已达到最优解<<<

Process finished with exit code 0

 5、致谢

https://blog.csdn.net/m0_46927406/article/details/109903843