草庐IT

数学规划（Python cvxpy、scipy.optimize）

学习啊ZzZ 2024-05-10 原文

一、线性规划模型 

1.模型结构

①决策变量，x =（x1,x2,x3…,xn）。

②目标函数，f(x) 

③可行域，，常用一组不等式(约束条件)表示：

当目标函数和约束条件对于决策变量而言都是线性的时，称为线性规划

2.模型特征

①比例性，决策变量对目标函数和约束条件的“贡献”，与决策变量的取值成正比

②可加性，决策变量对目标函数和约束条件的“贡献”，与决策变量的取值无关

③连续性，决策变量的取值是连续的

二、线性规划模型求解

（以数学模型第86页模型为例）

1.基本模型

(1) 

(2) 

(3) 

(4) 

(5) 

2.代码求解

import cvxpy as cp
import numpy as np

coef = np.array([72, 64])#输入目标函数系数
left = np.array([[1, 1], [12, 8], [3, 0]])#输入约束条件系数
right = np.array([50, 480, 100])#输入约束条件上限值
x = cp.Variable(2)#构造决策变量
obj = cp.Maximize(coef @ x)#构造目标函数
cons = [x >= 0, left @ x <= right]#构造约束条件
prob = cp.Problem(obj, cons)#构建模型
prob.solve(solver='CPLEX')#模型求解
print("最优值:", prob.value)
print("最优解:", x.value)
print("剩余牛奶：", right[0] - sum(left[0] * x.value))
print("剩余劳动时间：", right[1] - sum(left[1] * x.value))
print("A1剩余加工能力：", right[2] - sum(left[2] * x.value))

3.输出结果

最优值: 3360.0
最优解: [20. 30.]
剩余牛奶： 0.0
剩余劳动时间： 0.0
A1剩余加工能力： 40.0

三、结果分析

1.有效约束

约束条件的右端不妨称为“资源”，目标函数不妨称为“效益”。

若“资源”在最优解下的剩余为0，则称该约束为有效约束。

2.影子价格

当成为有效约束的“资源”增加，“效益”必然会增长。

“资源”增加带来的“效益”的增量可以看作“资源”的潜在价值，称为影子价格

四、复杂线性规划模型举例

（以数学模型第96页模型为例）

1.基本模型

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

2.代码求解

import cvxpy as cp
import numpy as np
#输入目标函数系数
coef = np.array([160, 130, 220, 170,
                 140, 130, 190, 150,
                 190, 200, 230])
#输入约束条件系数
left = np.array([[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
                 [0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0],
                 [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1],
                 [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]])
right_min = np.array([30, 70, 10, 10])#输入约束条件下限值
right_max = np.array([80, 140, 30, 50])#输入约束条件上限值
x = cp.Variable(11)#构造决策变量
obj = cp.Minimize(coef @ x)#构造目标函数
#构造约束条件
cons = [x >= 0,
        left @ x <= right_max,
        left @ x >= right_min,
        cp.sum(x[0:4]) == 50,
        cp.sum(x[4:8]) == 60,
        cp.sum(x[8:11]) == 50]
prob = cp.Problem(obj, cons)#构建模型
prob.solve(solver="CPLEX")#模型求解
print("管理费用最小值为：", prob.value)
print("最优分配方案为：", x.value)

3.输出结果

管理费用最小值为： 24400.0
最优分配方案为： [ 0. 50.  0.  0.  0. 50.  0. 10. 40.  0. 10.]

五、整数规划模型

在线性规划模型中，增加约束条件使得决策变量均为整数，这样得到的新模型就称为整数规划

例一

（以数学模型第101页为例）

1.基本模型

(1) 

(2) 

(3) 

(4) 

(5)  均为整数

2.代码求解

import cvxpy as cp
import numpy as np

coef = np.array([2, 3, 4])#输入目标函数系数
left = np.array([[1.5, 3, 5], [280, 250, 400]])#输入约束条件系数
right = np.array([600, 60000])#输入输入约束条件上限值
x = cp.Variable(3, integer=True)#创建决策变量，并且为整数
obj = cp.Maximize(coef @ x)#构造目标函数
cons = [x >= 0, left @ x <= right]#构造约束条件
prob = cp.Problem(obj, cons)#构建模型
prob.solve(solver="CPLEX")#模型求解
print("最优值：", prob.value)
print("最优解：", x.value)
print("钢材剩余量：", right[0] - sum(left[0] * x.value))
print("劳动时间剩余量：", right[1] - sum(left[1] * x.value))

3.输出结果

最优值： 632.0
最优解： [ 64. 168.   0.]
钢材剩余量： 0.0
劳动时间剩余量： 80.0

例二

（以数学模型第103页解法二为例）

在例一的基础上，增加新的约束条件——“如果生成某一类型汽车，则至少要生产80辆”，选择引入0-1变量，化为整数规划

1.基本模型

(1) 

(2) 

(3) 

(4) 

(5)  均为整数

(6) 

(7) 

(8) 

对于约束条件(6)-(8)，0-1变量y的存在就是为了控制决策变量x的范围，约束条件的上限就是为了保证y等于0时，决策变量x也一定为0。因此M的取值无需对决策变量x产生影响，取值时只需保证足够大即可。

2.代码求解

import cvxpy as cp
import numpy as np

coef = np.array([2, 3, 4])
left = np.array([[1.5, 3, 5], [280, 250, 400]])
right = np.array([600, 60000])
x = cp.Variable(3, integer=True)
y = cp.Variable(3, integer=True)
obj = cp.Maximize(coef @ x)
cons = [x >= 0, left @ x <= right,
        y >= 0, y <= 1,
        x[0] >= 80 * y[0], x[0] <= 1000 * y[0],
        x[1] >= 80 * y[1], x[1] <= 1000 * y[1],
        x[2] >= 80 * y[2], x[2] <= 1000 * y[2], ]
prob = cp.Problem(obj, cons)
prob.solve(solver="CPLEX")
print("最优值：", prob.value)
print("最优解：", x.value)
print("钢材剩余量：", right[0] - sum(left[0] * x.value))
print("劳动时间剩余量：", right[1] - sum(left[1] * x.value))

3.输出结果

最优值： 610.0
最优解： [ 80. 150.   0.]
钢材剩余量： 30.0
劳动时间剩余量： 100.0

六、复杂整数规划模型求解

（以数学模型第107页解法三为例）

该模型的难点在于对分段函数进行变换，使得其可以使用线性规划、整数规划模型加以处理。

1.基本模型

(1)  

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 

(11) 

(12) 

(13) 

约束条件(9)-(11)是对分段函数(1)转换得到的；每个z都是对分段函数每个分点的表示，两个相邻的z可以确定对应段x的取值；每个0-1变量y指的是是否在该分段；条件(12),(13)可用变量z表示x和c(x)。

2.代码求解

import cvxpy as cp
import numpy as np

coef_x = np.array([4.8, 5.6, 4.8, 5.6])#输入目标函数x对应系数
coef_cx = np.array([0, 5000, 9000, 12000])#输入用z表示cx的系数
coef_buy_x = np.array([0, 500, 1000, 1500])#输入用z表示x的系数
left = np.array([[0, 0, 1, 1], [-1, 0, 1, 0], [0, -2, 0, 3]])#输入约束条件系数
right = np.array([1000, 0, 0])#输入约束条件上限值
x = cp.Variable(4)#创建决策变量x
y = cp.Variable(3, integer=True)#创建0-1变量y
z = cp.Variable(4)#创建变量z
obj = cp.Maximize(coef_x @ x - coef_cx @ z)#构造目标函数
#构造约束条件
cons = [cp.sum(x[0:2]) <= 500 + cp.sum(coef_buy_x @ z),
        left @ x <= right,
        sum(coef_buy_x @ z) <= 1500,
        x >= 0,
        z[0] <= y[0], z[1] <= y[0] + y[1], z[2] <= y[1] + y[2], z[3] <= y[2],
        cp.sum(z[:]) == 1, z >= 0,
        cp.sum(y[:]) == 1,
        y >= 0, y <= 1]
prob = cp.Problem(obj, cons)#构造模型
prob.solve(solver="CPLEX")#求解模型
print("最优值：", prob.value)
print("最优解：", x.value)
print("购买原油A：", sum(coef_buy_x * z.value), "t")

3.输出结果

最优值： 5000.0
最优解： [   0. 1500.    0. 1000.]
购买原油A： 1000.0 t

七、0-1规划模型

在线性规划模型中，增加决策变量只能为0或1的约束条件，这样得到的模型就称为0-1规划模型

 （以数学模型第110页为例）

1.基本模型

(1) 

(2) 

(3) 

(4) 

2.代码求解

import cvxpy as cp
import numpy as np
#输入目标函数系数
coef = np.array([66.8, 75.6, 87, 58.6,
                 57.2, 66, 66.4, 53,
                 78, 67.8, 84.6, 59.4,
                 70, 74.2, 69.6, 57.2,
                 67.4, 71, 83.8, 62.4])
x = cp.Variable(20, integer=True)#构造决策变量
#构造目标函数
obj = cp.Minimize(coef @ x)
#输入约束条件
cons = [x >= 0, x <= 1,
        cp.sum(x[0:4]) <= 1,
        cp.sum(x[4:8]) <= 1,
        cp.sum(x[8:12]) <= 1,
        cp.sum(x[12:16]) <= 1,
        cp.sum(x[16:20]) <= 1,
        cp.sum(x[0:20:4]) == 1,
        cp.sum(x[1:20:4]) == 1,
        cp.sum(x[2:20:4]) == 1,
        cp.sum(x[3:20:4]) == 1]
prob = cp.Problem(obj, cons)#构造模型
prob.solve(solver="CPLEX")#模型求解
print("最优值：", prob.value)
print("最优解：", x.value)

3.输出结果

最优值： 253.2
最优解： [ 0. -0. -0.  1.  1. -0. -0.  0. -0.  1. -0. -0. -0. -0.  1.  0. -0. -0. -0. -0.]

八、多目标规划模型

只有一个优化目标的规划问题称为单目标规划，而将多于一个目标的规划问题称为多目标规划。

多目标规划的目标函数可以表示为一个向量：  ,"V-min"是向量最小化的意思，

对于任意最大化目标a，只需在a前加负号即可统一为最小化问题

（以数学模型第113页同学丙为例）

1.基本模型

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 

(11) 

(12) 

(13) 

2.代码求解

import cvxpy as cp
import numpy as np
#输入目标函数的系数
coef_obj = np.array([-0.8, -0.5, -0.5, -0.2, -0.5, -0.2, 0.1, 0.1, -0.2])
coef_credits = np.array([5, 4, 4, 3, 4, 3, 2, 2, 3])#输入课程学分系数
x = cp.Variable(9, integer=True)#构造决策变量
obj = cp.Minimize(coef_obj @ x)#构造目标函数
#输入约束条件
cons = [cp.sum(x[0:5]) >= 2,
        x[2] + [4] + x[5] + x[7] + x[8] >= 3,
        x[3] + x[5] + x[6] + x[8] >= 2,
        2 * x[2] - x[0] - x[1] <= 0,
        x[3] - x[6] <= 0,
        2 * x[4] - x[0] - x[1] <= 0,
        x[5] - x[6] <= 0,
        x[7] - x[4] <= 0,
        2 * x[8] - x[0] - x[2] <= 0,
        x >= 0, x <= 1]
prob = cp.Problem(obj, cons)#模型构建
prob.solve(solver="CPLEX")#模型求解
print("选课结果：", x.value)
print("学分总和：", sum(coef_credits * x.value))

3.输出结果

选课结果： [ 1.  1.  1.  1.  1.  1.  1. -0.  1.]
学分总和： 28.0

九、非线性规划模型

若数学规划模型的目标函数或约束条件对于决策变量而言是非线性的，则称为非线性规划。 

由于CVXPY库处理非线性函数需要区分凸函数、非凸函数等，比较复杂；而SCIPY库可以处理任意非线性规划模型，因此决定用SCIPY求解非线性规划模型

例一

1.基本模型

(1) 

(2) 

(3) 

2.代码求解

from scipy.optimize import minimize
import numpy as np

#输入目标函数
def obj(x):
    return x[0] ** 3 - x[1] ** 3 + x[0] * x[1] + 2 * x[0] ** 2

#输入第一个约束条件
def cons1(x):
    return -x[0] ** 2 - x[1] ** 2 + 6

#输入第二个目标函数
def cons2(x):
    return x[0] * x[1] - 2

#初始化
x0 = np.zeros(2)
x0[0] = 1
x0[1] = 2

bound = (0, 3)#输入决策变量范围
bounds = (bound, bound)#构造决策变量界限
#构造约束条件
cons1 = {"type": "ineq", "fun": cons1}
cons2 = {"type": "eq", "fun": cons2}
cons = ([cons1, cons2])
#构造模型
prob = minimize(obj, x0, method="SLSQP", bounds=bounds, constraints=cons)
#模型求解
x = prob.x
print("最优值：", obj(x))
print("最优解：", x[0], x[1])

 需要注意的是

① 对于“ineq”（即大于等于，小于等于等），约束条件默认的格式是 ，因此需要对模型中约束条件进行标准化再输入

② 初始化这一步骤不可省略

3.输出结果

最优值： -7.785844454002188
最优解： 0.8740320488968545 2.2882456112715115

例二

1.基本模型

(1) 

(2) 

(3) 

(4) 

(5) 

2.代码求解

from scipy.optimize import minimize
import numpy as np

#输入目标函数
def obj(x):
    return 0.5 * x[0] ** 2 + x[1] ** 2 - x[0] * x[1] - 2 * x[0] - 6 * x[1]

#输入第一个约束条件
def cons1(x):
    return 2 - x[0] - x[1]

#输入第二个约束条件
def cons2(x):
    return 2 + x[0] - 2 * x[1]

#输入第三个约束条件
def cons3(x):
    return 3 - 2 * x[0] - x[1]

#初始化
x0 = np.zeros(2)
x0[0] = 0
x0[1] = 1
#输入决策变量范围
bound = (0, None)
#构造决策变量界限
bounds = (bound, bound)
#构造约束条件
cons1 = {"type": "ineq", "fun": cons1}
cons2 = {"type": "ineq", "fun": cons2}
cons3 = {"type": "ineq", "fun": cons3}
cons = ([cons1, cons2, cons3])
#构造模型
prob = minimize(obj, x0, method="SLSQP", bounds=bounds, constraints=cons)
#模型求解
x = prob.x
print("最优值：", obj(x))
print("最优解：", x[0], x[1])

3.输出结果

最优值： -8.222222222222175
最优解： 0.6666666666666612 1.333333333333325

optimizePythonlatex7D7B大数据numpyscipy

有关数学规划（Python cvxpy、scipy.optimize）的更多相关文章

  1. ruby - 我可以在 Ruby 中动态调用数学运算符吗？ - 2

    ruby中有这样的东西吗？send(+,1,2)我想让这段代码看起来不那么冗余ifop=="+"returnarg1+arg2elsifop=="-"returnarg1-arg2elsifop=="*"returnarg1*arg2elsifop=="/"returnarg1/arg2 最佳答案 是的，只需像这样使用send(或者更好的是public_send):arg1.public_send(op,arg2)这是可行的，因为Ruby中的大多数运算符(包括+、-、*、/、andmore)只需调用方法。所以1+2与1.+(2)相同

  2. ruby |设计数学？ - 2

    情况:我正在编写一个程序来求解素数。我需要解决4x^2+y^2=n的问题，其中n是一个已知变量。是的，必须是Ruby。我愿意在这个项目上花费大量时间。我最好自己编写方程式的求解算法，并将其作为该项目的一部分。我真正喜欢的是:如果任何人都可以向我提供指南、网站的链接，或者关于与求解代数方程特别相关的形式算法的构造的歧义消除，或者向我提供似乎你是读者它会帮助我完成任务。请不要建议我使用其他语言。如果您在回答之前接受我真的非常想这样做，我将不胜感激。该项目没有范围或时间限制，也不以营利为目的。这是为了我自己的教育。注意:我并不直接反对为Ruby实现和使用现存的数学库/模块/其他东西，但我更喜

  3. ruby - 我在哪里可以找到 Ruby 中的数学密集型应用程序 - 2

    我发现许多Rails应用程序主要针对企业、社交网络类型的Web应用程序。我看到有人将Ruby与一些出色的OOPS语言(如Java和C#)进行了比较，但我确实发现很难获得一些数学密集型应用程序。非常感谢任何知识渊博的输入(指向示例程序的链接等)，其中轻松显示了语言的用法，就像快速启动或显示该语言如何用于各种数学问题一样。 最佳答案 不幸的是，Ruby并没有在数学和科学计算领域涉足太多。目前，有一个名为SciRuby的pre-alpha库它试图为Ruby带来更多面向数学的功能。他们正试图构建一个NumPy/SciPy等价物。SciRub

  4. 论文解读OTA: Optimal Transport Assignment for Object Detection - 2

    CSDN优秀解读：https://blog.csdn.net/jiaoyangwm/article/details/1266387752021https://arxiv.org/pdf/2103.14259.pdf关键解读在目标检测中标签分配的最新进展主要寻求为每个GT对象独立定义正/负训练样本。在本文中，我们创新性地从全局的角度重新审视标签分配，并提出将分配程序制定为一个最优传输（OT）问题——优化理论中一个被充分研究的课题。具体来说，我们将每个需求方（锚框）和供应商（GT标签）的单位传输成本定义为他们的分类和回归损失加权之和。在公式化后，找到最好的分配方案即为最小传播成本解决最优传输方案，

  5. 【动态规划】背包问题（详细总结，很全） - 2

    【动态规划】一、背包问题1.背包问题总结1）动规四部曲：2）递推公式总结：3）遍历顺序总结：2.01背包1）二维dp数组代码实现2）一维dp数组代码实现3.完全背包代码实现4.多重背包代码实现一、背包问题1.背包问题总结暴力的解法是指数级别的时间复杂度。进而才需要动态规划的解法来进行优化！背包问题是动态规划（DynamicPlanning）里的非常重要的一部分,关于几种常见的背包，其关系如下：在解决背包问题的时候，我们通常都是按照如下五部来逐步分析，把这五部都搞透了，算是对动规来理解深入了。1）动规四部曲：（1）确定dp数组及其下标的含义（2）确定递推公式（3）dp数组的初始化（4）确定遍历顺

  6. 科研中论文常见数学符号及其含义（科研必备，建议收藏） - 2

    论文常见数学符号及其含义（科研必备）返回论文和资料目录数学符号在数学领域是非常重要的。在论文中，使用数学符号可以使得论文更加简洁明了，同时也能够准确地描述各种概念和理论。在本篇博客中，我将介绍一些常见的数学符号及其含义（省去特别简单的符号），希望能够帮助读者更好地理解数学论文。高等数学∑i=1nxi\sum_{i=1}^nx_i∑i=1n​xi​（求和符号）：表示将x1,x2,…,xnx_1,x_2,\dots,x_nx1​,x2​,…,xn​中的所有数相加，例如∑i=1nxi\sum_{i=1}^nx_i∑i=1n​xi​表示将x1,x2,…,xnx_1,x_2,\dots,x_nx1​,x

  7. 中国民用飞机制造行业市场现状规模及发展战略规划报告2021-2027年 - 2

    中国民用飞机制造行业市场现状规模及发展战略规划报告2021-2027年详情内容请咨询鸿晟信合研究院！【全新修订】：2022年2月【撰写单位】：鸿晟信合研究研究【报告目录】第1章：中国民用飞机制造行业发展综述1.1民用飞机制造行业概述1.1.1民用飞机的概念1.1.2飞机制造的概念1.1.3民用飞机的分类1.2民机制造行业周期特性1.2.1影响行业周期的因素（1）GDP增速分析（2）运量增量分析（3）飞机更替分析（4）航空公司获利水平1.2.2行业现阶段周期分析1.2.3行业现阶段景气分析1.3民机制造信息化分析1.3.1信息化技术应用状况分析（1）MDO技术应用分析（2）供应链协同研发分析（3

  8. ruby 数学 gem - 2

    有人知道是否有适合数学/几何运算的好工具吗？我需要一些这样的东西:点积叉积向量之间的角度矩阵表面法线向量计算(加、减...)奇异值分解谢谢，晚礼服 最佳答案 我强烈推荐Ruby/GSLgem。您可能需要引用GSLmanual计算出您需要调用的函数的名称，如果它没有记录在Ruby/GSL文档中的话。 关于ruby数学gem，我们在StackOverflow上找到一个类似的问题： https://stackoverflow.com/questions/541665

  9. 牛客竞赛每日俩题 - 动态规划3 - 2

    目录类01背包问题，选or不选变种走方格类01背包问题，选or不选不同的子序列_牛客题霸_牛客网问题翻译：        S有多少个不同的子串与T相同        S[1:m]中的子串与T[1:n]相同的个数        由S的前m个字符组成的子串与T的前n个字符相同的个数状态：        子状态：由S的前1,2,...,m个字符组成的子串与T的前1,2,...,n个字符相同的个数        F(i,j):S[1:i]中的子串与T[1:j]相同的个数状态递推：        在F(i,j)处需要考虑S[i]=T[j]和S[i]!=T[j]两种情况        当S[i]=T[j]

  10. 数学建模之马尔可夫链模型详解（附详细Matlab程序） - 2

    🔗运行环境：Matlab🚩作者：左手の明天🥇精选专栏：《python》🔥推荐专栏：《算法研究》📚选自专栏：《数学建模》🧿优秀专栏：《Matlab神经网络案例分析》目前持续更新的专栏：🥇专栏：MatlabGUI编程技巧🔥专栏：Matlab从无到有系列大家好，我是左手の明天！今天和大家分享数学建模重要模型——马尔可夫链模型。在对数学建模之马尔可夫链模型进行介绍时，首先需要明确两个问题：马氏链模型用来干什么马尔可夫预测法是应用概率论中马尔可夫链（Markovchain）的理论和方法来研究分析时间序列的变化规律，并由此预测其未来变化趋势的一种预测技术。马氏链模型什么时候用应用马尔可夫链的计算方法进行马

随机推荐