数学建模(1)线性规划模型

此处分享数模学习笔记。Anaconda 下载地址。安装时，记得勾选环境变量。

线性规划模型^[1]

矩阵的相关运算

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6]])
b = np.array([[1, 2], [3, 4], [5, 6]])
c = np.array([[1, 2, 3]])
d = np.array([[9, 8, 7], [3, 2, 1]])
e = np.array([[1, 2], [3, 4]])

s = a+d  # 加法
f = 3*a  # 数乘
g = np.dot(a, b)  # 矩阵乘
h = a*d  # 元素乘
i = c.T  # 转置
j = np.linalg.inv(e)  # 逆矩阵
k = np.linalg.det(e)  # 行列式
m = np.linalg.matrix_rank(d)  # 秩

---

求一次方程组的解

主要通过 scipy sympy numpy 这三个库就能实现各种各样的一次方程组求解。sympy 主要用于符号解，numpy 和 scipy 主要用于数值解。

from sympy import symbols, Eq, solve
import numpy as np

A = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]])
b = np.array([72, 83, 42])

inv_A = np.linalg.inv(A)
x = inv_A.dot(b)

y = np.linalg.solve(A, b)

print(x)  # 7-8行与10行等价，x与y输出相同

x, y, z = symbols('x y z')
eqs = [Eq(10*x-y-2*z, 72),
       Eq(-x+10*y-2*z, 83),
       Eq(-x-y+5*z, 42)]

print(solve(eqs, [x, y, z]))

• 输出：

[11. 12. 13.]
{x: 11, y: 12, z: 13}

---

求线性规划的解

定义：在线性等式和不等式约束下，最小化/最大化线性目标函数。

\[ \begin{gather} \min_{\mathbf x}{\mathbf c^\intercal\mathbf x} \\ \text{s.t.} \begin{cases} \mathbf A\mathbf x \leq \mathbf b \\ \mathbf A_{eq}\mathbf x = \mathbf b_{eq}\\ \mathbf{l}_{\mathbf b} \leq\mathbf x\leq \mathbf u_{\mathbf b} \end{cases} \end{gather} \]

其中 \(\mathbf c\) 为价值向量，\(\mathbf b\) 为资源向量，\(\mathbf x\) 为决策变量的向量，\(\mathbf A\)，\(\mathbf A_{eq}\) 分别为不等式约束和等式约束对应的矩阵。

• 例 1：求解

\[ \begin{gather} \max_{\mathbf x}z=4x_1+3x_2 \\ \text{s.t.} \begin{cases} 2x_1+x_2 \leq 10 \\ x_1+x_2\leq 8 \\ x_2\leq 7 \\ x_1,x_2\geq0 \end{cases} \end{gather} \]

%基于求解器求解
clear;
c = [4; 3];
A = [2, 1; 1, 1; 0, 1];
b = [10; 8; 7];
lb = zeros(2, 1);
[x, fval] = linprog(-c, A, b, [], [], lb); %有等号约束，则[]依次为Aeq和beq
x, fval = -fval %最大化则需要把fval变负

%基于问题求解
clear;
prob = optimproblem('ObjectiveSense', 'max'); %最小化则去除整个小括号，只留optimproblem
c = [4; 3];
A = [2, 1; 1, 1; 0, 1];
b = [10; 8; 7];
lb = zeros(2, 1);
x = optimvar('x', 2, 'LowerBound', 0); %2为决策维度
prob.Objective = c' * x;
prob.Constraints.con = A * x <= b;
[sol, fval, flage, out] = solve(prob) %fval返回最优值
sol.x %显示决策变量的值

其他代码：

clear;
prob = optimproblem('ObjectiveSense', 'max');
x = optimvar('x', 3, 'LowerBound', 0);
prob.Objective = 2 * x(1) + 3 * x(2) - 5 * x(3);
prob.Constraints.con1 = x(1) + x(2) + x(3) == 7;
prob.Constraints.con2 = 2 * x(1) - 5 * x(2) + x(3) >= 10;
prob.Constraints.con3 = x(1) + 3 * x(2) + x(3) <= 12;
[sol, fval, flag, out] = solve(prob);
sol.x
fval

clear;
prob = optimproblem;
x = optimvar('x', 4, 4, 'LowerBound', 0);
prob.Objective = 2800 * sum(x(:, 1)) + 4500 * sum(x(1:3, 2)) + ...
    6000 * sum(x(1:2, 3)) + 7300 * x(1, 4);
prob.Constraints.con = [sum(x(1, :)) >= 15,
                        sum(x(1, 2:4)) + sum(x(2, 1:3)) >= 10,
                        x(1, 3) + x(1, 4) + x(2, 2) + x(2, 3) + x(3, 1) + x(3, 2) >= 20,
                        x(1, 4) + x(2, 3) + x(3, 2) + x(4, 1) >= 12];
[sol, fval, flag, out] = solve(prob);
sol.x
fval

clear;
a = load('data1LP5.txt');
c = a(1:end - 1, 1:end - 1);
d = a(end, 1:end - 1);
e = a(1:end - 1, end);
prob = optimproblem;
x = optimvar('x', 6, 8, 'LowerBound', 0);
prob.Objective = sum(sum(c .* x));
prob.Constraints.con1 = sum(x, 1) == d;
prob.Constraints.con2 = sum(x, 2) <= e;
[sol, fval, flag, out] = solve(prob);
xx = sol.x
fval
writematrix(xx, 'out1LP5.xlsx')

---

Python 方面

我们运用到：

scipy.optimize.linprog(c,A=None,b=None,A_eq=None,b_eq=None,bounds=None,method='',callback=None,options=None,x0=None)​

其中

• c：【一维数组】线性目标函数的系数
• A(可选)：【二维数组】不等式约束矩阵，\(\mathbf A_{ub}\) 的每一行指定\(x\)上的线性不等式约束的系数
• b(可选)：【一维数组】不等式约束向量，每个元素代表\(\mathbf A_{ub}x\) 的上限
• A*eq(可选)：【二维数组】等式约束矩阵，\(\mathbf A*{eq}\) 的每一行指定\(x\)上的线性等式约束的系数
• b*eq(可选)：【一维数组】等式约束向量，\(\mathbf A*{eq}x\) 的每个元素必须等于\(\mathbf b_{eq}\) 的对应元素
• bounds(可选)：【\((min, max)\)序列对】定义决策变量\(x\)的最小值和最大值
• None：使用 None 表示没有界限，默认情况下，界限为\((0，None)\)(所有决策变量均为非负数)。
• 如果提供一个元组\((min, max)\)，则最小值和最大值将用作所有决策变量的界限。
• method(可选)：【字符串】{'interior-point', 'revised simplex', 'simplex'}，三种算法可选
• callback(可选)：调用回调函数，等待被调用的参数，如果提供了回调函数，则算法的每次迭代将至少调用一次。回调函数必须接受单个 scipy.optimize.OptimizeResult 由以下字段组成：
• https://blog.csdn.net/weixin_45288557/article/details/109139303

解

\[ \begin{gather} \max \ z=2x_1+3x_2-5x_3 \\ \text{s.t.} \begin{cases} x_1+x_2+x_3 & = 7 \\ 2x_1-5x_2+x_3 & \geq 10 \\ x_1+3x_2+x_3 & \leq 12 \\ x_1,x_2,x_3 & \geq 0 \end{cases} \end{gather} \]

• 注意：求解 max 则取负的-c。约束处小于等于取正，大于等于取负。

from scipy import optimize
import numpy as np

c = np.array([2, 3, -5])
A = np.array([[-2, 5, -1], [1, 3, 1]])
b = np.array([-10, 12])
A_eq = np.array([[1, 1, 1]])
b_eq = np.array([7])
x1, x2, x3 = (0, None), (0, None), (0, None)

res = optimize.linprog(-c, A, b, A_eq, b_eq, bounds=(x1, x2, x3))
print(res)

• 输出：

    con: array([1.80713222e-09])
    fun: -14.57142856564506
message: 'Optimization terminated successfully.'
    nit: 5
  slack: array([-2.24583019e-10,  3.85714286e+00])
 status: 0
success: True
      x: array([6.42857143e+00, 5.71428571e-01, 2.35900788e-10])

• con：【一维数组】等式约束的残差(名义上为 0)\(\mathbf b_{eq}-\mathbf A_{eq}\mathbf x\)
• fun：【浮点数】目标函数的最佳值(\(\mathbf c^T\mathbf x\))
• message：【字符串】算法退出状态的字符串描述符
• nit：【整数】在所有阶段中执行的迭代总数
• slack：【一维数组】不等式约束的松弛值(名义上为正值)\(\mathbf b-\mathbf{Ax}\)
• status：【整数】算法退出状态，0优化成功终止，1达到迭代限制，2问题似乎不可行，3问题似乎不收敛，4遇到数值困难
• success：【布尔值】当算法成功找到最佳解决方案时为True
• x：【一维数组】在满足约束的情况下将目标函数最小化的决策变量的值。2.35900788e-10 由于精度问题，可视为0。

---

可以转化为线性规划的问题

解

\[ \begin{gather} \min \ z=\left|x_1\right|+\left|x_2\right|+\dots+\left|x_n\right| \\ \text{s.t.}\ \mathbf A\mathbf x\leq\mathbf b \end{gather} \]

令 \(u_i=\frac{\left|x_i\right|+x_i}{2}\geq0\)，\(v_i=\frac{\left|x_i\right|-x_i}{2}\geq0\)，则将问题转化为

\[ \begin{gather} \min \ \sum_{i=1}^n\left(u_i+v_i\right) \\ \text{s.t.} \begin{cases} \mathbf A\left(\mathbf u-\mathbf v\right)\leq \mathbf b \\ \mathbf u, \mathbf v\geq0 \end{cases} \end{gather} \]

例：

\[ \begin{gather} \min \ z=\left|x_1\right|+2\left|x_2\right|+3\left|x_3\right|+4\left|x_4\right| \\ \text{s.t.} \begin{cases} x_1-x_2-x_3+x_4\leq -2 \\ x_1-x_2+x_3-3x_4\leq -1 \\ x_1-x_2-2x_3+3x_4\leq -\frac{1}{2} \end{cases} \end{gather} \]

化为

\[ \begin{gather} \min \ \begin{pmatrix} 1 & 2 & 3 & 4 \end{pmatrix}\left(\mathbf u+\mathbf v\right) \\ \text{s.t.} \begin{cases} \begin{pmatrix} 1 & -1 & -1 & 1\\ 1 & -1 & 1 & -3\\ 1 & -1 & -2 & 3 \end{pmatrix}\left(\mathbf u-\mathbf v\right)\leq \begin{pmatrix} -2\\ -1\\ -\frac{1}{2} \end{pmatrix} \\ \mathbf u, \mathbf v\geq0 \end{cases} \end{gather} \]

clear, close all;
c = [1:4]';
A = [1, -1, -1, 1; 1, -1, 1, -3; 1, -1, -2, 3];
b = [-2; -1; -1/2];
prob = optimproblem;
u = optimvar('u', 4, 'LowerBound', 0);
v = optimvar('v', 4, 'LowerBound', 0);
prob.Objective = sum(c' * (u + v));
prob.Constraints.con = A * (u - v) <= b;
[sol, fval, flag, out] = solve(prob);
x = sol.u - sol.v
fval

---

补充

【数学建模 之 投资的收益与风险(线性规划)】1998 年全国大学生数学建模竞赛 A 题

参考文献

1. 司守奎，孙玺菁. 数学建模算法与应用(第 3 版). ↩︎