Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

线性规划算法

线性规划定义

线性规划问题求解向量 x,它在线性约束下最小化线性函数 fTx:

minxfTx

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

A·x ≤ b
Aeq·x = beq
l ≤ x ≤ u。

内点 linprog 算法

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

  1. 预求解,也就是将问题简化和转换为标准形式。

  2. 生成初始点。初始点的选择对于高效求解内点算法尤其重要,此步骤可能相当耗时。

  3. 使用预测-校正迭代求解 KKT 方程。此步骤通常最耗时。

预求解

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

  • 检查是否有任何变量具有相等的上界和下界。如果有,请检查可行性,然后固定变量的值以去除变量。

  • 检查是否有任何线性不等式约束只涉及一个变量。如果是,请检查可行性,然后将线性约束更改为一个边界。

  • 检查是否有任何线性等式约束只涉及一个变量。如果是,请检查可行性,然后固定该变量的值以去除该变量。

  • 检查是否有任何线性约束矩阵具有全零行。如果是,请检查可行性,然后删除这些行。

  • 确定边界和线性约束是否一致。

  • 检查是否有任何变量仅作为线性项出现在目标函数中,而不出现在任何线性约束中。如果是这样,请检查可行性和有界性,然后将变量固定为其适当的边界值。

  • 通过添加松弛变量,将任何线性不等式约束更改为线性等式约束。

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

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

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

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

生成初始点

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

  1. x0 初始化为 ones(n,1),其中 n 是目标函数向量 f 的元素数。

  2. 转换所有有界分量,使其下界为 0。如果分量 i 具有有限上界 u(i),则 x0(i) = u/2

  3. 对于只有一个边界的分量,如有必要,请修改分量,使其严格位于边界内。

  4. 要使 x0 接近中心路径,请执行一个预测-校正步骤,然后修改结果位置和松弛变量,使其完全处于任何边界内。有关中心路径的详细信息,请参阅 Nocedal 和 Wright 的论著 [7](第 397 页)。

预测-校正

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

minxfTx subject to {A¯x=b¯x+t=ux,t0.

  • 现在假设所有变量都有至少一个有限边界。此假设根据需要对分量进行平移和求反,使 x 的所有分量均有下界 0。

  • A¯ 是同时包含线性不等式和线性等式的扩展线性矩阵。b¯ 是对应的线性等式向量。A¯ 还包括用于扩展向量 x 的项,这些项使用松弛变量 s 将不等式约束转换为等式约束:

    A¯x=(Aeq0AI)(x0s),

    其中,x0 表示原始 x 向量。

  • t 是将上界转换为等式的松弛变量组成的向量。

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

  • y,与线性等式相关联的拉格朗日乘数

  • v,与下界相关联的拉格朗日乘数(正性约束)

  • w,与上界相关联的拉格朗日乘数

此拉格朗日函数写作

L=fTxyT(A¯xb¯)vTxwT(uxt).

因此,此方程组的 KKT 条件是

fA¯Tyv+w=0A¯x=b¯x+t=uvixi=0witi=0(x,v,w,t)0.

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

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

(0A¯T0IIA¯0000I0I00V00X000W0T)(ΔxΔyΔtΔvΔw)=(fA¯Tyv+wA¯xb¯uxtVXWT)=(rdrprubrvxrwt),(1)

其中,X、V、W、T 是分别对应于向量 x、v、w 和 t 的对角矩阵。方程最右边的各个残差向量分别是:

  • rd,对偶残差

  • rp,原始残差

  • rub,上界残差

  • rvx,下界互补残差

  • rwt,上界互补残差

迭代输出将报告这些量:

Primal infeasibility=rp1+rub1Dual infeasibility=rd.

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

(DA¯TA¯0)(ΔxΔy)=(Rrp),(2)

其中

D=X1V+T1WR=rdX1rvx+T1rwt+T1Wrub.

D 和 R 定义中的所有矩阵的逆都很容易计算,因为矩阵是对角矩阵。

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

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

公式 2 的第二组行是

A¯Δx=rp

第一组行是

DΔx+A¯TΔy=R.

对下式进行代入

Δx=D1A¯TΔy+D1R

可得

A¯D1A¯TΔy=A¯D1Rrp.(3)

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

有关更多算法细节,请参阅 Mehrotra 的论著 [6]

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

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

停止条件

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

ρ=max(1,A¯,f,b¯).

当满足所有下列条件时,算法停止:

rp1+rub1ρTolConrdρTolFunrcTolFun,

其中

rc=maxi(min(|xivi|,|xi|,|vi|),min(|tiwi|,|ti|,|wi|)).

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

内点传统线性规划

简介

默认内点传统方法基于 LIPSOL ([52]),它是 Mehrotra 预测-校正算法 ([47]) 的变体,该算法是一种原始-对偶内点方法。

主算法

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

minxfTx such that {Ax=b0xu.(4)

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

minxfTx such that {Ax=bx+s=ux0, s0.(5)

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

maxbTyuTw  such that  {ATyw+z=fz0, w0,(6)

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

F(x,y,z,s,w)=(Axbx+suATyw+zfxizisiwi)=0,                 x0, z0, s0, w0,(7)

其中,xizi 和 siwi 表示按分量乘法。

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

线性规划的二次方程 xizi = 0siwi = 0 称为互补条件;其他(线性)方程称为可行性条件。量

xTz + sTw

对偶间隙,它度量当 (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 首先计算所谓的预测方向

Δvp=(FT(v))1F(v),

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

Δvc=(FT(v))1F(v+Δvp)μe^,

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

v+=v+α(Δvp+Δvc),

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

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

满足

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

在求解上述预测/校正方向的过程中,该算法根据修正的 Cholesky 因子 A·AT 来计算(稀疏)直接分解。如果 A 有稠密列,则算法改用 Sherman-Morrison 公式。如果该解不充分(残差太大),它将对步方程的增广方程组形式执行 LDL 分解来求得解。(请参阅MATLAB® ldl 函数参考页中的示例 4 - D 的结构。)

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

max(rbmax(1,b),rfmax(1,f),rumax(1,u),|fTxbTy+uTw|max(1,|fTx|,|bTyuTw|))tol,

其中

rb=Axbrf=ATyw+zfru={x}+su

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

fTxbTy+uTw

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

原始不可行性的度量是 ||rb||,对偶不可行性的度量是 ||rf||,其中范数是欧几里德范数。

预处理

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

  • 检查是否有任何变量具有相等的上界和下界。如果有,请检查可行性,然后固定变量的值以去除变量。

  • 检查是否有任何线性不等式约束只涉及一个变量。如果是,请检查可行性,然后将线性约束更改为一个边界。

  • 检查是否有任何线性等式约束只涉及一个变量。如果是,请检查可行性,然后固定该变量的值以去除该变量。

  • 检查是否有任何线性约束矩阵具有全零行。如果是,请检查可行性,然后删除这些行。

  • 确定边界和线性约束是否一致。

  • 检查是否有任何变量仅作为线性项出现在目标函数中,而不出现在任何线性约束中。如果是这样,请检查可行性和有界性,然后将变量固定为其适当的边界值。

  • 通过添加松弛变量,将任何线性不等式约束更改为线性等式约束。

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

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

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

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

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

对偶单纯形算法

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

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

minxfTx such that {Ax=b0xu.

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

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

x+={xif x>00if x0.

原始不可行性的度量是

Primal infeasibility=((lbx)+)2+((xub)+)2+((Axb)+)2+|Aeqxbeq|2.

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

maxbTyuTw  such that  {ATyw+z=fz0, w0.

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

对偶不可行性的度量是

Dual infeasibility=ATy+zwf2.

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

为了帮助减轻退化(请参阅 Nocedal 和 Wright 的论著 [7],第 366 页),对偶单纯形算法从扰动目标函数开始。

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

 阶段 1 概述

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

 阶段 2 概述

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

基变量和非基变量

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

minxfTx such that {Ax=b,lbxub.

(请注意,A 和 b 不是原始问题中定义不等式的矩阵和向量。)假设 A 是 m×n 矩阵,秩 m < n,包含列 {a1, a2, ..., an}。假设 {ai1,ai2,...,aim} 是 A 的列空间的基,索引集 B = {i1, i2, ..., im},N = {1, 2, ..., n}\B 是 B 的互补矩阵。子矩阵 AB 称为,互补子矩阵 AN 称为非基。由基变量组成的向量是 xB,由非基变量组成的向量是 xN。在阶段 2 的每次迭代中,该算法用非基的一列替换当前基的一列,并相应地更新变量 xB 和 xN

如果 x 是 A·x = b 的解,并且 xN 中的所有非基本变量都等于其下界或上界,则 x 称为基本解。此外,如果 xB 中的基变量满足其下界和上界,使得 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 https://www.maths.ed.ac.uk/~gondzio/software/correctors.ps.

[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.