线性规划算法 -凯发k8网页登录

线性规划算法

线性规划定义

线性规划问题求解向量 x，它在线性约束下最小化线性函数 f^tx：

$\min_{x} f^{t} x$

使得以下一个或多个条件成立：

a·x ≤ b
aeq·x = beq
l ≤ x ≤ u。

内点 `linprog` 算法

linprog 'interior-point' 算法与interior-point-convex quadprog 算法非常相似。它还有许多与 linprog 'interior-point-legacy' 算法相同的特征。这些算法具有相同的整体框架：

预求解，也就是将问题简化和转换为标准形式。
生成初始点。初始点的选择对于高效求解内点算法尤其重要，此步骤可能相当耗时。
使用预测-校正迭代求解 kkt 方程。此步骤通常最耗时。

预求解

该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容：

检查是否有任何变量具有相等的上界和下界。如果有，请检查可行性，然后固定变量的值以去除变量。
检查是否有任何线性不等式约束只涉及一个变量。如果是，请检查可行性，然后将线性约束更改为一个边界。
检查是否有任何线性等式约束只涉及一个变量。如果是，请检查可行性，然后固定该变量的值以去除该变量。
检查是否有任何线性约束矩阵具有全零行。如果是，请检查可行性，然后删除这些行。
确定边界和线性约束是否一致。
检查是否有任何变量仅作为线性项出现在目标函数中，而不出现在任何线性约束中。如果是这样，请检查可行性和有界性，然后将变量固定为其适当的边界值。
通过添加松弛变量，将任何线性不等式约束更改为线性等式约束。

如果算法检测到不可行或无界的问题，它会停止并发出相应的退出消息。

该算法可能到达一个表示解的单个可行点。

如果算法在预求解步骤中没有检测到不可行或无界的问题，并且预求解没有产生解，则算法会继续其后续步骤。到达停止条件后，算法会重新构造原始问题，这将撤消任何预求解变换。这最后一步是后求解步骤。

为简单起见，如果问题在预求解步骤中没有求得解，该算法会将所有有限下界移至零。

生成初始点

要设置初始点 x0，该算法执行以下操作。

将 x0 初始化为 ones(n,1)，其中 n 是目标函数向量 f 的元素数。
转换所有有界分量，使其下界为 0。如果分量 i 具有有限上界 u(i)，则 x0(i) = u/2。
对于只有一个边界的分量，如有必要，请修改分量，使其严格位于边界内。
要使 x0 接近中心路径，请执行一个预测-校正步骤，然后修改结果位置和松弛变量，使其完全处于任何边界内。有关中心路径的详细信息，请参阅 nocedal 和 wright 的论著 [7]（第 397 页）。

预测-校正

与 fmincon 内点算法相似，interior-point 算法尝试找到条件成立的点。要在线性规划问题中描述这些方程，请考虑使用经预处理的线性规划问题的标准形式：

$\min_{x} f^{t} x subject to {\begin{matrix} \bar{a} x = \bar{b} \\ x t = u \\ x, t \geq 0. \end{matrix}$

现在假设所有变量都有至少一个有限边界。此假设根据需要对分量进行平移和求反，使 x 的所有分量均有下界 0。
$\bar{a}$ 是同时包含线性不等式和线性等式的扩展线性矩阵。 $\bar{b}$ 是对应的线性等式向量。 $\bar{a}$ 还包括用于扩展向量 x 的项，这些项使用松弛变量 s 将不等式约束转换为等式约束：
$\bar{a} x = (\begin{matrix} a_{e q} & 0 \\ a & i \end{matrix}) (\begin{matrix} x_{0} \\ s \end{matrix}),$
其中，x₀ 表示原始 x 向量。
t 是将上界转换为等式的松弛变量组成的向量。

此方程组的拉格朗日函数包括以下向量：

y，与线性等式相关联的拉格朗日乘数
v，与下界相关联的拉格朗日乘数（正性约束）
w，与上界相关联的拉格朗日乘数

此拉格朗日函数写作

$l = f^{t} x - y^{t} (\bar{a} x - \bar{b}) - v^{t} x - w^{t} (u - x - t) .$

因此，此方程组的 kkt 条件是

$\begin{matrix} f - {\bar{a}}^{t} y - v w = 0 \\ \bar{a} x = \bar{b} \\ x t = u \\ v_{i} x_{i} = 0 \\ w_{i} t_{i} = 0 \\ (x, v, w, t) \geq 0. \end{matrix}$

linprog 算法对返回的拉格朗日乘数使用与本文讨论内容不同的符号约定。本讨论内容使用与大多数文献相同的符号。请参阅。

该算法首先根据 newton-raphson 公式计算一个预测步，然后计算一个校正步。校正器尝试减小非线性互补方程 s_iz_i = 0 中的残差。newton-raphson 步是

(\begin{matrix} 0 & - {\bar{a}}^{t} & 0 & - i & i \\ \bar{a} & 0 & 0 & 0 & 0 \\ - i & 0 & - i & 0 & 0 \\ v & 0 & 0 & x & 0 \\ 0 & 0 & w & 0 & t \end{matrix}) (\begin{matrix} δ x \\ δ y \\ δ t \\ \begin{array}{l} δ v \\ δ w \end{array} \end{matrix}) = - (\begin{matrix} f - {\bar{a}}^{t} y - v w \\ \bar{a} x - \bar{b} \\ u - x - t \\ \begin{array}{l} v x \\ w t \end{array} \end{matrix}) = - (\begin{matrix} r_{d} \\ r_{p} \\ r_{u b} \\ \begin{array}{l} r_{v x} \\ r_{w t} \end{array} \end{matrix}),

(1)

其中，x、v、w、t 是分别对应于向量 x、v、w 和 t 的对角矩阵。方程最右边的各个残差向量分别是：

r_d，对偶残差
r_p，原始残差
r_ub，上界残差
r_vx，下界互补残差
r_wt，上界互补残差

迭代输出将报告这些量：

$\begin{array}{l} primal infeasibility = {‖ r_{p} ‖}_{1} {‖ r_{u b} ‖}_{1} \\ dual infeasibility = {‖ r_{d} ‖}_{\infty} . \end{array}$

要求解公式 1，首先将其转换为对称矩阵形式

(\begin{matrix} - d & {\bar{a}}^{t} \\ \bar{a} & 0 \end{matrix}) (\begin{matrix} δ x \\ δ y \end{matrix}) = - (\begin{matrix} r \\ r_{p} \end{matrix}),

(2)

其中

$\begin{array}{l} d = x^{- 1} v t^{- 1} w \\ r = - r_{d} - x^{- 1} r_{v x} t^{- 1} r_{w t} t^{- 1} w r_{u b} . \end{array}$

d 和 r 定义中的所有矩阵的逆都很容易计算，因为矩阵是对角矩阵。

要从公式 1 推导公式 2，请注意公式 2 的第二行与公式 1 的第二个矩阵行相同。公式 2 的第一行来自于先求解公式 1 的最后两行以得到 δv 和 δw，再求解以得到 δt。

公式 2 是对称的，但由于 –d 项，它不是正定方程。因此，不能使用 cholesky 分解对其求解。多执行几步可得到另一个方程，它是正定方程，因此可以通过 cholesky 分解高效求解。

公式 2 的第二组行是

$\bar{a} δ x = - r_{p}$

第一组行是

$- d δ x {\bar{a}}^{t} δ y = - r .$

对下式进行代入

$δ x = d^{- 1} {\bar{a}}^{t} δ y d^{- 1} r$

可得

\bar{a} d^{- 1} {\bar{a}}^{t} δ y = - \bar{a} d^{- 1} r - r_{p} .

(3)

通常，找到牛顿步的最高效方法是使用 cholesky 分解求解公式 3 以获得 δy。之所以可以使用 cholesky 分解，是因为乘以 δy 的矩阵明显对称，并且在没有退化的情况下是正定矩阵。然后，要找到牛顿步，请回代以求得 δx、δt、δv 和 δw。但是，当 $\bar{a}$ 有稠密列时，求解公式 2 更为高效。linprog 内点算法根据列的稠密度选择求解算法。

有关更多算法细节，请参阅 mehrotra 的论著 [6]。

在计算校正的牛顿步后，算法执行更多计算，以获得更长的当前步，并准备好计算更优的后续步。这些多重校正计算可以提高性能和稳健性。关详细信息，请参阅 gondzio 的论著 [4]。

除了二次项，预测-校正算法与 quadprog 'interior-point-convex' 满算法基本相同。请参阅满预测-校正算法。

停止条件

预测-校正算法会持续迭代，直至达到一个可行（在容差范围内满足约束）且相对步长较小的点。具体来说，定义

$ρ = \max (1, ‖ \bar{a} ‖, ‖ f ‖, ‖ \bar{b} ‖) .$

当满足所有下列条件时，算法停止：

$\begin{matrix} {‖ r_{p} ‖}_{1} {‖ r_{u b} ‖}_{1} \leq ρ tolcon \\ {‖ r_{d} ‖}_{\infty} \leq ρ tolfun \\ r_{c} \leq tolfun, \end{matrix}$

其中

$r_{c} = \max_{i} (\min (| x_{i} v_{i} |, | x_{i} |, | v_{i} |), \min (| t_{i} w_{i} |, | t_{i} |, | w_{i} |)) .$

r_c 本质上是在度量互补残差 xv 和 tw 的大小，它们均为由零组成的向量，位于同一个解处。

内点传统线性规划

简介

内点传统方法基于 lipsol ()，它是 mehrotra 预测-校正算法 () 的变体，该算法是一种原始-对偶内点方法。

主算法

该算法首先应用一系列预处理步骤（请参阅预处理）。在预处理后，问题具有以下形式

\min_{x} f^{t} x such that {\begin{matrix} a \cdot x = b \\ 0 \leq x \leq u . \end{matrix}

(4)

上界约束隐式包含在约束矩阵 a 中。增加原始松弛变量 s 后，公式 4 变为

\min_{x} f^{t} x such that {\begin{matrix} a \cdot x = b \\ x s = u \\ x \geq 0, s \geq 0. \end{matrix}

(5)

以上为原始问题：x 由原始变量组成，s 由原始松弛变量组成。其对偶问题是

\max b^{t} y - u^{t} w such that {\begin{matrix} a^{t} \cdot y - w z = f \\ z \geq 0, w \geq 0, \end{matrix}

(6)

其中 y 和 w 由对偶变量组成， z 由对偶松弛变量组成。此线性规划的最优性条件，即原始公式 5 和对偶公式 6，是

\begin{array}{l} f (x, y, z, s, w) = (\begin{matrix} a \cdot x - b \\ x s - u \\ a^{t} \cdot y - w z - f \\ x_{i} z_{i} \\ s_{i} w_{i} \end{matrix}) = 0, \\ x \geq 0, z \geq 0, s \geq 0, w \geq 0, \end{array}

(7)

其中，x_iz_i 和 s_iw_i 表示按分量乘法。

linprog 算法对返回的拉格朗日乘数使用与本文讨论内容不同的符号约定。本讨论内容使用与大多数文献相同的符号。请参阅。

线性规划的二次方程 x_iz_i = 0 和 s_iw_i = 0 称为互补条件；其他（线性）方程称为可行性条件。量

x^tz s^tw

是对偶间隙，它度量当 (x,z,s,w) ≥ 0 时 f 的互补部分的残差。

该算法是一种原始-对偶算法，意味着同时求解原始规划和对偶规划。它可以被视为一种牛顿法，应用于公式 7 中的线性二次方程组 f(x,y,z,s,w) = 0，同时保持迭代 x、z、w 和 s 为正，因此称为内点方法。（迭代在由公式 5 中的不等式约束表示的严格内部区域中。）

该算法是 mehrotra 提出的预测-校正算法的变体。假设有迭代 v = [x;y;z;s;w]，其中 [x;z;s;w] > 0 首先计算所谓的预测方向

$δ v_{p} = - {(f^{t} (v))}^{- 1} f (v),$

这是牛顿方向；然后是所谓的校正方向

$δ v_{c} = - {(f^{t} (v))}^{- 1} f (v δ v_{p}) - μ \hat{e},$

其中，μ > 0 称为中心化参数，必须谨慎选择。 $\hat{e}$ 是由 0 和 1 组成的向量，其中 1 对应于 f(v) 中的二次方程，即扰动仅适用于互补条件，这些条件均为二次的；但不适用于可行性条件，这些条件均为线性的。这两个方向通过一个步长参数 α > 0 进行合并，并更新 v 以获得新迭代 v：

$v^{} = v α (δ v_{p} δ v_{c}),$

其中步长参数 α 的取值使得

v = [x;y;z;s;w]

满足

[x;z;s;w] >0。

在求解上述预测/校正方向的过程中，该算法根据修正的 cholesky 因子 a·a^t 来计算（稀疏）直接分解。如果 a 有稠密列，则算法改用 sherman-morrison 公式。如果该解不充分（残差太大），它将对步方程的增广方程组形式执行 ldl 分解来求得解。（请参阅matlab^® ldl 函数参考页中的。）

然后算法执行循环计算，直到迭代收敛。主终止条件是一个标准条件：

$\max (\frac{‖ r_{b} ‖}{\max (1, ‖ b ‖)}, \frac{‖ r_{f} ‖}{\max (1, ‖ f ‖)}, \frac{‖ r_{u} ‖}{\max (1, ‖ u ‖)}, \frac{| f^{t} x - b^{t} y u^{t} w |}{\max (1, | f^{t} x |, | b^{t} y - u^{t} w |)}) \leq t o l,$

其中

$\begin{matrix} r_{b} = a x - b \\ r_{f} = a^{t} y - w z - f \\ r_{u} = {x} s - u \end{matrix}$

分别是原始残差、对偶残差和上界可行性（{x} 是指具有有限上界的那些 x），并且

$f^{t} x - b^{t} y u^{t} w$

是原始目标值和对偶目标值之间的差，而 tol 是一定的容差。终止条件中的总和度量公式 7 中最优性条件中的总相对误差。

原始不可行性的度量是 ||r_b||，对偶不可行性的度量是 ||r_f||，其中范数是欧几里德范数。

预处理

该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容：

检查是否有任何变量具有相等的上界和下界。如果有，请检查可行性，然后固定变量的值以去除变量。
检查是否有任何线性不等式约束只涉及一个变量。如果是，请检查可行性，然后将线性约束更改为一个边界。
检查是否有任何线性等式约束只涉及一个变量。如果是，请检查可行性，然后固定该变量的值以去除该变量。
检查是否有任何线性约束矩阵具有全零行。如果是，请检查可行性，然后删除这些行。
确定边界和线性约束是否一致。
检查是否有任何变量仅作为线性项出现在目标函数中，而不出现在任何线性约束中。如果是这样，请检查可行性和有界性，然后将变量固定为其适当的边界值。
通过添加松弛变量，将任何线性不等式约束更改为线性等式约束。

如果算法检测到不可行或无界的问题，它会停止并发出相应的退出消息。

该算法可能到达一个表示解的单个可行点。

为简单起见，该算法将所有下界移至零。

虽然这些预处理步骤可以大大加快算法的迭代部分，但如果需要拉格朗日乘数，则必须撤消预处理步骤，因为在算法中计算的乘数适用于变换的问题，而不适用于原始问题。因此，如果不需要乘数，则不计算这种反向变换，因而可能节省一些计算时间。

对偶单纯形算法

就整体层面而言，linprog 'dual-simplex' 算法本质上是对对偶问题执行单纯形算法。

该算法首先执行预处理，如预处理中所述。有关详细信息，请参阅 andersen 和 andersen 的论著 [1] 以及 nocedal 和 wright 的论著 [7]，第 13 章。这种预处理将原始线性规划问题简化为公式 4 形式：

$\min_{x} f^{t} x such that {\begin{matrix} a \cdot x = b \\ 0 \leq x \leq u . \end{matrix}$

a 和 b 是原始约束矩阵的变换后的版本。这是原始问题。

原始可行性可以用函数来定义

$x^{} = {\begin{array}{l} x & if x > 0 \\ 0 & if x \leq 0. \end{array}$

原始不可行性的度量是

$primal infeasibility = \sqrt{{({(lb - x)}^{})}^{2} {({(x - ub)}^{})}^{2} {({(a x - b)}^{})}^{2} {| aeq x - beq |}^{2}} .$

正如公式 6 中所述，对偶问题是找到向量 y 和 w 以及松弛变量向量 z 来求解

$\max b^{t} y - u^{t} w such that {\begin{matrix} a^{t} \cdot y - w z = f \\ z \geq 0, w \geq 0. \end{matrix}$

linprog 算法对返回的拉格朗日乘数使用与本文讨论内容不同的符号约定。本讨论内容使用与大多数文献相同的符号。请参阅。

对偶不可行性的度量是

$dual infeasibility = {‖ a^{t} y z - w - f ‖}_{2} .$

众所周知（例如，请参阅 [7]），对偶问题的任何有限解都是原始问题的一个解，而原始问题的任何有限解都是对偶问题的一个解。而且，如果原始问题或对偶问题中有一个无界，则另一个问题不可行。并且，如果原始问题或对偶问题中有一个不可行，则另一个问题要么不可行，要么无界。因此，这两个问题在求有限解（如有）方面是等效的。由于原始问题和对偶问题在数学上是等效的，但计算步骤不同，因此有时更适合通过求解对偶问题来求解原始问题。

为了帮助减轻退化（请参阅 nocedal 和 wright 的论著 [7]，第 366 页），对偶单纯形算法从扰动目标函数开始。

对偶单纯形算法的阶段 1 是求得一个对偶可行点。该算法通过求解辅助线性规划问题来实现这一点。

在阶段 1 中，该算法通过求解辅助分段线性规划问题找到初始基可行解（请参阅基变量和非基变量了解定义）。辅助问题的目标函数是线性罚函数 $p = \sum_{j} p_{j} (x_{j})$ ，

其中 p_j(x_j) 的定义如下

$p_{j} (x_{j}) = {\begin{array}{l} x_{j} - u_{j} & if x_{j} > u_{j} \\ 0 & if l_{j} \leq x_{j} \leq u_{j} \\ l_{j} - x_{j} & if l_{j} > x_{j} . \end{array}$

p(x) 度量 x 点违反下界和上界条件的程度。辅助问题是

$\min_{x} \sum_{j} p_{j} subject to {\begin{matrix} a \cdot x \leq b \\ a e q \cdot x = b e q . \end{matrix}$

当且仅当辅助问题的最小值为 0 时，原始问题才有可行的基点。

算法通过启发式方法为辅助问题找到初始点，该方法根据需要添加松弛变量和人工变量。然后，该算法使用此初始点和单纯形算法来求解辅助问题。该解成为主算法阶段 2 的初始点。

在阶段 2 中，求解器反复选择一个入基变量和一个出基变量。该算法根据 forrest 和 goldfarb·[3] 提出的一种称为对偶最陡边定价的方法选择出基变量。该算法使用 koberstein·[5] 提出的 harris 比率检验变体来选择入基变量。为了帮助减轻退化，该算法可以在阶段 2 中引入额外的扰动。

求解器执行迭代，尝试在降低原始不可行性的同时保持对偶可行性，直到经过简化和扰动的问题的解既是原始可行的又是对偶可行的。该算法可消除它引入的扰动。如果（扰动后问题的）解对于未扰动（原始）问题是对偶不可行的，则求解器使用原始单纯形或阶段 1 算法还原对偶可行性。最后，求解器消除预处理步骤，以返回原始问题的解。

基变量和非基变量

本节定义线性规划问题的术语基、非基和基可行解。定义假设问题以如下标准形式给出：

$\min_{x} f^{t} x such that {\begin{matrix} a \cdot x = b, \\ l b \leq x \leq u b . \end{matrix}$

（请注意，a 和 b 不是原始问题中定义不等式的矩阵和向量。）假设 a 是 m×n 矩阵，秩 m < n，包含列 {a₁, a₂, ..., a_n}。假设 ${a_{i_{1}}, a_{i_{2}}, ..., a_{i_{m}}}$ 是 a 的列空间的基，索引集 b = {i₁、i₂、...、i_m}，n = {1, 2, ..., n}\b 是 b 的互补矩阵。子矩阵 a_b 称为基，互补子矩阵 a_n 称为非基。由基变量组成的向量是 x_b，由非基变量组成的向量是 x_n。在阶段 2 的每次迭代中，该算法用非基的一列替换当前基的一列，并相应地更新变量 x_b 和 x_n。

如果 x 是 a·x = b 的解，并且 x_n 中的所有非基本变量都等于其下界或上界，则 x 称为基本解。此外，如果 x_b 中的基变量满足其下界和上界，使得 x 是可行点，则 x 称为基可行解。

参考

[1] andersen, e. d., and k. d. andersen. presolving in linear programming. math. programming 71, 1995, pp. 221–245.

[2] applegate, d. l., r. e. bixby, v. chvátal and w. j. cook, the traveling salesman problem: a computational study, princeton university press, 2007.

[3] forrest, j. j., and d. goldfarb. steepest-edge simplex algorithms for linear programming. math. programming 57, 1992, pp. 341–374.

[4] gondzio, j. “multiple centrality corrections in a primal dual method for linear programming.” computational optimization and applications, volume 6, number 2, 1996, pp. 137–156. available at .

[5] koberstein, a. progress in the dual simplex algorithm for solving large scale lp problems: techniques for a fast and stable implementation. computational optim. and application 41, 2008, pp. 185–204.

[6] mehrotra, s. “on the implementation of a primal-dual interior point method.” siam journal on optimization, vol. 2, 1992, pp 575–601.

[7] nocedal, j., and s. j. wright. numerical optimization, second edition. springer series in operations research, springer-verlag, 2006.

线性规划算法 -凯发k8网页登录