碎念系列是对相关知识非结构化的记录，在这里不会考虑知识梯度和文章结构，而更像是基于笔者水平对相关知识的一些讨论的言语记录，或称碎碎念. 因本人水平有限，难免有错漏或不足之处，如有批评指正请联系penguinpi@163.com，同时也欢迎读者一同参与讨论，谢谢.
本篇讨论笔者对数学规划模型的理解，以及在一些数学建模竞赛中的应用心得，最后提供一些简单代码以便快速查阅.
方法思想
- 我们最原始的想法应该是：欲达到目标，发现此目标受某些变量影响，故寻找这些变量该如何取定，致使目标实现.
- 若采用数学的方法研究，则将目标写成依赖于变量的函数形式，将目标量化，求函数极小值以及变量取值. 由此衍生出一系列的专业名词：目标函数、决策变量、可行域、约束条件等. 更抽象地，我们可以说在一些约束条件下，找到目标函数的极小值和它的解，这便是数学规划研究的范畴.
- 所谓的数学规划模型，则是将现实世界的目标和变量合理地量化并写成数学形式，采用数学规划的手段优化目标函数，从而等价地实现目标（换句话说，也就是将极小化目标函数等价为目标实现，函数取极小值时，变量的取值映射为现实中的变量相应的取值）.
- 一些简单的规划问题很容易求出最小值，而一些复杂的则只能求出极小值，也就是局部最优解.
- 基于规划问题的具体形式，我们又可以将其简单地归纳为线性规划、整数规划、非线性规划、多目标规划.
- 关于非线性规划. 绝大部分的非线性规划我们都只能通过一些优化算法求得局部最优解. 深度神经网络的 l o s s loss loss优化就是一个非线性规划求解，不同的优化算法可以很好地求解对应的规划问题，但没有通用的优化方法.
- 关于多目标规划. 很多情况下我们都会用到多目标规划，如考虑成本降低同时收益增加，又如考虑利润增加同时加班时间减少. 具体多个目标同时优化求解，我们一般找不到最优解，取而代之的是有效解，或称Pareto最优解. 此类问题求解方法大致三种：一是给每个目标加权，整合到一个式子里优化，即优化它们的线性组合；二是分别求出各个目标的极小值，作为标签，像神经网络的优化那样，构造 l o s s loss loss函数，最终优化所有的 l o s s loss loss的线性组合，这也叫做理想点法；三是分级，一级一级地求解这些目标函数，每求出一级的极小值，就将其作为等式约束放到下一级的约束条件中，这也叫做序贯解法.
- 多目标规划的问题的具体求解手段一定是基于对实际问题背景的理解上的，具有解释性的，而非随机选一个方法.
- 对于数学建模而言，我们更关注其建模的过程，而求解的过程需要的是理论数学研究的，不必过于纠结. 不然随便一个优化算法都够学好久的了.
实际应用
- 谈到规划模型的实际应用，最关键的就是将实际问题抽象为数学形式.
- 这个时候就不像前面提及的线性或非线性的分类了，而是讨论如何将实际问题转化为数学形式.
- 其实在建模的过程中，我个人会更关注目标函数和约束条件的形式，例如：
- 目标函数该如何量化才更为合理，如果设置惩罚项又该如何设置.
- 对多目标规划的目标函数，其权重和优先级如何设置.
- 我们写目标函数可以写的非常丰富. 例如在供应链问题中，我们考虑众多因素给到目标函数优化，如最基础的收益和成本，我们可以进一步划分出各阶段所带来的收益和成本，不同阶段的收益和成本可能很难统一量纲，甚至难以量化；又例如在投资组合问题中，我们期待收益尽可能多，但同时又要评估风险，风险的评价指标是非常丰富的，最简单的办法是用股价波动的方差衡量，如此一来便有多目标规划模型，但如果给定的数据量非常多，如给出若干年内所有工作日的股价变化，我们对买和卖如何决策又需要考虑，因为灵活买卖可以更精细调节使得收益极大化，但同时又要考虑买卖需额外支付的佣金.
- 对目标函数的每一个变量，我们都应该考虑其具体意义，且对多个变量，我们将其放到一起又必须考虑量纲问题，若是加权，又要考虑权重设置，若是分级优化，我们又要考虑优先级问题，且是否存在一些变量是可以在同一级下优化.
- 对约束条件，例如简单的不等式和等式控制、比例约束，分段约束，我们都可以采用一些数学手段将其表达.
- 最关键之处还是从问题中提取这些约束条件.
- 不过在实际问题中，很多约束条件是需要挖掘的，并非像书本上的练习题那么理想，这就需要我们对实际业务场景的分析.
- 而且我们建模过程是需要发散思维的，我们可以融合多种思想构造目标函数和约束条件. 如利用差分构造变化率的约束条件，引入概率分析我们模型的参数敏感度.
- 对于求解规划模型的一些优化算法，其中的步长搜索利用到插值的办法找近似的最佳优化步长. 当然，这些都已经封装好，不需要我们考虑，我想借此表达的意思是对具体的应用场景应该充分挖掘并抽象出问题，不限于单一的知识点去解决问题，这才是真正将知识应用出来.
代码细节
- 避坑注意：
scipy.optimize的minimize的不等式约束默认是大于等于 ≥ \geq ≥，版本为1.8.1. - 下面给出各类规划模型的代码示例，省略具体问题，仅供快速查阅.
- 线性规划示例：
from scipy.optimize import linprog c = [-2, -3] # 转化为求极小值的最优化问题 A = [[1, 2]] b = [8] x0_bounds = (0, 4) x1_bounds = (0, 3) res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds]) print(res) - 二次规划示例：
from cvxopt import matrix, solvers """ Notes: solver.qp() 求解二次规划形如 minimize (1/2)*x'*P*x + q'*x subject to G*x <= h A*x = b. 每一个 matrix 必须为浮点型 """ P = matrix([ [8., 5, -20], [5, 72, -30], [-20, -30, 200] ]) q = matrix([0., 0, 0]) G = matrix([ [-5., -8, -10], [20, 25, 30], [-1, 0, 0], [0, -1, 0], [0, 0, -1] ]).T # 转置 h = matrix([-1000., 5000, 0, 0, 0]) sol = solvers.qp(P, q, G=G, h=h) print(f"最优解为 -------\n{sol['x']}最优值为 -------\n{sol['primal objective']}") - 非线性规划示例：
from scipy.optimize import minimize import numpy as np """ Notes: minimize 的约束部分 constraints 字典中 键 'ineq' 是指大于等于 >=, 不同于一般的 <= """ np.set_printoptions(precision=2, suppress=True) # 调整 numpy.ndarray 的输出格式 np.random.seed(6) # 固定随机种子 def obj(x): """ x shape is (16,) 一维向量 前 12 个变量为运送量 c_ij 后 4 个变量为料场坐标 x_1, y_1, x_2, y_2 """ c1 = x[0:6] c2 = x[6:12] x1, y1, x2, y2 = x[12:16] 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]) f = 0 for i in range(6): d1 = ((x1 - a[i]) ** 2 + (y1 - b[i]) ** 2) ** 0.5 d2 = ((x2 - a[i]) ** 2 + (y2 - b[i]) ** 2) ** 0.5 f = f + c1[i] * d1 + c2[i] * d2 return f def eq_fun(x): C = np.array([x[0:6], x[6:12]]) d = np.array([3, 5, 4, 7, 6, 11]) return C.T @ np.ones(C.shape[0]) - d def ineq_fun(x): C = np.array([x[0:6], x[6:12]]) e = np.array([20, 20]) return - (C @ np.ones(C.shape[1]) - e) # 不等条件转换为 >= cons = ({'type': 'eq', 'fun': eq_fun}, {'type': 'ineq', 'fun': ineq_fun}) bounds = ((0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0.5, 8.75), (0.75, 7.75), (0.5, 8.75), (0.75, 7.75)) x0 = np.append(np.random.rand(12) * 4, np.random.rand(4) * 7) # 随机初始化 res = minimize(obj, x0, constraints=cons, bounds=bounds, method='SLSQP') print(f"优化是否成功 -------\n{res.success}\n最优解为 -------\n{res.x}\n最优值为 -------\n{res.fun}") - 整数规划示例
import cvxpy as cp import numpy as np C = np.array([ [6, 7, 11, 2], [4, 5, 9, 8], [3, 1, 10, 4], [5, 9, 8, 2] ]) A0 = np.ones(C.T.shape) b0 = np.ones((C.T.shape[0], 1)) # 任务数量限制 A1 = np.ones(C.shape) b1 = np.ones((C.shape[0], 1)) # 人员数量限制 X = cp.Variable(C.shape, integer=True) expr = np.ones(C.shape[0]) @ cp.multiply(C, X) @ np.ones(C.shape[1]).T obj = cp.Minimize(expr) cons = [ A0 @ X == b0, A1 @ X.T == b1, X >= 0 ] prob = cp.Problem(obj, cons) prob.solve(solvers='GLPK_MI', verbose=False) print(f"最优解为 -------\n{prob.value}\n最优值为 -------\n{X.value}") - 多目标优化之序贯解法（分级优化）示例：
import cvxpy as cp import numpy as np """ Notes: X 优化变量，形状为 (10,) x1 甲产品产量 x2 乙产品产量 d1- 利润未达目标的量 d1+ 利润超出目标的量 d2- 甲乙比例低于目标比例的量 d2+ 甲乙比例高于目标比例的量 d3- 设备 A 的闲置时间 d3+ 设备 A 的加班时间 d4- 设备 B 的闲置时间 d4+ 设备 B 的加班时间 """ A_eq = np.array([ [2, 3, 1, -1, 0, 0, 0, 0, 0, 0], [1, -1, 0, 0, 1, -1, 0, 0, 0, 0], [2, 2, 0, 0, 0, 0, 1 ,-1, 0, 0], [1, 2, 0, 0, 0, 0, 0, 0, 1, -1] ]) b_eq = np.array([12, 0, 12, 8]) A_ub = np.array([ [4, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 4, 0, 0, 0, 0, 0, 0, 0, 0] ]) b_ub = np.array([16, 12]) X = cp.Variable((10,), integer=True) A1 = np.array([0, 0, 1, 0, 0, 0, 0, 0, 0, 0]) expr = A1 @ X obj = cp.Minimize(expr) cons = [ A_eq @ X == b_eq, A_ub @ X <= b_ub, X >= 0 ] prob = cp.Problem(obj, cons) prob.solve(solvers='GLPK_MI', verbose=False) print(f"最优解为 -------\n{prob.value}\n最优值为 -------\n{X.value}")# 第一级 expr 最优值为 0 # 将其作为约束代入下一级优化 cons.append(A1 @ X == 0) A2 = np.array([0, 0, 0, 0, 1, 1, 0, 0, 0, 0]) expr = A2 @ X obj = cp.Minimize(expr) prob = cp.Problem(obj, cons) prob.solve(solvers='GLPK_MI', verbose=False) print(f"最优解为 -------\n{prob.value}\n最优值为 -------\n{X.value}")# 第二级 expr 最优值为 0 # 将其作为约束代入下一级优化 cons.append(A2 @ X == 0) A3 = np.array([0, 0, 0, 0, 0, 0, 3, 3, 0, 1]) expr = A3 @ X obj = cp.Minimize(expr) prob = cp.Problem(obj, cons) prob.solve(solvers='GLPK_MI', verbose=False) print(f"最优解为 -------\n{prob.value}\n最优值为 -------\n{X.value}")
8199
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
