混合整数规划与混合整数二次规划
引言
在数值优化领域中，混合整数规划（Mixed Integer Programming, MIP）和混合整数二次规划（Mixed Integer Quadratic Programming, MIQP）是两类重要的优化模型。本文将对这两种规划的基本概念、常用求解方法进行介绍，并通过实例分析帮助大家更好地理解这些规划问题。
一、混合整数规划（MIP或MILP）基本概念
混合整数规划问题是指在优化过程中，部分优化变量必须取整数值，其他变量则可以取连续值的一类优化问题。
典型的MIP问题形式如下：
| 目标方程： | min c T x \min \ \mathbf{c} ^T \mathbf{x} min cTx |
|---|---|
| 约束： | A x ≤ b \mathbf{A} \mathbf{x} \le\mathbf{b} Ax≤b (线性约束) |
| x l ≤ x ≤ x h \mathbf{x} _l\le \mathbf{x} \le \mathbf{x} _h xl≤x≤xh(边界约束) | |
| x i ∈ Z \mathbf{x}_i∈Z xi∈Z, x j ∈ R \mathbf{x}_j∈R xj∈R (for some i, others j) （整数约束） |
这里， c \mathbf{c} c和 x \mathbf{x} x 是向量， A \mathbf{A} A是矩阵， b \mathbf{b} b是约束条件。换句话说，混合整数规划的目标函数与线性规划不同的是，优化变量的定义域可能是整数，这就给求解带来了极大的难度。
二、混合整数二次规划基本概念
混合整数二次规划是混合整数规划的扩展，其目标函数是一个二次函数，而约束条件通常为一次形式。
典型的MIQP问题形式如下：
| 目标方程： | min x T Q x + c T x \min \ \mathbf{x}^{T}\mathbf{Q}\mathbf{x}+ \mathbf{c} ^T \mathbf{x} min xTQx+cTx |
|---|---|
| 约束： | A x ≤ b \mathbf{A} \mathbf{x} \le\mathbf{b} Ax≤b (线性约束) |
| x l ≤ x ≤ x h \mathbf{x} _l\le \mathbf{x} \le \mathbf{x} _h xl≤x≤xh(边界约束) | |
| x i ∈ Z \mathbf{x}_i∈Z xi∈Z, x j ∈ R \mathbf{x}_j∈R xj∈R (for some i, others j) （整数约束） |
上述变量定义与MIP完全相同，可以看到MIQP其实就是再MIP的目标函数的基础上，加上了一个二次项。
三、常见求解方法
1. 分支定界法BnB（最常用）
分支定界法是求解MIP和MIQP问题的经典方法之一。其基本思想是将问题分解为一系列子问题，通过构建一个分支树来逐步逼近最优解。分支定界法主要包括三个步骤：分支、定界和剪枝。简单说一下基于LP的BnB方法：
首先，对初始的MILP删除所有整数约束，得到原MILP的线性规划松弛。接下来，求解这个线性规划松弛（LP）。如果得到的解恰好满足所有整数约束，那么这个解就是原始MILP的最优解，运算终止。如果解没有满足所有整数约束（这种情况较为常见），则需要选择某个整数约束实际上得到小数值的变量进行“分支”。为了便于说明，假设这个变量是x，其在LP解中的值是5.7。我们可以通过施加 x ≤ 5.0 x≤5.0 x≤5.0和 x ≥ 6.0 x≥6.0 x≥6.0的限制来排除该值5.7。
如果用 P P P表示原MILP问题，用 P 1 P_1 P1表示新的MILP（增加了 x ≤ 5.0 x≤5.0 x≤5.0的约束），用 P 2 P_2 P2表示新的MILP（增加了 x ≥ 6.0 x≥6.0 x≥6.0的约束）。变量 x x x被称为分支变量，我们在 x x x上进行了分支，产生了两个子MILP P 1 P_1 P1和 P 2 P_2 P2。显然，如果我们能够求解 P 1 P_1 P1和 P 2 P_2 P2的最优解，那么其中较好的解就是原始问题 P P P的最优解。
通过这种方式，我们用两个更简单（更少整数约束）的MILP取代了 P P P。我们现在将相同的思想应用于这两个子MILP，并在必要时选择新的分支变量。这样生成了所谓的搜索树。搜索过程生成的MILP称为树的节点， P P P为根节点。叶子节点是尚未分支的所有节点。如果所有叶节点的解都满足原MILP的整数约束，那么原始MILP问题就得到了求解。
2. 割平面法（Cutting Plane）
割平面法通过添加额外的约束条件（即割平面）来缩小可行解空间，从而更快地找到最优解。对于MIP和MIQP问题，可以通过不断添加有效的割平面来逼近整数解。
割平面通常认为是MIP中提升计算效率的最重要技巧。其基本思想为：思想是通过去除一些非整数解，就想MIP-based presolve一样，达到tighten formulation的目的。割平面不像branching，branching会产生两个子问题(两个分支)，cutting plane只是切掉一些边角(如下图)，不会产生2个MIP。
3. 启发式算法（Heuristic）
启发式算法是一类通过经验和启发信息来快速找到满意解的方法。常见的启发式算法包括遗传算法、模拟退火算法和粒子群优化算法。这些方法通常不能保证找到全局最优解，但可以在较短时间内找到质量较好的解。
4. 混合算法（Hybrid Algorithms）
混合算法结合了多种求解方法的优点，以期提高求解效率和解的质量。例如，可以将启发式算法与分支定界法结合，利用启发式算法快速找到初始解，然后通过分支定界法进行精细优化等等。
四、实例代码求解
本节我们举一个简单的生产实例，来说明如何使用Python的Pulp库进行MIP和MIQP问题求解，问题描述如下:
1. 生产计划问题（MILP问题）
假设某工厂生产两种产品产品 a a a和产品 b b b，每种产品都有其生产成本和销售利润，其中：
产品 a a a的生产成本为10元，销售价格为15元。
产品 b b b的生产成本为20元，销售价格为30元。
该工厂的生产能力和原材料有限，假设：
每月生产能力为70单位。
每月原材料限制为150单位，生产产品 a a a消耗1单位原材料，生产产品 b b b消耗3单位原材料。
问题一： 如何决定每种产品数量，以实现最大化利润？
问题分析：
从问题中，我们可以看出，待优化的决策变量为产品 a a a和产品 b b b的数量，用 x 1 x_1 x1和 x 2 x_2 x2表示，其中 x 1 x_1 x1和 x 2 x_2 x2必须取整数值，这个问题可以建模为一个MIP或MILP问题。
其中，我们的目标函数为：
max ( 15 − 10 ) x 1 + ( 30 − 20 ) x 2 \max {(15-10)x_1 + (30-20)x_2} max(15−10)x1+(30−20)x2
约束为：
x 1 + x 2 ≤ 70 x_1+x_2 \le 70 x1+x2≤70 x 1 + 3 x 2 ≤ 150 x_1+3x_2 \le 150 x1+3x2≤150 x 1 , x 2 ∈ Z + x_1, x_2 \in Z^{+} x1,x2∈Z+
Python求解代码（基于Pulp库）
import pulp as plimport numpy as npdef solver_main(): ProbLp=pl.LpProblem("ProbLp",sense=pl.LpMaximize) print(ProbLp.name) x1=pl.LpVariable('x_1',lowBound=0,upBound=None,cat='Integer') x2=pl.LpVariable('x_2',lowBound=0,upBound=None,cat='Integer') ProbLp+=(5*x1+10*x2) # Add Objective # Add Constraints ProbLp+=(x1+x2<=70) ProbLp+=(x1+3*x2<=150) ProbLp.solve() # Log for solver print("Shan Status:", pl.LpStatus[ProbLp.status]) print("Optimal objective value", pl.value(ProbLp.objective)) for v in ProbLp.variables(): print(v.name, "=", v.varValue)if __name__ == '__main__': solver_main()执行脚本，求得最优解： x 1 = 30 ， x 2 = 40 x_1=30，x_2=40 x1=30，x2=40，输出如下:
Welcome to the CBC MILP SolverVersion: 2.10.3Build Date: Dec 15 2019command line - D:\Anaconda3\envs\torch_env\lib\site-packages\pulp\solverdir\cbc\win\64\cbc.exe C:\Users\ADMINI~1\AppData\Local\Temp\7a689790a63c4d0bbd42e787d4d055db-pulp.mps -max -timeMode elapsed -branch -printingOptions all -solution C:\Users\ADMINI~1\AppData\Local\Temp\7a689790a63c4d0bbd42e787d4d055db-pulp.sol (default strategy 1)At line 2 NAME MODELAt line 3 ROWSAt line 7 COLUMNSAt line 18 RHSAt line 21 BOUNDSAt line 24 ENDATAProblem MODEL has 2 rows, 2 columns and 4 elementsCoin0008I MODEL read with 0 errorsOption for timeMode changed from cpu to elapsedContinuous objective value is 550 - 0.00 secondsCgl0004I processed model has 2 rows, 2 columns (2 integer (0 of which binary)) and 4 elementsCutoff increment increased from 1e-05 to 4.9999Cbc0012I Integer solution of -550 found by DiveCoefficient after 0 iterations and 0 nodes (0.00 seconds) Cbc0001I Search completed - best objective -550, took 0 iterations and 0 nodes (0.00 seconds)Cbc0035I Maximum depth 0, 0 variables fixed on reduced costCuts at root node changed objective from -550 to -550Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)Result - Optimal solution foundObjective value: 550.00000000Enumerated nodes: 0Total iterations: 0Time (CPU seconds): 0.00Time (Wallclock seconds): 0.00Option for printingOptions changed from normal to allTotal time (CPU seconds): 0.01 (Wallclock seconds): 0.01Shan Status: OptimalOptimal objective value 550.0x_1 = 30.0x_2 = 40.02. 考虑收益衰减的生产问题（MIQP问题）
在示例1中的生产问题是一种基本的生产模型，在实际生产过程中，产品的利润有可能会随着产品的数量增加而逐渐降低，呈现收益衰减的现象。我们在示例1的基础上，不妨假设：
产品 a a a的生产成本为10元，销售价格为 15 − 0.03 x 1 15-0.03{x_1} 15−0.03x1。
产品 b b b的生产成本为20元，销售价格为 30 − 0.02 x 2 30-0.02{x_2} 30−0.02x2元。
工厂生产能力和原材料与示例1相同，假设：
每月生产能力为70单位。
每月原材料限制为150单位，生产产品 a a a消耗1单位原材料，生产产品 b b b消耗3单位原材料。
问题二： 如何决定每种产品数量，以实现最大化利润？
问题分析：
从问题中，待优化的决策变量依然为产品 a a a和产品 b b b的数量，用 x 1 x_1 x1和 x 2 x_2 x2表示，其中 x 1 x_1 x1和 x 2 x_2 x2必须取整数值，但是收益模型变得更加复杂了。
其中，我们的目标函数为：
max ( 15 − 0.03 x 1 − 10 ) x 1 + ( 30 − 0.02 x 2 − 20 ) x 2 \max {(15-0.03x_1-10)x_1+(30-0.02x_2-20)x_2} max(15−0.03x1−10)x1+(30−0.02x2−20)x2 = max ( − 0.03 x 1 2 − 0.02 x 2 2 + 5 x 1 + 10 x 2 ) = \max {(-0.03{x_1}^{2}-0.02{x_2}^{2}+5x_1+10x_2)} =max(−0.03x12−0.02x22+5x1+10x2)
约束不变，仍为：
x 1 + x 2 ≤ 70 x_1+x_2 \le 70 x1+x2≤70 x 1 + 3 x 2 ≤ 150 x_1+3x_2 \le 150 x1+3x2≤150 x 1 , x 2 ∈ Z + x_1, x_2 \in Z^{+} x1,x2∈Z+
我们可以看出优化的目标函数为二次形式，因此此任务可以被建模为一个MIQP问题，我们将其转换为矩阵的表示形式：
max x T [ − 0.03 − 0.02 ] x + [ 5 10 ] T x \max {\mathbf{x} ^{T}\begin{bmatrix} -0.03 & \\ & -0.02 \end{bmatrix}\mathbf{x} + \begin{bmatrix} 5\\ 10 \end{bmatrix}^{T}\mathbf{x} } maxxT[−0.03−0.02]x+[510]Tx s . t . [ 1 1 1 3 ] x ≤ [ 70 150 ] , x 1 , x 2 ∈ Z + s.t.\ \begin{bmatrix} 1 & 1\\ 1 & 3 \end{bmatrix}\mathbf{x} \le \begin{bmatrix} 70\\ 150 \end{bmatrix}, x_1, x_2 \in Z^{+} s.t. [1113]x≤[70150],x1,x2∈Z+
Python求解代码（基于SCIP solver+cvxpy）
import cvxpy as cpimport numpy as npdef solver_main(): # Coefficient setting mat_Q = np.diag([-0.03, -0.02]) # print('Q= ', mat_Q) mat_c = np.array([5, 10]) # print('c= ', mat_c) mat_A = np.ones((2,2)) mat_A[1,1] = 3 mat_b = np.array([70, 150]) x = cp.Variable(2, integer=True) # Optimization setting objective = cp.Maximize(cp.QuadForm(x, mat_Q) + mat_c.T @ x) constraints = [x>=0, mat_A @ x <= mat_b] prob = cp.Problem(objective, constraints) prob.solve() print('The Optimal objective value: ', prob.value) # Optimal value print('[x_1, x_2]: ' , x.value)if __name__ == '__main__': solver_main()运行脚本，你将看到优化后的结果为：
The Optimal objective value: 491.0000000181662[x_1, x_2]: [30.00000001 40. ]注意：复杂的MIP和MIQP问题在数学领域中是极难求解的，一般求解器求得的都是近似最优解（取决于求解器的求解方法），求解整数规划的求解器主要有：Gurobi，SCIP，CBC，OSQP等等。
参考链接：
https://www.gurobi.com/resources/mixed-integer-programming-mip-a-primer-on-the-basics/
https://www.cvxpy.org/examples/basic/quadratic_program.html
https://blog.csdn.net/weixin_44077955/article/details/126607480
https://blog.csdn.net/weixin_46039719/article/details/121695205