14、单纯形法_两阶段法（目标最大，约束等于，包含对偶问题求解）

最新推荐文章于 2023-09-17 20:16:31 发布
chaoyuzhang
最新推荐文章于 2023-09-17 20:16:31 发布
阅读量1.1k
点赞数
分类专栏：运筹优化、机器学习基础算法实现
本文链接：https://blog.csdn.net/chaoyuzhang/article/details/88405244
版权
运筹优化、机器学习基础算法实现专栏收录该内容
42 篇文章 9 订阅
订阅专栏
本文介绍了一种基于单纯形法的线性规划求解算法，通过迭代计算寻找最优解。算法首先进行第一阶段迭代，检查是否存在可行解，然后进行第二阶段迭代，找到目标函数的最大值或最小值。文中提供了多个例子并展示了算法的具体实现。
摘要由CSDN通过智能技术生成
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 11 17:18:21 2019

@author: zhangchaoyu
"""
import copy
import math

def show_1(C, D, A, B, X, vc, vd):
    
    result = []
    result.append(copy.deepcopy(C+vc))
    result.append(copy.deepcopy(D+vd))
    for i in range(len(A)):
        result.append(copy.deepcopy(A[i]+[B[i]]))
    for r in result:
        print(r)
    print("当前解：")
    print("X="+str(X))
    print()

def find_new_1(D, A, B):
    idex_negative = -1
    for i in range(len(D)):
        if D[i] < 0 and (idex_negative == -1 or (idex_negative != -1 and D[i] < D[idex_negative])):
            idex_negative = i
            
    if idex_negative == -1:
        print("第一阶段检验数无小于0项，已至最优")
        return -1
    
    flag = -1
    row = -1
    rate = -1
    for i in range(len(B)):
        if A[i][idex_negative] <= 0:
            continue
        if rate == -1 or (rate != -1 and rate > B[i]/A[i][idex_negative]):
            flag = 1
            row = i
            rate = B[i]/A[i][idex_negative]
    if flag == 1:
        return [row, idex_negative]
    else:
        print("第一阶段检验数无小于0项，已至最优")
        return -1

def adjust_1(idex_row, idex_col, A, B, C, D, idexs, X, vc, vd):
    
    m = len(B)
    n = len(D)
    
    idexs[idex_row] = idex_col
    
    rate = A[idex_row][idex_col]
    for i in range(n):
        A[idex_row][i] /= rate
    B[idex_row] /= rate
    
    for i in range(m):
        if i == idex_row:
            continue
        rate = A[i][idex_col]/A[idex_row][idex_col]
        B[i] -= rate*B[idex_row]
        for j in range(n):
            A[i][j] -= rate*A[idex_row][j]
    
    rate = C[idex_col]/A[idex_row][idex_col]
    for i in range(n):
        C[i] -= rate*A[idex_row][i]
    vc[0] -= rate*B[idex_row]
    
    rate = D[idex_col]/A[idex_row][idex_col]
    for i in range(n):
        D[i] -= rate*A[idex_row][i]
    vd[0] -= rate*B[idex_row]
    
    for i in range(len(X)):
        if i not in idexs:
            X[i] = 0
        else:
            for j in range(len(idexs)):
                if i == idexs[j]:
                    X[i] = B[j]

def find_new_2(C, A, B):
    idex_positive = -1
    for i in range(len(C)):
        if C[i] > 0 and (idex_positive == -1 or (idex_positive != -1 and C[i] > C[idex_positive])):
            idex_positive = i
    
    if idex_positive == -1:
        print("第二阶段检验数无大于0项，已至最优")
        return -1
    
    flag = -1
    row = -1
    rate = -1
    flag_n = -1
    for i in range(len(B)):
        if A[i][idex_positive] < 0:
            flag_n = 1
        if A[i][idex_positive] <= 0:
            continue
        if rate == -1 or (rate != -1 and rate > B[i]/A[i][idex_positive]):
            flag = 1
            row = i
            rate = B[i]/A[i][idex_positive]
    if flag == 1:
        return [row, idex_positive]
    else:
        if flag_n == 1:
            print("第二阶段解无界")
        else:
            print("第二阶段检验数无大于0项，已至最优")
        return -1

def adjust_2(idex_row, idex_col, A, B, C, idexs, X, v):
    
    m = len(B)
    n = len(C)
    
    idexs[idex_row] = idex_col
    
    rate = A[idex_row][idex_col]
    for i in range(n):
        A[idex_row][i] /= rate
    B[idex_row] /= rate
    
    for i in range(m):
        if i == idex_row:
            continue
        rate = A[i][idex_col]/A[idex_row][idex_col]
        B[i] -= rate*B[idex_row]
        for j in range(n):
            A[i][j] -= rate*A[idex_row][j]
    
    rate = C[idex_col]/A[idex_row][idex_col]
    for i in range(n):
        C[i] -= rate*A[idex_row][i]
    v[0] -= rate*B[idex_row]
    
    for i in range(len(X)):
        if i not in idexs:
            X[i] = 0
        else:
            for j in range(len(idexs)):
                if i == idexs[j]:
                    X[i] = B[j]

def show_2(C, A, B, X, v):
    
    result = []
    result.append(copy.deepcopy(C+v))
    for i in range(len(A)):
        result.append(copy.deepcopy(A[i]+[B[i]]))
    for r in result:
        print(r)
    print("当前解：")
    print("X="+str(X))
    print()

def sub_data_integer(x):
    
    if math.ceil(100*x)/100 - x <= 1e-10:
        return math.ceil(100*x)/100
    
    if x - math.floor(100*x)/100 <= 1e-10:
        return math.floor(100*x)/100
    
    return x

def data_integer(A, B, C, D, X, vc, vd):
    
    for i in range(len(A)):
        for j in range(len(A[i])):
            A[i][j] = sub_data_integer(A[i][j])
    
    for i in range(len(B)):
        B[i] = sub_data_integer(B[i])
        
    for i in range(len(C)):
        C[i] = sub_data_integer(C[i])
    
    for i in range(len(D)):
        D[i] = sub_data_integer(D[i])
    
    for i in range(len(X)):
        X[i] = sub_data_integer(X[i])
        
    vc[0] = sub_data_integer(vc[0])
    vd[0] = sub_data_integer(vd[0])

def simplex(C0, A0, B0):
    
    C = copy.deepcopy(C0)
    A = copy.deepcopy(A0)
    B = copy.deepcopy(B0)

    for i in range(len(B)):
        if B[i] >= 0:
            continue
        B[i] = -B[i]
        for j in range(len(A[i])):
            A[i][j] = -A[i][j]

    C = C + [0 for i in range(len(B0))]
    for i in range(len(A)):
        temp = [0 for j in range(len(A))]
        temp[i] = 1
        A[i] = A[i] + copy.deepcopy(temp)
    D = [sum([-A[j][i] for j in range(len(B))]) for i in range(len(C0))] + [0 for i in range(len(B0))]
    vc = [0]
    vd = [sum([-B[i] for i in range(len(B))])]
    idexs = [len(C0)+i for i in range(len(A))]
    X = [0 for i in range(len(C0))] + copy.deepcopy(B)
    
    #开始第一阶段迭代计算
    loop = 0
    print("第一阶段第"+str(loop)+"次迭代后：")
    show_1(C, D, A, B, X, vc, vd)
    idex = find_new_1(D, A, B)
    
    while True:
        if idex == -1:
            break
        idex_row = idex[0]
        idex_col = idex[1]
        adjust_1(idex_row, idex_col, A, B, C, D, idexs, X, vc, vd)
        data_integer(A, B, C, D, X, vc, vd)
        loop += 1
        print("第一阶段第"+str(loop)+"次迭代后：")
        show_1(C, D, A, B, X, vc, vd)
        idex = find_new_1(D, A, B)
    
    #统计第一阶段结果，判断是否有可行解
    #planA
    for i in range(len(A)):
        if len(C0)+i in idexs:
            print("无可行解，当前近似解X="+str(X[:len(C0)]))
            return X[:len(C0)]
    """
    #planB
    if vd[0] != 0:
        print("无可行解，当前近似解X="+str(X[:len(C0)]))
        return X[:len(C0)]
    """
    #有可行解后，重新建立单纯形表
    C = C[:len(C0)]
    for i in range(len(A)):
        A[i] = A[i][:len(C0)]
    X = X[:len(C0)]
    v = vc
    
    #开始第二阶段迭代计算
    loop = 0
    print("第二阶段第"+str(loop)+"次迭代后：")
    show_2(C, A, B, X, v)
    idex = find_new_2(C, A, B)
    
    while True:
        if idex == -1:
            break
        idex_row = idex[0]
        idex_col = idex[1]
        adjust_2(idex_row, idex_col, A, B, C, idexs, X, v)
        data_integer(A, B, C, D, X, vc, vd)
        loop += 1
        print("第二阶段第"+str(loop)+"次迭代后：")
        show_2(C, A, B, X, v)
        idex = find_new_2(C,A,B)
    
    print("当前解：")
    print("X="+str(X))
    print("最优值：")
    print(sum([C0[i]*X[i] for i in range(len(X))]))
    return X 

def dual_simplex(C0, A0, B0):
    
    C = []
    A = []
    B = []
    
    for i in range(len(B0)):
        C.append(-B0[i])
        C.append(B0[i])
    for i in range(len(C0)):
        C.append(0)
    
    B = copy.deepcopy(C0)
    
    for i in range(len(A0[0])):
        temp = []
        for j in range(len(A0)):
            temp.append(A0[j][i])
            temp.append(-A0[j][i])
        for j in range(len(A0[0])):
            if i == j:
                temp.append(-1)
            else:
                temp.append(0)
        A.append(copy.deepcopy(temp))
    
    X = simplex(C, A, B)
    
    return X
    
 
examples = []
examples.append([[2,3,0,0,0],[[1,2,1,0,0],[4,0,0,1,0],[0,4,0,0,1]],[8,16,12]])
examples.append([[-1000,-800,0,0,0,0],[[1,0,-1,0,0,0],[0.8,1,0,-1,0,0],[1,0,0,0,1,0],[0,1,0,0,0,1]],[1,1.6,2,1.4]])
examples.append([[1,1,0,0],[[-2,1,1,0],[1,-1,0,1]],[4,2]])
examples.append([[0,-0.1,-0.2,-0.3,-0.8],
                 [[1,2,0,1,0],
                  [0,0,2,2,1],
                  [3,1,2,0,3]
                  ],
                 [100,100,100]])
examples.append([[-15,25,15,-30,10,0,-40,0,-10,0,0,0,0,0,0,0],
                 [[-0.5,0.5,0.5,0,0,0,0,0,0,1,0,0,0,0,0,0],
                  [-0.25,0.75,-0.25,0,0,0,0,0,0,0,1,0,0,0,0,0],
                  [0,0,0,-0.75,0.25,0.25,0,0,0,0,0,1,0,0,0,0],
                  [0,0,0,-0.5,0.5,-0.5,0,0,0,0,0,0,1,0,0,0],
                  [1,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0],
                  [0,1,0,0,1,0,0,1,0,0,0,0,0,0,1,0],
                  [0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,1]
                  ],
                 [0,0,0,0,100,100,60]])
examples.append([[-150,-150,-150,-80,-80,-80]+[0,0,0,0,0,0,0,0,0,0,0,0,0],
                 [[1,0,0,0,0,0]+[1,0,0,0,0,0,0,0,0,0,0,0,0],
                  [2,1,0,0,0,0]+[0,1,0,0,0,0,0,0,0,0,0,0,0],
                  [3,2,1,0,0,0]+[0,0,1,0,0,0,0,0,0,0,0,0,0],
                  [4,3,2,1,0,0]+[0,0,0,1,0,0,0,0,0,0,0,0,0],
                  [5,4,3,2,1,0]+[0,0,0,0,1,0,0,0,0,0,0,0,0],
                  [6,5,4,3,2,1]+[0,0,0,0,0,1,0,0,0,0,0,0,0],
                  [7,6,5,4,3,2]+[0,0,0,0,0,0,1,0,0,0,0,0,0],
                  [8,7,6,5,4,3]+[0,0,0,0,0,0,0,1,0,0,0,0,0],
                  [8,8,7,5,5,4]+[0,0,0,0,0,0,0,0,1,0,0,0,0],
                  [8,8,8,5,5,5]+[0,0,0,0,0,0,0,0,0,-1,0,0,0],
                  [4,3,2,1,0,0]+[0,0,0,0,0,0,0,0,0,0,-1,0,0],
                  [7,6,5,4,3,2]+[0,0,0,0,0,0,0,0,0,0,0,-1,0],
                  [1,1,1,1,1,1]+[0,0,0,0,0,0,0,0,0,0,0,0,1]
                  ],
                 [10,18,24,32,37,43,51,60,67,72,24,43,11]])
examples.append([[1,3,0,0,0],
                 [[5,10,1,0,0],
                  [1,1,0,-1,0],
                  [0,1,0,0,1]
                  ],
                 [50,1,4]])
    
examples.append([[-1,-1.5,0,0],
                 [[1,3,-1,0],
                  [1,1,0,-1]
                  ],
                 [3,2]])
examples.append([[0,0,0,1.15,1.25,1.4,0,0,0,0,1.06,0,0],
                 [[1,0,0,0,0,0,1,0,0,0,0,0,0],
                  [0,1,0,0,0,1,-1.06,1,0,0,0,0,0],
                  [-1.15,0,1,0,1,0,0,-1.06,1,0,0,0,0],
                  [0,-1.15,0,1,0,0,0,0,-1.06,1,0,0,0],
                  [0,0,-1.15,0,0,0,0,0,0,-1.06,1,0,0],
                  [0,0,0,0,1,0,0,0,0,0,0,1,0],
                  [0,0,0,0,0,1,0,0,0,0,0,0,1]
                  ],
                 [10,0,0,0,0,4,3]])

examples.append([[1,1,0,0],[[1,1,-1,0],[1,1,0,1]],[3,2]])

examples.append([[-1,-1,-1,-1,-1,-1,0,0,0,0,0,0],
                 [[1,0,0,0,0,1]+[-1,0,0,0,0,0],
                  [1,1,0,0,0,0]+[0,-1,0,0,0,0],
                  [0,1,1,0,0,0]+[0,0,-1,0,0,0],
                  [0,0,1,1,0,0]+[0,0,0,-1,0,0],
                  [0,0,0,1,1,0]+[0,0,0,0,-1,0],
                  [0,0,0,0,1,1]+[0,0,0,0,0,-1]
                  ],
                 [60,70,60,50,20,30]])
examples.append([[0.95,0.45,-0.05,1.4,0.95,0.45,1.9,1.45,0.95]+[0,0,0,0,0,0,0,0],
                 [[-0.4,0,0,0.6,0,0,0.6,0,0]+[1,0,0,0,0,0,0,0],
                  [0,-0.85,0,0.15,0,0,0.15,0,0]+[0,1,0,0,0,0,0,0],
                  [-0.2,0,0,-0.2,0,0,0.8,0,0]+[0,0,1,0,0,0,0,0],
                  [0,-0.6,0,0,-0.6,0,0,0.4,0]+[0,0,0,1,0,0,0,0],
                  [0,0,-0.5,0,0,-0.5,0,0,0.5]+[0,0,0,0,1,0,0,0],
                  [1,1,1,0,0,0,0,0,0]+[0,0,0,0,0,1,0,0],
                  [0,0,0,1,1,1,0,0,0]+[0,0,0,0,0,0,1,0],
                  [0,0,0,0,0,0,1,1,1]+[0,0,0,0,0,0,0,1]
                  ],
                 [0,0,0,0,0,2000,2500,1200]]) 

examples.append([[0.75,0.7753,-0.375,-0.44743,-0.35,-0.5,-0.2889,1.15,0.684371]+[0,0,0,0,0],
                 [[1,1,-1,-1,-1,0,0,0,0]+[0,0,0,0,0],
                  [0,0,0,0,0,1,1,-1,0]+[0,0,0,0,0],
                  [1,0,0,0,0,2,0,0,0]+[1,0,0,0,0],
                  [0,7,0,0,0,0,9,0,12]+[0,1,0,0,0],
                  [0,0,6,0,0,0,0,8,0]+[0,0,1,0,0],
                  [0,0,0,4,0,0,0,0,11]+[0,0,0,1,0],
                  [0,0,0,0,7,0,0,0,0]+[0,0,0,0,1]
                  ],
                 [0,0,1200,10000,4000,7000,4000]])
    
test = []
test.append(14)
test.append(-1640)
test.append("无界解")
test.append(-16)
test.append(500)
test.append(-1350)
test.append(14)
test.append(-2.25)
test.append(14.375)
test.append("无可行解")
test.append(-150)
test.append(6447)
test.append(1147)


C0 = examples[-1][0]
A0 = examples[-1][1]
B0 = examples[-1][2]
X = simplex(C0, A0, B0)
X_dual = dual_simplex(C0, A0, B0)