数学建模——运筹优化类
数学建模的过程，就是把题目“翻译”成数学语言的过程
一组公式，加上对这组公式含义的解释，就是一个数学模型
1 线性规划
1.1 基本概念
线性规划（Lincar programing，LP），是运筹学中研究较早、发展较快、应用广泛、方法较成熟的一个重要分支，是辅助人们进行科学管理的一种数学方法，是研究线性约束条件下线性目标函数的极值问题的数学理论和方法。
三要素：
- 决策变量：问题中要确定的未知量，用于表明规划问题中的用数量表示的方案、措施等，可由决策者决定和控制
- 目标函数：决策变量的函数，优化目标通常是求该函数的最大值或最小值
- 约束条件：决策变量的取值所受到的约束和限制条件，通常用含有决策变量的等式或不等式表示。
建模建立步骤：
- 根据影响所要达到目的的因素找到决策变量
- 由决策变量和所在达到目的之间的函数关系确定目标函数
- 由决策变量所受的限制条件确定决策变量所要满足的约束条件
1.2 模型原理
一般形式：
\max(或\min)\ z=c_1x_1+c_2x_2+\cdots+c_nx_n\\
\text{ s.t. }\left\{\begin{matrix}
a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n≤(或=,≥)\ b_1 \\
a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n≤(或=,≥)\ b_2 \\
\vdots \\
a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n≤(或=,≥)\ b_m\\
x_1,x_2,\cdots,x_n≥0
\end{matrix}\right.简写形式：
\max(或\min)\ z=\sum\limits_{i=1}^nc_jx_j\\
\text{ s.t. }\left\{\begin{matrix}
\sum\limits_{j=1}^na_{ij}x_j≤(或=,≥)\ b_i,\ i=1,2,\cdots,m\\
x_j≥0,\ j=1,2,\cdots,m
\end{matrix}\right.矩阵表现形式：
\max(或\min)\ z=c^\text{T}x\\
\text{ s.t. }\left\{\begin{matrix}
Ax≤(或=,≥)\ b\\
x≥0
\end{matrix}\right.其中
c=[c_1,c_2,\cdots,c_n]^\text{T}为目标函数的系数向量（价值向量）
x=[x_1,x_2,\cdots,x_n]^\text{T}为决策向量
A=(a_{ij})_{m\times n}为约束方程组的系数矩阵
b=[b_1,b_2,\cdots,b_m]^\text{T}为约束方程组的常数向量
模型特点：
（1）要解决的问题是优化类的（即在有限的资源条件下，获取最大的收益）
（2）目标函数和约束条件都是决策变量的线性函数（因此不可为x^2,\text{e}^x,\sin x,\log_2 x等）
（3）在一组线性约束条件下，求线性目标函数的最大值或最小值
手工求解可采用单纯形法，推荐采用MATLAB、Python的相关函数进行求解。
1.3 Python实现：linprog()函数
可使用SciPy库optimize模块中的linprog()函数快速求解如下标准形式的线性规划问题：
\min\limits_{x}cx\\
\text{s.t.}\left\{\begin{matrix}
A_\text{ub}\cdot x≤b_\text{ub}\\
A_\text{eq}\cdot x=b_\text{eq}\\
x \in \text{bounds}
\end{matrix}\right.函数语法：
result = linprog(c, A_ub, n_ub, A_eq, b_eq, bounds, method)c：目标函数的决策变量对应的系数向量（行、列向量皆可，下同）A_ub、b_ub：不等式约束条件的变量系数矩阵和常数项矩阵（必须为≤）A_eq、b_eq：等式约束条件的系数矩阵和常数项矩阵bounds：表示决策变量定义域的n * 2矩阵，None表示无穷（﹣∞或﹢∞）method：调用的求解方法，默认为highs- 返回值有多个参数，常用有如下几个：
x：最优解向量（最常用，通过result.x调用）fun：函数最小值nit：跌打次数
若不满足标准形式，可变号后继续求解。
1.4 典型例题
1.4.2 例1：简单最值问题
求解以下问题
\max y=20x_1+30x_2+45x_3\\
\text{s.t.}\left\{\begin{matrix}
4x_1+8x_2+15x_3≤100\\
x_1+x_2+x_3≤20\\
x_1,x_2,x_3≥0
\end{matrix}\right.对目标函数两边同乘-1，使之满足函数求解的标准形式，即
\min -y=-20x_1-30x_2-45x_3\\
\text{s.t.}\left\{\begin{matrix}
4x_1+8x_2+15x_3≤100\\
x_1+x_2+x_3≤20\\
x_1,x_2,x_3≥0
\end{matrix}\right.求解代码如下
import numpy as np
from scipy.optimize import linprog
c = [-20, -30, -45]
A_ub = [[4, 8, 15], [1, 1, 1]]
b_ub = [100, 20]
bounds = [[0, None], [0, None], [0, None]]
result = linprog(c, A_ub, b_ub, bounds=bounds)
print(result, '\n')
print(result.x, '\n')
y = -result.fun
print(f"max y = {y}") message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: -450.0
x: [ 1.500e+01 5.000e+00 0.000e+00]
nit: 2
lower: residual: [ 1.500e+01 5.000e+00 0.000e+00]
marginals: [ 0.000e+00 0.000e+00 2.500e+00]
upper: residual: [ inf inf inf]
marginals: [ 0.000e+00 0.000e+00 0.000e+00]
eqlin: residual: []
marginals: []
ineqlin: residual: [ 0.000e+00 0.000e+00]
marginals: [-2.500e+00 -1.000e+01]
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
[15. 5. 0.]
max y = 450.01.4.3 典型例题：1998年国赛A题
【例】市场上有n种资产（如股票、债券等）s_i\ (i=1,2,\cdots,n)供投资者选择，某公司有数额为M的一笔相当大的资金可用作一个时期的投资。公司财务分析人员对这n种资产进行了评估，估算出在这一时期内购买资产s_i的平均收益率为r_i，并预测出购买s_i的风险损失率为q_i。考虑到投资越分散，总的风险越小，公司确定，当用这笔资金购买若干种资产时，总体风险可用所投资的s_i中最大的一个风险来度量。
购买s_i要付交易费，费率为p_i，并且当购买额不超过给定值u_i时，交易费按购买u_i计算（不买当然无须付费）。另外，假定同期银行存款利率是r_0（r_0=5\%），且既无交易费又无风险。
已知n=4时的相关数据如表所示
投资收益问题：给上述公司设计投资组合方案，用给定资金M有选择地购买若干种资产或存银行生息，使净收益尽可能大，总体风险尽可能小。
1.4.3.1 问题分析
决策变量：投资不同项目s_i的为x_i\ (i=1,2,\cdots,n)
目标函数：净收益Q尽可能大，总风险尽可能小
约束条件：总资金M有限，每一笔投资都是非负数（隐藏信息）
由已知得，模板函数和约束条件都是决策变量的线性函数，因此可构建线性规划模型。
1.4.3.2 模型假设
可供投资的资金数额M相当大。
投资越分散，总的风险越小，总体风险可用所投资的s_i中最大的风险来度量。
可供选择的n+1种资产（含银行存款）之间是相互独立的。
每种资产可购买的数量为任意值。
在当前投资周期内，r_i,q_i,p_i,u_i\ (i=0,1,\cdots,n)固定不变。
不考虑在资产交易过程中产生的其他费用，如股票交易印花税等。
由于投资数额M相当大，而题目设定的定额u_i相对M很小，p_i,u_i更小，因此假设每一笔交易x_i都大于对应定额u_i。
1.4.3.3 模型建立
总体风险用所投资的s_i中最大的一个风险来衡量，即\max\{q_ix_i|i=1,2,\cdots,n\}。
购买s_i所付交易费本来是一个分段函数，但假设中已经假设每一笔交易x_i都大于u_i，因此交易费为p_ix_i，这样购买s_i的净收益可简化为(r_i-p_i)x_i。
由此可得目标函数（有两个目标）与约束条件为
\left\{\begin{matrix}
\max \sum\limits_{i=0}^n(r_i-p_i)x_i\\
\min\{\max\limits_{1≤i≤n}\{q_ix_i\}\}
\end{matrix}\right.\\
\text{s.t.}\left\{\begin{matrix}
\sum\limits_{i=0}^n(1+p_i)x_i=M\\
x_i≥0,\ i=0,1,\cdots,n
\end{matrix}\right.1.4.3.4 模型简化：减少目标
上述的目标函数属于多目标规划，一般求解策略：将其中一个目标转化为约束条件进行求解。
在实际投资中，投资者承受风险的程度不一样。给定一个风险界限a，使得最大的风险满足
\frac{q_ix_i}{M}≤a由此将第2个目标函数转化为了该约束条件（总体风险小于某个常数），将多目标变成了单目标，即
\max\sum\limits_{i=0}^n(r_i-p_i)x_i\\
\text{s.t.}\left\{\begin{matrix}
\frac{q_ix_i}{M}≤a\ i=1,2,\cdots,n\\
\sum\limits_{i=0}^n(1+p_i)x_i=M\\
x_i≥0,\ i=0,1,\cdots,n
\end{matrix}\right.1.4.3.5 代码实现
对1.3.4中的模型变形为函数标准形式，即
\min\sum\limits_{i=0}^n(p_i-r_i)x_i\\
\text{s.t.}\left\{\begin{matrix}
\frac{q_ix_i}{M}≤a\ i=1,2,\cdots,n\\
\sum\limits_{i=0}^n(1+p_i)x_i=M\\
x_i≥0,\ i=0,1,\cdots,n
\end{matrix}\right.将表中数据代入得
\min f=[0.05,0.27,0.19,0.185,0.185]\cdot[x_0,x_1,x_2,x_3,x_4]^\text{T}\\
\text{ s.t. }\left\{\begin{matrix}
x_0+1.01x_1+1.02x_2+1.045x_3+1.065x_4=M \\
0.025x_1≤aM \\
0.015x_2≤aM \\
0.055x_3≤aM \\
0.026x_4≤aM \\
x_i≥0,\ i=0,1,\cdots,4
\end{matrix}\right.不妨取M=1万元。由于a是任意给定的风险度，这里从a=0开始，以步长\delta a=0.001进行迭代搜索，迭代至a=5%（低风险者能够接受的风险）。
求解代码如下
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import linprog
%matplotlib inline
# 设置matplotlib的参数，使其支持LaTex文本和字体大小
plt.rc('text', usetex=True)
plt.rc('font', size=16)
# 目标函数系数
c = [-0.05, -0.27, -0.19, -0.185, -0.185]
# 线性不等式约束的系数矩阵
# np.c_[arr1, arr2]可以按列（左右）拼接数组arr1和arr2。同理np.r_[arr1, arr2]则是按行（上下）拼接
# np.diag()返回或提取一个对角矩阵
A_ub = np.c_[np.zeros(4), np.diag([0.025, 0.015, 0.055, 0.026])]
# 线性等式约束的系数矩阵和右侧的值
A_eq = [[1, 1.01, 1.02, 1.045, 1.065]]
b_eq = [1]
# 初始化参数和存储结果的列表
a = 0
list_a = []
list_x = []
list_Q = []
# 循环搜索
while a < 0.05:
# 创建线性不等式约束的右侧的值
b_ub = np.ones(4) * a * 1
# 求出该线性规划的解
res = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds=[(0, None), (0, None), (0, None), (0, None), (0, None)])
# 提取解向量和最优值并记录
x = res.x
Q = -res.fun
list_a.append(a)
list_x.append(x)
list_Q.append(Q)
# a自增0.001
a += 0.001
# 绘图
plt.plot(list_a, list_Q, 'r*') # 使用红色星号(*)描点
plt.xlabel('$a$')
plt.ylabel('$Q$', rotation=90) # 顺时针旋转90度
# 可将图像存储到文件中，设置分辨率为500dpi
# plt.savefig('figure_linprog.png', dpi=500)
plt.show()
风险a与收益Q的关系如图所示，可见
（1）风险不超过2.5%时，风险大，收益也大
（2）在a=0.006附近有一个转折点，在左侧的利润变化率大于在右侧的利润变化率。故对于风险和收益没有特殊的偏好者而言，应选择曲线的转折点作为最优投资组合，约为a=0.6%,Q=2000。
获取该点处对应的各资产投资额：
for i, a in enumerate(list_a):
if a == 0.006:
print(list_x[i])
break[0. 0.24 0.4 0.10909091 0.22122066]因此最优投资方案为：风险度a=0.006，收益Q=2019元；x_0=0元，x_1=2400元，x_2=4000元，x_3=1091元，x_4=2212元。
- 本题中做了很多模型假设，理论上来说不做也可以求解，且考虑更全面，但数模比赛时间很紧张，合理的假设也很重要
- 除了固定风险来简化模型，也可以固定收益，或者赋予风险收益相应的权重，来权衡二者的取舍
2 蒙特卡罗法
2.1 基本概念
蒙特卡罗法（Monte Carlo）又称统计模拟法，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。
原理：由大数定理可知，当样本容量足够大时，事件的发生频率即为其概率。
注意：蒙特卡洛不是一种算法，准确的来说只是一种思想，或者是一种方法，只要求解的问题与概率模型有关联就可以采用这种方法，从数学建模的角度来看，蒙特卡洛没有通用的代码，每个问题对应的代码都是不同的。
2.2 典型例题与Python实现
2.2.1 例1：估算圆周率π的值
【例1】用蒙特卡罗法估算圆周率\pi的值
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# 参数初始化
p = 10000 # 总共要投放的点数
r = 1 # 圆的半径
x0, y0 = 1, 1 # 圆心坐标
n = 0 # 圆内点数
# 设置绘图窗口
plt.figure()
plt.title("Monte Carlo Simulation for Estimating Pi")
plt.xlabel("x")
plt.ylabel("y")
# 投放p个点
for i in range(p):
# 随机生成要投放点的横、纵坐标
x_p = np.random.rand() * 2
y_p = np.random.rand() * 2
# 判断点是否在圆内
if (x_p - x0) ** 2 + (y_p - y0) ** 2 < r ** 2:
plt.plot(x_p, y_p, '.', color='b') # 圆内点用蓝色表示
n += 10000
else:
plt.plot(x_p, y_p, '.', color='r') # 圆外点用红色表示
plt.axis("equal")
plt.show()
# 计算pi的估计值
s = (n / p) * 4 # 计算圆面积占正方形面积的比例
pi_estimated = s
print(f"Estimated value of pi:{pi_estimated}")
Estimated value of pi: 31248.02.2.2 三门问题
【例2】你参加一档电视节目，节目组提供了A、B、C三扇门，主持人告诉你，其中一扇门后边有辆汽车，其他两扇门后面是一头山羊，你可以选择一扇门打开获得门后的东西。假如你选择了B们，这时，主持人打开了C门，让你看到C门后是只山羊，问你要不要改选A门？（你想要汽车）
法1：只考虑在成功条件下的概率
import numpy as np
n = 100000 # 模拟重复次数
a = 0 # 不改变主意时能赢得汽车的次数
b = 0 # 改变主意时能赢得汽车的次数
for i in range(n):
x = np.random.randint(1, 4) # 随机生成一个1-3之间的整数x表示汽车出现在第x扇门之后
y = np.random.randint(1, 4) # 随机生成一个1-3之间的整数y表示自己选的门
# 分两种情况讨论：x == y 和 x != y
if x == y: # 若x和y相同，则只有不改变主意时才能赢
a += 1
else: # 若x和y不同，则只有改变主意时才能赢
b += 1
print(f"不改变主意时能赢得汽车的概率: {a / n}")
print(f"改变主意时能赢得汽车的概率: {b / n}")不改变主意时能赢得汽车的概率: 0.33555
改变主意时能赢得汽车的概率: 0.66445法2：考虑失败情况下的概率（无条件概率）
import numpy as np
n = 100000 # 模拟重复次数
a = 0 # 不改变主意时能赢得汽车的次数
b = 0 # 改变主意时能赢得汽车的次数
c = 0 # 新增变量c，表示没有获奖的次数
for i in range(n):
x = np.random.randint(1, 4) # 随机生成一个1-3之间的整数x表示汽车出现在第x扇门之后
y = np.random.randint(1, 4) # 随机生成一个1-3之间的整数y表示自己选的门
change = np.random.randint(0, 2) # 是否改变主意也是随机的，0：不改变主意，1：改变主意
# 分两种情况讨论 x == y 和 x != y，
if x == y: # 若x和y相同，则只有不改变主意时才能赢
if change == 0: # 不改变主意
a += 1
else: # 改变主意
c += 1
else: # 若x和y不同，则只有改变主意时才能赢
if change == 0: # 不改变主意
c += 1
else: # 改变主意
b += 1
print(f"不改变主意时能赢得汽车的概率: {a / n}")
print(f"改变主意时能赢得汽车的概率: {b / n}")
print(f"没有获奖的概率: {c / n}")不改变主意时能赢得汽车的概率: 0.16644
改变主意时能赢得汽车的概率: 0.3331
没有获奖的概率: 0.50046分析原因：改变/不改变主意时的获奖率实际取决于第一次自己选中的事物——由于主持人排除的是山羊选项，要想在不改变主意时赢得汽车，需保证自己第一次选择时即选中汽车，概率为1/3；要想在改变主意时能赢得汽车，只需保证自己第一次选择时选中山羊，概率为2/3。
3 非线性规划
3.1 基本概念
非线性规划是一种求解目标函数或约束条件中有一个或几个非线性函数的最优化问题的方法，是运筹学的一个重要分支。20世纪50年代初，库哈（H.W.Kuhn）和托克（A.W. Tucker）提出了非线性规划的基本定理，为非线性规划奠定了理论基础。这一方法在工业、交通运输、经济管理和军事等方面有广泛的应用，特别是在“最优设计”方面，它提供了数学基础和计算方法，因此有重要的实用价值。
模型特点：
（1）至少一个变量是非线性，即包含x^2,\text{e}^x,\frac1x,\sin x, \log_2x等
（2）目前非线性规划还没有适合各种问题的一般解法，各种方法都有其特定的适用范围
3.2 模型原理
非线性规划模型的标准形式
\min f(x)\\
\text{ s.t. }\left\{\begin{matrix}
Ax≤b,\ A_\text{eq}\cdot x=b_\text{eq} &（线性） \\
c(x)≤0,\ C_\text{eq}(x)=0 & （非线性） \\
\text{lb}≤x≤\text{ub}
\end{matrix}\right.3.3 Python实现：minimizer()函数
可以使用scipy.optimizer.minimize函数求解
minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)- fun: 目标函数，即需要最小化的函数。该函数的参数应该包括待优化的变量。
- x0: 变量的初始猜测值。
- args: 传递给目标函数的额外参数，以元组的形式传递。
- method: 优化算法的选择。可以是字符串，表示使用的具体算法，例如 'BFGS'、'L-BFGS-B'、'COBYLA' 等。也可以是一个优化算法的对象。
- jac: 目标函数的雅可比矩阵（梯度），如果提供，可以加速某些优化算法的收敛。
- bounds: 变量的边界，以元组的形式传递，有n个参数就由多个元组组成((x0,y0), (x1, y1), ...)。
- constraints: 约束条件，可以是等式约束或不等式约束，以字典的形式传递。
- options: 一个字典，包含优化过程的其他控制参数。
3.4 典型例题
3.4.1 例1：优化二次函数
import numpy as np
from scipy.optimize import minimize
# 目标函数
def quadratic_function(x):
return x[0]**2 + x[1]**2 + x[0]*x[1]
# 初始猜测值
x0 = np.array([1.0, 1.0])
# 最小化目标函数
result = minimize(quadratic_function, x0, method='BFGS')
# 输出结果
print(result) message: Optimization terminated successfully.
success: True
status: 0
fun: 1.5869182079705773e-16
x: [-7.273e-09 -7.273e-09]
nit: 2
jac: [-6.918e-09 -6.918e-09]
hess_inv: [[ 6.667e-01 -3.333e-01]
[-3.333e-01 6.667e-01]]
nfev: 9
njev: 33.4.2 例2：选址问题
临时料场：A(5,1)、A(2,7)；日储量各20吨
（1）试制定每天的供应计划，即从两料场分别向各工地运送多少吨水泥，使总的吨千米数最小？
（2）为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各20吨，问应建在何处，节省的吨千米数力多大？
3.4.2.1 确定决策变量
设第i个工地的坐标(a_i,b_i)，水泥日用量d_i\ (i=1,2,\cdots.6)；料场位置(x_j,y_j)，日储量e_j\ (j=1,2)；从料场j向工地i的运送量为x_{ij}
3.4.2.2 确定约束条件
（1）料场水泥运输总量不超过其日储量：
\sum\limits_{i=1}^6 x_{ij}≤e_j,\ j=1,2（2）两个料场向某工地运输量之和等于该工地水泥日用量：
\sum\limits_{j=1}^2x_{ij}=d_i,\ i=1,2,\cdots,63.4.2.3 确定目标函数
求总吨干米数最小，即运送量乘运送距离求和最小
\min f=\sum\limits_{j=1}^2\sum\limits_{i=1}^6x_{ij}\sqrt{(x_j-a_i)^2+(y_j-b_i)^2}3.4.2.4 建立模型
（1）对于第1问，因料场位置已知，故决策变量仅为x_{ij}，为线性规划模型
（2）对于第2问，新料场位置未知，所以x_j和y_j均为变量，且不是线性的，故为非线性规划模型
综上，共有8个约束。故模型为
\min f=\sum\limits_{j=1}^2\sum\limits_{i=1}^6x_{ij}\sqrt{(x_j-a_i)^2+(y_j-b_i)^2}\\
\text{ s.t. }\left\{\begin{matrix}
\sum\limits_{i=1}^6x_{ij}≤e_j \\
\sum\limits_{j=1}^2x_{ij}=d_{i} \\
x_{ij}≥0
\end{matrix}\right.\ ,\ i=1,2,\cdots,6,\ j=1,23.4.2.5 第1问答案
# 线性规划
import numpy as np
from scipy.optimize import linprog
# 6个工地坐标
a = np.array([1.25, 8.75, 0.5, 5.75, 3, 7.25])
b = np.array([1.25, 0.75, 4.75, 5, 6.5, 7.75])
# 临时料场位置
x = np.array([5, 2])
y = np.array([1, 7])
# 6个工地水泥日用量
d = np.array([3, 5, 4, 7, 6, 11])
# 计算目标函数系数
l = np.zeros((6, 2))
for i in range(6):
for j in range(2):
l[i, j] = np.sqrt((x[j] - a[i]) ** 2 + (y[j] - b[i]) ** 2)
f = np.concatenate((l[:, 0], l[:, 1]))
# 不等式约束条件
A = np.array([[1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]])
b = np.array([20, 20])
# 等式约束条件
A_eq = np.hstack((np.eye(6), np.eye(6)))
b_eq = d
# 变量下界
V_lb = np.zeros(12)
# 求解线性规划问题
result = linprog(f, A_ub=A, b_ub=b, A_eq=A_eq, b_eq=b_eq, bounds=(0, None))
x = result.x
fval = result.fun
print("线性规划结果:", x)
print("最小目标函数值:", fval)TBC