数学建模笔记—— 非线性规划

非线性规划
- 1. 模型原理
- - 1.1 非线性规划的标准型
  - 1.2 非线性规划求解的Matlab函数
- 2. 典型例题
- 3. matlab代码求解
- - 3.1 例1 一个简单示例
  - 3.2 例2 选址问题
  - - 1. 第一问线性规划
    - 2. 第二问非线性规划

非线性规划

非线性规划是一种求解目标函数或约束条件中有一个或几个非线性函数的最优化问题的方法。运筹学的一个重要分支。20世纪50年代初,库哈(H.W.Kuhn)和托克(A.W.Tucker)提出了非线性规划的基本定理，为非线性规划奠定了理论基础。这一方法在工业、交通运输、经济管理和军事等方面有广泛的应用，特别是在“最优设计”方面，它提供了数学基础和计算方法，因此有重要的实用价值。

非线性规划模型特点：

模型中至少一个变量是非线性，即包含 $x^2,e^x,\frac1x,\sin x,\log_2x$ 等形式
线性规划有通用求准确解的方法(单纯形法),它的最优解只存在于可行域的边界上；非线性规划的最优解(若存在)可能在其可行域的任一点达到，目前非线性规划还没有适合各种问题的一般解法，各种方法都有其特定的适用范围

1. 模型原理

1.1 非线性规划的标准型

$min\quad f(x)\\\text{s.t.}\begin{cases}Ax\leq b, Aeq\cdot x=beq&\text{(线性)}\\c\big(x\big)\leq0, Ceq\big(x\big)=0&\text{(非线性)}\\lb\leq x\leq ub\end{cases}$

1.2 非线性规划求解的Matlab函数

$f min co n$ 函数： $[x\:fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,option)$

$f u n ：$ 把目标函数定义为一个单独的函数文件(min)
$x 0 :$ 决策变量的初始值
$A, b :$ 线性约束的不等式变量系数矩阵和常数项矩阵 $\le或<$
$A e q, b e q :$ 线性约束的等式变量系数矩阵和常数项矩阵
$l b, u b :$ 决策变量的最小取值和最大取值
$n o n l co n :$ 非线性约束，包括不等式和等式
$o pt i o n :$ 求解非线性规划使用的方法

注意：

非线性规划中对于初始值 $x 0$ 的选取非常重要，因为非线性规划的算法求解出来的是一个局部优化解。如果要求全局最优解，有两个思路：
- 给定不同初始值，在里面找到一个最优解；
- 先用蒙特卡罗模拟，得到一个蒙特卡罗解，然后将这个解作为初始值来求最优解。
$o pt i o n$ 选项可以给定求解的算法，一共有五种，interior-point(内点法)、trust-region-reflective(信赖域反射法)、sqp(序列二次规划法)、sqp-legacy(约束非线性优化算法)、active-set (有效集法)。不同的算法有其各自的优缺点和适用情况，我们可以改变求解的算法来对比求解的结果。
$ fun $表示目标函数，我们要编写一个独立的”m文件“储存目标函数
$n o n l co n$ 表示非线性部分的约束，也要编写一个独立的”m文件“存储非线性约束条件
决策变量的下表要改括号，比如 $x_1$ 要改为 $x (1)$ ,matlab才能识别
若不存在某种约束，可以用”[]“代替，若后面全为"[]“且option使用默认，后面的”[]"可以省略

2. 典型例题

选址问题：

临时料场： $A (5, 1)$ , $A (2, 7);$ 日储量各20吨

工地位置坐标及日需求量
Ⅰ Ⅱ Ⅲ Ⅳ Ⅴ Ⅵ
横坐标 1.25 8.75 0.5 5.75 3 7.25
纵坐标 1.25 0.75 4.75 5 6.5 7.25
日需求量 3 5 4 7 6 11

(1)试制定每天的供应计划，即从两料场分别向各工地运送多少吨水泥，使总的地千米数最小?

(2)为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各为20吨，问应建在何处，节省的吨千米数为多大?

工地位置坐标及日需求量
	Ⅰ	Ⅱ	Ⅲ	Ⅳ	Ⅴ	Ⅵ
横坐标	1.25	8.75	0.5	5.75	3	7.25
纵坐标	1.25	0.75	4.75	5	6.5	7.25
日需求量	3	5	4	7	6	11

确定决策变量

设第 $i$ 个工地的坐标 $a_i,b_i)$ ，水泥日用量 $d_i,i=1,2,\dots,6$ ，料场位置 $x_j,y_j)$ ，日储量 $e_j,j=1,2$ ；从料场 $j$ 向工地 $i$ 的运送量为 $x_{ij}$
确定约束条件
- 料场水泥运输总量不超过其日储量： $\sum_{i=1}^{6}x_{ij}\leq e_{j} ,j=1 ,2$
- 两个料场向某工地运输量之和等于该工地水泥日用量： $\sum_{j=1}^{2}x_{ij}=d_{i} ,i=1 ,2 ,\cdots,6$
确定目标函数

求总吨千米数最小，即运送量乘运送距离求和最小： $\min f=\sum_{j=1}^{2}\sum_{i=1}^{6}x_{ij}\sqrt{\left(x_{j}-a_{i}\right)^{2}+\left(y_{j}-b_{i}\right)^{2}}$
建立模型
$\begin{aligned}&min\quad f=\sum_{j=1}^{2}\sum_{i=1}^{6}x_{ij}\sqrt{\left(x_{j}-a_{i}\right)^{2}+\left(y_{j}-b_{i}\right)^{2}}\\&s.t.\begin{cases}\sum_{i=1}^{6}x_{ij}\leq e_{j},j=1,2\\\sum_{j=1}^{2}x_{ij}=d_{i},i=1,2,\ldots,6\\x_{ij}\geq0,i=1,2,\ldots,6;j=1,2\end{cases}\end{aligned}$
求解
- 对于第一问：因料场位置已知，故决策变量仅为 $x_{ij}$ ，为线性规划模型
- 对于第二问：新料场位置未知，所以 $x_j$ 和 $y_j$ 均为变量，且不是线性的，故为非线性规划模型
- 共有8个约束
  $\begin{aligned}&min\quad f=\sum_{j=1}^{2}\sum_{i=1}^{6}x_{ij}\sqrt{\left(x_{j}-a_{i}\right)^{2}+\left(y_{j}-b_{i}\right)^{2}}\\ &s.t.\begin{cases}\sum_{i=1}^{6}x_{ij}\leq e_{j},j=1,2\left(x_{11}+x_{21}+\ldots+x_{61}\leq e_{1},x_{12}+x_{22}+\ldots+x_{62}\leq e_{2}\right)\\ \sum_{j=1}^{2}x_{ij}=d_{i},i=1,2,\ldots,6\begin{pmatrix}x_{11}+x_{12}=d_1,\ldots, x_{61}+x_{62}=d_6\end{pmatrix}\\ x_{ij}\geq0,i=1,2,\ldots,6;j=1,2\end{cases}\end{aligned}$
  注意：在matlab里这些双角标的变量要改为单角标的变量，如 $x_{11}\to x_{1} ,\quad x_{21}\to x_{2} ,\quad\ldots ,\quad x_{62}\to x_{12}$

3. matlab代码求解

3.1 例1 一个简单示例

求解：
$\begin{aligned}&min\quad\mathrm{y}=x_{1}^{2}+x_{2}^{2}-x_{1}x_{2}-2x_{1}-5x_{2},\\&\mathrm{s.t.}\begin{cases}-\left(x_1-1\right)^2+x_2\geq0,\\2x_1-3x_2+6\geq0\end{cases}\end{aligned}$
非线性规划的目标函数fun1.m:

function f=fun1(x)
%FUN1 非线性规划的目标函数
%   这里的f实际上就是目标函数，函数的返回值也是f
%   输入值x实际上就是决策向量，由x1和x2组成的向量
% min f(x)=x1^2+x2^2-x1*x2-2x1-5x2f=x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2);
end

非线性规划中的非线性约束nonlfun1.m:

function [c ceq]=nonlfun1(x)
%NONLFUN1 非线性规划中的非线性约束，c为非线性不等式约束，ceq为非线性等式约束
%   输入值x为决策变量
%   返回值为 c（非线性不等式约束），ceq（非线性等式约束）
%   -(x1-1)^2+x2>=0c=(x(1)-1)^2-x(2);ceq=[];
end

给定任意初始值进行求解：

clear;
clc;
format long g %将matlab的计算结果显示为一般的长数字格式（默认保留两位小数或者使用科学计数法）
% min f(x)=x1^2+x2^2-x1*x2-2x1-5x2
% s.t. -(x1-1)^2+x2>=0; 2x1-3x2+6>=0
x0=[0 0];%任意给定一个初始值
A=[-2 3];
b=6;
disp("使用内点法求解：")
[x,fval]= fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)% 默认使用内点法
disp("使用SQP求解：")
option=optimoptions('fmincon','Algorithm','sqp');
[x,fval]= fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)

输出：

使用内点法求解：找到满足约束的局部最小值。优化已完成，因为目标函数沿
可行方向在最优性容差值范围内呈现非递减，
并且在约束容差值范围内满足约束。<停止条件详细信息>x =2.99941592955142          3.99922426270024fval =-12.9999995101786使用SQP求解：找到满足约束的局部最小值。优化已完成，因为目标函数沿
可行方向在最优性容差值范围内呈现非递减，
并且在约束容差值范围内满足约束。<停止条件详细信息>x =3.00000000090774          4.00000000060516fval =-13

使用蒙特卡罗的方法来找初始值在进行非线性规划求解：

%% 使用蒙特卡罗的方法来找初始值（推荐）
clc;
clear;
n=10000000;%生成的随机数组数
x1=unifrnd(-100,100,n,1); %生成在[-100,100]之间均匀分布的随机数组成n行1列的向量构成x1
x2=unifrnd(-100,100,n,1); %生成在[-100,100]之间均匀分布的随机数组成n行1列的向量构成x1
fmin=+inf; % 初始化函数f的最小值为正无穷
for i=1:nx=[x1(i),x2(i)];%构造x向量if ((x(1)-1)^2-x(2)<=0) && (-2*x1(i)-3*x2(i)-6<=0)result=x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2);if result<fminfmin=result;x0=x;endend
end
disp("蒙特卡罗选取的初始值为：")
disp(x0)
A=[-2,3];
b=6;
[x,fval]= fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)

输出：

蒙特卡罗选取的初始值为：3.00156691366464          4.03556138516905找到满足约束的局部最小值。优化已完成，因为目标函数沿
可行方向在最优性容差值范围内呈现非递减，
并且在约束容差值范围内满足约束。<停止条件详细信息>x =2.9992425325257          3.99899914772717fval =-12.9999991826508

3.2 例2 选址问题

模型：
$\begin{aligned}&min\quad f=\sum_{j=1}^{2}\sum_{i=1}^{6}x_{ij}\sqrt{\left(x_{j}-a_{i}\right)^{2}+\left(y_{j}-b_{i}\right)^{2}}\\&s.t.\begin{cases}\sum_{i=1}^{6}x_{ij}\leq e_{j},j=1,2\\\sum_{j=1}^{2}x_{ij}=d_{i},i=1,2,\ldots,6\\x_{ij}\geq0,i=1,2,\ldots,6;j=1,2\end{cases}\end{aligned}$

1. 第一问线性规划

代码：

%% 第一问：线性规划
clear
clc
% 6个工地坐标
a=[1.25 8.75 0.5 5.75 3 7.25];
b=[1.25 0.75 4.75 5 6.5 7.75];
% 临时料场位置
x=[5,2];
y=[1,7];
% 6个工地水泥日用量
d=[3 5 4 7 6 11];
% 计算目标函数系数，即六个工地与两个料场的距离，总共12个值
l=zeros([6,2]);
for i=1:6for j=1:2l(i,j)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);end
end
f=[l(:,1);l(:,2)]; % 目标函数系数向量，共12个值
% 不等式约束条件的变量系数和常数项
% 双下标转换成单下标：x11=x1,x21=x2,...,x62=x12
A=[ 1 1 1 1 1 1 0 0 0 0 0 00 0 0 0 0 0 1 1 1 1 1 1];
% 两个临时料场日储量
b=[20;20];% 矩阵的行数等于约束条件的个数，列是变量的个数
% 等式约束的变量系数和常数项
Aeq=[eye(6),eye(6)];
beq=[d(1);d(2);d(3);d(4);d(5);d(6)];
% 所有变量的下限全为0
Vlb=[0 0 0 0 0 0 0 0 0 0 0 0];
disp("第一问：")
[x,fval]=linprog(f,A,b,Aeq,beq,Vlb)

输出：

第一问：找到最优解。x =3507010040610fval =136.227519883182

2. 第二问非线性规划

非线性规划的目标函数fun2.m定义如下：

function f = fun2(x)
%FUN2 非线性规划的目标函数
%   这里的f实际上就是目标函数，函数的返回值也是f
%   输入值x实际上就是决策向量，由x1和x2组成的向量
% x前面12个是每个工地运输多少，后面四个为料场坐标
% 6个工地坐标
a=[1.25 8.75 0.5 5.75 3 7.25];
b=[1.25 0.75 4.75 5 6.5 7.75];
n=0;
f=0;
for j=13:2:16for i=1:6 n=n+1;f=f+x(n)*(sqrt((a(i)-x(j))^2+(b(i)-x(j+1))^2));end
end
end

求解代码：

%% 第二问 非线性规划
%注意，第二问中求新料场的位置，所以两个料场的横纵坐标也是变量，所以多了四个变量
% 对新坐标没有不等式约束，所以其不等式约束条件里面的系数为0
A2=[1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0];
B2=[20 ;20];
% 对新坐标也没有等式约束，所以相应项也为0
Aeq2=[eye(6),eye(6),zeros(6,4)];
beq2=[3 5 4 7 6 11]';
vlb2=[zeros(12,1);-inf;-inf;-inf;-inf];
% 非线性规划必须设置初始值，可以基于问题情况来设，设置rand()随机树等等
% 初始值设置为线性规划的计算结果，即临时料场的坐标
x0=[3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]';
disp("第二问")
[x2,fval2]=fmincon(@fun2,x0,A2,B2,Aeq2,beq2,vlb2)
% 注意，若约束条件里面有非线性函数，可在fmincon里使用nonlcon项

输出：

第二问可能存在局部最小值。满足约束。fmincon 已停止，因为当前步长小于
步长容差值并且在约束容差值范围内满足约束。<停止条件详细信息>x2 =2.940992189650064.840558656820213.87793737832666.943070403079231.303130777556460.02206785079521080.05900781034994370.1594413431797950.1220626216734020.05692959692077314.6968692224435410.97793214920485.729798455203994.975789925150467.249999954976637.74999993108167fval2 =90.4919073875194

第二问中可以使用蒙特卡罗方法求得近似值作为初始值：

求解过程中的不等式约束函数constraint.m如下

function [g,k] = constraint(x)
%CONSTRAINT 不等式约束条件
%   sum(x(:,1:6),2)是对矩阵前6列按行求和，即对前6个元素求和
%   对于6个工地接收第一个料场的总量。再减去20，即把不等式右边常数项移到左边
g=[sum(x(:,1:6),2)-20sum(x(:,7:12),2)-20];
%   等式约束条件，6个工地从两个料场收到总量分别为3，5，4，7，6，11
k=[x(1)+x(7)-3x(2)+x(8)-5x(3)+x(9)-4x(4)+x(10)-7x(5)+x(11)-6x(6)+x(12)-11];
end

求解过程：

%% 若有条件，可使用蒙特卡罗法求一个近似的解作为初始值
p0=inf;
n=10^6;
ticfor i =1:n% 前12个数是6个工地从料场接收的量，不会超过日需求量，为了加速计算取整数% 后四个变量是料场的横纵坐标，根据题目工地的坐标都在0-9，这里也取该范围x_m=[randi(4)-1,randi(6)-1,randi(5)-1,randi(8)-1,randi(7)-1,randi(12)-1,...randi(4)-1,randi(6)-1,randi(5)-1,randi(8)-1,randi(7)-1,randi(12)-1,...9*rand(1,4)];% 约束条件[g,k]=constraint(x_m);if all(g<=0) % 等式约束难以满足，此处相差不大即可算近似if all(abs(k)<=1)ff=fun2(x_m); %目标函数if ff<p0x_m0=x_m;p0=ff;endendend
end
x_m0,p0,toc
disp("以蒙特卡罗求得近似值作为初始值的线性规划结果")
[x3,fval3]=fmincon(@fun2,x_m0,A2,B2,Aeq2,beq2,vlb2)

输出：

x_m0 =列 1 至 40                         0                         0                         5列 5 至 83                         9                         2                         4列 9 至 123                         1                         2                         1列 13 至 166.85179793730359          7.45156987818458          5.78450172159294          4.84001343131759p0 =87.3564817958772历时 2.518048 秒。
以蒙特卡罗求得近似值作为初始值的线性规划结果可能存在局部最小值。满足约束。fmincon 已停止，因为当前步长小于
步长容差值并且在约束容差值范围内满足约束。<停止条件详细信息>x3 =列 1 至 40.0267009295509488          4.82102444621973        0.0235678116431545         0.426803450498299列 5 至 80.0304197766741595           10.979350696337          2.97329907044905         0.178975553780268列 9 至 123.97643218835685           6.5731965495017          5.96958022332584         0.020649303662989列 13 至 167.2500000010943          7.74999998555883          3.22063993810178          5.66691666664995fval3 =85.9490103544715