在上一篇文章中，我们简要介绍了单纯形法的原理，但仍存在一个问题，就是没法获得初始的可行基。 这可以通过求解另一个线性规划来得到，这就是二阶段单纯形法。

1. 二阶段单纯形法

1.1. Phase I 寻找可行基

在上文中，不知道你有没有发现一个问题：第一步我们计算初始可行解 $\hat{x}_B=A_B^{-1}b$，但如果它并不满足 $\hat{x}_B\ge 0$ 呢？

我们将原始问题记为 $\eqref{eq:p1}$

\[\begin{aligned} \min_{x\in \mathbb{R}^n} &\quad c^{\top} x \\ \text { s.t. } &\quad Ax=b \\ &\quad x \ge 0 \end{aligned} \tag{P} \label{eq:p1}\]

那么可以通过添加人工变量(artificial variables)，将其转化为另一个辅助问题 $\eqref{eq:p2}$

\[\begin{aligned} \min_{x\in \mathbb{R}^n, w\in \mathbb{R}^m} &\quad e^{\top} w \\ \text { s.t. } &\quad Ax+w=b \\ &\quad x,w \ge 0 \end{aligned} \tag{P2} \label{eq:p2}\]

很明显，$(0,b)^{\top}$ 是 $\eqref{eq:p2}$ 的一个基本可行解. 并且如果原问题(P)的可行域 $\mathcal{F}_P\neq \emptyset$，那么(P2)的最优解一定能让 $w$ 取到0. 这时会出现两种情况：

  1. $w$ 中都是非基变量，此时得到的解 $x$ 就可作为(P)的基本可行解.
  2. $w$ 中存在基变量，这时需要将基变量出基，再从 $x$ 中寻找变量进基.

1.2. Phase II 代入求解

在得到了可行基之后，我们就能使用之前的修正单纯形法进行求解了。

1.3. 测试

对于优化问题

\[\begin{aligned} \min_{x} &\quad -x_{1}-x_{2} \\ \text { s.t. } &\quad x_{1}+2 x_{3}+x_{4} \leq 1 \\ &\quad x_{2}-x_{3}+x_{5} \leq 1 \\ &\quad x_{1}+x_{2}+x_{3} \leq 2 \\ &\quad x_{\nu} \geq 0 \quad(1 \leq \nu \leq 5) \end{aligned}\]

1.3.1. Phase I

我们先求解该问题获得可行基：

\[\begin{aligned} \min_{x} &\quad x_{6} \\ \text { s.t. } &\quad x_{1}+2 x_{3}+x_{4} \leq 1 \\ &\quad x_{2}-x_{3}+x_{5} \leq 1 \\ &\quad x_{1}+x_{2}+x_{3}+x_6 \leq 2 \\ &\quad x_{\nu} \geq 0 \quad(1 \leq \nu \leq 6) \end{aligned}\]

很显然，该问题有一个显然的可行基 $B=(4,5,6)$，解得可行基为 $\hat{B}=(1,2,5)$，最优解为 $\hat{x}=(1,1,0,0,1,0)^T$.

1.3.2. Phase 2

可行基中没有人工变量，于是原问题的可行基依然为 $B=(1,2,5)$，更换目标函数后继续代入单纯形法，得到最终解.

import numpy as np


class Data(object):
    pass


def LP_simplex(c, A, b, J):
    """
    solve linearly constrained optimization problems
        min  c^T x
        s.t. Ax=b
             x>=0
    B: feasible basis
    """
    (m, n) = A.shape
    Sn = np.arange(n)
    data = Data()
    Max_iter_num = 300
    for _ in range(Max_iter_num):
        K = np.setdiff1d(Sn, J)
        A_J = A[:, J]
        A_K = A[:, K]
        A_J_inv = np.linalg.inv(A_J)
        x = np.zeros(n)
        x[J] = A_J_inv@b
        data.x = x
        data.z = np.dot(c, x)
        data.J = J
        c_J = c[J]
        y = A_J_inv.T@c_J
        c_bar = c[K]-A_K.T@y
        S = np.where(c_bar < 0)
        S = S[0]
        if len(S) == 0:
            data.info = "Success"
            return data
        sigma = np.argmin(c_bar)
        s = K[sigma]
        A_s_bar = A_J_inv@A[:, s]
        d_J = -A_s_bar
        d_K = np.zeros(n-m)
        d_K[sigma] = 1
        Theta = x[J]/d_J
        Theta[Theta >= 0] = -1e+9
        rho = np.argmax(Theta)
        J[rho] = s
        J.sort()
    data.info = "Maximum number of iterations exceeded."
    return data


if __name__ == "__main__":
    A = np.array([[1, 0, 2, 1, 0], [0, 1, -1, 0, 1], [1, 1, 1, 0, 0]])
    b = np.array([1, 2, 2])
    c = np.array([-1, -1, 0, 0, 0])

    # Phase 1
    A1 = np.hstack((A, np.array([[0], [0], [1]])))
    c1 = np.array([0, 0, 0, 0, 0, 1])
    B1 = np.array([4, 5, 6])-1
    data = LP_simplex(c1, A1, b, B1)
    print(data.info)
    print(data.x)
    print(data.J)

    # Phase 2
    data = LP_simplex(c, A, b, data.J)
    print(data.info)
    print(data.x)
    print(data.J)

其结果如下：

Success
[1. 1. 0. 0. 1. 0.]
[0 1 4]
Success
[1. 1. 0. 0. 1.]
[0 1 4]

2. 有界线性规划

以下内容参考于 Optimization—Theory and Practice.

Let the following linear program be given:

\[\begin{aligned} &c^{T} x \longrightarrow \min \\ &A x=b, \quad \ell \leq x \leq u \end{aligned}\]

In the following we will use the notations from section $4.1$ (cf. p. $152 \mathrm{ff}$ ).

Show that this problem can be solved via a modification of the revised simplex method:

(1) Determine a feasible basic point $\bar{x}$ to a basis $J$ with

\[\begin{aligned} &\ell_{j} \leq \bar{x}_{j} \leq u_{j} \text { for } j \in J \text { and } \\ &\bar{x}_{j}=\ell_{j} \text { or } \bar{x}_{j}=u_{j} \text { for } j \in K . \end{aligned}\]

(2) Set \(\bar{c}^{T}:=c_{K}^{T}-y^{T} A_{K}\) with \(y:=A_{J}^{-T} c_{J}\). For all \(s=k_{\sigma} \in K\) let \(\bar{c}_{\sigma} \geq 0\) if \(\bar{x}_{s}=\ell_{s}\) and \(\bar{c}_{\sigma} \leq 0\) if \(\bar{x}_{s}=u_{s}:\) STOP

(3) Otherwise there exists an index \(s=k_{\sigma} \in K\) such that \(\bar{x}_{s}=\ell_{s}\) and \(\bar{c}_{\sigma}<0\) (Case A) or \(\bar{x}_{s}=u_{s}\) and \(\bar{c}_{\sigma}>0\) (Case B).

(4) Define \(d \in \mathbb{R}^{n}\) with \(d_{J}:=-\bar{a}_{s}, d_{K}:=e_{\sigma}\) and \(\bar{a}_{s}:=A_{J}^{-1} a_{s}\).

  • Case A: Determine the largest possible \(\tau \geq 0\) such that
\[\begin{aligned} \tau & \leq \frac{\bar{x}_{j_{i}}-u_{j_{i}}}{\bar{a}_{i, s}}, \text { if } \bar{a}_{i, s}<0 \\ \tau & \leq \frac{\bar{x}_{j_{i}}-\ell_{j_{i}}}{\bar{a}_{i, s}}, \text { if } \bar{a}_{i, s}>0 \end{aligned}\]

Assume that the maximum is attained for $i=\varrho .$ If $\tau \leq u_{s}-$ $\ell_{s}: x_{r}$ with $r=j_{\varrho}$ leaves the basis, $x_{s}$ enters. If $\tau>u_{s}-\ell_{s}$ : Set $\tau:=u_{s}-\ell_{s} ; J$ remains unchanged.

  • Case B: Determine the smallest possible $\tau \leq 0$ such that
\[\begin{aligned} &\tau \geq \frac{\bar{x}_{j_{i}}-\ell_{j_{i}}}{\bar{a}_{i, s}}, \text { if } \bar{a}_{i, s}<0 \\ &\tau \geq \frac{\bar{x}_{j_{i}}-u_{j_{i}}}{\bar{a}_{i, s}}, \text { if } \bar{a}_{i, s}>0 \end{aligned}\]

Assume that the minimum is reached for $i=\varrho .$ If $\tau \geq \ell_{s}-u_{s}$ : $x_{r}$ with $r=j_{\varrho}$ leaves the basis, $x_{s}$ enters. If $\tau<\ell_{s}-u_{s}:$ Set $\tau:=\ell_{s}-u_{s} ; J$ remains unchanged. $\bar{x}:=\bar{x}+\tau d ;$ update of $J$ as described above; goto (2).

3. 实验

Solve the following linear optimization problem

\[\begin{gathered} 2 x_{1}+x_{2}+3 x_{3}-2 x_{4}+10 x_{5} \quad \longrightarrow \min > \\ x_{1}+x_{3}-x_{4}+2 x_{5}=5 \\ x_{2}+2 x_{3}+2 x_{4}+x_{5}=9 \\ 0 \leq x \leq(7,10,1,5,3)^{T} . \end{gathered}\]

Firstly, verify that \(J:=(1,2), K:=(3,4,5)\) and \(\bar{x}_{K}:=(0,0,> 0)^{T}\) or \(\bar{x}_{K}:=(1,0,0)^{T}\) yield a feasible basic > solution.

废话就不多说了，直接摆代码：

import numpy as np


class Data(object):
    pass


def simplex_bound(c, A, b, L, U, J, x0):
    """
    solve linearly constrained optimization problems
        min  c^T x
        s.t. Ax=b
             L<=x<=U
    B: feasible basis
    x0: feasible basic solution
    """
    (m, n) = A.shape
    Sn = np.arange(n)
    data = Data()
    Max_iter_num = 300
    for _ in range(Max_iter_num):
        K = np.setdiff1d(Sn, J)
        data.x = x0
        data.z = np.dot(c, x0)
        data.J = J
        A_J = A[:, J]
        A_K = A[:, K]
        A_J_inv = np.linalg.inv(A_J)
        c_J = c[J]
        y = A_J_inv.T@c_J
        c_bar = c[K]-A_K.T@y
        Case = None
        for sigma in range(n-m):
            s = K[sigma]
            if c_bar[sigma] < 0 and x0[s] == L[s]:
                Case = 'A'
                break
            elif c_bar[sigma] > 0 and x0[s] == U[s]:
                Case = 'B'
                break
        # print(f'Case {Case}: ')
        A_s_bar = A_J_inv@A[:, s]
        d_J = -A_s_bar
        d_K = np.zeros(n-m)
        d_K[sigma] = 1
        if Case is None:
            data.info = "Success"
            return data
        elif Case == 'A':
            u = U[J]
            l = L[J]
            u[A_s_bar > 0] = l[A_s_bar > 0]
            Theta = (x0[J]-u)/A_s_bar
            rho = np.argmin(Theta)
            tau = Theta[rho]
            r = J[rho]
            if tau < (U[s]-L[s]):
                x0[J] += tau*d_J
                x0[K] += tau*d_K
                J[rho] = s
            else:
                tau = U[s]-L[s]
                x0[J] += tau*d_J
                x0[K] += tau*d_K
        elif Case == 'B':
            u = U[J]
            l = L[J]
            u[A_s_bar < 0] = l[A_s_bar < 0]
            Theta = (x0[J]-u)/A_s_bar
            rho = np.argmax(Theta)
            tau = Theta[rho]
            r = J[rho]
            if tau < (L[s]-U[s]):
                x0[J] += tau*d_J
                x0[K] += tau*d_K
                J[rho] = s
            else:
                tau = L[s]-U[s]
                x0[J] += tau*d_J
                x0[K] += tau*d_K
        J.sort()
    data.info = "Maximum number of iterations exceeded."
    return data


if __name__ == "__main__":
    A = np.array([[1, 0, 1, -1, 2], [0, 1, 2, 2, 1]])
    b = np.array([5, 9])
    c = np.array([2, 1, 3, -2, 10])
    J = np.array([0, 1])
    L = np.zeros(5)
    U = np.array([7, 10, 1, 5, 3])
    x0 = np.array([5.0, 9, 0, 0, 0])
    data = simplex_bound(c, A, b, L, U, J, x0)
    print(data.x)
    print(data.info)

输出结果为:

[7. 1. 1. 3. 0.]
Success

4. 参考资料

  1. W. Forst and D. Hoffmann, Optimization—Theory and Practice. New York, NY: Springer New York, 2010. https://doi.org/10.1007/978-0-387-78977-4.
  2. Math 340: Linear Programming - math340notes.pdf