数学建模(3)非线性规划模型

非线性规划模型^[1]

\[ \begin{gather} \min z=x_1^2+x_2^2-4x_1+4 \\ \text{s.t.} \begin{cases} -x_1+x_2-2\leq0 \\ x_1^2-x_2+1\leq0 \\ x_1,x_2\geq0 \end{cases} \end{gather} \]

clear;
prob = optimproblem;
x = optimvar('x', 2, 'LowerBound', 0);
prob.Objective = sum(x .^ 2) - 4 * x(1) + 4;
prob.Constraints.con = [-x(1) + x(2) - 2 <= 0,
                        x(1) ^ 2 - x(2) + 1 <= 0];
x0.x = rand(2, 1); %非线性规划必须赋初值
[sol, fval, flag, out] = solve(prob, x0);
sol.x

---

彩电问题

clear, close all;
format long g
syms x1 x2 %定义符号变量
f = (339 - 0.01 * x1 - 0.003 * x2) * x1 + (399 - 0.004 * x1 - 0.01 * x2) * x2 - (400000 + 195 * x1 + 225 * x2);
f = simplify(f); %化简目标函数
f1 = diff(f, x1);
f2 = diff(f, x2); %求目标函数关于x1, x2的偏导数
[x10, x20] = solve(f1, f2); %解代数方程求驻点
x10 = round(double(x10)) %取整
x20 = round(double(x20)) %取整
f0 = subs(f, {x1, x2}, {x10, x20}); %求目标函数的取值
f0 = double(f0)
subplot(121), fmesh(f, [0, 10000, 0, 10000]), title('') %画三维图形
xlabel('$x_1$', 'Interpreter', 'latex')
ylabel('$x_2$', 'Interpreter', 'latex')
subplot(122)
c = fcontour(f, [0, 10000, 0, 10000]);
contour(c.XData, c.YData, c.ZData, 'ShowText', 'on') %等高线标注
xlabel('$x_1$', 'Interpreter', 'latex')
ylabel('$x_2$', 'Interpreter', 'latex')
p1 = 339 - 0.01 * x10 - 0.003 * x20 %计算19英寸的平均售价
p2 = 399 - 0.004 * x10 - 0.01 * x20 %计算21英寸的平均售价
c = 400000 + 195 * x10 + 225 * x20 %计算总支出
rate = f0 / c %计算利润率

• 灵敏度分析

clear;
format long g
syms x1 x2 a %定义符号变量
f = (339 - a * x1 - 0.003 * x2) * x1 + (399 - 0.004 * x1 - 0.01 * x2) * x2 - (400000 + 195 * x1 + 225 * x2);
f = simplify(f) %化简目标函数
f1 = diff(f, x1);
f2 = diff(f, x2); %求目标函数关于x1,x2的偏导数
[x10, x20] = solve(f1, f2); %求驻点
pretty(x10)
pretty(x20) %以书写习惯的方式显示
subplot(121), fplot(x10, [0.002, 0.02]), title('') %画x1关于a的曲线
xlabel('$a$', 'Interpreter', 'latex')
ylabel('$x_1$', 'Interpreter', 'latex', 'Rotation', 0)
subplot(122), fplot(x20, [0.002, 0.02]), title('') %画x2关于a的曲线
xlabel('$a$', 'Interpreter', 'latex')
ylabel('$x_2$', 'Interpreter', 'latex', 'Rotation', 0)
dx1 = diff(x10, a);
dx10 = subs(dx1, a, 0.01);
dx10 = double(dx10)
sx1a = dx10 * 0.01/4735
dx2 = diff(x20, a);
dx20 = subs(dx2, a, 0.01);
dx20 = double(dx20)
sx2a = dx20 * 0.01/7043
F = subs(f, {x1, x2}, {x10, x20}); %求关于a的目标函数
F = simplify(F) %对目标函数进行化简
figure, fplot(F, [0.002, 0.02]), title('')
xlabel('$a$', 'Interpreter', 'latex')
ylabel('$y$', 'Interpreter', 'latex', 'Rotation', 0)
Sya = -4735 ^ 2 * 0.01/553641
f3 = subs(f, {x1, x2, a}, {4735, 7043, 0.011});
f3 = double(f3) %计算近似最优利润
f4 = subs(F, a, 0.011);
f4 = double(f4) %计算最优利润
delta = (f4 - f3) / f4 %计算利润的相对误差

---

二次规划

\[ \begin{gather} \min -x_1^2-0.3x_1x_2-2x_2^2+98x_1+277x_2 \\ \text{s.t.} \begin{cases} x_1+x_2\leq100 \\ x_1-x_2\leq0 \\ x_1,x_2\geq0 \end{cases} \end{gather} \]

矩阵化转化为

\[ \begin{gather} \min \begin{pmatrix} x_1 & x_2 \end{pmatrix} \begin{pmatrix} -1 & -0.15 \\ -0.15 & -2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} + \begin{pmatrix} 98 & 277 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \\ \text{s.t.} \begin{cases} \begin{pmatrix} 1 & 1 \\ 1 & -2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \leq \begin{pmatrix} 100 \\ 0 \end{pmatrix} \\ x_1,x_2\geq0 \end{cases} \end{gather} \]

clear;
H = [-1, -0.15; -0.15, -2];
f = [98, 277];
A = [1, 1; 1, -2];
b = [100; 0];
prob = optimproblem;
x = optimvar('x', 2, 'LowerBound', 0);
prob.Objective = x' * H * x + f * x;
prob.Constraints = A * x <= b;
[sol, fval, flag, out] = solve(prob);
sol.x
%下面使用fmincon求解，效果较好
fx = @(x)x' * H * x + f * x; %目标函数的匿名函数，x为列向量
[x, y] = fmincon(fx, rand(2, 1), A, b, [], [], zeros(2, 1))

---

极值与极值点

求 \(f(x,y)=x^3-y^3+3x^2+3y^2-9x\) 在 \((0,0)\) 附近的极值。

clear;
f = @(x) x(1) ^ 3 - x(2) ^ 3 + 3 * x(1) ^ 2 + 3 * x(2) ^ 2 - 9 * x(1); %定义匿名函数
g = @(x) -f(x);
[xy1, z1] = fminunc(f, rand(2, 1)) %求极小值点
[xy2, z2] = fminsearch(g, rand(2, 1)); %求极大值点
xy2, z2 = -z2 %显示极大点及对应的极大值

求 \(f(x)=100(x_2-x_1^2)^2+(1-x_1)^2\) 的极小点。

clear;
prob = optimproblem;
x = optimvar('x', 2);
prob.Objective = 100 * (x(2) - x(1) ^ 2) ^ 2 + (1 - x(1)) ^ 2;
x0.x = rand(2, 1); %初始值
[sol, fval, flag, out] = solve(prob, x0);
sol.x

---

fmincon

\[ \begin{gather} \min z=x_1^2+x_2^2+4x_3^2+8 \\ \text{s.t.} \begin{cases} x_1^2-x_2+x_3^2\geq0 \\ x_1+x_2^2+x_3^2\leq20 \\ -x_1-x_2^2+2=0 \\ x_2+2x_3^2=3 \\ x_1,x_2,x_3\geq0 \end{cases} \end{gather} \]

clear;
fun1 = @(x) sum(x .^ 2) + 8;
[x, y] = fmincon(fun1, rand(3, 1), [], [], [], [], zeros(3, 1), [], @fun2)

function [c, ceq] = fun2(x)
    c = [-x(1) ^ 2 + x(2) - x(3) ^ 2,
         x(1) + x(2) ^ 2 + x(3) ^ 3 - 20]; %非线性不等式约束
    ceq = [-x(1) - x(2) ^ 2 + 2,
           x(2) + 2 * x(3) ^ 2 - 3]; %非线性等式约束
end

clear;
prob = optimproblem;
x = optimvar('x', 3, 'LowerBound', 0);
prob.Objective = sum(x .^ 2) + 8;
prob.Constraints.con1 = [-x(1) ^ 2 + x(2) - x(3) ^ 2 <= 0,
                         x(1) + x(2) ^ 2 + x(3) ^ 3 <= 20]; %非线性不等式约束
prob.Constraints.con2 = [-x(1) - x(2) ^ 2 + 2 == 0,
                         x(2) + 2 * x(3) ^ 2 == 3]; %非线性等式约束
x0.x = rand(3, 1); %非线性规划必须赋初值
[sol, fval, flag, out] = solve(prob, x0);
sol.x, fval

---

Lagrange 乘数法与 KKT 条件

参见：Lagrange 乘数法与 KKT 条件。

---

求约束极值函数

https://github.com/guofei9987/scikit-opt

• 安装 sko.GA：

pip install scikit-opt

【例 1】三台火电机组的运行成本(单位：t/h)与出力限制(单位：MW)分别如下：

\[ \begin{gather} F_{G1}=4+0.3P_{G1}+0.0007P_{G1}^2, & 100\leq P_{G1}\leq 200 \\ F_{G2}=3+0.32P_{G2}+0.0004P_{G2}^2, & 120\leq P_{G2}\leq 250 \\ F_{G3}=3.5+0.3P_{G3}+0.00045P_{G3}^2, & 150\leq P_{G3}\leq 300 \end{gather} \]

当负荷为\(700\)MW 时，求经济调度的结果。

1. 方法一：scipy

from scipy.optimize import minimize
import numpy as np
# 目标函数即min(FG1+FG2+FG3)

def func(x):
    return 4+0.3*x[0]+0.0007*x[0]**2+3+0.32*x[1]+0.0004*x[1]**2+3.5+0.3*x[2]+0.00045*x[2]**2

def con():
    # 约束条件 分为eq和ineq
    # eq表示函数结果等于0；ineq表示表达式大于等于。
    cons = ({'type': 'eq', 'func': lambda x: x[0]+x[1]+x[2]-700})
    # {'type': 'ineq', 'func': 1ambda x:-x[2] + x2max}#如果有不等式约束
    # cons=([con1, con2, con3, con4, con5, con6, con7, con8]) #如果有多个约束, 则最后返回结果是这个
    # x[0]其中的0必须是具体数字, 不能是t等参数
    return cons

# 上下限约束
b1, b2, b3 = (100, 200), (120, 250), (150, 300)
bnds = (b1, b2, b3)  # 边界约束
if __name__ == "__main__":
    cons = con()  # 约束
    # 设置x初始猜测值
    x0 = np.array((100, 200, 400))
    res = minimize(func, x0, method='BFGS', constraints=cons, bounds=bnds)
    print(res)
    # 代价res.fun
    # 解res.x

• 输出：

     fun: -135.64285713335553
hess_inv: array([[ 695.56014419,  -27.36205954,    1.88143735],
      [ -27.36205954, 1188.1171721 ,    2.12363692],
      [   1.88143735,    2.12363692, 1113.34590189]])
     jac: array([ 0.00000000e+00, -3.81469727e-06, -1.90734863e-06])
 message: 'Optimization terminated successfully.'
    nfev: 100
     nit: 14
    njev: 25
  status: 0
 success: True
       x: array([-214.28690489, -400.00384852, -333.33573003])

2. 方法二：遗传算法

from sko.GA import GA

def object_func(x):
    return 4+0.3 * x[0]+0.0007 * x[0]**2+3+0.32*x[1]+0.0004 * x[1]**2+3.5+0.3*x[2]+0.00045*x[2]**2

# 等式约束
def cons(x): return x[0]+x[1]+x[2]-700

ga = GA(func=object_func, n_dim=3, size_pop=500, max_iter=500,
        lb=[100, 120, 150], ub=[200, 250, 300], constraint_eq=[cons])
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)

• 输出：

best_x: [169.92779809 232.75054415 297.32165777]
 best_y: [306.81703248]

---

【例 2】某公司有 6 个建筑工地要开工，每个工地的位置(用平面坐标系 \(a,b\) 表示，距离单位：千米)及水泥日用量 \(d\)(吨)由下表给出。规划设立两个料场位于 A、B，日储量各为 \(20\) 吨。假设从料场到工地之间均有直线道路相连，试确定料场的位置，并制定每天的供应计划，即从 A、B 两料场分别向各工地运送多少吨水泥，使总的吨千米数最小。

表 1 工地位置 \((a,b)\) 及水泥日用量 \(d\)

	1	2	3	4	5	6
a	1.25	8.75	0.5	5.75	3	7.25
b	1.25	0.75	4.75	5	6.5	7.25
d	3	5	4	7	6	11

\[ \begin{gather} \min f=\sum_{i=1}^2\sum_{j=1}^6w_{ij}\sqrt{(x_i-a_j)^2+(y_i-b_j)^2} \\ \text{s.t.} \begin{cases} \sum_{i=1}^2x_{ij}\geq d_j \\ \sum_{j=1}^6x_{ij}\leq e_i \end{cases} \end{gather} \]

参考文献

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