单纯形法python实现

'''
max Z=2*x_1+3*x_2
        x_1+2*x_2<=8
      4*x_1<=16
            4*x_2<=12
        x_1,x_2>=0
'''
#add slack variables and transform the problem into standard form
'''
max Z=2*x_1+3*x_2+0*x_3+0*x_4+0*x_5
        x_1+2*x_2+x_3               ==8
      4*x_1              +x_4       ==16
            4*x_2              +x_5 ==12
        x_1,x_2,x_3,x_4,x_5>=0
'''
import numpy as np
Basic=[2,3,4]
Nonbasic=[0,1]
c=np.array([2,3,0,0,0]).astype(float)
c_B=np.array([0,0,0]).astype(float)
c_N=np.array([2,3]).astype(float)
A=np.array([[1,2,1,0,0]
            ,[4,0,0,1,0]
            ,[0,4,0,0,1]]).astype(float)
A_N=np.array([[1,2]
              ,[4,0]
              ,[0,4]])
b=np.array([8,16,12]).astype(float)
B_inv=np.array([[1,0,0]
                ,[0,1,0]
                ,[0,0,1]]).astype(float)
x_opt=np.array([0,0,0,0,0]).astype(float)
z_opt=0
solutionStatus=None

row_num=len(A)
column_num=len(A[0])
reductedCost=c_N-np.dot(np.dot(c_B,B_inv),A_N)
print(reductedCost)
max_sigma=max(reductedCost)
#print(np.argmax(reductedCost))
eps=0.001
iterNum=1

while(max_sigma>=eps):
    #determine the entering basic  variable
    enter_var_index=Nonbasic[np.argmax(reductedCost)]
    print('enter_var_index:',enter_var_index)

    #determine the leaving basic variable :minmum ratio test
    min_ratio=100
    leave_var_index=enter_var_index
    for i in range(row_num):#3
        #print('b:',b[i],'\t A:',A[i][enter_var_index],'\t  ratio:',b[i]/A[i][enter_var_index])
        if(A[i][enter_var_index]==0):
            #solutionStatus
           continue
            #return solutionStatus
        elif(b[i]/A[i][enter_var_index]<min_ratio and b[i]/A[i][enter_var_index]>0):
            min_ratio=b[i]/A[i][enter_var_index]
            leave_var_index=i
    print('min_ratio',min_ratio)
    # process entering basis and leaving basis
    print('leave_var_index:',leave_var_index)
    leave_var=Basic[leave_var_index]
    Basic[leave_var_index]=enter_var_index
    Nonbasic.remove(enter_var_index)
    Nonbasic.append(leave_var)
    #Nonbasic.sort()
    #Gaussian elimination
    #update pivot row
    pivot_number=A[leave_var_index][enter_var_index]
    print('pivot_number',pivot_number)#4
    for col in range(column_num):
        A[leave_var_index][col]=A[leave_var_index][col]/pivot_number
    b[leave_var_index]=b[leave_var_index]/pivot_number
    #update other rows
    for row in range(row_num):
        if(row!=leave_var_index):
            factor=-A[row][enter_var_index]/1.0
            for col in range(column_num):
                A[row][col]=A[row][col]+factor*A[leave_var_index][col]
            b[row]=b[row]+factor*b[leave_var_index]
    # update c_N,c_B,A_N,B_inv
    for i in range(len(Nonbasic)):
        var_index=Nonbasic[i]
        c_N[i]=c[var_index]
    for i in range(len(Basic)):
        var_index=Basic[i]
        c_B[i]=c[var_index]
    for i  in range(row_num):
        for j  in  range(len(Nonbasic)):
            var_index=Nonbasic[j]
            A_N[i][j]=A[i][var_index]
    for i in range(len(Basic)):
        col=Basic[i]
        for row in range(row_num):
            B_inv[row][i]=A[row][col]

    #update reduced cost
    reductedCost=c_N-np.dot(np.dot(c_B,B_inv),A_N)
    max_sigma=max(reductedCost)
    iterNum +=1
#check the solution status
for i in range(len(reductedCost)):
    if(reductedCost[i]==0):
        solutionStatus='Alternative optimal solution'
        break
    else:
        solutionStatus='optimal'
#get the solution
x_basic=np.dot(B_inv,b)
x_opt=np.array([0.0]*column_num).astype(float)
for i in range(len(Basic)):
    basic_var_index=Basic[i]
    x_opt[basic_var_index]=x_basic[i]
z_opt=np.dot(np.dot(c_B,B_inv),b)
print('simplex iteration:',iterNum)
print('objective:',z_opt)
print('optimal solution:',x_opt)

来源于运筹优化常用模型（刘兴禄）