混合整数规划与混合整数二次规划

混合整数规划与混合整数二次规划

引言

在数值优化领域中，混合整数规划（Mixed Integer Programming, MIP）和混合整数二次规划（Mixed Integer Quadratic Programming, MIQP）是两类重要的优化模型。本文将对这两种规划的基本概念、常用求解方法进行介绍，并通过实例分析帮助大家更好地理解这些规划问题。

一、混合整数规划（MIP或MILP）基本概念

混合整数规划问题是指在优化过程中，部分优化变量必须取整数值，其他变量则可以取连续值的一类优化问题。

典型的MIP问题形式如下：

目标方程：	$\min {\mathbf{c}}^T\mathbf{x}$
约束：	$\mathbf{Ax}\le \mathbf{b}$ (线性约束)
	${\mathbf{x}}_l\le \mathbf{x}\le {\mathbf{x}}_h$(边界约束)
	${\mathbf{x}}_i\in Z$, ${\mathbf{x}}_j\in R$ (for some i, others j) （整数约束）

这里，$\mathbf{c}$和$\mathbf{x}$ 是向量，$\mathbf{A}$是矩阵，$\mathbf{b}$是约束条件。换句话说，混合整数规划的目标函数与线性规划不同的是，优化变量的定义域可能是整数，这就给求解带来了极大的难度。

二、混合整数二次规划基本概念

混合整数二次规划是混合整数规划的扩展，其目标函数是一个二次函数，而约束条件通常为一次形式。

典型的MIQP问题形式如下：

目标方程：	$\min {\mathbf{x}}^T\mathbf{Qx}+{\mathbf{c}}^T\mathbf{x}$
约束：	$\mathbf{Ax}\le \mathbf{b}$ (线性约束)
	${\mathbf{x}}_l\le \mathbf{x}\le {\mathbf{x}}_h$(边界约束)
	${\mathbf{x}}_i\in Z$, ${\mathbf{x}}_j\in R$ (for some i, others j) （整数约束）

上述变量定义与MIP完全相同，可以看到MIQP其实就是再MIP的目标函数的基础上，加上了一个二次项。

三、常见求解方法

1. 分支定界法BnB（最常用）

分支定界法是求解MIP和MIQP问题的经典方法之一。其基本思想是将问题分解为一系列子问题，通过构建一个分支树来逐步逼近最优解。分支定界法主要包括三个步骤：分支、定界和剪枝。简单说一下基于LP的BnB方法：

首先，对初始的MILP删除所有整数约束，得到原MILP的线性规划松弛。接下来，求解这个线性规划松弛（LP）。如果得到的解恰好满足所有整数约束，那么这个解就是原始MILP的最优解，运算终止。如果解没有满足所有整数约束（这种情况较为常见），则需要选择某个整数约束实际上得到小数值的变量进行“分支”。为了便于说明，假设这个变量是x，其在LP解中的值是5.7。我们可以通过施加$x\le 5.0$和$x\ge 6.0$的限制来排除该值5.7。

如果用$P$表示原MILP问题，用$P_1$表示新的MILP（增加了$x\le 5.0$的约束），用$P_2$表示新的MILP（增加了$x\ge 6.0$的约束）。变量$x$被称为分支变量，我们在$x$上进行了分支，产生了两个子MILP $P_1$和 $P_2$。显然，如果我们能够求解$P_1$和$P_2$的最优解，那么其中较好的解就是原始问题$P$的最优解。

通过这种方式，我们用两个更简单（更少整数约束）的MILP取代了$P$。我们现在将相同的思想应用于这两个子MILP，并在必要时选择新的分支变量。这样生成了所谓的搜索树。搜索过程生成的MILP称为树的节点，$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$和产品$b$，每种产品都有其生产成本和销售利润，其中：

• 产品$a$的生产成本为10元，销售价格为15元。

• 产品$b$的生产成本为20元，销售价格为30元。

  该工厂的生产能力和原材料有限，假设：

• 每月生产能力为70单位。

• 每月原材料限制为150单位，生产产品$a$消耗1单位原材料，生产产品$b$消耗3单位原材料。

问题一： 如何决定每种产品数量，以实现最大化利润？

• 问题分析：

  从问题中，我们可以看出，待优化的决策变量为产品$a$和产品$b$的数量，用$x_1$和$x_2$表示，其中$x_1$和$x_2$必须取整数值，这个问题可以建模为一个MIP或MILP问题。

  其中，我们的目标函数为：

  $$\max (15-10)x_1+(30-20)x_2$$

  约束为：

  $$x_1+x_2\le 70$$$$x_1+3x_2\le 150$$$$x_1,x_2\in Z^{+}$$

Python求解代码（基于Pulp库）

import pulp as pl
import numpy as np

def 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$，输出如下:

Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019

command 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          MODEL
At line 3 ROWS
At line 7 COLUMNS
At line 18 RHS
At line 21 BOUNDS
At line 24 ENDATA
Problem MODEL has 2 rows, 2 columns and 4 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 550 - 0.00 seconds
Cgl0004I processed model has 2 rows, 2 columns (2 integer (0 of which binary)) and 4 elements
Cutoff increment increased from 1e-05 to 4.9999
Cbc0012I 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 cost
Cuts at root node changed objective from -550 to -550
Probing 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 found

Objective value:                550.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

Shan Status: Optimal
Optimal objective value 550.0
x_1 = 30.0
x_2 = 40.0

2. 考虑收益衰减的生产问题（MIQP问题）

在示例1中的生产问题是一种基本的生产模型，在实际生产过程中，产品的利润有可能会随着产品的数量增加而逐渐降低，呈现收益衰减的现象。我们在示例1的基础上，不妨假设：

• 产品$a$的生产成本为10元，销售价格为$15-0.03x_1$。

• 产品$b$的生产成本为20元，销售价格为$30-0.02x_2$元。

  工厂生产能力和原材料与示例1相同，假设：

• 每月生产能力为70单位。

• 每月原材料限制为150单位，生产产品$a$消耗1单位原材料，生产产品$b$消耗3单位原材料。

问题二： 如何决定每种产品数量，以实现最大化利润？

• 问题分析：

  从问题中，待优化的决策变量依然为产品$a$和产品$b$的数量，用$x_1$和$x_2$表示，其中$x_1$和$x_2$必须取整数值，但是收益模型变得更加复杂了。

  其中，我们的目标函数为：

  $$\max (15-0.03x_1-10)x_1+(30-0.02x_2-20)x_2$$$$=\max (-0.03{x_1}^2-0.02{x_2}^2+5x_1+10x_2)$$

  约束不变，仍为：

  $$x_1+x_2\le 70$$$$x_1+3x_2\le 150$$$$x_1,x_2\in Z^{+}$$

  我们可以看出优化的目标函数为二次形式，因此此任务可以被建模为一个MIQP问题，我们将其转换为矩阵的表示形式：

  $$\max {\mathbf{x}}^T\left[\begin{array}{cc}-0.03& \\ & -0.02\end{array}\right]\mathbf{x}+{\left[\begin{array}{c}5\\ 10\end{array}\right]}^T\mathbf{x}$$$$s.t.\left[\begin{array}{cc}1& 1\\ 1& 3\end{array}\right]\mathbf{x}\le \left[\begin{array}{c}70\\ 150\end{array}\right],x_1,x_2\in Z^{+}$$

Python求解代码（基于SCIP solver+cvxpy）

import cvxpy as cp
import numpy as np

def 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