python求解线性规划问题———单纯形法（二）

两阶段法解一般的标准形式的线性规划问题

1、两阶段法是什么

单纯形法的关键在与如何找到初始基可行解，而两阶段法通过添加人工变量构造第一阶段的另一个线性规划问题，使得基可行解易求得，然后通过迭代得到第一阶段的最优解，再讨论原问题的最优解情况。

2、例题

原LP：（已经化为标准形式）
$$\begin{gathered} maxz=x_1+3x_2-x_3 \\ s.t.\{\begin{array}{l}{x}_1+x_2+2x_3+x_4=4\\ -x_1+2x_2+x_3+x_4=4\\ 3x_1+3x_3+x_4=4\\ x_i\ge 0i=1,2,3,4\end{array} \end{gathered}$$
LP1:
$$\begin{gathered} max{z}_1=-(x_5+x_6+x_7) \\ s.t.\{\begin{array}{l}{x}_1+x_2+2x_3+x_4+x_5=4\\ -x_1+2x_2+x_3+x_4+x_6=4\\ 3x_1+3x_3+x_4+x_7=4\\ x_i\ge 0i=1,2,3,4,5,6,7\end{array} \end{gathered}$$
易知LP1若有最优解，且最优解中$X_M=(x_5,x_6,x_7)=(0,0,0)$那么原问题就有基本可行解。

3、python实现

a、获取输入

def getinput():
	global m,n
	m = int(input('输入约束的个数'))# m个约束条件，n个变量
	n = int(input('输入变量的个数'))
    a = []   #约束系数矩阵
    for i in range(m):
        a.append(
            [eval(input('a({}{})='.format(i+1, j+1))) for j in range(n)]
            )
    e = [list(i) for i in np.diag([1]*m)]   #人工变量系数矩阵（m阶单位阵)
    b = [eval(input('b({})='.format(i+1))) for i in range(m)] #约束条件右端常数
    r_1 = []  #LP_1检验数
    for i in range(n):
        r_1.append(sum([a[j][i] for j in range(m)]))
    r_1 += [0]*m
    r_1.append(sum(b))
    r = [eval(input('r({})='.format(i+1))) for i in range(n)]+[0]*(m+1)
    vect = [n+i+1 for i in range(m)]  # 基变量足标
    a = [i+j+[k] for i, j, k in zip(a, e, b)]+[r_1]+[r]
    return a,vect

b、判断是否需要转轴

def judge(matrix):
    if max(matrix[-2][:-1]) <= 0:  # 先判断是否需要转轴（对于LP1）
        flag = False
    else:
        flag = True
    return flag

c、转轴

def trans():  #转轴
	global matrix,a,vect
    maxi = max(matrix[-2][:-1])
    index = a[-2].index(maxi)
    d = {}   #用来寻找满足min的变量
    for i in a[:-2]:
        if i[index] >0:
            d[i[-1]/i[index]] = a.index(i)
    pivot = d[min(d)]
    b[pivot] = b[pivot]/b[pivot][index]
    for i in range(m+2):
        if i != pivot:
            b[i] = b[i] - b[i][index]*b[pivot]
    vect[pivot] = index+1   #基变量的变动
    matrix = [list(i) for i in b]
    a = [list(i) for i in matrix]

d、格式化输出单纯性表

def pr(matrix,vect):   #输出单纯形表
    print('*'*20)   #单纯为了好看
    print('XB',end='\t')
    for i in range(n+m):
        print('X_{}'.format(i+1),end='\t')
    print('b')
    for i in range(m+2):
        if i <= m-1:
            if vect[i]>m:
                print('(x_{})'.format(vect[i]),end='\t')
            else:
                print('x_{}'.format(vect[i]),end='\t')
        elif i == m:
            print('r_1',end='\t')
        elif i == m+1:
            print('r',end='\t')
        for j in matrix[i]:
            print(f(str(j)).limit_denominator(),end='\t')   #输出分数形式
        print(end='\n')

e、输出最优解

def print_solution(matrix):
    print('*'*20)
    for i in range(n+m):
        print('x{}*={}'.format(i+1,matrix[-2][i]),end='，')
    print('\nz*={}'.format(-matrix[-2][-1]))

f、主函数

def main():
	global a,matrix,vect
    a,vect = getinput() 
    matrix = np.array(a, dtype=np.float64)   #array可以方便地进行整行的操作，而列表可以方便索引（其实就是我不会array的索引23333）
    pr(matrix,vect)
    while judge(matrix):
        trans()
        pr(matrix,vect)
    print_solution(matrix)

源代码

import numpy as np
from fractions import Fraction as f
def getinput():
    global m,n
    m = int(input('输入约束的个数'))# m个约束条件，n个变量
    n = int(input('输入变量的个数'))
    a = []   #约束系数矩阵
    for i in range(m):
        a.append(
            [eval(input('a({}{})='.format(i+1, j+1))) for j in range(n)]
            )
    e = [list(i) for i in np.diag([1]*m)]   #人工变量系数矩阵（m阶单位阵)
    b = [eval(input('b({})='.format(i+1))) for i in range(m)] #约束条件右端常数
    r_1 = []  #LP_1检验数
    for i in range(n):
        r_1.append(sum([a[j][i] for j in range(m)]))
    r_1 += [0]*m
    r_1.append(sum(b))
    r = [eval(input('r({})='.format(i+1))) for i in range(n)]+[0]*(m+1)
    vect = [n+i+1 for i in range(m)]  # 基变量足标
    a = [i+j+[k] for i, j, k in zip(a, e, b)]+[r_1]+[r]
    return a,vect
def judge(matrix):
    if max(matrix[-2][:-1]) <= 0:  # 先判断是否需要转轴（对于LP1）
        flag = False
    else:
        flag = True
    return flag
def trans():  #转轴
    global a,matrix,vect
    maxi = max(matrix[-2][:-1])
    index = a[-2][:-1].index(maxi)
    d = {}   #用来寻找满足min的变量
    for i in a[:-2]:
        if i[index] >0:
            d[i[-1]/i[index]] = a.index(i)
    pivot = d[min(d)]
    matrix[pivot] = matrix[pivot]/matrix[pivot][index]
    for i in range(m+2):
        if i != pivot:
            matrix[i] = matrix[i] - matrix[i][index]*matrix[pivot]
    vect[pivot] = index+1   #基变量的变动
    a = [list(i) for i in matrix]
def pr(matrix,vect):   #输出单纯形表
    print('*'*20)   #单纯为了好看
    print('XB',end='\t')
    for i in range(n+m):
        print('X_{}'.format(i+1),end='\t')
    print('b')
    for i in range(m+2):
        if i <= m-1:
            if vect[i]>m:
                print('(x_{})'.format(vect[i]),end='\t')
            else:
                print('x_{}'.format(vect[i]),end='\t')
        elif i == m:
            print('r_1',end='\t')
        elif i == m+1:
            print('r',end='\t')
        for j in matrix[i]:
            print(f(str(j)).limit_denominator(),end='\t')   #输出分数形式
        print(end='\n')
def print_solution(matrix):
    print('*'*20)
    for i in range(n+m):
        print('x{}*={}'.format(i+1,matrix[-2][i]),end='，')
    print('\nz*={}'.format(-matrix[-2][-1]))
def main():
    global a,matrix,vect
    a,vect = getinput() 
    matrix = np.array(a, dtype=np.float64)   #array可以方便地进行整行的操作，而列表可以方便索引（其实就是我不会array的索引23333）
    pr(matrix,vect)
    while judge(matrix):
        trans()
        pr(matrix,vect)
    print_solution(matrix)
main()

输入输出示例：

4、后续对原问题的讨论（略）