混合整数线性规划 (MILP) 算法
混合整数线性规划定义
混合整数线性规划 (MILP) 问题具有以下要素：
线性目标函数 fTx，其中 f 是由常数组成的列向量，x 是由未知数组成的列向量
边界和线性约束，但没有非线性约束（有关定义，请参阅编写约束）
对 x 的某些分量的限制，使其必须具有整数值
以数学语言表达，即根据向量 f、lb 和 ub，矩阵 A 和 Aeq，对应的向量 b 和 beq，以及索引集 intcon，求解向量 x 使下式成立
HiGHS MILP 算法
HiGHS 概述
intlinprog 算法基于 HiGHS 开源软件。intlinprog 将 MATLAB® 格式的输入和选项转换为等效的 HiGHS 参量，并将返回的解转换为标准 MATLAB 格式。
算法概要
算法执行以下步骤。
获取分支定界节点（如果没有，则停止）
重复，直到达到停止条件：
搜索，直到达到停止条件：
执行 下潜/潜水
传播域
剪除节点，更新边界，如果检测到不可行性或达到
ObjectiveCutOff选项值则退出如果满足重启条件，则返回步骤 2
安装下一个节点：
选择并计算节点
如果计算后剪除了此节点，则返回步骤 5
为节点生成切割
如果域不可行，则切掉该节点，并将开放节点放入节点队列
更新基变量
转至步骤 4
预求解
通常，我们可以减少问题中的变量数量（x 的分量数量），并减少线性约束的数量。虽然执行这些缩减可能会花费求解器的时间，但它们通常会缩短求解的总时间，并使更大型的问题可求解。这些算法可以使解在数值上更加稳定。此外，这些算法有时可以检测不可行问题。
预求解步骤旨在消除冗余变量和约束，改善模型的尺度和约束矩阵的稀疏性，加强变量的边界，检测模型的原始和对偶不可行性。有关背景，请参阅 Andersen 和 Andersen 的论著 [3] 以及 Mészáros 和 Suhl 的论著 [10]。
在混合整数规划预处理期间，intlinprog 分析线性不等式 A*x ≤ b 以及完整性约束，以确定是否有以下情形：
此问题不可行。
有些边界可以收紧。
有些不等式是冗余的，因此可以忽略或删除。
一些不等式可以得到加强。
一些整数变量可以固定。
有关整数预处理的背景，请参阅 Achterberg 等人的论著 [1]。
计算根节点
为了计算根节点，算法执行以下步骤。
检测对称性并简化问题。
评估根 LP。
生成并添加 LP 切割（请参阅切割生成）。
应用随机舍入。
生成并添加更多 LP 切割。
循环执行切割生成和启发式方法。
检查重启条件，如果有必要，则从步骤 2 重启。
在根节点计算完成后，算法将继续执行分支定界。
切割生成
切割是 intlinprog 为问题增加的额外线性不等式约束。这些不等式尝试限制 LP 松弛的可行域，使其解更接近整数。有关切割生成算法（亦称割平面方法）的背景信息，请参阅 Cornuéjols 的论著 [6] 以及 Atamtürk、Nemhauser 和 Savelsbergh 的论著 [4]。
下潜/潜水
为了找到整数可行点，intlinprog 使用类似于分支定界步骤的启发式方法，但只沿树的一个分支向下，而不创建其他分支。此单分支使得算法能够快速下“潜”到树的片段，因此命名为“潜水”。
潜水启发式方法通常选择一个应该取整数值、但当前解为小数的变量。然后，启发式方法引入一个强制变量取整数值的边界，并再次求解相关联的松弛 LP。各种潜水启发式方法之间的主要区别在于如何选择要定界的变量。请参阅 Berthold 的论著·[5]，第 3.1 节。
随机舍入、RINS 和 RENS
为了找到新的整数可行点，intlinprog 搜索当前最佳整数可行解点（如有）的邻域，以求得更好的新解。请参阅 Danna、Rothberg 和 Le Pape·的论著 [8]。同样，为了找到新的整数可行点，intlinprog 取节点处松弛问题的 LP 解，并以尝试保持可行性的方式舍入整数分量。通过采取随机舍入步骤，intlinprog 有时可以找到新的可行点。RENS 表示松弛强制邻域搜索，它是寻找整数可行点的另一种搜索方法。请参阅 Berthold 的论著 [5]。
分支定界
分支定界方法构造一系列尝试收敛到 MILP 解的子问题。子问题给出关于解 fTx 的一系列上界和下界。这些边界称为原始边界和对偶边界。对于最小化问题，第一个上界（原始）是任何可行解，第一个下界（对偶）是松弛问题的解。对于最大化问题，原始边界是下界，对偶边界是上界。线性规划松弛问题的任何解都比 MILP 解具有更低的目标函数值。此外，任何可行点 xfeas 都满足
| fTxfeas ≥ fTx， | (1) |
因为 fTx 是所有可行点中最小的。
在这种情况下，节点是具有以下特征的 LP：它具有与原始问题相同的目标函数、边界和线性约束，但没有整数约束，并对线性约束或边界作出特定改变。根节点是没有整数约束且线性约束或边界也不发生变化的原始问题，这意味着根节点是初始松弛 LP。
从起始边界开始，分支定界方法通过从根节点产生分支来构造新子问题。分支步骤采用启发式方法，遵循以下规则之一。每个规则都基于这样的理念：约束一个变量，使之小于或等于整数 J，或者大于或等于 J+1，从而拆分问题。当 xLP 中对应于 intcon 所指定整数的条目不是整数时，会出现这两个子问题。在此处，xLP 是松弛问题的解。以 J 为变量的下限（向下舍入），以 J+1 为上限（向上舍入）。由此产生的两个问题的解大于或等于 fTxLP，因为这两个问题具有更多约束。因此，对于最小化问题来说，此过程可能会提高下界。
在算法分支后，有两个新节点需要探查。intlinprog 通过将其目标函数值与当前全局边界进行比较，跳过对某些子问题的分析。
分支定界过程继续按部就班地生成子问题以进行分析，并丢弃那些不会改进目标函数上界或下界的子问题，直到满足以下停止条件之一：
问题被证明不可行。
目标函数值达到
ObjectiveCutOff限制。算法超过
MaxTime选项指定的运行时间。目标函数的上下界之差小于
AbsoluteGapTolerance或RelativeGapTolerance容差。探查的节点数量超过
MaxNodes选项指定的数量。
迭代输出
当您通过将 Display 选项设置为默认值 "iter" 来选择迭代输出时，求解器会显示它的一些步骤。HiGHS 迭代输出比其他求解器的迭代输出更加广泛和复杂。此外，HiGHS 算法可以重新开始其分支定界搜索，从而导致迭代输出也重新开始。
要选择迭代输出，请使用以下代码：
options = optimoptions("intlinprog",Display="iter"); [x,fval,exitflag,output] = intlinprog(f,intcon,A,b,Aeq,beq,... lb,ub,options)
预备步. 迭代输出从显示“预求解”的结果开始。预求解算法通过识别和删除线性约束矩阵中的冗余行和列并执行问题的相关简化来降低原始问题的复杂性。例如，
Presolving model 18018 rows, 26027 cols, 248579 nonzeros 15092 rows, 24343 cols, 217277 nonzeros Objective function is integral with scale 1
根节点计算. 根节点是问题的线性规划解，不考虑任何整数约束。对于最小化问题，根节点的目标函数值是包含整数约束的问题解的目标函数值的下界。对于最小化问题，上界（如果有）来自符合所有约束的一个可行点。如果还未找到可行点，则上界为 Inf。
迭代输出显示预求解后问题的大小。
Solving MIP model with: 15092 rows 24343 cols (24343 binary, 0 integer, 0 implied int., 0 continuous) 217277 nonzeros
binary是二元变量的个数。integer是整数变量的个数。implied int是暗示为整数的变量的个数。例如，如果x(1)为整数且x(1) + x(2) = 5，则暗示x(2)为整数。continuous是连续变量的个数。
变量总数是模型中的列数。
动态约束创建. 软件首先创建“动态约束”，在迭代输出中有三个标题：
切割 - 活动切割的数量
inLp - LP 矩阵中非模型行的数量
Confl. - 冲突的数量
软件可以通过重新启动动态约束创建来进一步扩展约束，这可以通过从新状态开始创建过程来创建更多约束。
在开始分支定界过程之前的最后一步中，软件会报告对称性检测的结果以及找到的生成器和轨道体 (orbitope) 的数量。例如，以下是出现在预备步和动态约束创建阶段中的部分迭代输出：
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 0 inf inf 0 0 0 0 1.5s
0 0 0 0.00% 0 inf inf 0 0 5 4934 3.6s
…
0 0 0 0.00% 0 inf inf 4630 553 289 140296 159.6s
0.0% inactive integer columns, restarting
Model after restart has 15090 rows, 24338 cols (24338 bin., 0 int., 0 impl., 0 cont.), and 217075 nonzeros
0 0 0 0.00% 0 inf inf 550 0 0 140296 160.5s
…
0 0 0 0.00% 0 inf inf 5602 524 260 277323 318.0s
6.2% inactive integer columns, restarting
Model after restart has 14185 rows, 22816 cols (22816 bin., 0 int., 0 impl., 0 cont.), and 200423 nonzeros
0 0 0 0.00% 0 inf inf 524 0 0 277323 318.6s
…
0 0 0 0.00% 0 inf inf 4683 535 525 408393 505.8s
Symmetry detection completed in 3.1s
Found 215 generators and 12 full orbitope(s) acting on 730 columns有关对称性检测和轨道体的信息，请参阅 Hojny、Pfetsch 和 Schmitt 的论著 [7] 以及 Pfetsch 和 Rehn 的论著 [12]。软件在进行下文描述的分支定界步骤时继续添加动态约束。
分支定界. 分支定界方法构造一系列尝试收敛到 MILP 解的子问题。子问题给出关于解 fTx 的一系列上界和下界。对于最小化问题，第一个上界是任何可行解，第一个下界是松弛问题的解。
在分支定界过程中，intlinprog 给出以下迭代输出。
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
72 0 2 1.56% 0 inf inf 4695 535 738 667009 686.9s
…
T 271 107 51 1.56% 0 449 100.00% 6051 271 1767 792786 776.8s
T 279 107 53 1.56% 0 439 100.00% 6053 271 1794 793241 777.5s
…
L 1223 538 295 1.98% 0 434 100.00% 6984 243 7628 1689k 1580.6s
…
1321 539 333 1.98% 0 434 100.00% 7029 243 8628 1898k 1650.6s
Restarting search from the root node
Model after restart has 13947 rows, 22426 cols (22426 bin., 0 int., 0 impl., 0 cont.), and 194633 nonzeros
1323 0 0 0.00% 0 434 100.00% 243 0 0 1902k 1653.5s
1323 0 0 0.00% 0 434 100.00% 243 76 10 1905k 1655.7s
…
1694 173 52 0.00% 0 434 100.00% 9411 318 2220 3205k 2584.3s
B 1710 167 55 0.00% 0 433 100.00% 9415 318 2263 3207k 2586.2s
1726 224 56 0.00% 0 433 100.00% 9999 378 2307 3237k 2608.7s最左边的一列显示代码，指示如何为输出中的该行找到新可行点：
L- 在原始启发式方法中求解子 MIP 问题时T- 在树搜索期间，计算节点时B- 在分支期间H- 通过启发式方法P- 在启动期间，求解 MIP 之前C- 通过中心舍入R- 通过随机舍入S- 求解 LP 时F- 通过可行性泵U- 基于无界 LP
完成. 迭代输出结束时会显示算法停止的原因和结果摘要。
Solving report
Status Time limit reached
Primal bound 14
Dual bound 0
Gap 100% (tolerance: 0.01%)
Solution status feasible
14 (objective)
0 (bound viol.)
1.7763568394e-15 (int. viol.)
0 (row viol.)
Timing 7200.02 (total)
3.08 (presolve)
0.00 (postsolve)
Nodes 5830
LP iterations 9963775 (total)
2657448 (strong br.)
324908 (separation)
2140011 (heuristics)
Solver stopped prematurely. Integer feasible point found.
Intlinprog stopped because it exceeded the time limit, options.MaxTime = 7200 . The intcon variables are integer within tolerance,
options.ConstraintTolerance = 1e-06.Status- 迭代停止的原因Primal bound- 对于最小化问题，显示目标函数值的上界；对于最大化问题，显示目标函数值的下界Dual bound- 对于最大化问题，显示目标函数值的下界；对于最大化问题，显示最小化问题的上界Gap- 原始边界和对偶边界之间的相对间隙Solution status解的状态
目标函数值，标注为
(objective)解对变量边界的最大违反值，标注为
(bound viol.)整数类型变量与整数值的最大违反值，标注为
(int. viol.)解对线性约束的最大违反值，标注为
(row viol.)
Timing- 各个求解阶段的计时，以秒为单位Nodes- 探查的节点数LP iterations- 各阶段及总共的线性规划迭代次数
参考
[1] Achterberg, T.,Bixby, R., Gu, Z., Rothberg, E., and Weninger, D. Presolve Reductions in Mixed Integer Programming. INFORMS J. Computing 32(2), 2019, pp. 473–506. Available at https://doi.org/10.1287/ijoc.2018.0857.
[2] Achterberg, T., T. Koch and A. Martin. Branching rules revisited. Operations Research Letters 33, 2005, pp. 42–54. Available at https://opus4.kobv.de/opus4-zib/files/788/ZR-04-13.pdf.
[3] Andersen, E. D., and Andersen, K. D. Presolving in linear programming. Mathematical Programming 71, pp. 221–245, 1995.
[4] Atamtürk, A., G. L. Nemhauser, M. W. P. Savelsbergh. Conflict graphs in solving integer programming problems. European Journal of Operational Research 121, 2000, pp. 40–55.
[5] Berthold, T. Primal Heuristics for Mixed Integer Programs. Technischen Universität Berlin, September 2006. Available at https://opus4.kobv.de/opus4-zib/files/1029/Berthold_Primal_Heuristics_For_Mixed_Integer_Programs.pdf.
[6] Cornuéjols, G. Valid inequalities for mixed integer linear programs. Mathematical Programming B, Vol. 112, pp. 3–44, 2008.
[7] Hojny, C., Pfetsch, M. E., Schmitt, A. Extended formulations for column-constrained orbitopes. TU Darmstadt, Department of Mathematics, Research Group Optimization, Dolivostr. 15, 64293 Darmstadt, June 2017. Available at https://optimization-online.org/2017/06/6092/.
[8] Danna, E., Rothberg, E., Le Pape, C. Exploring relaxation induced neighborhoods to improve MIP solutions. Mathematical Programming, Vol. 102, issue 1, pp. 71–90, 2005.
[9] Hendel, G. New Rounding and Propagation Heuristics for Mixed Integer Programming. Bachelor's thesis at Technische Universität Berlin, 2011. PDF available at https://opus4.kobv.de/opus4-zib/files/1332/bachelor_thesis_main.pdf.
[10] Mészáros C., and Suhl, U. H. Advanced preprocessing techniques for linear and quadratic programming. OR Spectrum, 25(4), pp. 575–595, 2003.
[11] Nemhauser, G. L. and Wolsey, L. A. Integer and Combinatorial Optimization. Wiley-Interscience, New York, 1999.
[12] Pfetsch, M. E. and Rehn, T. A computational comparison of symmetry handling methods for mixed integer programs. Mathematical Programming Computation (2019) 11:37–93. Available at https://optimization-online.org/wp-content/uploads/2015/11/5209.pdf.
[13] Wolsey, L. A. Integer Programming. Wiley-Interscience, New York, 1998.
