Simplex算法

文章分类:	Optimization
创建时间:	2013年4月

算法初始

单纯型算法是求解线性规划模型最有效的算法。根据单纯型算法，如果存在 $m$ 个约束条件，则最优解不为 0 的变量必然不多于 $m$ 个。现在求解 $m$ 约束条件的模型，假设目标函数存在可行解且最优解有限。可以选择 $m$ 个变量作为 基变量 用于求解，其他变量作为 非基变量 设定为 0。基变量有时也称基。

单纯型算法的基本步骤如下：

选择 $m$ 个基变量，剩余 $n-m$ 个变量为非基变量。令非基变量等于 0；
求取基变量的最优解，并用求解得到的基变量和非基变量计算目标函数值；
判断目标函数值是否最优。如不是最优解，则返回第 1 步。

单纯型算法基于标准模型，如果模型不是标准型，则用模型变换方法变换。

选择基变量后，变量和变量的系数就被分为两组：

基变量 $\vect{x}_B$ 和非基变量 $\vect{x}_N$ ，即 $\vect{x}=\begin{bmatrix}\vect{x}_B\\\vect{x}_N\end{bmatrix}$
基变量和非基变量在约束条件中对应的系数： $\vect{B}$ 和 $\vect{N}$ ，即 $\vect{A} = \begin{bmatrix} \vect{B} & \vect{N} \end{bmatrix}$
基变量和非基变量在目标函数中对应的系数： $\vect{c}_B$ 和 $\vect{c}_N$ ，即 $\vect{c}=\begin{bmatrix} \vect{c}_B & \vect{c}_N \end{bmatrix}$
约束条件中不等式或等式的右边的常数不变，即 $\vect{b} = \vect{b}$

$\begin{array}{lll} \vect{c}\vect{x} & = \begin{bmatrix} \vect{c}_B & \vect{c}_N \\ \end{bmatrix} \begin{bmatrix} \vect{x}_B \\ \vect{x}_N \\ \end{bmatrix} & = \vect{c}_B \vect{x}_B + \vect{c}_N \vect{x}_N \\[1.5em] \vect{A}\vect{x} & = \begin{bmatrix} \vect{B} & \vect{N} \\ \end{bmatrix} \begin{bmatrix} \vect{x}_B \\ \vect{x}_N \\ \end{bmatrix} & = \vect{B}\vect{x}_B + \vect{N}\vect{x}_N \end{array}$

算法过程

基本可行解

首先需要根据选定的基变量求解模型。因为 $\vect{x}_N=\vect{0}$ ，所以：

$\vect{A}\vect{x}=\vect{b}\quad\Rightarrow\quad \vect{B}\vect{x}_B=\vect{b}$

则基变量 $\vect{x}_B$ 的解：

$\vect{x}_B = \vect{B}^{-1}\vect{b} = \begin{bmatrix} d_1 \\ d_2 \\ \vdots \\ d_m \\ \end{bmatrix} \quad\Rightarrow\quad \vect{x}_d = \begin{bmatrix} \vect{B}^{-1}\vect{b} \\ \vect{0} \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ \vdots \\ d_m \\ 0 \\ \vdots \\ 0 \end{bmatrix}$

$\vect{x}_d$ 即线性优化模型的一个基本解。对于标准型，因为 $\vect{x}\ge 0$ ，所以必须满足 $\vect{x}_d\ge 0$ 。 $\vect{x}_d$ 是优化问题的一个 基本可行解，矩阵 $\vect{B}$ 为 可行基。

基本可行解与可行域多面体的顶点顶点对应，就像要判断多面体的顶点是否是最优解一样，需要判断基本可行解是否是最优解。如果不是，则需要调整。

最优解判定

用非基变量表示基变量：

$\vect{A}\vect{x} =\vect{B}\vect{x}_B + \vect{N}\vect{x}_N = \vect{b} \quad\Rightarrow\quad \vect{x}_B=\vect{B}^{-1}\vect{b}-\vect{B}^{-1}\vect{N}\vect{x}_N$

对应的目标函数用非基变量表示：

$z = \vect{c}\vect{x} = \vect{c}_B \vect{x}_B + \vect{c}_N \vect{x}_N = \vect{c}_B\vect{B}^{-1}\vect{b} + \left(\vect{c}_N-\vect{c}_B\vect{B}^{-1}\vect{N}\right)\vect{x}_N$

若非基变量 $\vect{x}_N=\vect{0}$ ，即基本可行解 $\vect{x}_d=\vect{B}^{-1}\vect{b}$ 时，目标函数值 $z=\vect{c}_B\vect{B}^{-1}\vect{b}$ 。如果这个基本可行解是最优解，那么该目标函数值就应该是最大可能值。因此即使某个或全部非基变量 $\vect{x}_N$ 不取 0，目标函数值也不会增加。因此对于 $\vect{x}\ge 0$ 必须有 $\vect{c}_N-\vect{c}_B\vect{B}^{-1}\vect{N}\le 0$ 。

令 $\sigma_i=c_i-\vect{c}_B\vect{B}^{-1}\vect{a}_i$ ：

$\vect{c}_N-\vect{c}_B\vect{B}^{-1}\vect{N} = \begin{bmatrix} \sigma_{\color{red}{1}} \\ \sigma_{\color{red}{2}} \\ \vdots \\ \sigma_{\color{red}{n-m}} \end{bmatrix}^T = \begin{bmatrix} c_{\color{red}{m+1}}-\vect{c}_B\vect{B}^{-1}\vect{a}_{\color{red}{m+1}} \\ c_{\color{red}{m+2}}-\vect{c}_B\vect{B}^{-1}\vect{a}_{\color{red}{m+2}} \\ \vdots \\ c_{\color{red}{n}}-\vect{c}_B\vect{B}^{-1}\vect{a}_{\color{red}{n}} \\ \end{bmatrix}^T$

则基本可行解是最优解的判定条件为：

所有 $\sigma_i\le 0$ ，且如所有 $\sigma_i<0$ ，则线性规划问题有唯一最优解。
如果所有 $\sigma_i\le 0$ 时，存在某个（几个） $\sigma_k=0$ ，则线性规划问题有多个最优解。
如果某个 $\sigma_k>0$ 且其对应的 $\vect{B}^{-1}\vect{a}_{m+k}\le \vect{0}$ ，则原优化问题有任意最优解。

基变量调整

如果基本可行解不是最优解，那么需要调整基变量，对应的 $\vect{B}$ 和 $\vect{N}$ 等都会发生变化。调整步骤如下：

从值为正的 $\sigma_i$ 中选择最大的；
在 $\vect{N}$ 中找到 $\sigma_k$ 对应的第 $k$ 列 $\vect{a}_{m+k}$ ，该列为 进基列；
计算 $\vect{B}^{-1}\vect{a}_{m+k}$ 和 $\vect{B}^{-1}\vect{b}$ ，求取两个向量的元素比 $\vect{\beta}$ ，即 $\beta_{j,m+k}=\frac{d_j}{\alpha_{j,m+k}},~j=1,2,\cdots,m$ 。在所有正的 $\alpha_{j,m+k}$ 中选择最小的 $\beta_{j,m+k}$ 。
在 $\vect{B}$ 中找到 $\beta_{l,m+k}$ 对应的第 $l$ 列 $\vect{a}_l$ ，该列为 出基列；
用进基列代替出基列，并对应地调整和：
- $\vect{x}_B$ 和 $\vect{c}_B$ 增加第 $m+k$ 个元素，去掉第 $l$ 个元素；
- $\vect{x}_N$ 和 $\vect{c}_N$ 增加第 $l$ 个元素，去掉第 $m+k$ 个元素。
利用调整后的基变量重新迭代。

计算程序

GLPK 是最好的优化软件。当然，我用 Go 语言实现的 Simplex 程序也不差。

$ glpsol -m mymodel.mod -o mymodel.sol

.mod 是模型设定文件，如下述模型：

$\begin{gather*} \text{maximize\quad } 3x_1+2x_2 \\[2mm] \begin{aligned} 2x_1 & + x_2 & \le 100 \\ x_1 & + x_2 & \le 80 \\ x_1 & & \le 40 \\ \end{aligned} \\ x_1\ge 0 \text{\ and\ } x_2\ge 0 \\ \end{gather*}$

模型设定代码为：

/* Decision variables */
var x1 >=0;  /* product 1 */
var x2 >=0;  /* product 2 */

/* Objective function */
maximize z: 3*x1 + 2*x2;

/* Constraints */
s.t. Labor    : 2*x1 + x2 <= 100;
s.t. Material : x1 + x2 <= 80;
s.t. Demand   : x1 <= 40;

end;