当前位置:首页 » 《休闲阅读》 » 正文

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

13 人参与  2024年09月12日 13:21  分类 : 《休闲阅读》  评论

点击全文阅读


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

引言

在数值优化领域中,混合整数规划(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.0

2. 考虑收益衰减的生产问题(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.03x1​2−0.02x2​2+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. [11​13​]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问题在数学领域中是极难求解的,一般求解器求得的都是近似最优解(取决于求解器的求解方法),求解整数规划的求解器主要有:GurobiSCIPCBC,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


点击全文阅读


本文链接:http://m.zhangshiyu.com/post/158757.html

<< 上一篇 下一篇 >>

  • 评论(0)
  • 赞助本站

◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。

关于我们 | 我要投稿 | 免责申明

Copyright © 2020-2022 ZhangShiYu.com Rights Reserved.豫ICP备2022013469号-1