【最优化】线性规划的MathBeauty

凸集 | Convex Set

定义

设$S\subset\Bbb R^n$ 是$n$ 维欧氏空间中的一个点集，若连接$S$ 中任意两个点$\boldsymbol x_1,\boldsymbol x_2\in S$ 的线段仍属于$S$；换言之，对每个实数$\lambda\in[0,1]$ 都有：

$$\lambda\boldsymbol x_1+(1-\lambda)\boldsymbol x_2\in S$$

则称集合$S$ 为凸集，称$\lambda\boldsymbol x_1+(1-\lambda)\boldsymbol x_2$ 是$\boldsymbol x_1,\boldsymbol x_2$ 的凸组合。

定义

集合 $C$  的凸包是  $C$ 中所有点的所有凸组合的集合，记作 $\text{conv }C$.

Note: 当$\lambda\in\Bbb R$ 时，集合为仿射集合（Affine sets），对应的线性组合也叫仿射组合，另外还有仿射包。

e.g. 上图中，$S_1$ 为凸集；$S_2$ 为非凸集.

🦄凸优化问题开源课程推荐：https://www.stat.cmu.edu/~ryantibs/convexopt/

常见的凸集

超平面 | Hyperplane

我们知道，以$xyz$ 坐标系为例，形如$ax+by+cz=d$ 的式子在三维空间中表示一个平面。更确切来说，应该是这样一个点集：$\{(x,y,z)|ax+by+cz=d\}$ 形成了一个三维平面。

将这种概念拓展到多维空间中时，我们可以定义满足$\boldsymbol p^T\boldsymbol x=\alpha$ 的$n$ 维点$\boldsymbol x\in\Bbb R^n$ 组成的集合为超平面：$H=\{\boldsymbol x|\boldsymbol p^T\boldsymbol x=\alpha\}$.

可以证明超平面是一个凸集.

半空间 & Others

不难想象，一个超平面可以把一个空间一分为二，而超平面往下的部分或者另一边的部分就是一个半空间（Half space）。用数学形式表示：$H^-=\{\boldsymbol x|\boldsymbol p^T\boldsymbol x\leq\alpha\}$.

可以证明半空间是一个凸集.

---

【多面体】有限个半空间的 交 我们称为多面体(Polyhedron)，其约束可以写为：

$$\begin{aligned} \begin{pmatrix}\boldsymbol p_1^T\\\boldsymbol p_2^T\\\vdots\\\boldsymbol p_n^T\end{pmatrix}\boldsymbol x \leq\begin{pmatrix}b_1\\b_2\\\vdots\\b_n\end{pmatrix}\Leftrightarrow\boldsymbol {Ax\leq b} \end{aligned}$$

多面体也是凸集。

【射线】以某一定点$\boldsymbol x^{(0)}$ 为起点，沿着方向$\boldsymbol d$ 延长出去的射线：$L=\{\boldsymbol x|\boldsymbol x=\boldsymbol x^{(0)}+\lambda\boldsymbol d,\;\lambda\geq0\}$ 是凸集。要求方向向量$\boldsymbol d$ 为非零向量。

【凸锥】若集合$C\subset\Bbb R^n$ 中每一个点$\boldsymbol x$，对任意非负$\lambda$ 都有$\lambda\boldsymbol x\in C$，则称$C$ 为锥（Cone）。当$C$ 为凸集时，称其为凸锥（convex cone）。

在多面集的表达式中，若$\boldsymbol {b=0}$ ，则它也是凸锥，称其为多面锥。

单纯形 | Simplex

单纯形是一种特殊的多面体。它之所以叫“单纯形”，是因为它是该维度上“最简单”的多面体。几何上，单纯形是三角形或四面体在任意维度的推广。

定义

设 $k+1$ 个点 $\boldsymbol x_0,...,\boldsymbol x_k$ 仿射独立，即 $\boldsymbol x_1−\boldsymbol x_0,  ...  ,\boldsymbol x_k−\boldsymbol x_0$ 线性无关。则这些点的凸组合构成了一个单纯形：

$$C=\text{conv}⁡\{\boldsymbol x_0,\cdots,\boldsymbol x_k\}=\{\lambda_0\boldsymbol x_0+\cdots+\lambda_k\boldsymbol x_k\mid\lambda\geq0,\sum_{i=0}^k\lambda_i=1\}$$

这个单纯形的仿射维度为 $k$ ，因此被称为 $\Bbb R^n$ 空间的$k$ 维单纯形 (k-simplex)。

仿射维度$k$	几何形态
0	点
1	线段
2	三角形
3	四面体
4	5-cell
…	…

单纯形与单纯形方法的关系

单纯形方法中并没有用到单纯形。那为什么单纯形方法会如此命名呢？
一个解释是该方法是在一些单纯锥上进行的，加上一些额外的约束，这些单纯锥就成了单纯形。

欧氏球 | Euclidean balls

$B(\boldsymbol x_{c},r)=\{\boldsymbol x \mid ||\boldsymbol x-\boldsymbol x_{c} ||_{2}{\;\leq\;} r \}=\{\boldsymbol x_{c}+ru \mid || u ||_{2}{\;\leq\;}1\}$，其中$\boldsymbol x_c$ 为球心，$r$ 为半径。
因为使用的是$L2$ 范数定义，所以就称之为欧几里得球体。欧氏球是凸集。

事实上这可以扩展到任意范数，都能证明其为凸集，推导原理是三角不等式。类似地，还可以得出各种范数里的椭球也是凸集。

常见性质

设$S_1,S_2\subset\Bbb R^n$ 为两个凸集，$\beta$ 为一实数，则：

1. $\beta S_1=\{\beta\boldsymbol x|\boldsymbol x\in S_1\}$ 为凸集；
2. $S_1\cap S_2$ 为凸集；
3. $S_1+S_2=\{\boldsymbol x^{(1)}+\boldsymbol x^{(2)}|\boldsymbol x^{(1)}\in S_1,\boldsymbol x^{(2)}\in S_2\}$ 为凸集；
4. $S_1-S_2$ 为凸集.

事实上，以上限制进一步推广就是保凸运算。此处不过多展开，可参考前面提到的开源课程或者下面的博客。

凸优化 - 凸集 - Infinity-SEU
凸集 - 水论文的程序猿

极点与极方向

定义

设 $S\subset\Bbb R^n$ ，则$x\in S$  称为 $S$ 的极点，当且仅当 $\not\exists x_1,x_2\in S\;(x_1≠x_2),\;\lambda∈(0,1)$  使得 $x=\lambda x_1+(1−\lambda)x_2$ 。

即极点不能用凸集中两个不同的点的严格凸组合表示。（注意$\lambda\neq0,1$）

Note: 为了书写方便，此处(甚至以后)的向量和矩阵在不会引起歧义的情况下都不再通过加粗的方式板书。

【Krein-Milman定理】：设 $S⊆\Bbb R^n$ 是紧凸集，则 $S=\text{conv(ext}(S))$. 即紧凸集可由其极点的所有凸组合表示出来(即极点的凸包)。其中 $\text{ext}(S)$ 表示$S$ 的所有极点构成的集合。

---

以图(a)为例，$x_6$ 在$x_5,x_7$ 构成的线段上，换言之它能被$x_5,x_7$ 的凸组合表示，所以它不是极点。类似地，$x_8$ 也不是极点，而其余点都不能被其他点的凸组合表示，所以它们都是极点。而图(b)中，圆周上所有的点都是极点。

图(a)中的凸集$S$ 有$\text{ext}(S)=\{x_1,x_2,x_3,x_4,x_5,x_7\}$，并且$S=\text{conv(ext}(S))$。

---

Krein-Milman定理只适用于紧凸集，但对于无界集来说就不成立了。为此我们需要引出极方向的概念。

定义

设 $S⊆\Bbb R^n$ 是闭凸集，向量$d\in\Bbb R^n,d\neq0$ ，若$\forall x\in S$ 都有：

$$\{ x+\lambda d\mid\lambda\geq0\}\subset S$$

则称$d$ 为$S$ 的方向。也就是说，无论$x$ 在$S$ 中的哪个位置，从$x$ 出发沿着$d$ 的方向延长出去，这条射线上的点也都被凸集$S$ 所包含，那么$d$ 就是$S$ 的一个方向。

若$d$ 不能由$S$ 中其他两个不同方向的线性组合表示，则$d$ 为$S$ 的一个极方向。

以下图为例，这是一个二维空间的多面体，其中$0$ 点是它的一个极点，$d^{(1)},d^{(2)}$ 是它唯二的两个极方向。可以发现区域中的其他方向都能由$d^{(1)},d^{(2)}$ 线性表出。

更近一步对多面体进行展开，我们可以得到其表示定理：
设$S=\{x\mid Ax=b,x\geq0\}$ 为非空多面体，则有：

1. $\text{ext}(S)\neq\varnothing$，极点集非空，且存在有限个极点；
2. 极方向集合为空的充要条件是$S$ 有界，否则存在有限个极方向；
3. $\boldsymbol x\in S$ 的充要条件是：

$$\begin{aligned} \boldsymbol x&=\sum_{j=1}^k\lambda_j\boldsymbol x^{(j)}+\sum_{j=1}^l\mu_j\boldsymbol d^{(j)}\\ &\sum_{j=1}^k\lambda_j=1\\ &\lambda_j\geq0,\quad j=1,...,k\\ &\mu_j\geq0,\quad j=1,...,l \end{aligned}$$

Note: 应当注意闭集与紧集的区别。下面补充部分集合的基础知识：

凸函数 | Convex Function

定义

设 $S⊆\Bbb R^n$ 是非空凸集，$f$ 是定义在$S$ 上的实函数，若对$S$ 中任意两个点$\boldsymbol x_1,\boldsymbol x_2\in S$ 及$\forall\lambda\in(0,1)$ 都有：

$$f(\lambda\boldsymbol x_1+(1-\lambda)\boldsymbol x_2)\leq\lambda f(\boldsymbol x_1)+(1-\lambda)f(\boldsymbol x_2)$$

则称$f$ 是凸集$S$ 上的凸函数，$-f$ 为偶函数。不等号取$\lt$ 时为严格凸函数。

注意：在凸优化领域，甚至是数学领域提到凸函数都是上图这种形式，系从坐标数值从低往高的方向看去时的形状是凸的。与初高中乃至高等数学中学的从上往下看的定义正好相反！

重要性质

1. 有限个凸函数的仿射组合也是凸函数；
2. 凸函数$f:S\to\Bbb R$ 的水平集$H_c(f)=\{x|x\in S,\;f(x)\leq c\}$ 也是凸集；
3. 凸函数$f$ 在$x_0\in S$ 处沿着方向$d$ 的方向导数$\text{D}f(x_0;d)$ 存在；
4. 凸函数$f$ 在凸集$S$ 上的局部极小点就是全局极小点（反证法可证明），极小点集合为凸集（根据性质2）

凸函数判别

【一阶条件】当$f(\boldsymbol x)$ 在$S$ 上可微时，$f$ 为凸函数的充要条件为：

$$f(\boldsymbol x^{(2)})\geq f(\boldsymbol x^{(1)})+\nabla f(\boldsymbol x^{(1)})^\top(\boldsymbol x^{(2)}-\boldsymbol x^{(1)})$$

几何意义：可微函数为凸函数的充要条件是在其定义域凸集中任一点处的切平面 (切线)都不在曲面(曲线)的上方。也就是说任意点处的切线增量不超过函数的增量。

【二阶条件】当$f(\boldsymbol x)$ 是定义在$S$ 上的二次可微函数时，$f$ 为凸函数的充要条件是对$\forall x\in S$，其 Hesse矩阵$\nabla^2f$ 是半正定矩阵。（是正定矩阵时则为严格凸函数）

凸优化的定义

设$f(\boldsymbol x),g_i(\boldsymbol x)$ 为$\Bbb R^n$ 上的凸函数，$h_j(\boldsymbol x)$ 为$\Bbb R^n$ 上的线性函数，则称最优化问题：

$$\begin{aligned} \min \quad&f(\boldsymbol x)\\ \text{s. t. }\quad&g_i(\boldsymbol x)\leq0,\quad i=1,...,m\\ &h_j(\boldsymbol x)=0,\quad j=1,...,l \end{aligned}$$

为凸规划。

也就是说，求凸函数在凸集上的极小点问题就是凸规划。
这是因为，$g_i(\boldsymbol x)\leq0$ 根据凸函数性质（水平集）可知满足此约束的点集是凸集，而$h_j(\boldsymbol x)$ 作为线性函数即是凸函数也是凹函数，满足$h_j(\boldsymbol x)=0$ 的点集自然构成凸集，而凸集之间的交集依然是凸集。

如果目标函数$f(\boldsymbol x)$ 是严格凸函数，且存在极小点，则其极小点唯一。

线性规划 | Linear  Programming

一般线性规划问题的标准形式如下：

$$\begin{aligned} \min\quad&\sum_{j=1}^nc_jx_j\\ \text{s. t. }\quad&\sum_{j=1}^na_{ij}x_j=b_i,\quad i=1,...,m\\ &x_j\geq0,\quad j=1,...,n \end{aligned}$$

我们更常用其矩阵形式：

$$\begin{aligned} \min\quad&\boldsymbol{cx}\\ \text{s. t. }\quad&\boldsymbol{Ax=b}\\ &\boldsymbol{x\geq0} \end{aligned}$$

为了计算方便，我们一般假定$\boldsymbol b\geq0$. 否则可以对方程两端添负号以统一格式。

形式转换

最大化

目标函数需要最大化时，对目标函数取反化为求最小值问题。

非负性

对于没有限制$\boldsymbol x\geq0$ 的问题，可用变量替换法。
如$x_j$ 没有非负限制，可以引入$x_j',x_j''\geq0$，令$x_j=x_j'-x_j''$。

上下界

当存在$\boldsymbol x\in[a,b]$ 的限制时，可用变量替换法。
令$x_j'=x_j-a$ 使得$x_j'\geq0$ 且$x_j'\leq b-a$，令$x_j''=(b-a)-x_j'$ 使得$x_j''\geq0$。

绝对值

在目标函数中存在$|x_j|$ 时，可用变量替换法。
令$x_j'=(x_j+|x_j|)/2,x_j''=(|x_j|-x_j)/2$。
如此一来，表示变量本身时用：$x_j=x_j'-x_j''$，表示绝对值时用：$|x_j|=x_j'+x_j''$。

不等式

当约束条件$Ax=b$ 不是严格等式时，引入 松弛变量（Slack Variable）或 剩余变量（Surplus Variable）：

$$\begin{aligned} a_{11}x_1+\cdots+a_{1n}x_n\leq b_i\quad\to\quad a_{11}x_1+\cdots+a_{1n}x_n+{\color{blue}{x_{n+1}}}= b_i \end{aligned}$$

图解法

方法略

最优极点

显然 LP 的标准形式中，约束项是一种多面体，利用表示定理可以得到其目标函数如下：

$$\sum_{j=1}^k(\boldsymbol c\boldsymbol x^{(j)})\lambda_j+\sum_{j=1}^l(\boldsymbol c\boldsymbol d^{(j)})\mu_j$$

并且对$\mu$ 的限制只有非负，因此当某一项$\boldsymbol c\boldsymbol d^{(j)}\lt0$ 时，随着$\mu_j$ 的增大，最后目标函数趋于负无穷。从几何角度来说，就是目标函数沿着没有边界的方向移动时函数值可以不断减小，这导致最终不存在有限最优值。我们也称这样的问题无界。

在此后的讨论中，我们默认这种无界情况就属于“不存在最优解”。

另外，若令$\boldsymbol c\boldsymbol x^{(p)}=\min\limits_{1\leq j\leq k}\boldsymbol c\boldsymbol x^{(j)}$ ，则有：

$$\begin{aligned} \boldsymbol c\boldsymbol x&=\sum_{j=1}^k(\boldsymbol c\boldsymbol x^{(j)})\lambda_j+\sum_{j=1}^l(\boldsymbol c\boldsymbol d^{(j)})\mu_j\\ &\geq\sum_{j=1}^k(\boldsymbol c\boldsymbol x^{(j)})\lambda_j\qquad(\text{$\mu_j可以都取0$})\\ &\geq\sum_{j=1}^k(\boldsymbol c\boldsymbol x^{(p)})\lambda_j\qquad(\text{$x^{(j)}可以都取x^{(p)}$})\\ &=\boldsymbol c\boldsymbol x^{(p)} \end{aligned}$$

从而说明目标函数的最小值可以在某个极点上取到。

---

以上两方面的推导告诉我们，如果$\boldsymbol c\boldsymbol d^{(j)}\lt 0$ 则无界；如果LP存在最优解，那么一定能在某个极点上取得。

🔔最优基本可行解

在 LP 模型中，我们将系数矩阵$A\in\Bbb R^{m\times n}$ 按列分块：$A=(p_1,p_2,...,p_n),\;p_i\in\Bbb R^m$。
从而有：

$$Ax=b\Rightarrow p_1x_1+p_2x_2+\cdots+p_nx_n=b$$

又假设$r(A)=m$ ，那么$A$ 中必定存在$m$ 个列向量$p_i$ ，它们共同组成$m$ 阶满秩方阵，记为$B$ 。那么经过有限次列变换，有$A=[B,N]$。同理，对应有$x=[x_B,x_N]$.
即：

$$\begin{aligned} &\begin{pmatrix}B&N\end{pmatrix}\begin{pmatrix}x_B\\x_N\end{pmatrix}=b\\ \Leftrightarrow\quad& Bx_B+Nx_N=b\\ \Rightarrow\quad& x_B=B^{-1}b-B^{-1}Nx_N \end{aligned}$$

而$x_N$ 就是线性代数中的自由未知量，它们取不同的值就能得到方程组不同的解。

特别地，我们都取为零，得到的$x$ 定义为$Ax=b$ 的基本解：

$$x=\begin{pmatrix}x_B\\x_N\end{pmatrix}=\begin{pmatrix}B^{-1}b\\0\end{pmatrix}$$

其中，$B$ 称为基矩阵，简称为基，$x_B$ 的各分量称为基变量。

如果基本解满足非负约束，即$B^{-1}b\geq0$ ，则称这样的基本解为基本可行解。
若$B^{-1}b$ 中至少有一个分量为$0$，则称这个基本解是退化的；
若$B^{-1}b\gt0$ 则称为非退化的。

---

不难发现，基本解的个数与基矩阵的个数一致，而基矩阵的个数等价于每次从$n$ 个列向量中选出$m$ 个。因此基本解的个数为：

$$C_n^m=\frac{n!}{m!(n-m)!}$$

此外还可以证明，可行域$K=\{x\mid Ax=b,\;x\geq0\}$ 的极点集$\text{ext}(K)$ 与基本可行解等价。而在上一节我们指出如果LP存在最优解，那么一定能在某个极点上取得。

于是，LP问题的求解就相当于求最优基本可行解。

单纯形法 | Simplex Method

上一章节我们已经论述了LP问题的求解等价于求最优基本可行解。我们能想到的一个最简单的求解方法就是将所有可能的基本可行解一一枚举，然后哪一个解使得目标函数值最小（标准形式的目标函数是最小化的），那么那个解就是最优解。

显然这个方法在变量过大时，计算是低效的（并且还要反复计算逆矩阵）。于是，单纯形法就提出了这样的思想：
先找出可行域的一个顶点（即一个基本可行解），然后根据一定规则判断其是否最优；否则对基本可行解进行一些改动，使得目标函数值比之前更优；如此下去，直到找到最优解为止。

最优解判别准则

为了推理更加清晰，我们下面都使用标准LP的矩阵形式，并且沿用上一章对基本可行解的等价形式转换，如下图所示。

假设我们目前得到了一个基本可行解：

$$\boldsymbol x^{(0)}=\begin{bmatrix}B^{-1}\bf b\\\bf0\end{bmatrix}$$

此时的目标函数值：

$$f_0=\boldsymbol{cx}^{(0)}=\boldsymbol c_BB^{-1}\boldsymbol b$$

又设任意一个基本可行解：

$$\boldsymbol x=\begin{bmatrix}\boldsymbol x_B\\\boldsymbol x_N\end{bmatrix}$$

由$A\boldsymbol x=\boldsymbol b$ 可得：

$$\boldsymbol x_B=B^{-1}\boldsymbol b-B^{-1}N\boldsymbol x_N$$

从而在$\boldsymbol x$ 处的目标函数值为：

$$\begin{aligned} f&=\boldsymbol {cx}=\boldsymbol c_B\boldsymbol x_B+\boldsymbol c_N\boldsymbol x_N\\ &=\boldsymbol c_B(B^{-1}\boldsymbol b-B^{-1}N\boldsymbol x_N)+\boldsymbol c_N\boldsymbol x_N\\ &=\boldsymbol c_BB^{-1}\boldsymbol b -(\boldsymbol c_BB^{-1}N-\boldsymbol c_N)\boldsymbol x_N\\ &=f_0-(\boldsymbol c_BB^{-1}N-\boldsymbol c_N)\boldsymbol x_N\quad//\because f_0=\boldsymbol c_BB^{-1}\boldsymbol b\\ &=f_0-\sum_{j\in R}(\boldsymbol c_BB^{-1}\boldsymbol p_j-c_j)x_j\quad//R\text{ 为非基变量的下标集合}\\ &=f_0-\sum_{j\in R}(z_j-c_j)x_j \end{aligned}$$

其中，定义$\boldsymbol c_BB^{-1}\boldsymbol p_j=z_j$.
显然，如果能在非基变量中（也就是可选的自由未知量）$x_j$ 中选取合适的值，使得

$$\sum_{j\in R}(z_j-c_j)x_j\gt0$$

那么，目标函数值就能得到改善。

由于非基变量也需要满足$\boldsymbol x_N\geq\boldsymbol 0$，所以要让上面的不等式成立（也就是让目标函数值得到改善），就必须要求所有的$z_j-c_j\gt0$。若不然，就说明目标函数值无法再继续下降，所以当前的基本可行解就是最优解。

所以，我们称$z_j-c_j$ 为判别数或检验数。

事实上，目标函数对非基变量的梯度$\nabla f(\boldsymbol x)|_{\boldsymbol x_N}=-(z_1-c_1,...,z_k-c_k,...)^T$，这对我们快速得出是否已经得到最优解很有帮助，即：

进基变量的选择

考虑还没有到达最优解的情况，即$\exists j\in R$ 使得$z_j-c_j\gt0$ 时，如果要让函数值下降最快，我们就需要最大化：

$$\max_{\boldsymbol x_N}\sum_{j\in R}(z_j-c_j)x_j$$

在调整基本可行解时，单纯形法的思想是“逐步”调整的。也就是从非基变量中选取一个变量，将其划归到基变量内，称为进基，选一个基变量中的变量划归到非基变量内，称为离基。所以这里不失一般性地，我们令$n-m$ 个自由未知量中的其中一个变量$x_k$ 为大于零的数，其余都为零，目的是将其作为进基变量。

在这样的情况下，就有：

$$\max_{\boldsymbol x_N}\sum_{j\in R}(z_j-c_j)x_j=(z_k-c_k)x_k\quad\text{where }k=\arg\max_{j\in R}(z_j-c_j)$$

因为$x_k$ 是自由未知量，所以选到最大的判别数$z_k-c_k$ 之后，只要$x_k$ 不断增大，目标函数值$f$ 就会下降越快。但是$x_k$ 不一定能无限增大，因为$\boldsymbol x_B=B^{-1}\boldsymbol b-B^{-1}N\boldsymbol x_N$，根据前面的假设，有$\boldsymbol x_B=B^{-1}\boldsymbol b-B^{-1}\boldsymbol p_kx_k:=\bar{\boldsymbol b}-\boldsymbol y_kx_k$。它作为基本可行解的一部分，要求$\boldsymbol x_B\geq \boldsymbol 0$。

首先，$B^{-1}\boldsymbol b=\bar{\boldsymbol b}\geq\boldsymbol 0$ 是一定的，因为它取自我们的初始可行解$\boldsymbol x^{(0)}$。但是，如果$\boldsymbol y_k\leq\boldsymbol 0$，那么$\boldsymbol x_B$ 就一定满足大于等于零的条件，此时$x_k$ 可以无限增大。这种情况就是原问题无界了。

接下来我们进一步考虑$\boldsymbol y_k\not\leq\boldsymbol 0$ 的情况。为了满足$\boldsymbol x_B\geq \boldsymbol 0$，我们需要$x_k$ 满足：

$$\forall i,\quad x_k\leq\frac{\bar b_i}{y_{ik}}$$

其中，$i$ 为基变量下标。

当然我们希望尽可能选最大的$x_k$ ，所以直接取：

$$x_k=\min_i\left\{\frac{\bar{b}_i}{y_{ik}}\bigg|y_{ik}\gt0\right\}=\frac{\bar{b}_r}{y_{rk}}$$

离基变量的选择

基变量中的$x_r$ 可以作为离基变量。其中，$r$ 是由上一节中$x_k=\bar{b}_r/y_{rk}$ 得到的。下面我们来论证为什么它可以作为出基变量。

$\boldsymbol p_k$ 是非基矩阵$N$ 的第$k$ 列，由于$\boldsymbol y_k=B^{-1}\boldsymbol p_k$ 所以：

$$\boldsymbol p_k=B\boldsymbol y_k=\sum_{i=1}^my_{ik}\boldsymbol p_{B_i}$$

它是基矩阵$B$ 的$m$ 个列向量$\boldsymbol p_{B_i}$ 的线性组合。
所以，原来$m$ 个线性无关的向量构成的基矩阵$B$ 将$\boldsymbol p_{B_r}$ 替换成$\boldsymbol p_k$ 之后依然线性无关。所以新的基本解$\boldsymbol x$ 仍然是一个基本可行解。

替换得到新的可行解后，在不进行新的调整时（取新的$\boldsymbol x_N=\boldsymbol 0$），则

$$f=\boldsymbol c_B\boldsymbol x_B=f_0-c_{B_r }x_r+c_kx_k$$

而$z_k=c_{B_r}·y_{rk},\;x_k=\bar{b}_r/y_{rk}=x_r/y_{rk}$ ，所以$f=f_0-(z_k-c_k)x_k$. 新的目标函数值比旧的值下降了$(z_k-c_k)x_k$。

🔔算法流程

综合以上论述，我们得到单纯形法的一般算法流程：

单纯形表 | 表格法

用单纯形法求解LP的过程，其实就是一个解线性方程组的过程。只是在求解过程中需要按照某种规则挑选自由未知量，但是其他变换过程就是简单的行变换。所以，我们可以通过预先的某种变化，将单纯形法计算过程中需要用到的变量（如$z_j,c_j,y_{ik}$ 等）组成一个矩阵/一张表。

下面我们在以前对原问题的转换中进一步修改，进行等价转换：

此时，整个约束条件就是一个关于$f(\boldsymbol x),\boldsymbol x_B$ 和$\boldsymbol x_N$ 的线性方程组。写出其增广矩阵：

快速填充方法

在实际利用表格法时，我们可以不用像上图那样严格按照基变量和非基变量的顺序填表。特别地，当选取的初始基矩阵$B$ 正好是$m$ 阶单位矩阵$I_m$ 时（请一定要记住这个前提，否则是不适用的！！！）我们可以很快地进行填表：

1. 填充约束$A\boldsymbol x=\boldsymbol b$ 的增广矩阵$[A,b]$（含不等式时引入松弛变量）；
2. 将$A$ 中可以组成$m$ 阶单位矩阵的元素作为基变量，按照 1 出现的次序依次填在最左边（表格外最左边一列的索引）；
3. 基变量对应的列，最后一行都填 0，非基变量 计算$\boldsymbol {c_By_j}-c_j$ 填在最后一行；
4. 计算初始目标函数值$f_0=\boldsymbol c_B\bar{\boldsymbol b}$ 填在最后（若$\boldsymbol c_B=\mathbf 0$ 则$f_0=0$）

以下面这个问题为例，我们给出一个快速填表的示例：

基本判别规则

填完表之后，因为最后一行是[判别数, 0, 目标函数值]。所以，

第一步：看是不是所有判别数都小于等于0，如果是，说明已经得到最优解。否则，
第二步：找到判别数最大（$z_k-c_k$）的那一列，这一列被称为主列。
第三步：查看主列的除最后一行外的所有值（其实就是$y_{ik}$），如果这些值全小于等于0，说明该问题无界。否则，
第四步：验证最后一列的值（其实就是$\bar{b}_i$）除以主列（$y_{ik}\gt0$ 的部分）后得到的所有结果，找到最小的那一行$r$ ，这一行被称为主行。

仍然以上面那么例子做演示，上面的四步可以可视化如下图：

主元消去法

由主行和主列得到的元素我们称为主元。
现在我们已经确定要将主列对应的变量$x_k$ 作为进基变量，$x_r$ 作为离基变量了。

整个可行解的调整过程很简单，就是利用主元消去法：对单纯形表作有限次行变换使得主元所在的位置调整为 1，主列其他位置为 0。

最后表格外最左边的基变量对应位置替换成$x_k$ 即可。

为什么主元消去法正好对应进基和离基过程呢？
首先将主元位置调整为 1 ，主列其他位置为 0。不考虑最后一行时，相当于线性方程组的等价变换，所以不改变问题。

接下来我们来思考主元消去法对最后一行的影响：

• 要使主元位置为 1，需要将第$r$ 行的元素除以$y_{rk}$ ；要使主列的最后一行为 0，需要在前面的基础上将最后一行减去第$r$ 行的$z_k-z_c$ 倍。

从而，非主列的其他最后一行的元素就变为：

$$(z_j-c_j)'=(z_j-c_j)-(y_{rj}/y_{rk})(z_k-c_k)$$

这是新的判别数。

$$(\boldsymbol c_BB^{-1}\boldsymbol b)'=\boldsymbol c_BB^{-1}\boldsymbol b-(\bar{b}_r/y_{rk})(z_k-c_k)$$

这就是下降后的新目标函数值。

此时，最后一列（不含最后一个）的值$\bar{\boldsymbol b}=\boldsymbol x_B$ 就是当前可行解基变量的值，再取非基变量值为0就能得到当前的可行解了。

如此往复进行迭代，直到满足最优解条件时，整个单纯形表的最后一列就对应了最优解和目标函数最优值。

根据新的判别数和新的目标函数值的更新公式，我们可以看出，当存在判别数$z_j-c_j=0$ 时，这一个非基变量$x_j$ 的新判别数小于或等于0。

• 新判别数小于0时，它以后不会再被考虑加入基变量；
• 新判别数等于0时，等到它被视为进基变量（因为只要还有判别数不小于0，就说明目标函数还有下降的可能，算法会继续找进基变量），此时新的目标函数值不会再变化。

从而我们可以得出一个结论：当$\exists j\in R,-\nabla f(\boldsymbol x)|_{\boldsymbol x_{N_j}}=0$ 时，LP有多个解（但不是无限个）。

上面那个例子的几次迭代与最终结果如下：

引入人工变量

前面我们介绍快速填充时，有一个强假设：系数矩阵中可以直接找出满秩单位矩阵$I_m$ 作为基矩阵。即便不找出单位阵，单纯形法的算法也要求我们至少能先找出一个基本可行解来。这个基本可行解满足$A\boldsymbol x=\boldsymbol b$。但现实中不一定能够这么幸运和容易的。

为了解决这个问题，我们可以人为引入变量$\boldsymbol x_\alpha$，使得约束中强行存在单位阵。这种变量被称为人工变量。引入方法如下：

与松弛变量不同，人工变量的引入会改变原来的约束！ 可以说它是一种“不合法”的变量。所以引入它之后，我们希望调整新的基本可行解使得在这个解下面$\boldsymbol x_\alpha=\boldsymbol 0$. 如果是这样，那么就等价于$A\boldsymbol x+\boldsymbol x_\alpha=A\boldsymbol x=\boldsymbol b$ 成立，我们也就找到了一个原问题的基本可行解了。

🔔两阶段法

两阶段法提出，在求解LP的第一阶段，先用单纯形法求解一个辅助LP：

$$\begin{aligned} \min\quad&\boldsymbol{e^Tx}_\alpha\\ \text{s. t. }\quad&\boldsymbol{Ax+x_\alpha=b}\\ &\boldsymbol{x\geq0,x_\alpha\geq0} \end{aligned}$$

其中，$\boldsymbol e=(1,1,...,1)^T$ 是分量全1的$m$ 维列向量。

显然，由于$\boldsymbol x_\alpha\geq0$ ，所以理想情况下目标函数$\boldsymbol e^T\boldsymbol x_\alpha\geq0$，当且仅当$\boldsymbol x_\alpha=\boldsymbol 0$ 时得到它的最小化。因此，两阶段法的第一阶段就是这样简单的数学原理和单纯形法的机制找到$\boldsymbol x_\alpha=\boldsymbol 0$ 的可能，以得到原问题的基本可行解。

假如求解这个辅助LP得到最优基本可行解$(\overline{\boldsymbol x}^\top,\overline{\boldsymbol x}^\top_\alpha)$，那么可能存在下列三种情况：

1. $\overline{\boldsymbol x}_\alpha\neq\bf0$ ：原始LP没有可行解（约束条件与目标函数不兼容或没有非负解）；
2. $\overline{\boldsymbol x}_\alpha=\bf0$ 且分量都是非基变量：此时$\overline{\boldsymbol x}$ 是原始LP的一个可行解；
3. $\overline{\boldsymbol x}_\alpha=\bf0$ 但是它有分量是基变量：利用主元消去法把基变量中含有的人工变量人为地替换出去。

Note：第三种情况提到的主元消去法与表格法略有不同，可以不遵循单纯形法限定的离基和进基规则，我们的目的主要是希望把人工变量转换成非基变量，从而得到原始LP的只与$\boldsymbol x$ 有关联的基本可行解。

当一阶段得到基本可行解之后，两阶段法的第二阶段就和单纯形法一样了。并且在一阶段时留下的单纯形表还可以继续使用（需要把最后一行的判别数与函数值按照原始LP的目标函数进行修正）。

🔔大M法

大$M$ 法给出了结合两个阶段的LP来计算的策略：合并目标函数，增加一个惩罚因子$M$ ，使得合并得到的新目标函数为$\boldsymbol{cx}+M\boldsymbol{e^Tx}_\alpha$ 。

当$M$ 足够大时，为了迫使目标函数最小化，理想情况下有$\boldsymbol{e^Tx}_\alpha=\bf0$ ，$\boldsymbol{cx}$ 达到最小。也就是原LP找到了最优解。

但是和两阶段法一样，求解这个新LP可能存在下列几种情况：

1. 得到最优解$(\overline{\boldsymbol x}^\top,\overline{\boldsymbol x}^\top_\alpha)$，且$\overline{\boldsymbol x}_\alpha=\bf0$：此时的$\overline{\boldsymbol x}$ 就是原始LP的最优解；
2. 得到最优解$(\overline{\boldsymbol x}^\top,\overline{\boldsymbol x}^\top_\alpha)$，但$\boldsymbol{e^Tx}_\alpha\gt0$：此时原始LP没有可行解；
3. 新LP无界/不存在有限最优值，且$\overline{\boldsymbol x}_\alpha=\bf0$：此时原始LP也无界；
4. 新LP无界/不存在有限最优值，但$\boldsymbol{e^Tx}_\alpha\gt0$：此时原始LP没有可行解。

单个人工变量技巧

之前我们探讨了我们为了保险起见——找到基本可行解，引入了$m$ 个人工变量，及其构造辅助LP的方法。但事实上，我们还可以只引入一个人工变量来解决问题。

首先我们同样将系数矩阵$A$ 拆分为$[B,N]$ 。但是这次只要求$B$ 是满秩的，不需要选出来的变量一定是可行的基变量。首先同样有：

$$B\boldsymbol x_B+N\boldsymbol x_N=\boldsymbol b$$

两边同时左乘$B^{-1}$，令$\bar{\boldsymbol b}=B^{-1}\boldsymbol b$，就有：

$$\boldsymbol x_B+B^{-1}N\boldsymbol x_N=\bar{\boldsymbol b}$$

取$N$ 对应的那些变量$\boldsymbol x_N=\bf0$，就能得到一个解：$\boldsymbol x=\begin{bmatrix}\boldsymbol x_B\\\boldsymbol x_N\end{bmatrix}=\begin{bmatrix}\bar{\boldsymbol b}\\\bf0\end{bmatrix}$.

当然，我们需要让这个解满足非负约束$\boldsymbol x\geq\bf0$，也就是$\bar{\boldsymbol b}\geq\bf0$。如果确实如此，那么这个解就完全是一个基本可行解，我们“瞎猫碰到死耗子”——随便找就找到符合条件的解了，此时直接用单纯形表来做就OK了。

通常我们不是先找任意一个$B$ 再通过它的逆去判断$\bar{\boldsymbol b}\geq\bf0?$ ，而是对约束作行变换找出一个单位矩阵$I_m$ 出来，再看最右边那一列是否有小于0的情况。

不过现在我们更想关注$\bar{\boldsymbol b}\not\geq\bf0$ 时该怎么处理？

我们引入一个人工变量$x_a\geq0$，构建新的约束：

$$\boldsymbol x_B+B^{-1}N\boldsymbol x_N-x_a\boldsymbol e=\bar{\boldsymbol b}$$

记$\bar{\boldsymbol b}=(\bar{b}_1,...,\bar{b}_m)^T$。现在我们想把人工变量通过”进基“的方式变成基变量，使得新的$\bar{\boldsymbol b}\geq0$.
首先我们找到所有$\bar{\boldsymbol b}$ 的负值中最小的那一个：

$$\bar{b}_r=\min_{i=1,..,m}\{\bar{b}_i\}\lt0$$

从而以第$r$ 行为主行，$x_a$ 所在的那一列作为主列，通过主元消去法把$x_a$ 引进基。
这样一来，新的$\bar{\boldsymbol b}$ 就能保证一定大于等于0（证明略），包含有人工变量的基本可行解也就找到了。之后，再用两阶段法或者大$M$ 法处理即可。

退化情形与摄动法

待更

修正单纯形法 | Correctional Simplex Method

待更

LP的对偶理论

每一个LP都存在另一个与它密切相关的LP，我们称其中之一为原问题，另一个称为它的对偶问题。它们之间的内在联系为进一步深入研究LP的理论与算法提供了理论依据。

下面我们以最经典的食谱问题切入：
假如针对维生素的配比和蛋奶的选购，有对应的“卖家”和“买家”。食谱的关系表如下。

对于“卖家”^[注1]来说，他希望对维生素$V_c,V_b$ 分别定价为$x,y$，使得自己的利润最大。但是，他的定价不能比市场上的奶和蛋的流行价更高，否则就会失去他的顾客。

对于“买家”来说，他希望能购买$m$ 件奶和$n$ 件蛋，花最少的钱买到能达到自身每日需求的维生素。

于是，从他们二者的角度出发，我们可以写出两个线性规划^[注2]：

注1：这里所谓的“卖家”，是指“卖维生素的商家”，与奶和蛋的商家是竞争关系。
注2：分别求解这两个LP可以得出一个有意思的结论——它们的最优目标函数值相等。

对偶问题的转换

线性规划中的对偶可以概括为三种形式。

对称形式的对偶

原问题

$$\begin{aligned} \min\quad&\boldsymbol{cx}\\ \text{s. t. }\quad&\boldsymbol{Ax\geq b}\\ &\boldsymbol{x\geq0} \end{aligned}$$

对偶问题

$$\begin{aligned} \max\quad&\boldsymbol{wb}\\ \text{s. t. }\quad&\boldsymbol{wA\leq c}\\ &\boldsymbol{w\geq0} \end{aligned}$$

其中，$\boldsymbol A\in\Bbb R^{m\times n}$ 是原问题约束的系数矩阵，而对偶问题的变量$\boldsymbol{w}$ 的长度是$m$ ——原问题的约束个数（除非负约束）。对应地，对偶问题的约束个数（除非负约束）正好是原问题变量$\boldsymbol x$ 的长度。

可以证明：对偶问题的对偶为原问题.

非对称形式的对偶

考虑具有等式约束的LP：

$$\begin{aligned} \min\quad&\boldsymbol{cx}\\ \text{s. t. }\quad&\boldsymbol{Ax=b}\\ &\boldsymbol{x\geq0} \end{aligned}$$

可以将其利用不等式规则：$(x\geq a)\land(x\leq a)\to x=a$，进行等价：

$$\begin{aligned} \min\quad&\boldsymbol{cx}\\ \text{s. t. }\quad&\boldsymbol{Ax\geq b}\\ &\boldsymbol{-Ax\geq -b}\\ &\boldsymbol{x\geq0} \end{aligned}$$

从而其对偶问题可以写作：

$$\begin{aligned} \min\quad&\boldsymbol{ub-vb}\\ \text{s. t. }\quad&\boldsymbol{uA-vA\leq c}\\ &\boldsymbol{u,v\geq0} \end{aligned}$$

令$\boldsymbol{w=u-v}$ 可得：

$$\begin{aligned} \min\quad&\boldsymbol{wb}\\ \text{s. t. }\quad&\boldsymbol{wA\leq c} \end{aligned}$$

此时，$\boldsymbol w$ 无非负约束！

一般情形下的对偶

当原问题约束存在多种不等式约束时，我们可以引入松弛变量$\boldsymbol x_s,\boldsymbol x_t$ 将其化为非对称形式的对偶。整个推导过程如下：

经过这样的推导，我们还可以总结出以下规律，便于我们以后快速得出对偶问题：

对偶定理及其推论

下面研究对偶的基本性质. 此处仅叙述对称形式的对偶，但对其他形式仍是成立的.

定理

设$\boldsymbol x^{(0)},\boldsymbol w^{(0)}$ 分别是原问题$(L)$ 和其对偶问题$(D)$ 的可行解，则$\boldsymbol{cx}^{(0)}\geq\boldsymbol w^{(0)}\boldsymbol b$.

证明

这个定理指出，极小化问题给出了极大化问题的目标函数上界；极大化问题给出了极小化问题的目标函数下界。对应的，如果对偶问题无上/下界则说明原问题无可行解。

根据该定理还可以得出以下推论：

• 若$\boldsymbol x^{(0)},\boldsymbol w^{(0)}$ 分别是$(L)$ 和$(D)$ 的可行解，且$\boldsymbol{cx}^{(0)}=\boldsymbol w^{(0)}\boldsymbol b$，则它们分别是对应问题的最优解；
• $(L)$ 和$(D)$ 有最优解的充要条件是它们同时有可行解。

---

定理

若原问题$(L)$ 和其对偶问题$(D)$ 中有一个问题存在最优解，则另一个一定也存在最优解，并且它们的最优值相等.

证明

• 推论：若$(L)$ 存在一个对应基$B$ 的最优基本可行解，则单纯形乘子$\boldsymbol w=\boldsymbol c_BB^{-1}$ 是$(D)$ 的一个最优解。

互补松弛性质

设$\boldsymbol x^{(0)},\boldsymbol w^{(0)}$ 分别是原问题$(L)$ 和其对偶问题$(D)$ 的可行解，那么$\boldsymbol x^{(0)},\boldsymbol w^{(0)}$ 都是最优解的充要条件是，$\forall i,j\quad 1\leq i\leq m;\;1\leq j\leq n$ 有：

1. 如果$x_j^{(0)}\gt0$，就有$\boldsymbol w^{(0)}\boldsymbol p_j=c_j$；
2. 如果$\boldsymbol w^{(0)}\boldsymbol p_j\lt c_j$，就有$x_j^{(0)}=0$；
3. 如果$w_i^{(0)}\gt0$，就有$\boldsymbol A_i\boldsymbol x^{(0)}=b_i$；
4. 如果$\boldsymbol A_i\boldsymbol x^{(0)}\gt b_i$，就有$w_i^{(0)}=0$。

对于非对称形式，则满足前两条即可。

这个规则可以用文字描述如下：
问题$(L)$ 的可行解$\boldsymbol x^{(0)}$ 的第$j$ 个变量$x_j^{(0)}$ 如果大于0，那么对偶问题$(D)$ 的第$j$ 行不等式约束取等号；相反，如果代入$\boldsymbol w^{(0)}$ 到第$j$ 行不等式约束中发现不能取等号，则说明$x_j^{(0)}=0$。
同理，$(D)$ 的可行解$\boldsymbol w^{(0)}$ 第$i$ 个变量$w_i^{(0)}$ 如果大于0 则$(L)$ 的第$i$ 行不等式约束取等号。

---

通过互补松弛性质，我们可以快速地根据对偶问题的最优解直接求出其原问题的最优解——根据对偶问题的最优解与0的大小关系，列出原问题最优解应该满足的等式线性方程组，并求解之。

比如下面这个例题：

对偶单纯形法

考虑我们之前一直研究的标准型 LP（记作$L$）：

$$\begin{aligned} \min\quad&\boldsymbol{cx}\\ \text{s. t. }\quad&\boldsymbol{Ax=b}\\ &\boldsymbol{x\geq0} \end{aligned}$$

之前我们为了得到一个基本可行解以便使用单纯形法，往往引入人工变量并通过两阶段法来处理。而利用对偶性质，我们可以给出一种不用人工变量的方法。为此，我们先引入对偶可行的基本解的概念：

设$\boldsymbol x^{(0)}$ 是标准 LP 问题$(L)$ 的一个基本解（可以“不可行”，即不一定满足约束条件），它对应的基矩阵为$\boldsymbol B$，单纯形乘子$\boldsymbol {w=c_B B}^{-1}$。如果$\boldsymbol w$ 是$(L)$ 的对偶问题$(D)$ 的可行解，即满足对偶问题$(D)$ 的所有约束条件，亦即$\forall j,\;\boldsymbol {wp}_j-c_j\leq0$，亦即其所有判别数都小于等于零，那么就称$\boldsymbol x^{(0)}$ 是$(L)$ 的对偶可行的基本解。

基本思想与计算步骤

原版单纯形法是保持原问题基本解可行的情况下（其实就是保持右端列$\bar{\boldsymbol b}$ 非负），不断调整进基变量和出基变量使得判别数均小于或等于零（此时目标函数值再无可下降的可能性）。

而对偶单纯形法的思路则是不断改进对偶问题的基本可行解，当对偶基本可行解改进后也能使得原问题的基本解变得可行时（其实就是右端列$\bar{\boldsymbol b}$ 变得非负），那么根据对偶性质可知，此时两个问题均达到了最值。

为了尽可能让右端列$\bar{\boldsymbol b}$ 变得非负，对偶单纯形法的步骤是先选择离基变量再选择出基变量。
其中，离基变量我们选择的是$\boldsymbol {\bar b}$ 负数中最小的那一行，因为这样可以在主元消去时尽可能让更多的$b_i$ 非负。即：

$$\bar{b}_r=\min_i\{\bar{b}_i\}\lt0$$

然后是进基变量的选择。为了保证主元消去后判别数仍然小于或等于零（即保证$\boldsymbol w$ 始终是对偶问题的可行解），我们需要找到判别数除以主行元素后最小的那一列，即：

$$\frac{z_k-c_k}{y_{rk}}=\min_{j}\left\{\frac{z_j-c_j}{y_{rj}}\bigg|y_{rj}\lt0\right\}$$

• 当主行无负元，即$\forall j,\;y_{rj}\geq0$ 时，原问题中的变量$\boldsymbol x$ 取任意非负的值都不能满足第$r$ 个约束条件，此时原问题无可行解。
• 对偶问题的最优解就是原问题最终单纯形表中松弛变量的判别数的相反数。（可根据对偶性质的推论得来）

初始对偶可行的基本解

对偶单纯形法同样面临一个问题，那就是要先找到一个对偶可行的基本解，也就是说要找到一个可以使得判别数非正的基本解。为此，类比两阶段法，我们需要先解一个扩充问题。

我们将原问题的约束条件$\boldsymbol{(B,N)·\begin{pmatrix}\boldsymbol x_B\\\boldsymbol x_N\end{pmatrix}=b}$ 两边同时乘以$B^{-1}$，得到：

$$\begin{aligned} B^{-1}B\boldsymbol x_B+B^{-1}N\boldsymbol x_N&=B^{-1}\boldsymbol b\\\\ i.e.\quad\; \boldsymbol x_B+\sum_{j\in R}\boldsymbol y_jx_j&=\bar{\boldsymbol b} \end{aligned}$$

其中，$R$ 为非基变量下标集，$\boldsymbol y_j=B^{-1}\boldsymbol p_j$，$\boldsymbol x$ 是我们任意给定的基本解（可以不是可行解，也可以不是对偶可行的基本解），它可以按线性无关性分为$\boldsymbol x_B,\boldsymbol x_N$。

此外，我们再增加一个变量$x_{n+1}$ 和约束条件：

$$\sum_{j\in R}x_j+x_{n+1}=M$$

其中，$M$ 是一个充分大的正数。从而我们得到扩充问题：

$$\begin{aligned} \min\quad&\boldsymbol{cx}\\ \text{s. t. }\quad&\boldsymbol x_B+\sum_{j\in R}\boldsymbol y_jx_j=\bar{\boldsymbol b}\\ &\sum_{j\in R}x_j+x_{n+1}=M\\ &\boldsymbol{x\geq0},\;x_{n+1}\geq0 \end{aligned}$$

通过扩充问题，我们得到如下求解策略：

1. 利用"原始" Simplex 进行主元消去，直到$-\nabla f\boldsymbol{(x)\leq0}$，即判别数全部非正；
2. 利用对偶 Simplex 进行主元消去，直到$\bar{\boldsymbol b}\geq0$。

注：上述策略中，“原始” Simplex 是打双引号的。
这主要是主元消去时的第一步我们选择完主列之后，要无条件选择$m+1$ 的这一行作为主行，这一步可以证明仅一次操作就能使得$-\nabla f\boldsymbol{(x)\leq0}$。

此外还有，可以证明：

1. 扩充问题没有可行解时，原问题也没有可行解；
2. 扩充问题的最优值与$M$ 无关时，即$x_{n+1}$ 对应判别数为0时，原问题的最优解就是前$n$ 个变量在扩充问题中的取值；
3. 扩充问题的最优值与$M$ 有关时，即$x_{n+1}$ 对应判别数小于0时，原问题无下界。

原始-对偶算法

待更

C++实现

建立单纯形法类

Simplex.h 文件

#ifndef SIMPLEX_H
#define SIMPLEX_H

#include <cstdlib>

struct solution{
    double f; //目标函数最优值
    double *discri; //判别数向量（单纯形表最后一行）
    double *x; //最优基本可行解
    int tag=0; //tag=0有基本可行解,-1无界
};

template<std::size_t M, std::size_t N>
class Simplex{
    private:
        int m = M, n = N; //约束条件个数m 和自变量个数 n
        double (*A)[N]; //约束系数矩阵
        double *b; //约束常数项向量
        double b_bar[M]; //=B^(-1)b
        double *c; //目标函数系数向量
        int base_list[M]; //基变量的索引

        solution res; //优化结果

    public:
        Simplex<M,N>(double A[M][N], 
                double b[M], 
                double c[N]):
                A(A),b(b),c(c){}

        ~Simplex(){}

        void UpdateBaseVars_Im(); //更新基变量
        void UpdateDIscri(); //更新判别数
        void Solve_Im(); //处理基矩阵是单位矩阵的问题
        void TwoStepMethod(); //两阶段法
        solution GetSolution(); //获取最优解

};

#endif

第一次迭代

接下来是针对第一次迭代进行初始化。

第一步是先找出可以组成单位矩阵$I_m$ 的基变量。这里用 base_list[M] 存储，内容是基变量对应的下标。通过类方法 Simplex::UpdateBaseVars_Im() 实现。

第二步是根据公式计算判别数。这里用 discri[N] 存储，基变量对应的判别数置0，对应着单纯形表的最后一行。通过类方法 Simplex::UpdateDIscri() 实现。

Simplex.cpp 文件

#include "Simplex.h"
#include "string.h"
#include <iostream>

using namespace std;

template<std::size_t M, std::size_t N>
void Simplex<M,N>::UpdateBaseVars_Im(){

    memset(base_list, 0, sizeof(int)*m);

    int cnt, idx, k=0;
    for(int j = 0; j < n; j++){
        cnt = 0; idx = -1;
        for(int i = 0; i < m; i++){
            if(A[i][j] != 0 && A[i][j] != 1) break; //第j列不满足单位矩阵要求
            if(A[i][j] == 0) cnt++;
            if(A[i][j] == 1) idx=i;
        }
        if(cnt == m-1 && idx != -1 && A[idx][j] == 1){
            base_list[idx] = j; //基变量的顺序应该与1所在的行数对应
            // 测试：cout << j << " ";
        }
    }    
}

template<std::size_t M, std::size_t N>
void Simplex<M,N>::UpdateDIscri(){

    res.discri = new double[n];
    memset(res.discri, 0, sizeof(double)*n);

    for(int j = 0; j < n; j++){
        for(int i = 0; i < m; i++){
            if(base_list[i] == j){
                res.discri[j] = 0; //基变量对应的最后一行元素为0
                break;
            }else{
                res.discri[j]+=A[i][j]*c[base_list[i]];
            }
        }
        res.discri[j] -= c[j];
        // 测试：cout << res.discri[j] << " ";
    }//现行判别数
}

多次迭代求解

接下来是从第二次迭代开始直到算法结束的整个迭代过程。

首先计算现行目标函数值 f，然后根据 discri[N] 找到主列 main_col，如果判别数全小于0则说明达到最优解，否则如果 A[main_col] 所有元素都小于0则说明原问题无界。

如果没有以上两种问题，则根据 b_bar[M]/A[main_col] 寻找主行 main_row，最后利用主元消去法进行矩阵 A 的更新、判别数 discri[N] 和基变量 base_list[M] 的更新。

为了应对可能出现循环的那一类 LP 问题，通过设置最大迭代次数 iter_nums 规避无限循环。

最终的最优解由 b_bar[M] 对应的基变量取特定的值得到，而非基变量的取值为0。通过类方法 Simplex::GetSolution() 得到结果。

Simplex.cpp 文件

template<std::size_t M, std::size_t N>
void Simplex<M,N>::Solve_Im(){

    UpdateBaseVars_Im(); //找到现行基变量

    res.f = 0;
    for(int i = 0; i < m; i++){
        b_bar[i] = b[i]; //现行 b_bar = B^(-1)b
        res.f += c[base_list[i]]*b_bar[i]; //现行目标函数值
    }

    UpdateDIscri(); //更新判别数

    int cnt,T=1;
    int iter_nums=100;//迭代次数，防止出现死循环
    int main_col;
    int main_row;
    double A_tmp[M][N];
    
    while(iter_nums--){ //不断迭代，直到找到最优解

        main_col = 0;
        main_row = 0;
        for(int j = 0; j < n; j++){
            if(res.discri[main_col] < res.discri[j]){
                main_col = j;
            }
        }
        if(res.discri[main_col] <= 0) return; //得到最优解, 否则主元消去法

        cnt = 0;
        for(int i = 0; i < m; i++){
            if(A[i][main_col] < 0) cnt++;
            if(A[main_row][main_col] < 0) main_row = i;
            if(b_bar[main_row]/A[main_row][main_col] > b_bar[i]/A[i][main_col] && A[i][main_col] >= 0){
                main_row = i;
            }
        }
        if(cnt == m){res.tag = -1; return;} //LP无界
        
        //主元消去法
        res.f -= b_bar[main_row]/A[main_row][main_col]*res.discri[main_col];

        for(int i = 0; i < m; i++){
            if(i == main_row) b[i] = b_bar[i]/A[main_row][main_col];
            else b[i] -= b_bar[main_row]/A[main_row][main_col]*A[i][main_col];
        }
        for(int j = 0; j < n; j++){
            if(j == main_col) continue;
            res.discri[j] -= A[main_row][j]/A[main_row][main_col]*res.discri[main_col];
        }
        res.discri[main_col] = 0;
        
        T++;
        cout << "＠第" << T << "次迭代：" << endl;
        cout << "--------------------------------------------------------" << endl;

        for(int i = 0; i < m; i++){
            for(int j = 0; j < n; j++){
                if(i == main_row) A_tmp[i][j] = A[i][j]/A[main_row][main_col];
                else A_tmp[i][j] = A[i][j]-A[main_row][j]/A[main_row][main_col]*A[i][main_col];
                cout << A_tmp[i][j] << "\t";
            }
            b_bar[i] = b[i];
            cout << "|" << b[i] << endl;
        }
        cout << "--------------------------------------------------------" << endl;

        for(int i = 0; i < m; i++){
            for(int j = 0; j < n; j++){
                A[i][j] =A_tmp[i][j];
            }
        }
        for(int j = 0; j < n; j++){
            cout << res.discri[j] << "\t";
        }
        cout << "f=" << res.f << endl;
        cout << "---------------------------------------------------------" << endl << endl;

        base_list[main_row] = main_col; //进离基更换
    }

}

template<std::size_t M, std::size_t N>
solution Simplex<M,N>::GetSolution(){
    res.x = new double[n];
    memset(res.x, 0, sizeof(double)*n);
    for(int i = 0; i < m; i++){
        res.x[base_list[i]] = b_bar[i];
    }
    return res;
}

求解测试

main.cpp 文件

#include <iostream>
#include <cstdlib>
#include "Simplex.h"
 // #include "Simplex.cpp"

using namespace std;

int main(void){
    
    // double A[3][6] = {{1,1,-2,1,0,0},
    //                   {2,-1,4,0,1,0},
    //                   {-1,2,-4,0,0,1}};
    // double b[3] = {10,8,4};
    // double c[6] = {1,-2,1,0,0,0};

    // Simplex<3,6> myQuestion(A,b,c);

    double A[3][5] = {{1,-2,1,0,0},
                      {0,1,-3,1,0},
                      {0,1,-1,0,1}};
    double b[3] = {2,1,2};
    double c[5] = {0,-1,2,0,0};

    Simplex<3,5> myQuestion(A,b,c);

    myQuestion.Solve_Im();
    solution res = myQuestion.GetSolution();
    
    cout << "最优解: [ ";
    for(int i = 0; i < 5; i++) cout << res.x[i] << ", ";
    cout << "];" << endl;
    cout << "最优值: f_min = " << res.f << endl;

    return 0;
}

Python

scipy.optimize模块也能实现，但不及下述两种库，此处略过

PuLP库

from pulp import *
import pandas as pd
import numpy as np
import math

# Read CSV files and Process them
WB = pd.read_csv('BCHAIN-MKPRU.csv')
WG = pd.read_csv('LBMA-GOLD.csv')
WBWG = pd.merge(WB, WG, how='outer', on='Date')
WBWG.to_csv('MergeData.csv',index=False)

# The day that market is not open
NanIdx = WBWG[WBWG['USD (PM)'].isnull()].index

# Setting LpProblem and Variable
prob = LpProblem("The_Ideal_situation",LpMaximize)
xGs = [f'xG{i}' for i in range(0,1826)]
xBs = [f'xB{i}' for i in range(0,1826)]
iGs = [f'iG{i}' for i in range(0,1826)]
iBs = [f'iB{i}' for i in range(0,1826)]
xG = LpVariable.dicts("xG",xGs,0,None,LpContinuous)
xB = LpVariable.dicts("xB",xBs,0,None,LpContinuous)
iG = LpVariable.dicts("iG",iGs,-1,1,LpInteger)
iB = LpVariable.dicts("iB",iBs,-1,1,LpInteger)

# Objective function
C = 1000
G = 0
B = 0
for i in range(10):
    C += -(iB[f'iB{i}']+0.0002)*xB[f'xB{i}']*WBWG.loc[i,'Value']
    G += xG[f'xG{i}']
    B += xB[f'xB{i}']
    if not math.isnan(WBWG.loc[i,'USD (PM)']):
        C += -(iG[f'iG{i}']+0.0001)*xG[f'xG{i}']*WBWG.loc[i,'USD (PM)']

prob += C+WBWG.loc[10,'Value']*B+WBWG.loc[10,'USD (PM)']*G

# Subject to...
C = 1000
G = 0
B = 0
for i in range(10):
    C += -(iB[f'iB{i}']+0.0002)*xB[f'xB{i}']*WBWG.loc[i,'Value']
    G += xG[f'xG{i}']
    B += xB[f'xB{i}']
    if not math.isnan(WBWG.loc[i,'USD (PM)']):
        C += -(iG[f'iG{i}']+0.0001)*xG[f'xG{i}']*WBWG.loc[i,'USD (PM)']
    prob += C >= 0
    prob += G >= 0
    prob += B >= 0

for i in NanIdx:
    if i > 9:
        break
    prob += var[f'xG{i}'] == 0

# Get the solution
prob.to_json("Ideal_solution.json")
prob.solve()

print("\n" ,"status: ",LpStatus [prob.status] ,"\n")
Name = []
Value = []
for v in prob.variables() :
    Name.append(v.name)
    Value.append(v.varValue)
pd.DataFrame({"VarName":Name,"VarValue":Value}).to_csv("IdealResult.csv",index=False)
print("Maximun Advantage =", "Rs", value(prob.objective))

cvxpy库

$$\max_{x_G^k,x_B^k,1\leq k\leq K} C^K+G^KW_G^K+B^KW_B^K\\ s.t.\begin{cases} C^0=1000,G^0=0,B^0=0\\ C^{k+1}=C^k-x_G^kW_G^k-x_B^kW_B^k-\alpha_\mathbf{bitcoin}\%|x_G^k|W_G^k-\alpha_\mathbf{gold}\%|x_B^k|W_B^k\\ G^{k+1}=G^k+x_G^k\\ B^{k+1}=B^k+x_B^k\\ C^k,G^k,B^k \geq 0 \\ k=1,2,...,K \end{cases}$$

以上述模型为例，

import math
import pandas as pd
import numpy as np
import cvxpy as cp

# Read CSV files and Process them
WB = pd.read_csv('BCHAIN-MKPRU.csv')
WG = pd.read_csv('LBMA-GOLD.csv')
WBWG = pd.merge(WB, WG, how='outer', on='Date')
WBWG.to_csv('MergeData.csv',index=False)

# The day that market is not open
NanIdx = WBWG[WBWG['USD (PM)'].isnull()].index

# Setting LpProblem and Variable
C = 1000;G = 0;B = 0

# Set Iterative
cons = []
xG = cp.Variable(1826)
xB = cp.Variable(1826)
for i in range(1826):
    C -= 0.0002*WBWG.loc[i,'Value']*cp.abs(xB[i])+xB[i]*WBWG.loc[i,'Value']
    G += xG[i]
    B += xB[i]
    if not math.isnan(WBWG.loc[i,'USD (PM)']):
        C -= 0.0001*WBWG.loc[i,'USD (PM)']*cp.abs(xG[i])+xG[i]*WBWG.loc[i,'USD (PM)']
    cons.append(C >= 0)
    cons.append(G >= 0)
    cons.append(B >= 0) # Subject to...

# Objective function
func = cp.Maximize(C+G*WBWG.loc[1825,'USD (PM)']+B*WBWG.loc[1825,'Value'])

# Subject to...
for i in NanIdx:
    cons.append(xG[i] == 0)

# Get the solution
prob = cp.Problem(func,cons)
prob.solve()
print("\nstatus:", prob.status)
print("\nThe Maximun Advantage is: ", prob.value)
pd.DataFrame({"Gold":np.around(xG.value, 2),"Bitcoin":np.around(xB.value, 2)}).to_csv("IdealResult.csv",index=False)

参考

1. 陈宝林 编著. 《最优化理论与算法》. 清华大学出版社. ISBN : 978-7-320-11376-8
2. 凸优化高质量开源课程 - 美国卡内基梅隆大学