【数学建模】《实战数学建模：例题与讲解》第三讲-非线性规划（含Matlab代码）
如果这篇文章对你有帮助，欢迎点赞与收藏~
非线性规划介绍
非线性规划（Nonlinear Programming，简称NLP）是一类涉及非线性目标函数和/或非线性约束的数学优化问题的解决方法。在数学建模的过程中，我们常常面临实际问题中的非线性规划，需要通过优化算法来寻找目标函数的最小值或最大值，同时满足一定的约束条件。本文将深入探讨非线性规划的基本概念、解决方法以及常见算法。
如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。
一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不像线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。
在非线性规划问题中，目标函数或约束中至少有一个是非线性的，这使得问题的求解相对复杂。常见的非线性规划算法包括梯度下降法、牛顿法、拟牛顿法等。
基本概念
非线性规划问题主要体现在目标函数或约束条件中至少包含一个非线性函数。与线性规划相比，非线性规划的求解难度更大，因为非线性函数引入了问题的复杂性。线性规划通常可以通过单纯形法等通用方法解决，而非线性规划则缺乏一种适用于所有问题的通用算法，各个方法往往具有特定的适用范围。
解决方法
在面对非线性规划问题时，寻找最优解的过程变得相对复杂。为了克服这一挑战，研究者们提出了多种解决方法，其中一些常见的包括梯度下降法、牛顿法、拟牛顿法等。
-
梯度下降法
梯度下降法是一种基于目标函数梯度信息的迭代优化算法。该方法通过不断沿着梯度的反方向更新参数，逐步接近最优解。梯度下降法在解决非线性规划问题时表现出色，但对初始值敏感，可能陷入局部最小值。 -
牛顿法
牛顿法是一种二阶优化算法，通过使用目标函数的二阶导数信息来加速收敛。相比于梯度下降法，牛顿法具有更快的收敛速度，但计算复杂度较高，尤其是在大规模问题上。 -
拟牛顿法
拟牛顿法是对牛顿法的改进，通过估计目标函数的海森矩阵来避免计算二阶导数。这样一来，拟牛顿法在一定程度上综合了梯度下降法和牛顿法的优点，具有较好的收敛性能和计算效率。
应用领域
非线性规划广泛应用于科学、工程和经济等领域。例如，在工程设计中，非线性规划可用于优化材料的使用；在经济学中，它可以帮助制定最优的资源分配方案。随着计算机算力的提升，非线性规划在实际问题中的应用变得更加普遍。
注意点
对于一个实际问题，在把它归结成非线性规划问题时， 一般要注意如下几点：
(1) 确定供选方案： 首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认问题的可供选择的方案，并用一组变量来表示它们。
(2) 提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。
(3) 给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好” 或“坏＂的价值标准，并用某种数量形式来描述它。
(4) 寻求限制条件： 由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。
注： 线性规划与非线性规划的区别在千，如果线性规划的最优解存在，则其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到），而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。
习题3.1
1. 题目要求
2.解题过程
解：
设 第一季度生产x1台，第二季度生产x2台，第三季度生产x3台。
第一季度：费用全部为生产费用
f
1
=
0.2
x
1
2
+
50
x
1
f_1=0.2x_1^2+50x_1
f1=0.2x12+50x1
第二季度：费用为生产费用和存储费用
f
2
=
0.2
x
2
2
+
50
x
2
+
4
(
x
1
−
40
)
f_2=0.2x_2^2+50x_2+4(x_1-40)
f2=0.2x22+50x2+4(x1−40)
第三季度：费用为生产费用和存储费用
f
3
=
0.2
x
2
3
+
50
x
3
+
4
(
x
1
+
x
2
−
40
−
60
)
f_3=0.2x_2^3+50x_3+4(x_1+x_2-40-60)
f3=0.2x23+50x3+4(x1+x2−40−60)
由题目所给数据可建立如下非线性规划模型：
min
f
=
0.2
x
1
2
+
50
x
1
+
0.2
x
2
2
+
50
x
2
+
4
(
x
1
−
40
)
+
0.2
x
2
3
+
50
x
3
+
4
(
x
1
+
x
2
−
40
−
60
)
s. t.
{
x
1
⩾
40
x
1
+
x
2
⩾
100
x
1
+
x
2
+
x
3
⩾
180
x
i
⩾
0
,
i
=
1
,
2
,
3
\min f=0.2x_1^2+50x_1+0.2x_2^2+50x_2+4(x_1-40)+0.2x_2^3+50x_3+4(x_1+x_2-40-60) \\ \text { s. t. }\left\{ \begin{array}{l} x_{1} \geqslant 40 \\ x_{1}+x_{2} \geqslant 100 \\ x_{1}+x_{2}+x_{3} \geqslant 180 \\ x_{i} \geqslant 0, i=1,2,3 \end{array}\right. \\
minf=0.2x12+50x1+0.2x22+50x2+4(x1−40)+0.2x23+50x3+4(x1+x2−40−60) s. t. ⎩
⎨
⎧x1⩾40x1+x2⩾100x1+x2+x3⩾180xi⩾0,i=1,2,3
化成MATLAB标准型，即：
min
f
=
1
2
[
x
1
x
2
x
3
]
[
0.4
0
0
0
0.4
0
0
0
0.4
]
[
x
1
x
2
x
3
]
+
[
58
54
50
]
[
x
1
x
2
x
3
]
−
560
s. t.
{
[
1
0
0
1
1
0
1
1
1
]
[
x
1
x
2
x
3
]
⩾
[
40
100
180
]
x
i
⩾
0
,
i
=
1
,
2
,
3
\min f= \frac{1}{2} \begin{bmatrix} x_1 & x_2 & x_3 \end{bmatrix} \begin{bmatrix} 0.4 & 0 & 0 \\ 0 & 0.4 & 0 \\ 0 & 0 & 0.4 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 58 & 54 & 50 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} -560 \\ \text { s. t. }\left\{ \begin{array}{l} \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} \geqslant \begin{bmatrix} 40 \\ 100 \\ 180 \end{bmatrix} \\ x_{i} \geqslant 0, i=1,2,3 \end{array}\right. \\
minf=21[x1x2x3]
0.40000.40000.4
x1x2x3
+[585450]
x1x2x3
−560 s. t. ⎩
⎨
⎧
111011001
x1x2x3
⩾
40100180
xi⩾0,i=1,2,3
3.程序
求解的MATLAB程序如下：
% 二次规划
clc , clear
H = 0.4 * eye(3); % 实对称矩阵
f = [58, 54, 50]; % 一次项
% 线性约束
A = [-1, 0, 0; ...
-1, -1, 0; ...
-1, -1, -1];
b = [-40, -100, -180];
%求解
[x, value] = quadprog(H, f, A, b, [], [], zeros(3, 1));
% 展示结果
x
value - 560
4.结果
求得的最优解是：
第一季度生产50台，第二季度生产60台，第三季度生产70台。
最大利润：11280（元）
习题3.2
1. 题目要求
2.解题过程
解：
模型二中，建立如下非线性规划模型：
min
∑
i
=
1
6
(
Δ
θ
i
)
2
,
s. t.
{
Δ
i
j
<
0
,
1
⩽
i
⩽
5
,
i
+
1
⩽
j
⩽
6
,
∣
Δ
θ
i
∣
⩽
π
6
,
i
=
1
,
2
,
⋯
,
6
。
\min \sum_{i=1}^{6}\left(\Delta \theta_{i}\right)^{2}, \\ \text { s. t. }\left\{\begin{array}{l} \Delta_{i j}<0,1 \leqslant i \leqslant 5, i+1 \leqslant j \leqslant 6, \\ \left|\Delta \theta_{i}\right| \leqslant \frac{\pi}{6}, i=1,2, \cdots, 6 \text { 。 } \end{array}\right.
mini=1∑6(Δθi)2, s. t. {Δij<0,1⩽i⩽5,i+1⩽j⩽6,∣Δθi∣⩽6π,i=1,2,⋯,6 。
其中：
Δ
i
j
=
b
~
i
j
2
−
4
a
~
i
j
c
~
i
j
<
0
。
a
~
i
j
=
4
sin
2
θ
i
−
θ
j
2
,
b
~
i
j
=
2
{
[
x
i
(
0
)
−
x
j
(
0
)
]
(
cos
θ
i
−
cos
θ
j
)
+
[
y
i
(
0
)
−
y
j
(
0
)
]
(
sin
θ
i
−
sin
θ
j
)
}
,
c
~
i
j
=
[
x
i
(
0
)
−
x
j
(
0
)
]
2
+
[
y
i
(
0
)
−
y
j
(
0
)
]
2
−
64
。
\begin{array}{l} \Delta_{i j}=\tilde{b}_{i j}^{2}-4 \tilde{a}_{i j} \tilde{c}_{i j}<0 \text { 。 } \\ \tilde{a}_{i j}=4 \sin ^{2} \frac{\theta_{i}-\theta_{j}}{2}, \\ \tilde{b}_{i j}=2 \left\{\left[x_{i}(0)-x_{j}(0)\right]\left(\cos \theta_{i}-\cos \theta_{j}\right)+\left[y_{i}(0)-y_{j}(0)\right]\left(\sin \theta_{i}-\sin \theta_{j}\right)\right\}, \\ \tilde{c}_{i j}=\left[x_{i}(0)-x_{j}(0)\right]^{2}+\left[y_{i}(0)-y_{j}(0)\right]^{2}-64 。 \end{array}
Δij=b~ij2−4a~ijc~ij<0 。 a~ij=4sin22θi−θj,b~ij=2{[xi(0)−xj(0)](cosθi−cosθj)+[yi(0)−yj(0)](sinθi−sinθj)},c~ij=[xi(0)−xj(0)]2+[yi(0)−yj(0)]2−64。
3.程序
编写非线性约束函数
function [A, Aeq] = fun_A(x) % 变量是Theta的变化量
% 记录飞机初始数据
x0 = [150, 85, 150, 145, 130, 0]';
y0 = [140, 85, 155, 50, 150, 0]';
angle0 = [243, 236, 220.5, 159, 230, 52]';
angle = angle0 + x; % 初始角度加上角度的变化量：变化后的方向角
% 从现在开始，创建不等式约束，有ij嵌套循环
num = 1; % 计数
for i = 1:5
for j = i + 1:6
aij = 4 * (sind((angle(i) - angle(j))/2))^2;
bij = (x0(i) - x0(j)) * (cosd(angle(i)) - cosd(angle(j))) + (y0(i) - y0(j)) * (sind(angle(i)) - sind(angle(j)));
bij = bij * 2;
cij = (x0(i) - x0(j))^2 + (y0(i) - y0(j))^2 - 64;
A(num) = bij^2 - 4 * aij * cij;
num = num + 1;
end
end
Aeq = []; % 没有等式约束
end
编写目标函数
function f = fun_f(x) % 变量是Theta的变化量
f = sum(x.^2);
end
主程序
clc, clear
[theta, w] = fmincon('fun_f', rand(6, 1), [], [], [], [], -30*ones(6, 1), 30*ones(6, 1), 'fun_A');
theta
w
4.结果
求得的最优解是：
Δ θ 1 = − 6.6683 , Δ θ 2 = 0.3315 , Δ θ 3 = 2.0589 , Δ θ 4 = − 0.5004 , Δ θ 5 = 6.3313 , Δ θ 6 = 1.5709 \Delta\theta_1= -6.6683, \Delta\theta_2= 0.3315, \Delta\theta_3= 2.0589, \Delta\theta_4= -0.5004, \Delta\theta_5= 6.3313, \Delta\theta_6= 1.5709 Δθ1=−6.6683,Δθ2=0.3315,Δθ3=2.0589,Δθ4=−0.5004,Δθ5=6.3313,Δθ6=1.5709
习题3.3
1. 题目要求
2.解题过程
解：
模型二中，建立如下非线性规划模型：
min
∑
i
=
1
6
(
Δ
θ
i
)
2
,
s. t.
{
Δ
i
j
<
0
,
1
⩽
i
⩽
5
,
i
+
1
⩽
j
⩽
6
,
∣
Δ
θ
i
∣
⩽
π
6
,
i
=
1
,
2
,
⋯
,
6
。
\min \sum_{i=1}^{6}\left(\Delta \theta_{i}\right)^{2}, \\ \text { s. t. }\left\{\begin{array}{l} \Delta_{i j}<0,1 \leqslant i \leqslant 5, i+1 \leqslant j \leqslant 6, \\ \left|\Delta \theta_{i}\right| \leqslant \frac{\pi}{6}, i=1,2, \cdots, 6 \text { 。 } \end{array}\right.
mini=1∑6(Δθi)2, s. t. {Δij<0,1⩽i⩽5,i+1⩽j⩽6,∣Δθi∣⩽6π,i=1,2,⋯,6 。
其中：
Δ
i
j
=
b
~
i
j
2
−
4
a
~
i
j
c
~
i
j
<
0
。
a
~
i
j
=
4
sin
2
θ
i
−
θ
j
2
,
b
~
i
j
=
2
{
[
x
i
(
0
)
−
x
j
(
0
)
]
(
cos
θ
i
−
cos
θ
j
)
+
[
y
i
(
0
)
−
y
j
(
0
)
]
(
sin
θ
i
−
sin
θ
j
)
}
,
c
~
i
j
=
[
x
i
(
0
)
−
x
j
(
0
)
]
2
+
[
y
i
(
0
)
−
y
j
(
0
)
]
2
−
64
。
\begin{array}{l} \Delta_{i j}=\tilde{b}_{i j}^{2}-4 \tilde{a}_{i j} \tilde{c}_{i j}<0 \text { 。 } \\ \tilde{a}_{i j}=4 \sin ^{2} \frac{\theta_{i}-\theta_{j}}{2}, \\ \tilde{b}_{i j}=2 \left\{\left[x_{i}(0)-x_{j}(0)\right]\left(\cos \theta_{i}-\cos \theta_{j}\right)+\left[y_{i}(0)-y_{j}(0)\right]\left(\sin \theta_{i}-\sin \theta_{j}\right)\right\}, \\ \tilde{c}_{i j}=\left[x_{i}(0)-x_{j}(0)\right]^{2}+\left[y_{i}(0)-y_{j}(0)\right]^{2}-64 。 \end{array}
Δij=b~ij2−4a~ijc~ij<0 。 a~ij=4sin22θi−θj,b~ij=2{[xi(0)−xj(0)](cosθi−cosθj)+[yi(0)−yj(0)](sinθi−sinθj)},c~ij=[xi(0)−xj(0)]2+[yi(0)−yj(0)]2−64。
罚函数为：
P
(
θ
,
M
)
=
∑
i
=
1
6
(
Δ
θ
i
)
2
+
M
∑
i
=
1
r
max
(
Δ
i
j
,
0
)
P(\boldsymbol{\theta}, M)=\sum_{i=1}^{6}\left(\Delta \theta_{i}\right)^{2}+M \sum_{i=1}^{r} \max \left(\Delta_{i j}, 0\right)
P(θ,M)=i=1∑6(Δθi)2+Mi=1∑rmax(Δij,0)
3.程序
编写罚函数
function P = penalty_function(x)
M = 1000000;
f = sum(x.^2);
% 记录飞机初始数据
x0 = [150, 85, 150, 145, 130, 0]';
y0 = [140, 85, 155, 50, 150, 0]';
angle0 = [243, 236, 220.5, 159, 230, 52]';
angle = angle0 + x; % 初始角度加上角度的变化量：变化后的方向角
% 从现在开始，创建g，有ij嵌套循环
num = 1; % 计数
for i = 1:5
for j = i + 1:6
aij = 4 * (sind((angle(i) - angle(j))/2))^2;
bij = (x0(i) - x0(j)) * (cosd(angle(i)) - cosd(angle(j))) + (y0(i) - y0(j)) * (sind(angle(i)) - sind(angle(j)));
bij = bij * 2;
cij = (x0(i) - x0(j))^2 + (y0(i) - y0(j))^2 - 64;
g(num) = bij^2 - 4 * aij * cij;
num = num + 1;
end
end
P = f + M * max([g, 0]);
end
主程序
clc, clear
[theta,w]=fminsearch('penalty_function',rand(6,1));
theta
w
4.结果
求得的最优解是：
Δ θ 1 = − 6.2039 , Δ θ 2 = 0.4905 , Δ θ 3 = 0.9688 , Δ θ 4 = − 2.3849 , Δ θ 5 = 6.7960 , Δ θ 6 = 3.4565 \Delta\theta_1= -6.2039, \Delta\theta_2= 0.4905, \Delta\theta_3= 0.9688, \Delta\theta_4= -2.3849, \Delta\theta_5= 6.7960, \Delta\theta_6= 3.4565 Δθ1=−6.2039,Δθ2=0.4905,Δθ3=0.9688,Δθ4=−2.3849,Δθ5=6.7960,Δθ6=3.4565
习题3.4
1.题目要求
2.解题过程
解：
化成MATLAB标准型，即：
min
w
=
(
−
1
)
∗
[
[
3
1
0
]
[
x
1
2
x
2
2
x
3
2
]
+
[
2
3
1
]
[
x
1
x
2
x
3
]
]
s. t.
{
[
2
2
0
]
[
x
1
2
x
2
2
x
3
2
]
+
[
1
1
1
]
[
x
1
x
2
x
3
]
−
10
⩽
0
[
1
1
0
]
[
x
1
2
x
2
2
x
3
2
]
+
[
1
1
−
1
]
[
x
1
x
2
x
3
]
−
50
⩽
0
[
1
0
0
]
[
x
1
2
x
2
2
x
3
2
]
+
[
2
2
1
]
[
x
1
x
2
x
3
]
−
40
⩽
0
[
1
0
0
]
[
x
1
2
x
2
2
x
3
2
]
+
[
0
0
1
]
[
x
1
x
2
x
3
]
−
2
=
0
[
−
1
−
2
0
−
1
1
0
]
[
x
1
x
2
]
⩾
[
−
1
0
]
x
1
⩾
0
\min w=(-1)*[ \begin{bmatrix}3 & 1 & 0 \end{bmatrix} \begin{bmatrix} x_1^2 \\ x_2^2 \\ x_3^2 \end{bmatrix} + \begin{bmatrix} 2 & 3 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} ]\\ \text { s. t. }\left\{ \begin{array}{l} \begin{bmatrix}2 & 2 & 0 \end{bmatrix} \begin{bmatrix} x_1^2 \\ x_2^2 \\ x_3^2 \end{bmatrix} + \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} -10 \leqslant 0 \\ \begin{bmatrix}1 & 1 & 0 \end{bmatrix} \begin{bmatrix} x_1^2 \\ x_2^2 \\ x_3^2 \end{bmatrix} + \begin{bmatrix} 1 & 1 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} -50 \leqslant 0 \\ \begin{bmatrix}1 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_1^2 \\ x_2^2 \\ x_3^2 \end{bmatrix} + \begin{bmatrix} 2 & 2 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} -40 \leqslant 0 \\ \begin{bmatrix}1 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_1^2 \\ x_2^2 \\ x_3^2 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} -2 = 0 \\ \begin{bmatrix} -1 & -2 & 0 \\ -1 & 1 & 0 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \geqslant \begin{bmatrix} -1 \\ 0 \end{bmatrix} \\ \\ x_{1} \geqslant 0 \end{array}\right. \\
minw=(−1)∗[[310]
x12x22x32
+[231]
x1x2x3
] s. t. ⎩
⎨
⎧[220]
x12x22x32
+[111]
x1x2x3
−10⩽0[110]
x12x22x32
+[11−1]
x1x2x3
−50⩽0[100]
x12x22x32
+[221]
x1x2x3
−40⩽0[100]
x12x22x32
+[001]
x1x2x3
−2=0[−1−1−2100][x1x2]⩾[−10]x1⩾0
3.程序
编写目标函数
function f = fun_f(x)
f = [3, 1, 0] * (x.^2) + [2, 3, 1] * x;
f = -f;
end
编写非线性约束函数
function [A, Aeq] = fun_A(x)
A = [[2, 2, 0] * (x.^2) + [1, 1, 1] * x - 10; ...
[1, 1, 0] * (x.^2) + [1, 1, -1] * x - 50; ...
[1, 0, 0] * (x.^2) + [2, 2, 1] * x - 40];
Aeq = [1, 0, 0] * (x.^2) + [0, 0, 1] * x - 2;
end
主程序
clc, clear
A = [-1, -2, 0; ...
-1, 0, 0];
b = [-1, 0];
[x, w] = fmincon('fun_f', rand(3, 1), A, b, [], [], [], [], 'fun_A');
x
w = -w
4.结果
求得的最优解是：
x 1 = 2.3333 , x 2 = 0.1667 , x 3 = − 3.4444 max f ( x ) = 18.0833 x_1=2.3333, x_2=0.1667, x_3=-3.4444 \\ \max f(x)=18.0833 x1=2.3333,x2=0.1667,x3=−3.4444maxf(x)=18.0833
如果这篇文章对你有帮助，欢迎点赞与收藏~
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
