线性规划的标准型

对于线性规划，先来看一个简单的数学规划模型，即：

max Z=70x1+30x2

s.t. {3x1+9x2≤5405x1+5x2≤4509x1+3x2≤720x1,x2≥0

显然这不是线性规划数学模型的标准形式，在线性规划求解方法中，模型的标准形式如下：

1. 目标函数求最大值；
2. 约束条件为等式约束；
3. 约束条件右边的常数项大于等于 0；
4. 所有变量大于或等于 0。

对于非标准形式的模型：

• 约束方程可以通过引入 松弛变量 使不等式约束转化为等式约束；
• 如果目标函数是求最小值，则两边乘以 −1 将求 min 转成求 max；
• 如果遇到约束方程右边常数项为负数，则将约束方程乘以 −1 使常数项非负；
• 如果变量 xk 没有约束，则既可以是正数也可以是负数， 另 xk=x′k−x″k，其中 x′kx″k≥0。

通过变换，上面模型的标准型如下：

max Z=70x1+30x2

s.t.{3x1+9x2+x3=5405x1+5x2+x4=4509x1+3x2+x5=720x1,x2,x3,x4,x5≥0

将线性规划模型准换成标准型后，就可以使用经典的线性规划方法求解了，包括单纯形法、内点法、列生成法等。

线性规划的求解方式

单纯形法

单纯形法是求解线性规划的经典方法，与多元消去法求解多元一次方程的原理类似，在具体实现上， 通过矩阵的变换对解空间进行搜索。

由于线性规划模型的目标函数和约束方程都是凸函数，所以单纯形法能够以很高的效率求解线性规划问题。 单纯形法是最优算法的基础算法，也是后续其他整数规划等算法的基础。

单纯形法原理

由于线性规划模型中目标函数和约束方程都是凸函数，因此从凸优化的角度来说， 线性规划的最优解在可行域的顶点上， 单纯形法的本质就是通过矩阵的线性变换来遍历这些顶点以计算最优解。

假设有一个规划问题，即：

max Z=CX s.t. AX=b

将变量 X 拆解为基变量 XB 和非基变量 XN 两部分， 即 X=[XB,XN]，同理，将 C 拆解为 C=[CB,CN]， 将 A 拆解为 A=[AB,AN] 两部分，则规划问题可以表示为：

max Z=[CB,CN]T⋅[XB,XN] s.t. [AB,AN]T⋅[XB,XN]=b

假设 C=(c1,c2,⋯,cm,cm+1,⋯,cn)； X=(x1,x2,⋯,xm,xm+1,⋯,xn)， 前 m 个变量 x1,x2,⋯,xm 为基变量， 后面的 xm+1,⋯,xn 为非基变量。b=(b1,b2,⋯,bm)， 则约束方程可以写为如下形式：

s.t.{x1+a1,m+1xm+1+⋯+a1,nxn=b1x2+a2,m+1xm+1+⋯+a2,nxn=b2⋯⋯xm+am,m+1xm+1+⋯+am,nxn=bm

则基变量的值为：

xi=bi−(ai,m+1xm+1+⋯+ai,nxn)=bi−n∑j=m+1ai,jxj, i=1,⋯,m

将基变量 xi 代入目标函数，并消去目标函数中的基变量 XB，则：

max Z=n∑j=1cjxj=m∑i=1cixi+n∑j=m+1cjxj=m∑i=1ci[bi−n∑j=m+1aijxj]+n∑j=m+1cjxj=m∑i=1cibi−m∑i=1n∑j=m+1ciaijxj+n∑j=m+1cjxj=m∑i=1cibi+n∑j=m+1[cj−m∑i=1ciaij]xj=Z0+n∑j=m+1[cj−Zj]xj=Z0+n∑j=m+1Rjxj

式中，目标函数 Z0 的值为：

Z0=m∑i=1cibi=(c1,c2,⋯,cm)(b1,b2,⋯,bm)T=CBb

式中，Rj 为非基变量检验数，即：

Rj=cj−m∑i=1ciaij=cj−(c1,c2,⋯,cm)(a1j,a2j,⋯,amj)T=cj−CBaj, j=m+1,⋯,n

根据上面的公式推导，可以计算出目标函数的值、非基变量的检验数， 同时也说明了单纯形表变换的内在数学原理。

max Z=CBb+n∑j=m+1(cj−CBaj)xj

单纯形法过程

单纯形法就是遍历可行域各顶点后判断最优解

要用单纯形法求解线性规划数学模型，还需要把模型转化成规范形，规范形的条件如下：

1. 数学模型已经是标准型。

2. 约束方程组系数矩阵中含有至少一个单位子矩阵，对应的变量称为 基变量， 基的作用是得到 初始基本可行解，这个初始基本可行解通常是原点。 在大部分的问题中，通常引入松弛变量得到单位子矩阵，即使约束条件是等式约束， 也可以引入 xN=0 的松弛变量。

   上述示例的约束方程组系数矩阵是：

   A=[391005501093001]=(a1,a2,a3,a4,a5)

   对应的单位子矩阵是：

   B=[100000001]=(a3,a4,a5)

如何理解线性规划中的基变量和非基变量呢？

线性规划的最优解只能在顶点处取到，所以单纯形法的思想就是从一个顶点出发， 连续访问不同的顶点，在每一个顶点处检查是否有相邻的其他顶点能取到更优的目标函数值。

线性规划里面的约束（等式或不等式）可以看作是超平面（Hyperplane）或半空间（Half space）。 可行域可以看作是被这组约束，或者超平面和半空间定义（围起来）的区域。 那么某一个顶点其实就是某组超平面的交点，这一组超平面对应的约束就是在某一个顶点取到 "=" 号的约束（也就是基）。 顶点对应的代数意义就是一组方程（取到等号的约束）的解。

3. 目标函数中不含基变量。在这里基变量 x3、x4 和 x5 是在约束方程中引进的变量， 所以目标函数中没有这些基变量。

单纯形法的计算过程可以表示成单纯形表：

单纯形法的具体计算过程如下：

1. 确定初始基本可行解。 利用规范型的数学模型，整理出目标函数和约束方程的形式。

   max Z=70x1+30x2 s.t.{3x1+9x2+x3=5405x1+5x2+x4=4509x1+3x2+x5=720x1,x2,x3,x4,x5≥0

   令非基变量 xj=0,j=1,2，这样可以直接得到基变量的取值，即 x3=540， x4=450，x5=720，将非基变量 xj=0,j=1,2 代入目标函数得到 Z=0， 初始基本可行解是：

   X=(x1,x2,x3,x4,x5)=(0,0,540,450,720)T Z=0

   此时顶点位置是原点 O。

2. 判断当前点 X 是否为最优解。 对于最大化问题， 目标函数中非基变量的系数 ci≤0 时为最优解。 而这里非基变量的系数 c1=70>0，c2=30>0， 意味着只要在可行域内随着非基变量 x1 和 x2 的增大， 目标函数就会继续增大，所以此时的解不是最优解。

   当然也可以利用梯度的知识来思考这个问题。对于最大化问题， 只需要沿着梯度方向搜索即可找到最大值。线性规划是一个典型的凸优化问题， 当梯度为零时，得到最大值。即当非基变量系数 cj≤0 时，可得到问题最优解。

3. 基变量出基与非基变量入基。 变量的入基和出基在几何图上表现为顶点的变化， 如从 a 顶点变换到 b 顶点。

   • 选择使目标函数 Z 变化最快的非基变量入基， 即选择目标函数系数 ci 最大且为正数的非基变量入基，所以选择 x1 入基(非基变量变为基变量)， 此时仍然 x2=0。从凸优化的角度来看，就是选择目标函数梯度最大的方向做下一步的计算。
   • 那该选择哪个基变量出基呢？可以利用计算约束方程中常数项（bj）与 x1 系数（a1j）的比值 θ， 选择最小的 θ 对应的约束方程的基变量出基，即 θ=bjai，根据下面的单纯形表，因此选择 x5 出基(基变量变为非基变量)。

   目标和基变量	x1	x2	x3	x4	x5	b	θ
   Z	c1=70	c2=30
   x3	a11=3	a21=9	a31=1			b1=540	θ1=540/3=180
   x4	a12=5	a22=5		a42=1		b2=450	θ2=450/5=90
   x5	a13=9	a23=3			a53=1	b3=720	θ3=720/9=80

   因为选择 x1 入基，x5 出基，因此单纯形表和约束方程的第三个式子变成：

   目标和基变量	x1	x2	x3	x4	x5	b
   Z	c1=70	c2=30
   x3	a11=3	a21=9	a31=1			b1=540
   x4	a12=5	a22=5		a42=1		b2=450
   x5	a13=9/9=1	a23=3/9=1/3			a53=1/9	b3=720/9=80

   s.t.3,new x1+13x2+19x5=80⇒x1=80−13x2−19x5

   约束方程第一个式子不变，即：

   s.t.1,old 3x1+9x2+x3=540

   根据多元消去法，将 s.t.3,new 中的 x1 代入 s.t.1,old 中：

   s.t.1,new 3(80−13x2−19x5)+9x2+x3=540⇒0x1+8x2+x3−13x5=300

   即：

   s.t.1,new=s.t.1,old−3×s.t.3,new

   同理，得到其他约束方程和目标函数的替代，即新的目标函数和新的约束表达式形式。

   下面重新计算新的单纯形表：

   目标和基变量	x1	x2	x3	x4	x5	b	θ
   Z		c2=203			c5=−703	5600
   x3		a21=8	a31=1		a51=−13	b1=300	θ1
   x4		a22=103		a42=1	a52=−59	b2=50	θ2
   x5	a13=1	a23=13			a53=19	b3=80	θ3

   新的问题形式为：

   max Z=203x2−703x5 s.t.{−13x5+8x2+x3=300−59x5+103x2+x4=5019x5+13x2+x1=80x1,x2,x3,x4,x5≥0

4. 计算新的解 X。 令非基变量 xj=0,j=2,5， 求出基变量 xi=bi,i=1,3,4，得到基变量的值：x1=80,x3=300,x4=50，即：

   X=(x1,x2,x3,x4,x5)=(80,0,300,50,0)T Z=5600

   变换后，x1 的值从 0 变成 80，称为入基；x5 的值从 720 变成 0，称为出基，此时对应顶点为 a 点。

5. 判断当前解 X 是否最优。 由于目标函数中 x2 的系数仍然大于 0，因此当前位置还不是最优， 因为在可行域内随着 x2 增大目标函数还会增大。

6. 基变量出基与非基变量入基。 在目标函数中，系数为正且最大的变量是 x2， 因此选择 x2 入基，并计算 θ 选择出基变量。

   目标和基变量	x1	x2	x3	x4	x5	b	θ
   Z		c2=203			c5=−703	5600
   x3		a21=8	a31=1		a51=−13	b1=300	θ1=37.5
   x4		a22=103		a42=1	a52=−59	b2=50	θ2=15
   x1	a13=1	a23=13			a53=19	b3=80	θ3=240

   经过计算，选择 x4 出基，重新计算单纯形表：

   目标和基变量	x1	x2	x3	x4	x5	b	θ
   Z		c2=−2			c5=−203	5600
   x3			a31=1	a41=125	a51=−53	b1=180	θ1
   x2		a22=1		a42=310	a52=−16	b2=15	θ2
   x1	a13=1			a43=110	a53=16	b3=75	θ3

   新的问题形式为：

   max Z=−2x2−203x5 s.t.{125x4−53x5+x3=180310x4−16x5+x2=15110x4+16x5+x1=75x1,x2,x3,x4,x5≥0

7. 确定新的解 X。 令非基变量 xj=0,j=4,5，求出基变量 xi=bi,i=1,2,3， 得到基变量的值 x1=75,x2=15,x3=180，即：

   X=(x1,x2,x3,x4,x5)=(75,15,180,0,0)T Z=5700

   变换后，x2 的值从 0 变成 15，称为入基；x4 的值从 50 变成 0，称为出基，此时对应顶点为 h 点。

8. 判断当前解 X 是否最优。 因为函数中所有变量的系数均小于 0，变量 x4 和 x5 变大只会使目标函数减小， 所以当前解是最优解，即：

   X=(x1,x2,x3,x4,x5)=(75,15,180,0,0)T Z=5700

   经过单纯形法的搜索，搜索路径是 O→a→， 而不是 O→b→k→h，最终得到最优解。

单纯形法示例

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd


# 定义线性规划求解函数
def lp_solver(matrix: pd.DataFrame):
    """
    输入线性规划的矩阵，根据单纯形法求解线性规划模型
    max cx
    s.t. ax <= b

    Args:
        matrix (pd.DataFrame): 
            b      x1    x2    x3   x4   x5
        obj 0.0    70.0  30.0  0.0  0.0  0.0
        x3  540.0  3.0   9.0   1.0  0.0  0.0
        x4  450.0  5.0   5.0   0.0  1.0  0.0
        x5  720.0  9.0   3.0   0.0  0.0  1.0

        - 第 1 行是目标函数的系数
        - 第 2~4 行是约束方程的系数
        - 第 1 列是约束方程的常数项
        - obj-b 交叉，即第 1 行第 1 列的元素是目标函数的负值
        - x3,x4,x5 既是松弛变量，也是初始可行解
    """
    # 检验数是否大于 0
    c = matrix.iloc[0, 1:]
    # TODO
    while c.max() > 0:
        # ------------------------------
        # 选择入基变量 
        # ------------------------------
        # 目标函数系数最大的变量入基
        c = matrix.iloc[0, 1:]
        print(c)
        in_x = c.idxmax()
        print(in_x)
        # in_x_v = c[in_x]  # 入基变量的系数
        # print(in_x_v)
        print("-" *40)
        # ------------------------------
        # 选择出基变量
        # ------------------------------
        # 选择正的最小比值对应的变量出基 min(b列/入基变量列)
        b = matrix.iloc[1:, 0]
        print(b)
        in_x_a = matrix.iloc[1:][in_x]  # 选择入基变量对应的列
        print(in_x_a)
        out_x = (b / in_x_a).idxmin()  # 得到出基变量
        print(out_x)
        # out_x_v = b[out_x]  # 出基变量的系数
        # print(out_x_v)
        print("-" * 40)
        # ------------------------------
        # 旋转操作
        # ------------------------------
        matrix.loc[out_x, :] = matrix.loc[out_x, :] / matrix.loc[out_x, in_x]
        # print(matrix)
        # print("-" * 40)
        for idx in matrix.index: 
            if idx != out_x:
                matrix.loc[idx, :] = matrix.loc[idx, :] - matrix.loc[out_x, :] * matrix.loc[idx, in_x]
        # print(matrix)
        # print("-" * 40)
        # 索引替换（入基与出基变量名称替换）
        matrix_index = matrix.index.tolist()
        i = matrix_index.index(out_x)
        print(matrix_index)
        print(in_x, out_x)
        print(i)
        matrix_index[i] = in_x
        print(matrix_index)
        print("-" * 40)
        matrix.index = matrix_index 
    # 打印结果
    print("最终的最优单纯形法是：")
    print(matrix)
    print(f"目标函数值是：{-matrix.iloc[0, 0]}")
    print("最优决策变量是：")
    x_count = (matrix.shape[1] - 1) - (matrix.shape[0] - 1)
    X = matrix.iloc[0, 1:].index.tolist()[:x_count]
    for xi in X:
        print(f"{xi} = {matrix.loc[xi, 'b']}")


# 测试代码 main 函数
def main():
    # 约束方程系数矩阵，包含常数项
    matrix = pd.DataFrame(
        np.array([
            [0.0, 70.0, 30.0, 0.0, 0.0, 0.0],
            [540.0, 3.0, 9.0, 1.0, 0.0, 0.0],
            [450.0, 5.0, 5.0, 0.0, 1.0, 0.0],
            [720.0, 9.0, 3.0, 0.0, 0.0, 1.0],
        ]),
        index = ["obj", "x3", "x4", "x5"],
        columns = ["b", "x1", "x2", "x3", "x4", "x5"]
    )
    print(matrix)
    print("-" * 40)
    # 调用前面定义的函数求解
    lp_solver(matrix)

if __name__ == "__main__":
    main()

最终的最优单纯形法是：
          b   x1   x2   x3   x4        x5
obj -5700.0  0.0  0.0  0.0 -2.0 -6.666667
x3    180.0  0.0  0.0  1.0 -2.4  1.000000
x2     15.0  0.0  1.0  0.0  0.3 -0.166667
x1     75.0  1.0  0.0  0.0 -0.1  0.166667
目标函数值是：5700.0
最优决策变量是：
x1 = 75.0
x2 = 15.0

内点法

内点法也是求解线性规划的一个方法，相比单纯形法， 内点法在 大规模线性优化、 二次优化、非线性规划 方面都有比较好的表现。 内点法是 多项式算法， 随着问题规模的增大，计算的复杂度却不会急剧增大， 因此在大规模问题上比单纯形法有更广泛的应用。

内点法原理

内点法的求解思路同 拉格朗日松弛法 的思路类似，将约束问题转化为无约束问题， 通过无约束函数的梯度下降进行迭代直至得到有效解。

内点法就是在梯度下降的过程中，如果当前迭代点是在可行域外，则会给损失函数一个非常大的值， 这样就能约束在可行域内求解。但是，内点法不能处理等式约束， 因为构造的内点惩罚函数是定义在可行域内的函数，而等式约束优化问题不存在可行域空间。 由此看来，内点法和单纯形法对优化问题的形式是不一样的。

下面介绍一下内点法是如何将约束问题转化为无约束问题的。

考虑一个最小化的线性规划问题：

minZ=cTx s.t.Ax≤b

借鉴拉格朗日松弛法的思路，这个线性规划问题可以表示成如下函数：

minf(x)=cTx+m∑i=1I(Ax−b)

其中，m 是约束方程的个数，I 是指示函数，一般定义如下：

I(u)={0,if u≤0∞,if u>0

通过指示函数可以将约束方程直接写到目标函数中，然后对目标函数求极小值。 但是这个指示函数 I(u) 是不可导的，需要用其他可导的函数近似替代， 常用的替代函数如下：

I−(u)=−1tlog(−u)

当 u>0 时，I−(u)=∞，参数 t 决定 I−(u) 对应 I(u) 的近似程度， 类似机器学习中损失函数的正则化参数的作用。因此，新的目标函数可以写成如下形式：

minf(x)=cTx−1tm∑i=1log(−Ax+b)=tcTx−m∑i=1log(−Ax+b)

由于指示函数 I−(u) 是凸函数，所以新的目标函数也是凸函数， 因此可以用凸优化中的方法求解该函数的极小值，例如梯度下降法、牛顿法、拟牛顿法、L-BFGS 等。 下面以经典的牛顿法讲解如何求函数最小化问题。

目标函数 f(x) 在 x0 做二阶泰勒公式展开时得到：

f(x)≈f(x0)+(x−x0)f′(x0)+12(x−x0)2f″(x0)

上式成立的条件是 f(x) 近似等于 f(x0)，对上面等式两边同时对 (x−x0) 求导， 并令导数为 0，可以得到下面的方程：

f′(x)=f′(x0)+(x−x0)f″(x0) x=x0−f′(x0)f″(x0)

这样就得到了下一点的位置，从 x0 走到 x1。 重复这个过程，直到到达导数为 0 的点，因此牛顿法的迭代公式是：

xn+1=xn−f′(xn)f″(xn)=xn−H−1∇f(xn)

式中，H−1 表示二阶导数矩阵（黑塞矩阵）的逆。因此，如果使用牛顿法求解目标函数最优值， 需要知道目标函数的一阶导数 f′(xn) 和二阶导数 f″(xn)。

内点法过程

再看前面的线性规划问题：

max Z=70x1+30x2 s.t.{3x1+9x2≤5405x1+5x2≤4509x1+3x2≤720x1,x2≥0

问题转换为内点法需要的形式：

min Z=−70x1−30x2 s.t.{3x1+9x2−540≤05x1+5x2−450≤09x1+3x2−720≤0−x1≤0−x2≤0

问题转换为无约束优化问题形式：

minZ=t(−70x1−30x2)=log(−3x1−9x2+540)−=log(−5x1−5x2+450)−=log(−9x1−3x2+720)−=log(x1)−=log(x2)

在问题中，目标函数的一阶导数是：

J=[∂f∂x1∂f∂x2]

∂f∂x1=−70t+3−3x1−9x2+540+5−5x1−5x2+450+9−9x1−3x2+720+1−x1 ∂f∂x2=−30t+9−3x1−9x2+540+5−5x1−5x2+450+3−9x1−3x2+720+1−x2

目标函数的二阶导数是：

H=[∂2f∂x1∂x1∂2f∂x1∂x2∂2f∂x2∂x1∂2f∂x2∂x2]

∂2f∂x1∂x1=9(3x1+x2−240)2+1(x1+3x2−180)2+1(x1+x2−90)2+1x21 ∂2f∂x1∂x2=3(3x1+x2−240)2+3(x1+3x2−180)2+1(x1+x2−90)2 ∂2f∂x2∂x1=3(3x1+x2−240)2+3(x1+3x2−180)2+1(x1+x2−90)2 ∂2f∂x2∂x2=1(3x1+x2−240)2+9(x1+3x2−180)2+1(x1+x2−90)2+1x22

选择一个恰当的初始解 x0 代入牛顿法迭代公式：

xn+1=xn−H−1∇f

不断迭代直至得到最优解。

总结：内点法和单纯形法的结果相差很大，只是因为内点法的搜索路径是在可行域内部， 而不在可行域的边界上，这也是内点法的局限性。并且通过前面的求解过程发现， 内点法不仅局限在线性规划上，二次规划等也是可以求解的，因为其本质是利用函数梯度求最优值， 这同机器学习算法的思路是一致的，真正的难点在于如何保证新的目标函数是否存在一阶导数和二阶导数， 以及如何得到一阶导数和二阶导数的信息。此外，初始迭代点的选择也是很重要的， 在线性规划问题中能够保证最后得到的是最优解，而非线性规划问题中，函数是非凸的， 因此很难保证最后的解是全局最优解。

内点法示例

import time
import numpy as np

def gradient(x1, x2, t):
    """
    计算目标函数在 x 处的一阶导数（雅可比矩阵 ）

    Args:
        x1 (_type_): _description_
        x2 (_type_): _description_
        t (_type_): _description_
    """
    j1 = -70*t + 3/(-3*x1 - 9*x2 + 540) \
               + 5/(-5*x1 - 5*x2 + 450) \
               + 9/(-9*x1 - 3*x2 + 720) \
               - 1/x1
    j2 = -30*t + 9/(-3*x1 - 9*x2 + 540) \
               + 5/(-5*x1 - 5*x2 + 450) \
               + 3/(-9*x1 - 3*x2 + 720) \
               - 1/x2

    return np.asmatrix([j1, j2]).T


def hessian(x1, x2):
    """
    计算目标函数在 x 处的二阶导数（黑塞矩阵）

    Args:
        x1 (_type_): _description_
        x2 (_type_): _description_
    """
    x1, x2 = float(x1), float(x2)
    h11 = 9/(3*x1 + x2 - 240)**2 + (x1 + 3*x2 - 180)**(-2) + (x1 + x2 - 90)**(-2) + x1**(-2)
    h12 = 3/(3*x1 + x2 - 240)**2 + 3/(x1 + 3*x2 - 180)**2 + (x1 + x2 - 90)**(-2)
    h21 = 3/(3*x1 + x2 - 240)**2 + 3/(x1 + 3*x2 - 180)**2 + (x1 + x2 - 90)**(-2)
    h22 = (3*x1 + x2 - 240)**(-2) + 9/(x1 + 3*x2 - 180)**2 + (x1 + x2 - 90)**(-2) + x2**(-2)
    
    return np.asmatrix([[h11, h12], [h21, h22]])


def invertible(H):
    """
    求黑塞矩阵的逆矩阵

    Args:
        H (_type_): _description_
    """
    H_inv = np.linalg.inv(H)

    return H_inv


def run():
    # 牛顿法的初始迭代值
    x = np.asmatrix(np.array([10, 10])).T
    # 指数函数中的 t
    t = 0.00001
    # 迭代停止的误差
    eps = 0.01
    # 记录迭代的次数
    iter_cnt = 0
    while iter_cnt < 20:
        iter_cnt += 1
        J = gradient(x[0, 0], x[1, 0], t)
        H = hessian(x[0, 0], x[1, 0])
        H_inv = invertible(H)
        # 牛顿法
        x_new = x - H_inv * J
        # 求二范数，判断迭代效果
        error = np.linalg.norm(x_new - x)
        print(f"迭代次数是：{iter_cnt}, x1={x_new[0, 0]:.2f}, x2={x_new[1, 0]:.2f}, 误差是: {error}")
        x = x_new
        if error < eps:
            break
        time.sleep(1)
    # 打印结果
    print(f"目标函数值是：{70*x[0, 0] + 30*x[1, 0]:.2f}")


# 测试代码 main 函数
def main():
    run()

if __name__ == "__main__":
    main()

迭代次数是：1, x1=15.91, x2=15.34, 误差是: 7.964345953591451
迭代次数是：2, x1=20.13, x2=18.19, 误差是: 5.09056325893654
迭代次数是：3, x1=21.00, x2=18.33, 误差是: 0.8795427059613212
迭代次数是：4, x1=21.02, x2=18.32, 误差是: 0.02676438103520062
迭代次数是：5, x1=21.02, x2=18.32, 误差是: 2.199053220249811e-05
目标函数值是：2021.17

列生成法

列生成法是一种用于求解大规模线性优化问题非常高效的算法。 本质上，列生成算法就是单纯形法的一种形式，它是用来求解线性规划问题的， 所不同的列生成法改善了大规模优化问题中单纯形法基变换计算效率低的问题， 列生成法在整数规划中已经得到了广泛应用。

列生成法原理

列生成法主要用于解决变量很多而约束相对较少的问题，特别是经常用于解决大规模整数规划问题。

单纯形法虽然能保证在数次迭代后找到最优解，但是其面对变量很多的线性规划问题就显得很弱了。 因为它需要在众多变量里进行基变换，所以这种遍历的计算量是很大的。

有人基于单纯形法提出了列生成法， 其思路是强制 原问题(Master Problem) 把一部分变量限定(Restrict)为基变量， 得到一个规模更小（即变量数比原问题少的）的 限制性主问题(Restrict Master Problem)， 在限制主问题上用单纯形法求解最优解，但是此时求得的最优解只是限制性主问题的解， 并不是原问题的最优解，就需要通过一个 子问题(Subproblem) 去检查在那些未被考虑的变量中， 是否有使限制主问题的 ReducedCost 小于 0，如果有，就把这个变量的相关系数列加入到限制主问题的系数矩阵中。

列生成法的形象化表示过程如下，考虑线性规划问题，即：

min c1x1+c2x2+⋯+cnxn s.t.{a11x1+a12x2+⋯+a1nxn=b1⋯am1x1+am2x2+⋯+amnxn=bm

把前面 j 个变量强制设定为基变量，后面 n−j 个变量设定为非基变量，由于非基变量等于 0， 用矩阵表示的线性规划形式是：

minc1x1+c2x2+⋯+cjxj s.t.[a11⋯a1j⋯⋯⋯am1⋯amj]x=b

此时的问题为限定主问题，通过求解限定主问题先得到对偶问题的最优解， 然后用子问题检查是否存在新的变量使目标变量继续朝优化方向变化（检验数小于 0）， 假设存在一个满足的变量 xj+1，那么原问题则变成：

minc1x1+c2x2+⋯+cjxj+cj+1xj+1 s.t.[a11⋯a1ja1j+1⋯⋯⋯⋯am1⋯amjamj+1]x=b

此时系数矩阵多了一列，这就是列生成法名称的由来。此时求解新限定主问题得到对偶问题的最优解， 通过子问题得到新的变量，如此循环往复直到无法添加新的变量。

通过这个过程可以看到，原问题的变量非常多，通过限定主问题和列生成后， 最终的线性规划模型的变量相比原问题少得多，这也是列生成法能够求解大规模线性规划的原因之一。

接下来讲解限定主问题和子问题的关系。

假设有如下线性规划问题：

min cTx s.t.{Ax=bx≥0

与单纯形的数学规范形的思路类似，令 x=[xB,xN]，其中 xB 表示基变量， xN 表示非基变量。类似地，A=[AB,AN]，cT=[cTB,cTN]， 约束方程变成：

Ax=b⇔BxB+NxN=b⇔xB=B−1b−B−1NxN

因为非基变量 xN=0，所以 xB=B−1b。

对目标函数有如下推导关系：

min cTx⇔min cTBxB+cTNxN⇔min cTB(B−1b−B−1NxN)+cTNxN⇔min (cTN−cTBB−1N)xN+cTBB−1b

其中：

• cTBB−1b 是常数
• xN 非负
• cTN−cTBB−1N 称为 ReducedCost， ReducedCost 即非基变量的检验数，如果 ReducedCost 小于 0， 则说明有非基变量 xN 使得目标函数可以更优。

在检验数 ReducedCost 的计算公式 cTN−cTBB−1N 中， cTBB−1 有以下两重含义：

1. 通过求解 RMP 问题 得到的 影子价格(Shadow Price)；
2. 通过求解 RMP 对偶问题 得到的 对偶变量(Dual Variable)。

因此，可以通过限制性主问题及其对偶问题来确定非基变量。

下面是 ReducedCost 和原问题的对偶问题的关系（非对称形式对偶），即：

Primal：

min cTx s.t.{Ax=bx≥0

⇔

Dual：

max λTb s.t. λTA≤cT

假设找到了原问题的最优解是 [xB,0]，那么 cTn−cTBB−1an≥0， 即 cTBB−1an≤cTn，cTBB−1A=cTBB−1[B,N]， 由于 B−1B=1，所以 cTBB−1A=cTBB−1[B,N]≤[cTB,cTN]=cT， 因此，cTBB−1A≤cT，cTBB−1 是对偶问题的可行解。