【病毒分析】Babyk加密器分析

单纯形法、大M法，及相关线性代数知识

单纯形法是求解一般线性规划问题的迭代算法。其基本思想是从一个初始的基可行解出发，通过迭代逐步逼近最优解。每一步迭代都选择一个新的基可行解，使得目标函数值逐步改善，直到达到最优解。

因为看其他资料可知，线性规划问题可以约定其他的标准形式，所以本文和《运筹学基础及应用（第六版）》一样，约定线性规划问题的标准形式为

 复制代码 隐藏代码
max CX
AX = b
X >= 0
C: 1 * n, X: n * 1, A: m * n, b: m * 1, 一般认为 m <= n

在我的代码中，C叫obj_func_coeff，A叫constraints_coeff，b叫the RHS of the constraints。

预备知识

• 凸集：如果一个集合中的任意两点的连线都完全包含在该集合内，那么这个集合就是一个凸集。换句话说，对于集合 ( C ) 中的任意两点x和y，以及任意\lambda满足0 \leq \lambda \leq 1，点\lambda x + (1 - \lambda) y也在集合 ( C ) 中。
• 顶点：在多面体（例如多边形或多面体）中，顶点是多面体的角点。对于线性规划问题，顶点是可行区域的边界上的点，这些点是由约束条件的交点形成的。代数定义：任意X1, X2 ∈ C，不存在X = a * X1 + (1-a) * X2, 0 < a < 1。

以上概念主要出现在讨论线性规划问题的可行域时。

• 基：在线性代数中，基是一个向量空间中的一组向量，这些向量线性无关且能够生成整个向量空间。在线性规划中，基是指约束条件矩阵中的一组线性无关的列向量。为什么基很重要？因为单位向量组是最简单的一组基，通过单位向量组我们立刻可以得到初始基可行解，从而开始运行单纯形法。值得注意的是，（1）m个m维向量线性无关/每个向量都不能由其他m-1个向量线性表出。（2）它们按列或按行组成的矩阵可逆/行列式不为0/满秩。都是对同一个事物的不同角度的描述。我们在系数矩阵的列向量组中选出一组基，则这些列向量对应的变量叫做基变量，其他变量叫做非基变量。
• 基解：令n - m个非基变量都为0，根据上一段的（1）和（2）（或者按《运筹学基础及应用（第六版）》的说法“根据克拉默规则”），基变量可得唯一解。于是一组基对应一个基解。当选择的基为单位向量组时，我们立刻可以得到初始基可行解，从而开始运行单纯形法。
• 可行解：满足所有约束条件的解就是可行解。所有可行解组成的集合叫可行域。
• 基可行解：基解和可行解的交集。
• 可行基：可行基是指对应于基可行解的基。

单纯形法的计算步骤

在单纯形法的计算步骤中，主要是需要计算检验数（Relative Profits）和theta，它们的公式的推导过程分别在《运筹学基础及应用（第六版）》第一章 3-5 最优性检验和解的判别和 3-4 从初始基可行解转换为另一基可行解中有详细说明。不理解也没事，把式子背下来就足够应付考试了。计算步骤总结如下：

1. 标准化问题：将线性规划问题转化为标准形式，即目标函数为最大化形式，约束条件为等式形式，且所有变量非负。这一步需要加上若干松弛变量（Slack Variables）。
2. 构造初始单纯形表：找到一个初始的基可行解，并构造相应的单纯形表。在大M法中再具体讨论。
3. 判断最优解：如果单纯形表中所有检验数都≤0，则当前单纯形表中的基可行解就是最优解；否则继续迭代。
4. 迭代计算：从一个基可行解转换到另一个更优秀的基可行解。（1）确定入基变量。首先算n个检验数：sigma[j] = C[j] - C_B * A[:, j]。显然基变量对应的检验数都是0。每个大于0的检验数都可以作为入基变量，但是我们一般会贪心地选择检验数最大的变量，因为这样似乎可以更快收敛。（2）确定出基变量。假设上一步确定入基变量是xk，那么theta = min(bi / A[i, k], A[i, k] > 0)。最小theta对应的每个i都可以作为出基变量。假设其中一个是p，A[p, k]就是下一步初等行变换的主元，我们会在单纯形表中把主元框起来。（3）更新单纯形表。我们要进行初等行变换，把主元列变成单位列向量。

文字说明看不懂也没事，可以结合下面的例题理解。

例题：《运筹学基础及应用（第六版）》第一章例5

 复制代码 隐藏代码
max z = 2*x1 + 3*x2
2*x1 + 2*x2 <= 12
4*x1 <= 16
5*x2 <= 15
x1, x2 >= 0

懒得画精美的markdown表格了，直接展示我代码输出的单纯形表吧！

• CB：基变量xj对应的C[j]，写出来是为了方便算内积，求出检验数。
• b：当前基可行解。
• σ：检验数。
• 被框起来的元素：初等行变换的主元。所在的行对应出基变量，所在的列对应入基变量。

 复制代码 隐藏代码
|    | cj    |    | 2 | 3   | 0 | 0 | 0 |
| CB | Basis | b  | x1| x2  | x3| x4| x5|
| 0  | x3    | 12 | 2 | 2   | 1 | 0 | 0 |
| 0  | x4    | 16 | 4 | 0   | 0 | 1 | 0 |
| 0  | x5    | 15 | 0 | [5] | 0 | 0 | 1 |
|    | sigma |    | 2 | 3   | 0 | 0 | 0 |

|    | cj    |    | 2   | 3 | 0 | 0 | 0      |
| CB | Basis | b  | x1  | x2| x3| x4| x5     |
| 0  | x3    | 6  | [2] | 0 | 1 | 0 | -2 / 5 |
| 0  | x4    | 16 | 4   | 0 | 0 | 1 | 0      |
| 3  | x2    | 3  | 0   | 1 | 0 | 0 | 1 / 5  |
|    | sigma |    | 2   | 0 | 0 | 0 | -3 / 5 |

|    | cj    |   | 2 | 3 | 0     | 0 | 0      |
| CB | Basis | b | x1| x2| x3    | x4| x5     |
| 2  | x1    | 3 | 1 | 0 | 1 / 2 | 0 | -1 / 5 |
| 0  | x4    | 4 | 0 | 0 | -2    | 1 | 0.8    |
| 3  | x2    | 3 | 0 | 1 | 0     | 0 | 0.2    |
|    | sigma |   | 0 | 0 | -1    | 0 | -1 / 5 |

步骤：

1. 选择x2为入基变量，计算theta = min(12 / 2, -, 15 / 5)确定x5为出基变量。初等行变换：(1) -= 2/5*(3), (2) = (2), (3) /= 5。
2. 选择x1为入基变量，计算theta = min(6 / 2, 16 / 4, -)确定x3为出基变量。初等行变换：(1) /= 2, (2) -= 4/2*(1), (3) = (3)。
3. 所有检验数都小于0，没有等于0的，所以有唯一最优解x1 = 3, x2 = 3, x4 = 3, x3 = x5 = 0，目标函数值15。

解的判别

解有几种情况：

1. 有唯一最优解
2. 有无穷多最优解
3. 有无界解
4. 无解

如何判断无穷多最优解：我们知道求得最优解时，所有变量的检验数都小于等于0。如果存在一个非基变量的检验数等于0，并且选择它为入基变量时可以正常求出theta，那么就能正常求出出基变量+进行初等行变换，于是就找到了另一个顶点也是最优解。于是这两个顶点连线上的所有点都是最优解。

如何判断无界解：求出基变量时，如果发现系数矩阵的每一列都小于等于0，导致theta = min(-, -, ...)，那么线性规划问题具有无界解。

人工变量法（大M法）

化为标准形式后，如果系数矩阵中不存在单位矩阵，那么我们需要先写代码找到一组基，找到基后得到的基解还不一定是可行解。所以我们有一个很朴素的想法，就是人工地添加所有缺失的单位列向量。为此，对应新增的人工变量（Artificial Variables）的系数要设为负无穷大，用-M表示，以保证解出的人工变量取值为0。建立初始单纯形表后，正常运行单纯形法求解即可。

比如《运筹学基础及应用（第六版）》第一章例6

 复制代码 隐藏代码
max z = -3 * x1 + x3
x1 + x2 + x3 + x4 = 4
-2 * x1 + x2 - x3 - x5 = 1
3 * x2 + x3 = 9
x1~x5 >= 0

其列向量

 复制代码 隐藏代码
P1 P2 P3 P4 P5
1  1  1  1  0
-2 1  -1 0  -1
0  3  1  0  0

可见只有一个单位列向量P4。那么我们再加两个单位列向量P6, P7

 复制代码 隐藏代码
P1 P2 P3 P4 P5 P6 P7
1  1  1  1  0  0  0
-2 1  -1 0  -1 1  0
0  3  1  0  0  0  1

目标函数相应地变为max z = -3 * x1 + x3 - M * x6 - M * x7。

核心代码讲解

项目传送门。目录结构：

 复制代码 隐藏代码
SIMPLEX-BIG-M
│  .gitignore
│  pytest.ini
│  README.md
│
├─outp
│      big_m_method.txt
│      main.txt
│      simplex.txt
│
├─src
│  │  big_m.py
│  │  big_m_method.py
│  │  consts.py
│  │  main.py
│  │  simplex.py
│  │  utils.py
│  └─ __init__.py
│
└─test
    │  test_big_m.py
    │  test_big_m_method.py
    │  test_simplex.py
    │  test_utils.py
    └─ __init__.py

开发过程主要分两阶段。

第一阶段：实现BigM类和Simplex类

BigM类就像复数一样，支持各种符号运算及与浮点数的混合运算。Simplex类支持输入一个列向量组含有所有单位向量的矩阵（行数m≤列数n），首先求出初始基可行解，然后实现单纯形法的各步。Simplex类应支持目标函数系数含有BigM类实例。

BigM类构造函数：

 复制代码 隐藏代码
class BigM:
    def __init__(self, a, b) -> None:
        self.a = a
        self.b = b

BigM类要支持的运算符不少，所以代码量略大，完整代码。为了减少代码量，一些运算符可以由其他运算符实现。比如只实现＜和=即可实现其他4个比较运算符≤、＞、≥、≠。又比如__radd__, __rmul__可以由__add__, __mul__实现。接下来看一些运算符的实现。＜：

 复制代码 隐藏代码
    def __lt__(self, other: object) -> bool:
        if isinstance(other, (int, float)):
            if self.a < 0:
                return True
            if self.a == 0:
                return self.b < other
            return False
        if not is_big_m_like(other):
            raise ValueError(f'Expect a BigM instance, but got type {type(other)}')
        if self.a < other.a:
            return True
        elif self.a == other.a:
            return self.b < other.b
        return False

这里用is_big_m_like(other)而非isinstance(other, BigM)的原因是单测文件导入的BigM位于sys.modules['src/big_m']，而src文件夹下其他文件导入的BigM位于sys.modules['big_m']，两者不一样，于是我写了is_big_m_like(other)暂时规避了这个问题。

乘法：可类比复数类的实现。

 复制代码 隐藏代码
    def __mul__(self, other: object):
        if isinstance(other, (int, float)):
            return BigM(self.a * other, self.b * other)
        if not is_big_m_like(other):
            raise ValueError(f'Expect a BigM instance, but got type {type(other)}')
        return BigM(self.a * other.a + self.b * other.a + self.a * other.b, self.b * other.b)

完整代码，BigM类单测。

接下来看最难的部分——单纯形法的实现。Simplex类完整代码，Simplex类单测。

 复制代码 隐藏代码
class Simplex:
    def __init__(
            self,
            obj_func_coeff: np.ndarray,
            constraints_coeff: np.ndarray,
            b: np.ndarray,
            should_print_table: bool = False,
            is_debug: bool = False) -> None:
        self.obj_func_coeff = obj_func_coeff  # 目标函数向量
        # 保证外部数组不会被内部代码修改
        self.tmp_constraints_coeff = constraints_coeff.copy()  # 会被修改的约束系数矩阵
        self.constraints_coeff = constraints_coeff.copy()  # 约束系数矩阵
        self.b = b  # 约束右侧向量，the RHS of the constraints
        self.should_print_table = should_print_table  # 是否打印各步的单纯形表
        self.is_debug = is_debug  # 单纯为了方便测试
        self.n = len(self.obj_func_coeff)
        self.m = len(self.constraints_coeff)
        self._check_problem_shape()

通过单位列向量求初始基解：

 复制代码 隐藏代码
    def _get_initial_basis(self) -> List[int]:
        m = self.m
        basis = [0] * m
        unit_vector_to_idx = get_unit_vector_to_idx(m)
        for i in range(self.n):
            col = tuple(self.tmp_constraints_coeff[:, i])
            if col not in unit_vector_to_idx:
                continue
            idx = unit_vector_to_idx[col]
            del unit_vector_to_idx[col]
            basis[idx] = i
        if unit_vector_to_idx:
            raise ValueError('Get initial basis failed')
        return basis

    def solve(self):
        m, n = self.m, self.n
        self.tmp_constraints_coeff = self.constraints_coeff.copy()
        basis = self._get_initial_basis()
        basis_sol = self.b.tolist()
        cb = [self.obj_func_coeff[v] for v in basis]

求检验数和入基变量：求内积即可。

 复制代码 隐藏代码
            relative_profits = [self.obj_func_coeff[i] - np.dot(self.tmp_constraints_coeff[:, i], cb) for i in range(n)]
            enter_basis = np.argmax(relative_profits)

成功条件：if relative_profits[enter_basis] <= 0。

求theta和出基变量：