运输问题

import numpy as np
from collections import deque
arr=input("输入各产地产量:\n")
a=[float(n)for n in arr.split()]
arr=input("输入各销地销量:\n")
b=[float(n)for n in arr.split()]
print("单位运价表:")
c=[]
for i in a:
    c.append(list(map(float, input().rstrip().split())))
str=[]
nbs=0
#产销不平衡：添加产地或销地
if sum(a) > sum(b):
    _add = np.zeros(len(a))
    c=np.column_stack((c,_add))#添加一个销售地
    b.append(sum(a)-sum(b))
elif sum(a) <sum(b):
    _add=np.zeros(len(b))
    c=np.row_stack((c,_add))#添加一个产地
    a.append(sum(b)-sum(a))
#确定初始可行基
_c=np.array(c)
c_test=np.full((_c.shape[0],_c.shape[1]),np.nan)
#西北角法
def ifb_1():
    for i in range(_c.shape[0]):
        for j in range(_c.shape[1]):
            if a[i]>0 and b[j]>0:
                if a[i] > b[j]:
                    c_test[i, j] = b[j]
                    a[i] = a[i] - c_test[i, j]
                    b[j] = b[j] - c_test[i, j]
                else:
                    c_test[i, j] = a[i]
                    a[i] = a[i] - c_test[i, j]
                    b[j] = b[j] - c_test[i, j]
#最小元素法（注意：适用于产销平衡问题）
def ifb_2():
    for k in range(_c.shape[0]+_c.shape[1]-1):
        min_index = np.argwhere(_c == np.min(_c))
        i=min_index[0,0]
        j=min_index[0,1]
        if a[i] > 0 and b[j] > 0:
            if a[i] > b[j]:
                c_test[i, j] = b[j]
                a[i] = a[i] - c_test[i, j]
                b[j] = b[j] - c_test[i, j]
                _c[:,j]=np.max(c)*_c.shape[0]
            else:
                c_test[i, j] = a[i]
                a[i] = a[i] - c_test[i, j]
                b[j] = b[j] - c_test[i, j]
                _c[i,:] = np.max(c) * _c.shape[1]
#函数：获取行差额和列差额
M=pow(10,9)
r_Dif=np.zeros(_c.shape[0])
c_Dif=np.zeros(_c.shape[1])
def get_dif():
    for row in range(_c.shape[0]):
            Row = sorted(_c[row, :])
            if Row[1]==M:
                r_min=-1
            else:r_min = Row[1] - Row[0]
            r_Dif[row] = r_min
    for column in range(_c.shape[1]):
            Column = sorted(_c[:, column])
            if Column[1]==M:
                c_min=-1
            else:c_min = Column[1] - Column[0]
            c_Dif[column] = c_min
#伏格尔法
def ifb_3():
    dgs=[]
    for k in range(_c.shape[0]+_c.shape[1]-2):
        get_dif()
        if max(r_Dif) > max(c_Dif):
            r_in=np.argwhere(r_Dif==max(r_Dif))[0][0]
            ren=np.argwhere(_c[r_in,:]==min(_c[r_in,:]))[0][0]
            min_rc=min(a[r_in],b[ren])
            c_test[r_in,ren]=min_rc
            a[r_in] -= min_rc
            b[ren] -= min_rc
            if a[r_in]==0 :
                _c[r_in, :] = [M] * _c.shape[1]
            else :
                _c[:,ren]=[M]*_c.shape[0]

        else :
            c_in=np.argwhere(c_Dif==max(c_Dif))[0][0]
            ren = np.argwhere(_c[:,c_in] == min(_c[:,c_in]))[0][0]
            min_cr=min(a[ren],b[c_in])
            c_test[ren,c_in]=min_cr
            a[ren]-=min_cr
            b[c_in]-=min_cr
            if b[c_in]==0 :
                _c[:, c_in] = [M] * _c.shape[0]
            else:
                _c[ren,:]=[M]*_c.shape[1]
    r_in = np.argwhere(_c != M)[0][0]
    c_in = np.argwhere(_c != M)[0][1]
    c_test[r_in, c_in] = min(a[r_in], b[c_in])
#位势法（potential method）
c_ = np.array(c)
c_pot = np.full((_c.shape[0] + 1, _c.shape[1] + 1), np.nan)
def potential_method():
    #新建矩阵，存放检验数
    c_pot[0, -1] = 0#u1=0
    #添加vj
    def add_v(r, c, c_, c_pot):
        c_pot[-1, c] = c_[r, c] - c_pot[r, -1]
        return c_pot
    #添加ui
    def add_u(r, c, c_, c_pot):
        c_pot[r, -1] = c_[r, c] - c_pot[-1, c]
        return c_pot
    #存放读取的初始可行基
    def add_stack(r, c, st):
        if r == 0:
            for i in range(len(ini)):
                if ini[i][1] == c and ini[i][0] != st:
                    stack.append(ini[i])
        elif c == 0:
            for i in range(len(ini)):
                if ini[i][0] == r and ini[i][1] != st:
                    stack.append(ini[i])
    #获取初始可行基的位置
    ini = np.argwhere(c_test >= 0).tolist()
    visit = []#存放已访问的初始可行基
    stack = []
    #初始化stack
    for i in range(len(ini)):
        if ini[i][0] == 0:
            stack.append(ini[i])
    #利用循环填写全部的ui，vj
    while len(visit) != (_c.shape[0] + _c.shape[1] - 1):
        if len(stack) > 0:
            r, c = stack.pop()
            c_pot[r, c] = 0
            if [r, c] not in visit:
                if np.isnan(c_pot[r, -1]) == True:
                    add_u(r, c, c_, c_pot)
                    add_stack(r, 0, c)
                else:
                    add_v(r, c, c_, c_pot)
                    add_stack(0, c, r)
                visit.append([r, c])
    str.clear()
    for i in range(c_.shape[0]):
        for j in range(c_.shape[1]):
            if [i, j] not in visit:
                c_pot[i, j] = c_[i,j]-c_pot[-1, j] - c_pot[i, -1]
                if c_pot[i,j]==0:
                    str.append([i,j])
#闭回路调整法
def Closed_circuit_adjustment():
    c_pot[-1, :] = np.nan* c_pot.shape[1]
    c_pot[:, -1] =np.nan* c_pot.shape[0]
    if bool(nbs)==True:
        r0,c0=str[0]
    else:
        adjust = np.argwhere(c_pot < 0)
        r0, c0 = adjust[0]
    c_pot[r0,-1]=0
    def add_stack(r, c, st):
        if r == -1:
            for i in range(len(ini)):
                if ini[i][1] == c and ini[i][0] != st:
                    stack.append(ini[i])
        elif c == -1:
            for i in range(len(ini)):
                if ini[i][0] == r and ini[i][1] != st:
                    stack.append(ini[i])
    #获取初始可行基的位置
    ini = np.argwhere(c_test >= 0).tolist()
    visit = []#存放已访问的初始可行基
    stack = []
    #初始化stack
    for i in range(len(ini)):
        if ini[i][0] == r0:
            stack.append(ini[i])
    #利用循环填写全部的ui，vj
    cmp = stack.copy()
    while 1:
        if len(stack) > 0:
            r, c = stack.pop()
            if len(visit) >= 1 and [r, c] in cmp:
                visit.clear()
            if [r, c] not in visit:
                if np.isnan(c_pot[r, -1]) == True:
                    add_stack(r, -1, c)
                    c_pot[r, -1] = 0
                else:
                    add_stack(-1, c, r)
                    c_pot[-1, c] = 0
                visit.append([r, c])
            if visit[-1][1]==c0:
                break
    visit=np.array(visit)
    way=[]
    way.append([r0,c0])
    way.append(visit[0])
    i=0
    k=1
    while 1:
        if k%2!=0:
            corner = np.argwhere(visit[:,1]==visit[i][1])
            i=corner[-1][0]
            way.append(visit[i])
        else:
            corner = np.argwhere(visit[:, 0] == visit[i][0])
            i = corner[-1][0]
            way.append(visit[i])
        k+=1
        if i==len(visit)-1:
            break
    visit=np.array(way)
    print("闭回路\n",visit)
    min_odd = []
    deg=[]
    min_c_odd=[]
    for i in visit[1::2]:
        r,c=i
        min_odd.append(c_test[r,c])
    min_ad=min(min_odd)
    r0,c0=visit[0]
    c_test[r0,c0]=0
    for i in range(len(visit)):
        r,c=visit[i]
        if i%2==0:
            c_test[r,c]+=min_ad
        else:
            c_test[r,c]-=min_ad
            if c_test[r,c]==0:
                c_test[r,c]=np.nan
                deg.append([r,c])
                min_c_odd.append(c_[r,c])
    if len(deg)>1:
        min_c=min(min_c_odd)
        for i in deg:
            r,c=i
            if c_[r,c]==min_c:
                c_test[r,c]=0
ifb_3()
print(c_test)
potential_method()
print(c_pot)
c_r=np.array(c)
def judge():
    c_pot1 = np.delete(c_pot, -1, axis=0)
    c_pot0 = np.delete(c_pot1, -1, axis=1)
    jud = np.argwhere(c_pot0 < 0)
    if len(jud) == 0:
        return 1
    else:
        return 0
result = 0
while bool(judge())==False:
    Closed_circuit_adjustment()
    print("运输方案:\n",c_test)
    c_pot = np.full((_c.shape[0] + 1, _c.shape[1] + 1), np.nan)
    potential_method()
    print("检验数:\n",c_pot)

basic_solution = np.argwhere(c_test >= 0)
for k in basic_solution:
    r, c = k
    result += c_test[r, c] * c_r[r, c]
print("最优解为：\n", result)
if len(str)!=0:
    print("无穷多最优解")
    nbs=1
    Closed_circuit_adjustment()
    print("运输方案2：\n",c_test)
    c_pot = np.full((_c.shape[0] + 1, _c.shape[1] + 1), np.nan)
    potential_method()
    print("检验数:\n", c_pot)