【运筹优化】单纯形法

1、处理问题

单纯形法，主要用于解决线性规划问题。

线性规划问题：指目标函数和约束条件皆为线性的最优化问题。

2、预备知识

1：线性规划的约束条件的所有点组成的集合即为可行域，可行域是一个凸集，且顶点有限。

2：如果线性规划的可行域有界，则可行域的顶点对应可行解，则最优解必然在可行域的顶点上。

3：线性规划问题局部最优一定是全局最优，因为线性规划是一个凸优化问题，可行域是一个凸集。

4：线性规划的标准形式特点：约束条件均为等于约束；决策变量取值为非负；右端常数非负。

5：线性规划标准型表达式：

6：检验数公式：

7：初始可行解定理：构建线性规划标准型、找到基本可行解。

8：最优性检验：针对min问题，所有检验数≥0时最优；针对max问题，所有验数≤0的时最优。

9：入基、出基规则：针对min问题，选择检验数<0的变量入基，最小的入基；针对max问题，选择检验数>0的变量入基，最小的入基。

3、处理流程

• 有选择地取基本可行解，而不是枚举所有的解。

1. 根据初始可行解定理，从可行域的一个顶点出发。
2. 根据最优性检验，判断当前解是不是最优的。如果是，则停止，如果不是，则3。
3. 根据入基、出基规则，沿着可行域的边界移到另一个相邻的顶点，回到2。

4、举个例子

5、代码实现

# Simplex Method Algorithm

'''
max 2 * x_1 + 3 * x_2
s.t.
        x_1 + 2 * x_2 <= 8
    4 * x_1           <= 16
              4 * x_2 <= 12
        x_1,      x_2 >= 0


standard form
max 2 * x_1 + 3 * x_2
s.t.
        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 pandas as pd
import numpy as np

# Construct the initial feasible solution and calculate:C、C_B、C_N；A、B、N；b；basic、nobasic；sigma、max_sigma；

nobasic = [0,1]
basic = [2,3,4]

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)
B = np.array([[1,0,0],
              [0,1,0],
              [0,0,1]]).astype(float)
N = np.array([[1,2],
              [4,0],
              [0,4]]).astype(float)

b = np.array([8,16,12]).astype(float)

x_opt = np.array([0,0,0,0,0]).astype(float)
z_opt = 0

row_num = len(A)
column_num = len(A[0])

sigma = C_N - np.dot(np.dot(C_B,B),N)
print('sigma:',sigma)
max_sigma = max(sigma)

eps = 0.001
iterNum = 1


while max_sigma >= eps:

    #calculate entering basic variable and outgoing variable
    enter_var_index = nobasic[np.argmax(sigma)]
    print('enter_var_index:',enter_var_index)

    min_ration = 1000
    leave_var_index = 0
    for i in range(row_num):
        print('b:',b[i],'\t A:',A[i][enter_var_index],'\t ration:',b[i]/A[i][enter_var_index])
        if A[i][enter_var_index] == 0:
            continue
        elif b[i]/A[i][enter_var_index] > 0 and b[i]/A[i][enter_var_index] < min_ration:
            min_ration = b[i]/A[i][enter_var_index] 
            leave_var_index = i
    
    leave_var = basic[leave_var_index]
    basic[leave_var_index] = enter_var_index
    nobasic.remove(enter_var_index)
    nobasic.append(leave_var)
    nobasic.sort()

    # Gaussian elimination iteration
    pivot_number = A[leave_var_index][enter_var_index]
    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
    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]
    for i in range(len(nobasic)):
        var_index = nobasic[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(nobasic)):
            var_index = nobasic[j]
            N[i][j] = A[i][var_index]
    for i in range(len(basic)):
        col = basic[i]
        for row in range(row_num):
            B[row][i] = A[row][col]
    
    # update sigma
    sigma = C_N - np.dot(np.dot(C_B,B),N)
    max_sigma = max(sigma)
    iterNum += 1

# Calculate the optimal solution
for i in range(len(sigma)):
    if sigma[i] == 0:
        solution_status = 'Alternative optimal solution'
        break
    else:
        solution_status = 'Optimal'

x_basic = np.dot(B,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),b)

print('Simplex iteration:',iterNum)
print('objective:',z_opt)
print('optimal solution:',x_opt)

#输出：
#b: 16.0          A: 0.0          ration: inf
#b: 12.0          A: 4.0          ration: 3.0
#enter_var_index: 0
#b: 2.0   A: 1.0          ration: 2.0
#b: 16.0          A: 4.0          ration: 4.0
#b: 3.0   A: 0.0          ration: inf
#enter_var_index: 4
#b: 2.0   A: -0.5         ration: -4.0
#b: 8.0   A: 2.0          ration: 4.0
#b: 3.0   A: 0.25         ration: 12.0
#Simplex iteration: 4
#objective: 14.0
#optimal solution: [4. 2. 0. 0. 4.]