线性规划：单纯形算法

最新推荐文章于 2023-03-08 16:49:35 发布

胡拉哥

最新推荐文章于 2023-03-08 16:49:35 发布

阅读量793

点赞数 4

分类专栏：优化算法文章标签：线性规划单纯形算法算法运筹优化

本文链接：https://blog.csdn.net/qx3501332/article/details/118460050

版权

优化算法专栏收录该内容

21 篇文章 55 订阅

订阅专栏

考虑线性规划的标准形式（更多介绍参考《线性规划的标准形》）：
$\begin{aligned} \min~ & c^T x\\ \text{s.t.}~ & Ax=b\\ & x\geq 0 \end{aligned}$
其中 $\in \mathbb{R}^n$ ， $A\in\mathbb{R}^{m\times n}$ ， $b\in\mathbb{R}^m \geq \mathbf{0}$ 。

单纯形算法的思路是从多面体的一个顶点出发，然后沿着降低目标的方向，迭代到另一个顶点，直到目标值无法降低，于是得到最优解。

基本可行解

把 $A$ 拆成两个部分：
$\quad N]$
其中 $B\in \mathbb{R}^{m\times m}$ ， $\in \mathbb{R}^{m \times (n-m)}$ 。

相应地，把 $c$ 和 $x$ 分别拆成两个部分，即
$\begin{bmatrix} c_B \\ c_N\end{bmatrix} \quad x = \begin{bmatrix} x_B \\ x_N \end{bmatrix}.$
注意： $c_B, x_B$ 的下标与 $B$ 的列分别对应。

假设 $A$ 满秩，我们可以交换列的位置保证 $B$ 可逆，所以 $B$ 的列向量是一组基。 $B$ 称为 基矩阵， $N$ 称为 非基矩阵。把 $x_B$ 对应的变量称为 基变量， $x_N$ 对应的变量称为 非基变量。

这样一来，约束条件可以写成
$\quad N] \begin{bmatrix} x_B\\ x_N \end{bmatrix} = b \Leftrightarrow Bx_B + N x_N = b.$
我们有
$x_B = B^{-1}b - B^{-1}Nx_N.$
令 $x_N = 0$ ，得到 基本解
$x_0 = \begin{bmatrix} B^{-1}b\\ 0 \end{bmatrix}.$
换句话说，基本解 $x_0$ 满足约束条件 $Ax_0 = b$ ，但有可能不可行。如果 $x_0 \geq \mathbf{0}$ ，那么称它为 基本可行解。换句话说，基本可行解是多面体的一个顶点（证明略）。

假设已知一个初始的基本可行解，接下来我们要想办法迭代到另一个基本可行解，并且使得目标函数值下降。

迭代条件

要想目标函数降低，我们自然想到梯度下降法。令 $f (x)$ 代表目标函数，梯度 $\nabla f(x)$ 就是函数值增加的方向。所以应该沿着负梯度的方向移动 $x$ 。

把 $x_B$ 代入目标函数 $f(x):=c^Tx$ 可得
$\begin{aligned} f(x) & = c^Tx \\ & = c_B^Tx_B + c_N^T x_N \\ & = c_B^T B^{-1}b + ( c_N^T - c_B^T B^{-1}N) x_N. \end{aligned}$
对 $f (x)$ 关于 $b$ 求导，得到 Shadow Price （也称为 对偶变量）：
$\lambda^T := c_B^T B^{-1}.$
对 $f (x)$ 关于 $x_N$ 求导（求梯度），得到 Reduced Cost ：
$\mu^T := c_N^T - c_B^TB^{-1}N = c_N^T - \lambda^T N.$
令 $J$ 代表非基变量的下标。 $\forall j\in J$ ， $x_j$ 每增加一个单位，目标函数就增加 $\mu_j$ 。

如果存在 $\mu_j < 0$ ，我们只要把 $x_j$ 从 0 增加 $\delta$ ，那么目标函数就可以降低 $-\mu_j \delta$ 。

当 $\mu_j \geq 0$ ， $\forall j\in J$ ，这意味着目标函数值无法降低，此时达到最优解（证明略）。

如何迭代

假设存在 $j\in J$ 使得 $\mu_j < 0$ ，于是我们想增加 $x_j$ 使得目标函数降低。但是问题来了，增加多少可以保证 $x$ 可行？

回顾基本可行解的定义
$\begin{bmatrix} x_B\\ x_N \end{bmatrix} = \begin{bmatrix} B^{-1}b - B^{-1}Nx_N \\ x_N \end{bmatrix}.$
注意到 $x_N = 0$ ，我们只要保证 $x_B \geq 0$ 。接下来把 $x_B$ 换一种写法，令 $\tilde{b} = B^{-1}b$ ， $\tilde{a}_j = B^{-1}a_j$ ，我们有
$x_B = \begin{bmatrix} x_{B_1}\\ x_{B_2}\\ \vdots\\ x_{B_m}\end{bmatrix} = \begin{bmatrix} \tilde{b}_1\\ \tilde{b}_1\\ \vdots\\ \tilde{b}_m \end{bmatrix} - \begin{bmatrix} \tilde{a}_{j_1}\\ \tilde{a}_{j_2}\\ \vdots\\ \tilde{a}_{j_m} \end{bmatrix}x_j = \begin{bmatrix} \tilde{b}_1 - \tilde{a}_{j_1}x_j\\ \tilde{b}_2 - \tilde{a}_{j_2}x_j\\ \vdots\\ \tilde{b}_m - \tilde{a}_{j_m}x_j \end{bmatrix}.$
现在要增加 $x_j$ ，那么只要保证 $x_B$ 的每个分量非负即可。

注意到 $x_0$ 是基本可行解，所以 $\tilde{b} = B^{-1}b \geq 0$ 。容易验证，下面的取值（称之为 Minimum Ratio Test）可以保证 $x_B \geq 0$ 。
$x_j := \min \left\{ \frac{\tilde{b}_i}{\tilde{a}_{j_i}} \text{ and } \tilde{a}_{j_i}>0, \quad i=1,2,\ldots, m\right\}.$
这样一来， $x_j \geq 0$ ，相应地，另一个变量 $x_{B_i} = \tilde{b}_i - \tilde{a}_{j_i} = 0$ ，其中 $i$ 是上式达到最小值的下标。换句话说，此时 $x_j$ 成为基变量（入基），而 $x_{B_i}$ 成为非基变量（出基）。可以证明，新得到的解仍然是一个基本可行解。

注意：如果 $\tilde{a}_{j_i} \leq 0$ ， $\forall i$ ，这意味着 $x_j$ 可以无限大，最优目标函数值为 $-\infty$ ，于是最优解不存在。

算法描述

第0步：输入基本可行解对应的基矩阵 $B$ 。

第1步：判断当前的解是否最优。如果 $\mu \geq 0$ ，当前是最优解，算法停止。

第2步：计算入基变量和出基变量。如果存在 $j\in J$ 使得 $\mu_j < 0$ ，那么 $x_j$ 是入基变量。执行 Minimum Ratio Test，找到出基变量 $x_i$ 。

第3步：判断问题是否无界。如果 $\tilde{a}_{j_i} \leq 0$ ， $\forall i$ ，则代表无界，算法停止。

第4步：执行出入基操作，更新基矩阵 $B$ ，然后执行第1步。

算法实现

下面我们用Python来实现单纯形算法。

先定义算法的输入和输出。

class SimplexA(object):
    """
    单纯形算法（基本版）。
    Note:
    	1、系数矩阵满秩。
        2、未处理退化情形。
        3、输入基本可行解（对应的列）。
    """
    def __init__(self, c, A, b, v0):
        """
        :param c: n * 1 vector
        :param A: m * n matrix
        :param b: m * 1 vector
        :param v0: basic variables, list of variable indices
        注意：v0是 B 的列下标。x0 = B^{-1}b 即为基本可行解（需要保证x0非负）。
        """
        # 输入
        self._c = np.array(c)
        self._A = np.array(A)
        self._b = np.array(b)
        self._basic_vars = v0  # basic variables
        self._m = len(A)
        self._n = len(c)
        self._non_basic_vars = self._init_non_basic_vars()  # non basic variables
        # 辅助变量
        self._iter_num = 0
        self._B_inv = None  # inverse of B
        self._lambda = None  # shadow price
        self._mu = None  # reduced cost
        # 输出
        self._obj = None  # objective function value
        self._sol = None  # solution
        self._status = None

接下来要实现单纯形算法 SimplexA.solve()，思路如下。

class SimplexA(object):
    
    # ...
    # 其它函数省略……
    
    def solve(self):
        self._iter_num = 0  # 记录迭代次数
        self._check_init_solution()  # 检查初始基本解是否可行
        self._update_reduced_cost()
        self._update_obj()
        self._update_solution()
        while not self._is_optimal():  # 判断是否最优或者无界
            if self._status == "UNBOUNDED":
                break
                self._pivot()  # 迭代（选主元入基，执行Minimum Ratio Test，然后出基）
                self._update_reduced_cost()  # 更新Reduced Cost: mu = c_N - lambda * N
                self._update_obj()
                self._update_solution()
                self._iter_num += 1
                if self._status != "UNBOUNDED":
                    self._status = 'OPTIMAL'
                    print("Done >> status: {}".format(self._status))
        return self

目标函数值和Reduced Cost的计算可以直接套公式。关键是实现每次迭代的出入基操作，即 SimplexA._pivot()。

class SimplexA(object):
    
    # ...
    # 其它函数省略……

    def _pivot(self):
        """ 选主元，入基和出基。
        """
        j_ind = np.argmin(self._mu)
        # 入基变量 x_j
        j = self._non_basic_vars[j_ind]
        # 出基变量 x_i
        i = self._minimum_ratio_test(j)
        if i is None:
            self._status = 'UNBOUNDED'
            return
        # update basic vars
        for k in range(self._m):
            if self._basic_vars[k] == i:
                self._basic_vars[k] = j
                break
        # update non basic vars
        self._non_basic_vars[j_ind] = i

    def _minimum_ratio_test(self, j):
        """ Minimum Ratio Test.
        给定入基的非基变量，返回出基的基变量。
        :param j: 入基变量 x_j 的下标 j
        :return: 出基变量 x_i 的下标 i
        """
        b_bar = np.dot(self._B_inv, self._b)
        a_in = np.dot(self._B_inv, self._A[:, j])
        ratios = list(map(lambda b, a: b/a if a > 1e-6 else np.infty, b_bar, a_in))
        i_ind = np.argmin(ratios)
        if ratios[i_ind] != np.infty:
            return self._basic_vars[i_ind]
        else:
            return None