线性规划算法
线性规划定义
线性规划问题求解向量 x,它在线性约束下最小化线性函数 fTx:
使得以下一个或多个条件成立:
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 的论著 [8](第 397 页)。
预测-校正
与 fmincon
内点算法相似,interior-point
算法尝试找到 卡鲁什-库恩-塔克 (KKT) 条件成立的点。要在线性规划问题中描述这些方程,请考虑使用经预处理的线性规划问题的标准形式:
现在假设所有变量都有至少一个有限边界。此假设根据需要对分量进行平移和求反,使 x 的所有分量均有下界 0。
是同时包含线性不等式和线性等式的扩展线性矩阵。 是对应的线性等式向量。 还包括用于扩展向量 x 的项,这些项使用松弛变量 s 将不等式约束转换为等式约束:
其中,x0 表示原始 x 向量。
t 是将上界转换为等式的松弛变量组成的向量。
此方程组的拉格朗日函数包括以下向量:
y,与线性等式相关联的拉格朗日乘数
v,与下界相关联的拉格朗日乘数(正性约束)
w,与上界相关联的拉格朗日乘数
此拉格朗日函数写作
因此,此方程组的 KKT 条件是
linprog
算法对返回的拉格朗日乘数使用与本文讨论内容不同的符号约定。本讨论内容使用与大多数文献相同的符号。请参阅 lambda
。
该算法首先根据牛顿-拉夫逊公式计算一个预测步,然后再计算一个校正步。校正器尝试减小非线性互补方程 sizi = 0 中的残差。牛顿-拉夫逊步是
(1) |
其中,X、V、W、T 是分别对应于向量 x、v、w 和 t 的对角矩阵。方程最右边的各个残差向量分别是:
rd,对偶残差
rp,原始残差
rub,上界残差
rvx,下界互补残差
rwt,上界互补残差
迭代输出将报告这些量:
要求解 公式 1,首先将其转换为对称矩阵形式
(2) |
其中
D 和 R 定义中的所有矩阵的逆都很容易计算,因为矩阵是对角矩阵。
要从公式 1 推导公式 2,请注意 公式 2 的第二行与 公式 1 的第二个矩阵行相同。公式 2 的第一行来自于先求解公式 1 的最后两行以得到 Δv 和 Δw,再求解以得到 Δt。
公式 2 是对称的,但由于 –D 项,它不是正定方程。因此,不能使用乔列斯基分解对其求解。多执行几步可得到另一个方程,它是正定方程,因此可以通过乔列斯基分解高效求解。
公式 2 的第二组行是
第一组行是
对下式进行代入
可得
(3) |
通常,找到牛顿步的最高效方法是使用乔列斯基分解求解公式 3 以获得 Δy。之所以可以使用乔列斯基分解,是因为乘以 Δy 的矩阵明显对称,并且在没有退化的情况下是正定矩阵。然后,要找到牛顿步,请回代以求得 Δx、Δt、Δv 和 Δw。但是,当 有稠密列时,求解公式 2 更为高效。linprog
内点算法根据列的稠密度选择求解算法。
有关更多算法细节,请参阅 Mehrotra 的论著 [7]。
在计算校正的牛顿步后,算法执行更多计算,以获得更长的当前步,并准备好计算更优的后续步。这些多重校正计算可以提高性能和稳健性。关详细信息,请参阅 Gondzio 的论著 [5]。
除了二次项,预测-校正算法与 quadprog
'interior-point-convex'
满算法基本相同。请参阅满预测-校正算法。
停止条件
预测-校正算法会持续迭代,直至达到一个可行(在容差范围内满足约束)且相对步长较小的点。具体来说,定义
当满足所有下列条件时,算法停止:
其中
rc 本质上是在测量互补残差 xv 和 tw 的大小,它们均为由零组成的向量,位于同一个解处。
内点传统线性规划
简介
内点传统方法基于 LIPSOL ([52]),它是 Mehrotra 预测-校正算法 ([47]) 的变体,该算法是一种原始-对偶内点方法。
主算法
该算法首先应用一系列预处理步骤(请参阅预处理)。在预处理后,问题具有以下形式
(4) |
上界约束隐式包含在约束矩阵 A 中。增加原始松弛变量 s 后,公式 4 变为
(5) |
以上为原始问题:x 由原始变量组成,s 由原始松弛变量组成。其对偶问题是
(6) |
其中 y 和 w 由对偶变量组成, z 由对偶松弛变量组成。此线性规划的最优性条件,即原始公式 5 和对偶公式 6,是
(7) |
其中,xizi 和 siwi 表示按分量乘法。
linprog
算法对返回的拉格朗日乘数使用与本文讨论内容不同的符号约定。本讨论内容使用与大多数文献相同的符号。请参阅 lambda
。
线性规划的二次方程 xizi = 0 和 siwi = 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 首先计算所谓的预测方向
这是牛顿方向;然后是所谓的校正方向
其中,μ > 0 称为中心化参数,必须谨慎选择。 是由 0 和 1 组成的向量,其中 1 对应于 F(v) 中的二次方程,即扰动仅适用于互补条件,这些条件均为二次的;但不适用于可行性条件,这些条件均为线性的。这两个方向通过一个步长参数 α > 0 进行合并,并更新 v 以获得新迭代 v+:
其中步长参数 α 的取值使得
v+ = [x+;y+;z+;s+;w+]
满足
[x+;z+;s+;w+] >0。
在求解上述预测/校正方向的过程中,该算法根据修正的乔列斯基因子 A·AT 来计算(稀疏)直接分解。如果 A 有稠密列,则算法改用谢尔曼-莫里森公式。如果该解不充分(残差太大),它将对步方程的增广方程组形式执行 LDL 分解来求得解;请参阅 ldl
。
然后算法执行循环计算,直到迭代收敛。主终止条件是一个标准条件:
其中
分别是原始残差、对偶残差和上界可行性({x} 是指具有有限上界的那些 x),并且
是原始目标值和对偶目标值之间的差,而 tol 是一定的容差。终止条件中的总和测度公式 7 中最优性条件中的总相对误差。
原始不可行性的测度是 ||rb||,对偶不可行性的测度是 ||rf||,其中范数是欧几里德范数。
预处理
该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容:
检查是否有任何变量具有相等的上界和下界。如果有,请检查可行性,然后固定变量的值以去除变量。
检查是否有任何线性不等式约束只涉及一个变量。如果是,请检查可行性,然后将线性约束更改为一个边界。
检查是否有任何线性等式约束只涉及一个变量。如果是,请检查可行性,然后固定该变量的值以去除该变量。
检查是否有任何线性约束矩阵具有全零行。如果是,请检查可行性,然后删除这些行。
确定边界和线性约束是否一致。
检查是否有任何变量仅作为线性项出现在目标函数中,而不出现在任何线性约束中。如果是这样,请检查可行性和有界性,然后将变量固定为其适当的边界值。
通过添加松弛变量,将任何线性不等式约束更改为线性等式约束。
如果算法检测到不可行或无界的问题,它会停止并发出相应的退出消息。
该算法可能到达一个表示解的单个可行点。
如果算法在预求解步骤中没有检测到不可行或无界的问题,并且预求解没有产生解,则算法会继续其后续步骤。到达停止条件后,算法会重新构造原始问题,这将撤消任何预求解变换。这最后一步是后求解步骤。
为简单起见,该算法将所有下界移至零。
虽然这些预处理步骤可以大大加快算法的迭代部分,但如果需要拉格朗日乘数,则必须撤消预处理步骤,因为在算法中计算的乘数适用于变换的问题,而不适用于原始问题。因此,如果不需要乘数,则不计算这种反向变换,因而可能节省一些计算时间。
对偶单纯形传统算法
就整体层面而言,linprog
'dual-simplex-legacy'
算法本质上是对对偶问题执行单纯形算法。
该算法首先执行预处理,如预处理中所述。有关详细信息,请参阅 Andersen 和 Andersen 的论著 [1] 以及 Nocedal 和 Wright 的论著 [8],第 13 章。这种预处理将原始线性规划问题简化为公式 4 形式:
A 和 b 是原始约束矩阵的变换后的版本。这是原始问题。
原始可行性可以用 + 函数来定义
原始不可行性的测度是
正如公式 6 中所述,对偶问题是找到向量 y 和 w 以及松弛变量向量 z 来求解
linprog
算法对返回的拉格朗日乘数使用与本文讨论内容不同的符号约定。本讨论内容使用与大多数文献相同的符号。请参阅 lambda
。
对偶不可行性的测度是
众所周知(例如,请参阅 [8]),对偶问题的任何有限解都是原始问题的一个解,而原始问题的任何有限解都是对偶问题的一个解。而且,如果原始问题或对偶问题中有一个无界,则另一个问题不可行。并且,如果原始问题或对偶问题中有一个不可行,则另一个问题要么不可行,要么无界。因此,这两个问题在求有限解(如有)方面是等效的。由于原始问题和对偶问题在数学上是等效的,但计算步骤不同,因此有时更适合通过求解对偶问题来求解原始问题。
为了帮助减轻退化(请参阅 Nocedal 和 Wright 的论著 [8],第 366 页),对偶单纯形算法从扰动目标函数开始。
对偶单纯形算法的阶段 1 是求得一个对偶可行点。该算法通过求解辅助线性规划问题来实现这一点。
在阶段 2 中,求解器反复选择一个入基变量和一个出基变量。该算法根据福雷斯特和戈德法布[3] 提出的一种称为对偶最陡边定价的方法选择出基变量。该算法使用科伯斯坦[6]提出的哈里斯比率检验变体来选择入基变量。为了帮助减轻退化,该算法可以在阶段 2 中引入额外的扰动。
求解器执行迭代,尝试在降低原始不可行性的同时保持对偶可行性,直到经过简化和扰动的问题的解既是原始可行的又是对偶可行的。该算法可消除它引入的扰动。如果(扰动后问题的)解对于未扰动(原始)问题是对偶不可行的,则求解器使用原始单纯形或阶段 1 算法还原对偶可行性。最后,求解器消除预处理步骤,以返回原始问题的解。
对偶单纯形 Highs 算法
'dual-simplex-highs'
算法基于 HiGHS 开源软件。就整体层面而言,linprog
'dual-simplex-highs'
算法本质上执行一种单纯形算法,该算法在朝着原始可行性迭代的同时保持对偶可行性。请参阅原始和对偶可行性。有关算法的详细信息,请参阅 Huangfu 和 Hall 的论著 [4],第 2 节。
该算法首先执行预处理,如预处理中所述。有关详细信息,请参阅 Andersen 和 Andersen 的论著 [1] 以及 Nocedal 和 Wright 的论著 [8],第 13 章。这种预处理通常会将原始线性规划问题简化为更容易求解的较小问题。
众所周知(例如,请参阅 [8]),对偶问题的任何有限解都是原始问题的一个解,而原始问题的任何有限解都是对偶问题的一个解。而且,如果原始问题或对偶问题中有一个无界,则另一个问题不可行。并且,如果原始问题或对偶问题中有一个不可行,则另一个问题要么不可行,要么无界。因此,这两个问题在求有限解(如有)方面是等效的。由于原始问题和对偶问题在数学上是等效的,但计算步骤不同,因此有时更适合通过求解对偶问题来求解原始问题。
为了帮助减轻退化(请参阅 Nocedal 和 Wright 的论著 [8],第 366 页),对偶单纯形算法从扰动目标函数系数开始。
对偶单纯形算法的阶段 1 旨在找到一个对偶可行点。算法通过对阶段 1 问题应用阶段 2 算法(见下文)来实现这一点。此问题具有经过下述修改的(扰动的)目标成本系数以及变量和约束的边界。单边上(下)界变量和约束的边界设置为 [0, 1] ([-1, 0]),固定变量和方程上的边界设置为零,自由变量和约束上的边界设置为 [–1000, 1000]。由于所有边界均为有限值,通过将原始变量设置在适当的边界上,会使任何基变量都对偶可行。如果一个(非基)变量的对偶值相对于原始边界不可行,则它对对偶目标函数值的贡献为负。注意,对偶目标函数的上界被限制为零。当阶段 1 问题求解到最优时,如果目标值为零,则当前点相对于原始边界是对偶可行的。负目标值意味着原始问题在扰动的目标成本系数下对偶不可行,并且强烈暗示未扰动的原始问题的对偶不可行。通过恢复为未受扰动的目标函数系数并尝试从当前对偶解继续求解阶段 1 问题来评估这种情况。
阶段 1 的解作为主算法的阶段 2 的初始点。
在阶段 2 中,求解器反复选择一个出基变量和一个入基变量。该算法根据福雷斯特和戈德法布[3] 提出的一种称为对偶最陡边定价的方法选择出基变量。该算法使用科伯斯坦[6]提出的哈里斯比率检验变体来选择入基变量。通过移动成本系数,可以消除由于舍入而导致的任何轻微的对偶不可行性。
在阶段 2 中,对偶单纯形算法,从阶段 1 的对偶可行点开始,用于求解原始问题。在每次迭代中,算法测试最优性条件(原始可行性),如果当前解是最优解,则算法停止执行。如果当前解不是最优解,则算法使用 基变量和非基变量 进行定义并执行以下操作:
从原始值不可行的基变量中选择一个变量(称为出基变量)。
从非基变量中选择一个变量(称为入基变量),并用它来替换基变量中的出基变量。
更新当前原始解和对偶解以及当前目标值。
在所谓的对偶比率测试中,尽可能增大出基变量的对偶值,同时保持对偶可行性。入基变量是非基变量,其对偶值在界限内被归零。计算对偶比率的数据需要解一个方程组。在基变量改变后,必须求解第二个方程组来更新原始解。如果对偶比率测试中没有界限,则扰动问题是对偶无界的(因此是原始不可行的)。
求解器执行迭代,在降低原始不可行性的同时保持对偶可行性,直到经过简化和扰动的问题的解既是原始的又是对偶可行的。然后,算法再消除目标成本扰动。如果(扰动问题的)解对于未扰动(原始)问题是对偶不可行的,则求解器将使用原始单纯形阶段 2 算法来恢复对偶可行性。最后,求解器消除预处理步骤,以返回原始问题的解。
基变量和非基变量
本节定义线性规划问题的术语基、非基和基可行解。定义假设问题以如下标准形式给出:
(请注意,A 和 b 不是原始问题中定义不等式的矩阵和向量。)假设 A 是 m×n 矩阵,秩 m < n,包含列 {a1, a2, ..., an}。假设 是 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 称为基可行解。
原始和对偶可行性
原始不可行性的度量报告为最大绝对约束违反值,用原始问题变量表示如下,
max(0, lb – x, x – ub, abs(Aeqx – beq), Ax – b)。 | (8) |
要定义对偶可行性,请首先将线性矩阵 A(来自 基变量和非基变量)分成两个模块,即一个包含线性无关列的基变量矩阵 B 和一个非基变量矩阵 N:
然后,方程 A x = b 的变量 x 自然地分成 xB 和 xN:
由于 B 由线性无关列组成,因此它是可逆的,即
同样,目标函数 fTx 可以写为
缩减成本是非基变量的调整成本系数。
如果对于处于其下界的所有非基变量,此缩减成本值 cj(cNT 中的第 j 个条目)大于或等于 0,并且处于其上界的所有非基变量的 cj 小于或等于 0,则当前解是对偶可行的。对偶单纯形算法保持此对偶可行性,并尝试使单纯形迭代原始可行。
参考
[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] Huangfu, Q. and J. A. J. Hall. Parallelizing the dual revised simplex method. Math. Prog. Comp. 10, pp. 119–142, 2018. Available at https://link.springer.com/article/10.1007/s12532-017-0130-5
.
[5] 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
.
[6] 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.
[7] Mehrotra, S. “On the Implementation of a Primal-Dual Interior Point Method.” SIAM Journal on Optimization, Vol. 2, 1992, pp 575–601.
[8] Nocedal, J., and S. J. Wright. Numerical Optimization, Second Edition. Springer Series in Operations Research, Springer-Verlag, 2006.