第三章 数学规划模型
确定优化模型要确定优化的目标和寻求的决策，可以用x来表示为决策变量，f(x)表示目标函数。实际问题一般对决策变量x的取值范围有限制，不妨记作：x\in\Omega，\Omega称为可行域。，优化问题的数学模型可表示为
\min(或\max)f(x),x\in\Omega
在第三章，x通常是1维或2维变量，\Omega通常是1维和2维的非负域。
实际中的优化问题通常有多个决策变量，用n维向量\mathbf x=(x_1,x_2,...,x_n)^T表示，目标函数f(\mathbf x)是多元函数，可行域\Omega比较复杂，常用一组不等式（也可以有等式）g_i(\mathbf x)\le0(i=1,2,...,m)来界定，称为约束条件。一般地，这类模型可表述成如下形式：
\min_xz=f(\mathbf x)\\s.t.~~g_i(\mathbf x)\le0,i=1,2,..,m
上述模型属于多元函数的条件极值问题的范围。然而许多实际问题归结出的这种形式的优化模型，其决策变量个数n和约束条件个数m一般较大，并且最优解往往在可行域的边界中取得，这样就不能简单的用微分法求解，数学规划就是解决这类问题的有效方法。
线性规划
目标函数和约束条件对于决策变量而言都是线性的，称为线性规划。
例题 奶制品的生产和销售
问题提出
奶制品加工厂用牛奶生产A_1,A_2两种奶制品，1桶牛奶可以在甲类设备上用12h加工成3kgA_1，或者在乙类设备上用8h加工成4kgA_2。根据市场需求，生产的A_1,A_2全部能售出，且每千克A_1，获利24元，每千克A_2，获利16元.现在加工厂每天能得到50桶牛奶的供应，每天正式工人总的劳动时间为480h,并且甲类设备每天至多能加工100kgA_1,乙类设备的加工能力没有限制.试为该厂制订一个生产计划，使每天获利最大，并进一步讨论以下3个附加问题：
- 若用35元可以买到1桶牛奶，应否作这项投资？若投资，每天最多购买多少桶牛奶？
- 若可以聘用临时工人以增加劳动时间，付给临时工人的工资最多是每小时几元？
- 由于市场需求变化，每千克A_1,的获利增加到30元，应否改变生产计划？
问题分析
- 决策变量：每天用x_1桶牛奶生产A_1，x_2桶牛奶生产A_2。
- 目标函数：设每天获利z元，总获利z=72x_1+64x_2。
- 约束条件：
- 原料供应：x_1+x_2\le50
- 劳动时间：12x_1+8x_2\le480
- 设备能力：3x_1\le 100
- 非负约束：x_1\ge0,x_2\ge0
综上可得：
\begin{align} &\max z = 72x_1+64x_2\\ s.t.~~~&x_1+x_2\le50\\ &12x_1+8x_2\le480\\ &3x_1\le 100\\ &x_1\ge0,x_2\ge0 \end{align}
模型求解
图解：
根据高中知识易知：取B点(20,30)时，z取得最大值。
软件实现：
import numpy as np
from scipy import optimize as op
np.set_printoptions(suppress=True)
# 给出变量取值范围
x1 = (0, None)
x2 = (0, None)
c = np.array([-72, -64]) # 目标函数系数，原题取最大值，目标函数两端加负号
A_ub = np.array([[1,1], [12,8], [3,0]]) # 不等式约束系数A
B_ub = np.array([50, 480, 100]) # 不等式约束系数b
res = op.linprog(c, A_ub, B_ub, bounds=(x1, x2)) # 调用函数求解
res
得到结果：
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: -3360.0
x: [ 2.000e+01 3.000e+01]
nit: 2
lower: residual: [ 2.000e+01 3.000e+01]
marginals: [ 0.000e+00 0.000e+00]
upper: residual: [ inf inf]
marginals: [ 0.000e+00 0.000e+00]
eqlin: residual: []
marginals: []
ineqlin: residual: [ 0.000e+00 0.000e+00 4.000e+01]
marginals: [-4.800e+01 -2.000e+00 -0.000e+00]
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
由fun可得-z最小值是-3360，故z最大值是3360，取x为[20，30]。
Lingo程序：
MODEL:
max = 72*x1+64*x2;
[milk] x1+x2<50;
[time] 12*x1+8*x2<480;
[cpct] 3*x1<100;
END
结果：
Global optimal solution found.
Objective value: 3360.000
Infeasibilities: 0.000000
Total solver iterations: 2
Elapsed runtime seconds: 0.09
Model Class: LP
Total variables: 2
Nonlinear variables: 0
Integer variables: 0
Total constraints: 4
Nonlinear constraints: 0
Total nonzeros: 7
Nonlinear nonzeros: 0
Variable Value Reduced Cost
X1 20.00000 0.000000
X2 30.00000 0.000000
Row Slack or Surplus Dual Price
1 3360.000 1.000000
MILK 0.000000 48.00000
TIME 0.000000 2.000000
CPCT 40.00000 0.000000
结论与上述相符。
整数规划
决策变量为整数的线性规划称为整数规划。
例题 汽车生产
问题
一汽车厂生产小、中、大三种型的汽车，已知各类型每辆车对钢材、劳动时间的需求，利润以及每月工厂钢材、劳动时间的现有量如表1所示.试制订月生产计划，使工厂的利润最大。
| 小型 | 中型 | 大型 | 现有量 | |
|---|---|---|---|---|
| 钢材/t | 1.5 | 3 | 5 | 600 |
| 劳动时间/h | 280 | 250 | 400 | 60000 |
| 利润/万元 | 2 | 3 | 4 |
很容易得到以下线性规划模型：
\begin{aligned} &\text{max} ~z =2x_{1}+3x_{2}+4x_{3} \\ \text{s.t.} ~~~&1.5x_{1}+3x_{2}+5x_{3}\leqslant600 \\ &280x_1+250x_2+400x_3\leqslant60000 \\ &x_{1},x_{2},x_{3}\geq0\\ &x_1,x_2,x_3均为整数 \end{aligned}
向这样得到的式子称为整数规划。
Python求解：
import numpy as np
import cvxpy as cp
n=3 #设置目标函数中变量个数
c=np.array([2,3,4]) #输入目标函数的系数
a=np.array([[1.5,3,5],[280,250,400]]) #输入约束条件的系数矩阵
b=np.array([600,60000]) #输入b值
x=cp.Variable(n,integer=True) #创建x，个数是3，类型是整数
objective=cp.Maximize(cp.sum(c*x)) #明确目标函数（此时c是3×1，x是3×1,但python里面可以相乘）
constriants=[0<=x,a*x<=b] #明确约束条件
prob=cp.Problem(objective,constriants) #明确问题
res=prob.solve(solver=cp.CPLEX) #求解问题
print(prob.value) #输出目标函数的值
print(x.value) #输出x的值
# 运行结果为：
# 632.0
# [ 64. 168. 0.]
可知：最优解x_1=64,x_2=168,x_3=0，最优值：z=632。
LINGO程序：
MODEL:
max = 2*x1 + 3*x2 + 4*x3;
[STELL] 1.5*x1 + 3*x2 + 5*x3 <= 600;
[TIME] 280*x1 + 250*x2 + 400*x3 <= 60000;
@GIN(x1);
@GIN(x2);
@GIN(x3);
END
Global optimal solution found.
Objective value: 632.0000
Objective bound: 632.0000
Infeasibilities: 0.000000
Extended solver steps: 0
Total solver iterations: 4
Elapsed runtime seconds: 0.12
Model Class: PILP
Total variables: 3
Nonlinear variables: 0
Integer variables: 3
Total constraints: 3
Nonlinear constraints: 0
Total nonzeros: 9
Nonlinear nonzeros: 0
Variable Value Reduced Cost
X1 64.00000 -2.000000
X2 168.0000 -3.000000
X3 0.000000 -4.000000
Row Slack or Surplus Dual Price
1 632.0000 1.000000
STELL 0.000000 0.000000
TIME 80.00000 0.000000
非线性规划
若规划问题的目标函数或约束条件中包含非线性函数，则称为非线性规划。目前非线性规划还没有适合各种问题的一般解法，各种方法都有其特定的适用范围。
例题 原油采购
问题某公司用两种原油(A和B)混合加工成两种汽油（甲和乙）.甲、乙两种汽油含原油A的最低比例分别为50%和60%，售价分别为4800元/t和5600元/t.该公司现有原油A和B的库存量分别为500t和1000t,还可以从市场上买到不超过1500t的原油A.原油A的市场价为：购买量不超过500t时的单价为10000元/：购买量超过500t但不超过1000t时，超过500t的部分8000元/t:购买量超过1000t时，超过 1000t的部分6000元/L.该公司应如何安排原油的采购和加工？
模型建立
设原油A的购买量为x，根据题目可以读出数据，采购的支出c(x)可以表示为如下的线性分段函数：
c\left(x\right)=\left\{\begin{array}{ll}10x,&\quad0\leqslant x\leqslant500\\1000+8x,&\quad500\leq x\leq1000\\3000+6x,&\quad1000\leq x\leq1500\end{array}\right.
总利润为（目标函数）：
\mathrm{max}z=4.8\left(x_{_{11}}+x_{_{21}}\right)+5.6\left(x_{_{12}}+x_{_{22}}\right)-c\left(x\right)
约束条件为：
\begin{aligned} &x_{11}+x_{12}\leq500+x \\ &x_{21}+x_{22}\leqslant1000 \\ &x\leqslant1500 \\ &\frac{x_{11}}{x_{_{11}}+x_{_{21}}}\geqslant0.5\\ &\frac{x_{12}}{x_{12}+x_{22}}\geqslant0.6\\&x_{11},x_{12},x_{21},x_{22},x\geqslant0 \end{aligned}
由于c(x)不是一个线性函数，则上述模型给出的是一个非线性规划。下面给出一种解法：
将原油A的采购分解为三个量，x_1,x_2,x_3。
此时目标函数变为：
\mathrm{max}z=4.8\left(x_{11}+x_{21}\right)+5.6\left(x_{12}+x_{22}\right)-\left(10x_{1}+8x_{2}+6x_{3}\right)
只有当
应该注意到，只有当以10千元/t：的价格购买x_1=500t时，才能以8千元/t的价格购买x_2(x_2>0),这个条件可以表示为：
\left(x_{1}-500\right)x_{2}=0
同理：
\left(x_{2}-500\right)x_{3}=0
此外还有非负的约束条件：
0\leqslant x_{1},x_{2},x_{3}\leqslant500
由于有非线性约束，则构成了非线性规划模型。
MODEL:
max = 4.8*(x11+x21) + 5.6 *(x12 + x22) - (10*x1 + 8*x2+ 6*x3);
x = x1 + x2 + x3;
x11 + x12 <= 500 + x;
x21 + x22 <= 1000;
x <= 1500;
x11 / (x11 + x21)>=0.5;
x12 / (x12 + x22)>=0.6;
x11 > 0;
x12 > 0;
x21 > 0;
x22 > 0;
x > 0;
(x1 - 500)*x2 = 0;
(x2 - 500)*x3 = 0;
x1 <= 500;
x2 <= 500;
x3 <= 500;
x1 >= 0;
x2 >= 0;
x3 >= 0;
END
Feasible solution found.
Objective value: 5000.000
Infeasibilities: 0.000000
Total solver iterations: 97
Elapsed runtime seconds: 0.13
Model Class: NLP
Total variables: 8
Nonlinear variables: 7
Integer variables: 0
Total constraints: 20
Nonlinear constraints: 4
Total nonzeros: 36
Nonlinear nonzeros: 8
Variable Value Reduced Cost
X11 0.000000 0.000000
X21 0.000000 1.600002
X12 1500.000 0.000000
X22 1000.000 0.000000
X1 500.0000 0.000000
X2 500.0000 2.400000
X3 0.2322321E-03 0.4000000
X 1000.000 0.000000
Row Slack or Surplus Dual Price
1 5000.000 1.000000
2 0.000000 5.600000
3 0.000000 5.600000
4 0.000000 5.600000
5 499.9998 0.000000
6 0.6484883E-06 0.000000
7 0.3715700E-07 0.000000
8 0.000000 0.000000
9 1500.000 0.000000
10 0.000000 0.000000
11 1000.000 0.000000
12 1000.000 0.000000
13 0.000000 -0.8800000E-02
14 0.000000 0.000000
15 0.000000 0.000000
16 0.000000 0.000000
17 499.9998 0.000000
18 500.0000 0.000000
19 500.0000 0.000000
20 0.2322321E-03 0.000000
购买1000tA与库存的500tA和1000tB一起生产2500t乙，利润为5000000元。
0-1 规划
0-1型整数规划是整数规划的一种特殊形式，它的变量x_i仅取值0或1。这种只能取0或1的变量称为0-1变量或二进制变量。
例题 选课策略
某学校规定，运筹学专业的学生毕业时必须至少学习过两门数学课、三门运筹学课和两门计算机课.这些课程的编号、名称、学分、所属类别和先修课要求如表所示.那么，毕业时学生最少可以学习这些课程中的哪些课程？
| 课程编号 | 课程名称 | 学分 | 所属类别 | 先修课要求 |
|---|---|---|---|---|
| 1 | 微积分 | 5 | 数学 | |
| 2 | 线性代数 | 4 | 数学 | |
| 3 | 最优化方法 | 4 | 数学;运筹学 | 微积分;线性代数 |
| 4 | 数据结构 | 3 | 数学;计算机 | 计算机编程 |
| 5 | 应用统计 | 4 | 数学;运筹学 | 微积分;线性代数 |
| 6 | 计算机模拟 | 3 | 计算机;运筹学 | 计算机编程 |
| 7 | 计算机编程 | 2 | 计算机 | |
| 8 | 预测理论 | 2 | 运筹学 | 应用统计 |
| 9 | 数学实验 | 3 | 运筹学;计算机 | 微积分;线性代数 |
模型建立
用x_1来表示表中按序号顺序的9们课程，问题的目标为选修的课程总数最少，即：
\min z=\sum_{i=1}^{9}x_{i}
约束条件包含两个方面：
- 每人最少要学习2门数学课、3门运筹学课和2门计算机课:
\begin{aligned} x_{1}+x_{2}+x_{3}+x_{4}+x_{5}\geq2 \\ x_{3}+x_{5}+x_{6}+x_{8}+x_{9}\geq3 \\ x_{4}+x_{6}+x_{7}+x_{9}\geq2 \end{aligned}
-
某些课程有先修课的要求。这个条件可以表示为：x_i\le x_j，表示在修完x_j后才能修x_i。
\begin{gathered} 2x_{3}-x_{1}-x_{2}\leq0 \\ x_{4}-x_{7}\leq0 \\ 2x_{5}-x_{1}-x_{2}\leq0 \\ x_{6}-x_{7}\leq0 \\ x_{8}-x_{5}\leq0 \\ 2x_{9}-x_{1}-x_{2}\leq0 \end{gathered}
这样就得到了需要求解的模型，使用Lingo求解得：
MODEL:
min = x1+x2+x3+x4+x5+x6+x7+x8+x9;
x1+x2+x3+x4+x5>=2;
x3+x5+x6+x8+x9>=3;
x4+x6+x7+x9>=2;
2*x3-x1-x2<=0;
x4-x7<=0;
2*x5-x1-x2<=0;
x6-x7<=0;
x8-x5<=0;
2*x9-x1-x2<=0;
@BIN(x1);@BIN(x2);@BIN(x3);@BIN(x4);@BIN(x5);@BIN(x6);@BIN(x7);@BIN(x8);@BIN(x9);
END
Variable Value Reduced Cost
X1 1.000000 1.000000
X2 1.000000 1.000000
X3 1.000000 1.000000
X4 0.000000 1.000000
X5 0.000000 1.000000
X6 1.000000 1.000000
X7 1.000000 1.000000
X8 0.000000 1.000000
X9 1.000000 1.000000
选编号为1、2、3、6、7、9的课程即可。
多目标规划
多目标规划是含有多个目标函数的规划问题。其要在多个具有相同约束条件优化的目标当中寻找合适的变量，使得多个目标集合的整体能够达到最优的优化问题。
在求解多目标，我们需要确定不同目标的权重，来确定新的目标。此时我们的约束条件可以不做改变。
例题 选课策略之二
例题如上述例题：选课策略，这里讨论：
讨论如果一个学生既希望选修课程数少，又希望所获得的学分数尽可能多，则除了前一个目标函数之外，还应根据表中的学分数写出另一个目标，即：
\max w=5x_{1}+4x_{2}+4x_{3}+3x_{4}+4x_{5}+3x_{6}+2x_{7}+2x_{8}+3x_{9}
对上述两个目标可以表示为对一个向量进行优化：
V-\min(z,-w)
要得到多目标规划问题的解，常常需要知道决策者对每个目标的重视程度，称为**偏好程度。**下面给出几个例子来讨论处理这类问题的方法。
-
同学甲只考虑获得尽可能多的学分，而不管所修课程的多少，那么他可以仅仅以上面的目标进行优化。这就变成了一个单目标优化问题.显然，这个问题不必计算就知道最优解是选修所有9门课程。
-
同学乙认为选修课程数最少是基本的前提，那么他可以只考虑之前目标而不管现在的目标，这就是前面得到的，最少为6门.如果这个解是唯一的，则他已别无选择，只能选修上面的6门课，总学分为21。但是LING0无法告诉我们一个优化问题的解是否唯一，所以他还可能在选修6门课的条件下，使总学分多于21.为探索这种可能，应在上面的规划问题中增加约束：
\sum_{j=1}^{9}x_{i}=6
然后再用刚才的目标进行求解计算，可以的得到最优解为：选取1、2、3、5、7、9。总学分是22。
-
同学丙不像甲、乙那样，只考虑学分最多或以课程最少为前提，而是觉得学分数和课程数这两个目标大致上应该三七开.这时可以将目标函数z和-w分别乘以0.7和0.3,组成一个新的目标函数y,有：
\begin{aligned}\min y&=0.7z-0.3w\\&=-0.8x_1-0.5x_2-0.5x_3-0.2x_4-0.5x_5-0.2x_6+0.1x_7+0.1x_8-0.2x_9\end{aligned}
将这个式子宗伟目标，再用LINGO进行求解，求解得到结果为：1,2,3,4,5,6,7,9。共28学分。
动态规划
动态规划（英语：Dynamic programming，简称 DP）是一种把原问题分解为相对简单的子问题的方式求解复杂问题的方法。
动态规划的4个组成部分：
- 最后一步和子问题
- 状态转移方程
- 初始条件和边界情况
- 计算顺序
下面会用一道简单的例题来演示动态规划的步骤：
例题 最佳硬币数
假设你有三种硬币，面值为2元、5元和7元，每种硬币都足够多。而买一本书需要27元，那么如何用最少的硬币数买到这本书且恰好付清不用找钱？
假设x是我们需要交付的钱数，f(x)是硬币数。
在解决这个问题前，我们需要确定最后一步的问题，假设我们最后一枚硬币的面值是x_k，且f(27)=k。表示最少我们需要k枚硬币来付款27元。那么最后一步就可以表示成：
f(x-x_k)=k-1
表示我们在交付最后一枚硬币前（k-1），已经能够交付价值为金额的钱。
接下来我们就能确定状态转移方程，由于我们想要的是最少的硬币数，而且我们的硬币面值只有2,5,7。那么状态转移方程为：
f(x)=\min\{f(x-2)+1,f(x-5)+1,f(x-7)+1\}
和线性规划一样，动态规划也需要边界条件。根据题意可以知道：
f(0)=0\\ f(x)=\infty,~~x<0
这些就是我们求解的边界条件。
接下来确定计算顺序，我们依次从x=1计算至x=27。最终可以求得最终解。
# 最少硬币支付，动态规划
def min_coin(k):
dp = [float('inf')] * (k+1) # 初始化dp数组
dp[0] = 0 # 边界条件
for i in range(1, k+1):
dp[i] = min(dp[i-2]+1, dp[i-5]+1, dp[i-7]+1) # 状态转移方程
return dp[k]
if __name__ == '__main__':
print(min_coin(27))
求得最终解为5，即用5，5，5，5，7支付。
例题 背包问题
一共有N件物品，第i（i从1开始）件物品的重量为w_i，价值为v_i。在总重量不超过背包承载上限W的情况下，能够装入背包的最大价值是多少？
设f(x,w)为装进前x件物品时重量为w时能装进背包的最大价值。最后一步分析为：
如果最后选择装东西：
f(x-1,w-w_k)=v-v_k
解释为，在装进最后一个东西前，背包里的最大价值为v-v_k。
如果最后没有选择装东西(装不下)：
f(x-1,w)=v=f(x,w)
我们想要寻求装载东西价值的最大值，于是我们可以写出状态转移方程：
f(x,w)=\max\{f(x-1,w),f(x-1,w-w_i)+v_i\}
下面寻找边界条件：
f(0,0)=0\\f(0,w)=0\\f(x,0)=0
接下来就是确定计算顺序，经过分析得知，x可以从1-N进行枚举，w需要从W到0逆向枚举。
此外，这道题还可以用0-1规划。
评论区