线性规划：单纯形算法之初始化

最新推荐文章于 2022-09-07 10:22:55 发布

胡拉哥

最新推荐文章于 2022-09-07 10:22:55 发布

阅读量553

点赞数

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

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

版权

优化算法专栏收录该内容

21 篇文章 57 订阅

订阅专栏

考虑线性规划的标准形式：
$\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}$ 。

回顾单纯形算法：从一个基本可行解开始，沿着目标函数下降的方向，迭代到新的基本可行解，直到目标函数达到最优。

那么如何找到一个初始的基本可行解？本文介绍一个方法：两阶段法。

两阶段法

两阶段法的思路如下：

1、构造辅助变量，把原问题转化成“新问题”，使得新问题的基本可行解是显而易见的。

2、用单纯形算法求解新问题，从而得到原问题的一个基本可行解。

3、再次使用单纯形算法，求解原问题。

下面介绍新问题的定义。定义辅助变量 $\in \mathbb{R}^m$ ， $\ldots, 1)^T \in \mathbb{R}^n$ ， $I_m$ 是 $m$ 维单位矩阵。

考虑如下问题（下文称它为新问题）：
$\begin{aligned} \min~ & e^T y\\ \text{s.t.}~ & Ax + I_m y =b\\ & x, y\geq 0 \end{aligned}$
原问题有可行解等价于新问题的最优目标函数值为零（证明略）。

注意到 $I_m$ 是新问题的基，用单纯形算法计算新问题，得到最优解 $x^*, y^*$ 。如果原问题可行，那么 $y^*=\mathbf{0}$ ，意味着 $A x^* = b$ 。换句话说， $x^*$ 是原问题的可行解。

不过问题来了， $x^*$ 虽然是原问题的可行解，但是不能直接当作单纯形算法的输入，因为新问题和原问题的基变量不一定相同。接下来我们讨论如何计算原问题的基变量。

基变量

为了简化记号，把新问题的最优解记作 $z$ ，即
$(x_1,\ldots,x_n,y_1,\ldots,y_m)^T.$
令 $z_B, z_N$ 分别代表基变量和非基变量。于是 $z$ 也可以写成
$\begin{bmatrix} z_B\\ z_N \end{bmatrix}.$
考虑到 $z_B, z_N$ 有可能包含 $x$ 和 $y$ 中的分量，令
$z_B = \begin{bmatrix} x_B\\ y_B \end{bmatrix}, \quad z_N = \begin{bmatrix} x_N\\ y_N \end{bmatrix},$
其中 $x_B, y_B$ 是 $x, y$ 对应的基变量， $x_N, y_N$ 是 $x, y$ 对应的非基变量。

于是， $z$ 可以进一步写成
$\begin{bmatrix} z_B \\ z_N \end{bmatrix} = \begin{bmatrix} x_B\\ y_B\\ x_N\\ y_N \end{bmatrix}.$
接下来考虑两种情形：

情形1：基变量 $z_B$ 不包含 $y$ 。
$\begin{bmatrix} x_B\\ x_N\\ y_N \end{bmatrix}$
此时，原问题的基变量为 $x_B$ 。

情形2：基变量 $z_B$ 包含 $y$ 的分量。
$\begin{bmatrix} x_B\\ y_B\\ x_N\\ y_N \end{bmatrix}$
这种情况下，不能直接把 $x_B, x_N)$ 作为原问题的初始基本可行解，因为 $x_B$ 的维度小于 $m$ 。

能不能把 $y_B$ 与 $x_N$ 对应的列交换？

答案是可以的。因为 $x_N = \mathbf{0}$ ，选择某一列入基不会降低目标函数值。

具体来说，交换行和列的顺序，先把最优解和对应的系数矩阵写成如下形式。

交换的思路是把 $I_{m-k}$ 与 $R_2$ 中非零元素所在的行列，按列进行交换。

考虑矩阵 $R_2$ 的第 $i$ 行，找到一个非零元所在的列 $j$ ，然后交换 $y_{B_i}$ 和 $x_{N_j}$ 。如果 $j$ 不存在，那么 $R_2$ 的第 $i$ 行是 $\mathbf{0}$ ，这意味着 $y_{B_i}$ 对应的约束是冗余的。

从 $i = 0$ 开始，执行上述步骤，最终得到的结果是：

（1） $R_2$ 被转换成 $I_{m-k}$ 或者

（2）执行到某一步时 $R_2$ 中的元素全为0导致无法继续交换。

第一种情况意味着 $y_B$ 全部被替换，此时 $z_B = x_B$ ，即原问题的基变量。

第二种情况意味着原问题的系数矩阵不满秩：
$\sim \begin{bmatrix} I_k & R_1 \\ \mathbf{0} & \mathbf{0} \end{bmatrix}.$
最后 $m - k$ 行是冗余的。删除它们，然后把 $x_B$ 作为原问题的基变量。

算法实现

import numpy as np

from simplex_ad import SimplexAD  # 之前实现的单纯形算法


class Simplex2P(object):

    def __init__(self, c, A, b):
        self._c = c
        self._A = A
        self._b = b
        self._m = len(self._A)
        self._n = len(self._c)
        self._status = None
        # The following variables are for phase 1
        self._ins1 = self._construct_phase1_instance()
        self._basic_vars = []
        self._basic_art_vars = []
        self._basic_art_vars_num = 0
        # 处理退化情形的中间变量
        # R = B^{-1}*N.
        # R12 对应 R 的前 k 列，其中 k 代表 basic non artificial variable 的个数
        self._R12 = None

    def _construct_phase1_instance(self):
        c1 = [0] * self._n + [1] * self._m
        I = np.eye(self._m)
        A = self._A
        A1 = np.array(np.bmat("A I"))
        b1 = self._b
        v0 = [i + self._n for i in range(self._m)]
        return c1, A1, b1, v0

    def _initialize(self):
        # solve phase1 problem
        sim = SimplexAD(*self._ins1).solve(print_info=False)

        # check feasibility
        if abs(sim.get_objective()) > 1e-6 or sim.get_status() != 'OPTIMAL':
            self._status = 'INFEASIBLE'
            return

        # If feasible, consider the following two cases.
        # Case1(normal): artificial vars not in basic vars.
        self._basic_vars = sim.get_basic_vars()

        # Case2(degenerate): artificial vars in basic vars.
        # Then exchange artificial columns with non-basic columns.
        self._basic_art_vars = [i for i in self._basic_vars if i >= self._n]
        self._basic_art_vars_num = len(self._basic_art_vars)
        if self._basic_art_vars_num > 0:
            self._resolve_degeneracy()

    def _swapping_columns_and_rows(self, B_inv_N):
        basic_indices = list(zip(range(self._m), self._basic_vars))
        basic_non_art_indices = [item for item in basic_indices if item[1] < self._n]
        basic_art_indices = [item for item in basic_indices if item[1] >= self._n]

        # column swapping
        # reorder basic indices w.r.t. non basic
        basic_indices2 = basic_non_art_indices + basic_art_indices
        self._basic_vars = [x[1] for x in basic_indices2]
        # row swapping
        B_inv_N1 = np.zeros(B_inv_N.shape)
        ind2 = [x[0] for x in basic_indices2]
        for i in range(self._m):
            B_inv_N1[i, :] = B_inv_N[ind2[i], :]

        return B_inv_N1

    def _format_R12(self):
        B = self._get_submatrix_of_A1(self._basic_vars)
        non_basic_vars = list(set(range(self._n + self._m)) - set(self._basic_vars))
        N = self._get_submatrix_of_A1(non_basic_vars)
        B_inv_N = np.linalg.inv(B) @ N
        B_inv_N1 = self._swapping_columns_and_rows(B_inv_N)
        # The number of non basic, non artificial variables
        k = self._n - (self._m - self._basic_art_vars_num)
        self._R12 = B_inv_N1[:, :k]

    def _resolve_degeneracy(self):
        self._format_R12()
        # replace artificial basic columns with non-basic columns.
        self._replace_artificial_basis()

        # remove redundant rows if there exist
        rr = self._get_redundant_rows()
        if len(rr) > 0:
            self._remove_redundant_rows(rr)
            # remove redundant basis w.r.t. redundant rows.
            reserved_indices = set(range(len(self._basic_vars))) - set(rr)
            self._basic_vars = [self._basic_vars[i] for i in reserved_indices]

    def _remove_redundant_rows(self, redundant_rows):
        # remove from original instance
        reserved_rows = set(range(self._m)) - set(redundant_rows)
        self._A = np.array([self._A[i] for i in reserved_rows])
        self._b = np.array([self._b[i] for i in reserved_rows])
        self._m = len(self._A)
        # remove from ins1
        c1, A1, b1, v0 = self._ins1
        A1 = np.array([A1[i] for i in reserved_rows])
        b1 = np.array([b1[i] for i in reserved_rows])
        self._ins1 = (c1, A1, b1, v0)

    def _get_redundant_rows(self):
        # return zero-vector indices w.r.t. R2.
        k = len(self._basic_vars) - self._basic_art_vars_num
        return [i for i in range(k, self._m) if sum(np.absolute(self._R12[i])) < 1e-6]

    def _get_submatrix_of_A1(self, cols):
        """ Return the sub matrix of A1 w.r.t. column indices.
        """
        mat = []
        A = self._ins1[1]
        for j in cols:
            mat.append(A[:, j])
        return np.array(mat).transpose()

    def _replace_artificial_basis(self):

        for i in range(self._basic_art_vars_num):
            # "in" variable from non basic variables
            in_var = self._get_in_var(self._basic_art_vars_num, i)
            if in_var is None:
                continue
            # "out" variable from basic variables (which is artificial)
            out_var = self._basic_art_vars[i]

            # replacement
            for j in range(self._m):
                if self._basic_vars[j] == out_var:
                    self._basic_vars[j] = in_var
                    break

            for j in range(self._basic_art_vars_num):
                if self._basic_art_vars[j] == out_var:
                    row_ind = self._m - self._basic_art_vars_num + j
                    self._R12[:, in_var] = np.array([1 if i == row_ind else 0 for i in range(self._m)])
                    break

    def _get_in_var(self, basic_art_vars_num, art_ind):
        """
        :param basic_art_vars_num: the total number of "basic" artificial variables
        :param art_ind: the index of the artificial variable in "basic_art_vars"
        :return: "in" variable from non basic variables
        """
        row_id = self._m - basic_art_vars_num + art_ind
        k = len(self._R12[row_id])
        for j in range(k):
            if abs(self._R12[row_id][j]) > 1e-6:
                return j
        return None

    def solve(self):
        self._initialize()
        if self._status == 'INFEASIBLE':
            print('>> INFEASIBLE.')
        else:
            print("Init: feasible.")
            sim = SimplexAD(self._c, self._A, self._b, self._basic_vars).solve()
            self._status = sim.get_status()

        return self