数学建模 非线性规划

一.非线性规划模型
1.概念:

如果目标函数或约束条件中包含非线性函数,就称该规划问题为非线性规划问题.一般来说,求解非线性规划问题要比求解线性规划问题困难得多,也不像线性规划问题那样有单纯形法这一通用方法,各个方法都有自己的适用范围

注意:如果线性规划的最优解存在,则其最优解只可能在其可行域的边界上(特别是顶点上);而如果非线性规划的最优解存在,其最优解可能在其可行域的任意位置

2.一般形式:

非线性规划问题的一般形式为:$$\begin{gathered} minf(x) \\ s.t.\{\begin{array}{l}{h}_j(x)\le 0(j=1,2...q)\\ g_i(x)=0(i=1,2...p)\end{array}(3.1) \end{gathered}$$其中$x=[x_1,x_2...x_n{]}^T$称为模型式$(3.1)$的决策变量,$f$为目标函数,$g_i,h_j$为约束函数,$g_i(x)=0$为等式约束,$h_j(x)\le 0$为不等式约束,目标函数和所有约束函数中至少有1个为非线性函数

某企业有$n$个项目可供选择投资,且至少要投资其中的1个项目.已知该企业的总资金为$A$,投资于第$i$个项目需花费的资金为$a_i$,并预计收益为$b_i$.请给出最佳的投资方案
设投资决策变量为$$x_i=\{\begin{array}{l}1,投资第i个项目\\ 0,不投资第i个项目\end{array}$$则投资总额为$\sum _{i=1}^n{a}_i{x}_i$,总收益为$\sum _{i=1}^n{b}_i{x}_i$.由于至少要投资1个项目,且总资金为$A$,故$$0<\sum _{i=1}^n{a}_i{x}_i\le A$$另外,由于$x_i$只能取$0$或$1$,故$$x_i(1-x_i)=0$$最佳方案应使投资额尽可能小而总收益尽可能大,故该问题的数学模型为$$\begin{gathered} maxQ=\frac{\sum _{i=1}^n{b}_i{x}_i}{\sum _{i=1}^n{a}_i{x}_i} \\ s.t.\{\begin{array}{l}0<\sum _{i=1}^n{a}_i{x}_i\le A\\ x_i(1-x_i)=0\end{array} \end{gathered}$$

3.MATLAB解法
(1)标准形式:

MATLAB中非线性规划的数学模型的标准形式为$$\begin{gathered} minf(x) \\ s.t.\{\begin{array}{l}Ax\le b\\ Aeq\cdot x=beq\\ c(x)\le 0\\ ceq(x)=0\\ lb\le x\le ub\end{array}(3.1) \end{gathered}$$其中$f(x)$为标量函数,$A,b,Aeq,beq,lb,ub$为相应维数的矩阵和向量,$c(x),ceq(x)$为非线性向量函数

(2)求解:

MATLAB中使用[x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,option)求解非线性规划问题.例如求解下述问题:$$\begin{gathered} minf(x)=x_1^2+x_2^2+x_3^2+8 \\ s.t.\{\begin{array}{l}{x}_1^2-x_2+x_3^2\ge 0\\ x_1+x_2^2+x_3^3\le 20\\ -x_1-x_2^2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3\ge 0\end{array} \end{gathered}$$计算代码为:

#func.mlx中:
function f=func(x);
f=sum(x.^2)+8;
end

#nonlcon.mlx中:
function [g,h]=nonlcon(x);
g=[-x(1)^2+x(2)-x(3)^2 x(1)+x(2)^2+x(3)^3-20]%非线性不等式约束
h=[-x(1)-x(2)^2+2 x(2)+2*x(3)^2-3]%非线性等式约束
end

#命令行中:
>> [x,y]=fmincon('func',rand(3,1),[],[],[],[],zeros(3,1),[],'nonlcon');
>> x,y

x =

    0.5522
    1.2033
    0.9478


y =

   10.6511

二.无约束问题
1.无约束极值问题
(1)符号解:

求多元函数$f(x,y)=x^3-y^3+3x^2+3y^2-9x$的极值.计算代码为:

>> syms x y
>> f=x^3-y^3+3*x^2+3*y^2-9*x;
>> df=jacobian(f);%求1阶偏导数
>> d2f=jacobian(df);%求2阶偏导数
>> [xx,yy]=solve(df);%求驻点
>> xx=double(xx);
>> yy=double(yy);
>> for i=1:length(xx)
a=subs(d2f,{x,y},{xx(i),yy(i)});
b=eig(a);%求矩阵的特征值
r=subs(f,{x,y},{xx(i),yy(i)});
r=double(r);
if all(b>0)
fprintf('极小值点(%f,%f),极小值%f\n',xx(i),yy(i),r);
elseif all(b<0)
fprintf('极大值点(%f,%f),极大值%f\n',xx(i),yy(i),r);
elseif any(b>0)&any(b<0)
fprintf('(%f,%f)非极值点\n',xx(i),yy(i));
else
fprintf('无法判断(%f,%f)\n',xx(i),yy(i))
end
end
极小值点(1.000000,0.000000),极小值-5.000000
(-3.000000,0.000000)非极值点
(1.000000,2.000000)非极值点
极大值点(-3.000000,2.000000),极大值31.000000

(2)数值解:

在MATLAB中,使用[x,fval]=fminunc(fun,x0,options)或[x,fval]=fminsearch(fun,x0,options)来求函数的极小值

求多元函数$f(x,y)=x^3-y^3+3x^2+3y^2-9x$的极值.计算代码为:

>> 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))

Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.

<stopping criteria details>


xy1 =

    1.0000
   -0.0000


z1 =

   -5.0000

>> [xy2,z2]=fminsearch(g,rand(2,1));
>> xy2,-z2

xy2 =

   -3.0000
    2.0000


ans =

   31.0000

求函数$f(x)=100(x_2-x_1^2{)}^2+(1-x_1{)}^2$的极小值.计算代码为:

#fun.mlx中:
function [f,g]=fun(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];%返回梯度向量
end
#命令行中:
>> options=optimset('GradObj','on');
>> [x,y]=fminunc('fun',rand(1,2),options)

Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.

<stopping criteria details>


x =

    1.0000    1.0000


y =

   7.3005e-17

求函数$f(x)=\sin (x)+3$的极小值.计算代码为:

>> f=@(x) sin(x)+3;
>> [x,y]=fminsearch(f,2)

x =

    4.7124


y =

    2.0000

2.函数的零点:

求多项式$f(x)=x^3-x^2+2x-3$的零点.计算代码为:

>> a=[1,-1,2,-3];
>> x0=roots(a)

x0 =

  -0.1378 + 1.5273i
  -0.1378 - 1.5273i
   1.2757 + 0.0000i

################################################################################################################################

>> syms x
>> x0=solve(x^3-x^2+2*x-3)%求零点的符号解
 
x0 =
 
 root(z^3 - z^2 + 2*z - 3, z, 1)
 root(z^3 - z^2 + 2*z - 3, z, 2)
 root(z^3 - z^2 + 2*z - 3, z, 3)
 
>> x0=vpa(x0,5)%转换为小鼠格式
 
x0 =
 
              1.2757
 - 0.13784 + 1.5273i
 - 0.13784 - 1.5273i

################################################################################################################################

>> y=@(x) x^3-x^2+2*x-3;
>> x=fsolve(y,rand)%求给定值附近的1个解

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>


x =

    1.2757

3.方程组的解:

如如下方程组的解:$$\{\begin{array}{l}{x}^2+y-6=0\\ y^2+x-6=0\end{array}$$计算代码为:

>> syms x y
>> [x,y]=solve(x^2+y-6,y^2+x-6)
 
x =
 
                2
               -3
 21^(1/2)/2 + 1/2
 1/2 - 21^(1/2)/2
 
 
y =
 
                2
               -3
 1/2 - 21^(1/2)/2
 21^(1/2)/2 + 1/2
 
################################################################################################################################

>> f=@(x) [x(1)^2+x(2)-6,x(2)^2+x(1)-6];
>> xy=fsolve(f,rand(2,1))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>


xy =

    2.0000
    2.0000

三.约束极值问题
1.概念:

带有约束条件的极值问题称为约束极值问题,也叫规划问题.求解约束极值问题要比求解无约束极值问题困难得多,通常采用将约束问题转换为无约束问题,非线性规划问题转换为线性规划问题的方法进行简化
库恩·塔克条件是非线性规划领域中的重要理论成果,是确定某点为最优点的必要条件(对凸规划,也是充分条件)

2.二次规划
(1)概念:

若某非线性规划的目标函数为二次函数,且约束条件均为线性约束条件,就称该规划为二次规划

(2)MATLAB中的标准形式:

MATLAB中二次规划问题的标准形式为:$$\begin{gathered} min\frac{1}{2}{x}^{T}H{x}_{+}{f}^{T}x \\ s.t.\{\begin{array}{l}Ax\le b\\ Aeq\cdot x=beq\\ lb\le x\le ub\end{array} \end{gathered}$$其中$H$为实对称矩阵,$f,b,beq,lb,ub$为列向量,$A,Aeq$为相应维数的矩阵

(3)求解:

MATLAB中使用[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)求解二次规划问题

求解下述二次规划问题:$$\begin{gathered} minf(x)=2x_1^2-4x_1{x}_2+4x_2^2-6x_1-3x_2 \\ s.t.\{\begin{array}{l}{x}_1+x_2\le 3\\ 4x_1+x_2\le 9\\ x_1,x_2\ge 0\end{array} \end{gathered}$$

>> h=[4,-4;-4,8];
>> f=[-6;-3];
>> a=[1,1;4,1];
>> b=[3;9];
>> [x,value]=quadprog(h,f,a,b,[],[],zeros(2,1))

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>


x =

    1.9500
    1.0500


value =

  -11.0250

3.罚函数法
(1)概念:

罚函数法又称序列无约束最小化技术(Sequential Unconstrained Minization Technique;SUMT).其基本思想是利用问题中的约束函数作出适当的罚函数,由此构造出带参数的增广目标函数,从而把问题转化为无约束非线性规划问题.该方法又分为外罚函数法和内罚函数法

(2)求解:

考虑问题$$\begin{gathered} minf(x) \\ s.t.\{\begin{array}{l}{g}_i(x)\le 0(i=1,2...r)\\ h_j(x)\ge 0(j=1,2...s)\\ k_m(x)=0(m=1,2...t)\end{array} \end{gathered}$$取1个充分大的数$M>0$,构造函数$$P(x,M)=f(x)+M\sum _{i=1}^{r}max(g_i(x),0)-M\sum _{j=1}^{s}min(h_j(x),0)+M\sum _{m=1}^t|k_m(x)|$$则以增广目标函数$P(x,M)$为目标函数的无约束极值问题$minP(x,M)$的最优解也是原问题的最优解

求解下述非线性规划问题$$\begin{gathered} minf(x)=x_1^2+x_2^2+8 \\ s.t.\{\begin{array}{l}{x}_1^2-x_2\ge 0\\ -x_1-x_2^2+2=0\\ x_1,x_2\ge 0\end{array} \end{gathered}$$

#test.mlx中:
function g=test(x);
M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+M*abs(-x(1)-x(2)^2+2);
end
#命令行中:
>> [x,y]=fminsearch('test',rand(2,1))

x =

    1.9991
    0.0296


y =

   11.9977

由于是非线性问题,很难求得全局最优解,只能求得1个局部最优解,且每次结果均不同

(3)适用:

①如果要求实时算法,可使用罚函数法,但精度较低
②如果不要求实时算法,但要求高精度,可使用Lingo或MATLAB中的fmincon()求解

3.其他方法
(1)有限区间上单变量非线性函数的最小值:

MATLAB使用[x,fval]=fminbnd(fun,x1,x2,options)来求解这类问题

求$f(x)=(x-3{)}^2-1(0\le x\le 5)$的最小值

#fun.mlx中:
function f=fun(x);
f=(x-3)^2-1
end
#命令行中:
>> [x,y]=fminbnd('fun',0,5)

f =

    0.1885


f =

   -0.9919


f =

   -0.3282


f =

    -1


f =

   -1.0000


f =

   -1.0000


x =

     3


y =

    -1

(2)存在非线性约束:

MATLAB中的标准形式为$$\begin{gathered} minf(x) \\ s.t.\{\begin{array}{l}Ax\le b\\ Aeq\cdot x=beq\\ lb\le x\le ub\\ c(x)\le 0\\ ceq(x)=0\\ K_i(x,w_i)\le 0(1\le i\le n)\end{array} \end{gathered}$$其中$x,b,beq,lb,ub$为向量,$A,Aeq$为矩阵,$c(x),ceq(x)$为向量函数,$K(x,w_i)$为标量函数,$w_i$为附加的变量

MATLAB中使用[x,fval]=fseminf(fun,x,ntheta,seminfcon,A,b,Aeq,beq,lb,ub)进行求解

求$f(x)=(x_1-0.5{)}^2+(x_2-0.5{)}^2+(x_3-0.5{)}^2$的最小值,约束为$$\begin{gathered} K_1(x,w_1)=\sin (w_1{x}_1)\cos (w_1{x}_2)-\frac{(w_1-50{)}^2}{1000}-\sin (w_1{x}_3)-x_3\le 1(1\le w_1\le 100) \\ {K}_2(x,w_2)=\sin (w_2{x}_2)\cos (w_2{x}_1)-\frac{(w_2-50{)}^2}{1000}-\sin (w_2{x}_3)-x_3\le 1(1\le w_2\le 100) \end{gathered}$$

#fun.mlx中:
function f=fun(x,s);
f=sum((x-0.5).^2);
end
#fun2.mlx中:
function [c,ceq,k1,k2,s]=fun2(x,s);
c=[];
ceq=[];
if isnan(s(1,1))
    s=[0.2,0;0.2,0];
end
w1=1:s(1,1):100;
w2=1:s(2,1):100;
k1=sin(w1*x(1)).*cos(w1*x(2))-((w1-50).^2)/1000-sin(w1*x(3))-x(3)-1;
k2=sin(w2*x(2)).*cos(w2*x(1))-((w2-50).^2)/1000-sin(w2*x(3))-x(3)-1;
plot(w1,k1,'-',w2,k2,'+');#见下图
end
#命令行中:
>> [x,y]=fseminf(@fun,x0,2,@fun2)

Local minimum possible. Constraints satisfied.

fseminf stopped because the size of the current search direction is less than
twice the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>


x =

    0.6675
    0.3012
    0.4022


y =

    0.0771

(3)双重最值:

求解$$\begin{gathered} \underset{x}{\min }\underset{i}{\max }{F}_i(x) \\ s.t.\{\begin{array}{l}Ax\le b\\ Aeq\cdot x=beq\\ c(x)\le 0\\ ceq(x)=0\\ lb\le x\le ub\end{array} \end{gathered}$$的MATLAB命令为[x,fval]=fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

求函数族 { f 1 ( x ) , f 2 ( x ) , f 3 ( x ) , f 4 ( x ) , f 5 ( x ) } \{f_1(x),f_2(x),f_3(x),f_4(x),f_5(x)\} {f1​(x),f2​(x),f3​(x),f4​(x),f5​(x)}取极小-极大值时的$x$值$$\{\begin{array}{l}{f}_1(x)=2x_1^2+x_2^2-48x_1-40x_2+304\\ f_2(x)=-x_1^2-3x_2^2\\ f_3(x)=x_1+3x_2-18\\ f_4(x)=-x_1-x_2\\ f_5(x)=x_1+x_2-8\end{array}$$

#fun.mlx中:
function f=fun(x);
f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
   -x(1)^2-3*x(2)^2
   x(1)+3*x(2)-18
   -x(1)-x(2)
   x(1)+x(2)-8];
end
#命令行中:
>> [x,y]=fminimax(@fun,rand(2,1))

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>


x =

    4.0000
    4.0000


y =

    0.0000
  -64.0000
   -2.0000
   -8.0000
   -0.0000

(4)利用梯度求解约束优化问题:

已知$f(x)=e^{x_1}(4x_1^2+2x_2^2+4x_1{x}_2+2x_2+1)$,求$$\begin{gathered} minf(x) \\ s.t.\{\begin{array}{l}{x}_1{x}_2-x_1-x_2\le -1.5\\ x_1{x}_2\ge -10\end{array} \end{gathered}$$

#fun.mlx中:
function [f,df]=fun(x);
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
df=[exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);
    exp(x(1))*(4*x(2)+4*x(1)+2)];%梯度
end
#fun2.mlx中:
function [c,ceq,dc,dceq]=fun2(x);
c=[x(1)*x(2)-x(1)-x(2)+1.5;
   -x(1)*x(2)-10];
dc=[x(2)-1,-x(2);
    x(1)-1,-x(1)];%不等式约束条件的梯度
ceq=[];
dceq=[];%等式约束条件的梯度
end
#命令行中:
>> options=optimset('GradObj','on','GradConstr','on');
>> [x,y]=fmincon(@fun,rand(2,1),[],[],[],[],[],[],@fun2,options)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>


x =

   -9.5474
    1.0474


y =

    0.0236

四.飞行管理问题
1.问题叙述:

在$10km$高空处某边长为$D=160km$的正方形区域(记为$\Omega$)中,经常有若干架飞机作水平飞行.当1架欲进入该区域的飞机到达区域边缘时,应立即记录其位置及速度,并判断其是否会与区域内其他飞机发生碰撞.如果会发生碰撞,则应计算如何调整各架飞机飞行的方向角(即飞行方向与x轴之间的夹角;记为${\theta }_i^0$)以避免碰撞.要求如下:
①不碰撞的标准为任意2架飞机的距离均大于$8km$
②飞机飞行的方向角调整的幅度(记为$\Delta {\theta }_i$)不应超过$30°$,新方向角记为${\theta }_i$
③所有飞机的速度均为$a=800km/h$
④进入该区域的飞机在到达区域边缘时,与区域内飞机的距离应大于$60km$
⑤最多考虑6架飞机
⑥不必考虑飞机离开此区域后的状况
⑦该区域的4个顶点的坐标分别为$(0,0),(160,0),(160,160),(0,160)$
⑧即第$i$架飞机的初始位置为$(x_i^0,y_i^0)$,在$t$时刻的位置为$(x_i(t),y_i(t))$
⑨初始飞行数据如下表:

飞机编号	横坐标 $x_i^0$	纵坐标 $y_i^0$	方向角 / °
1	150	140	243
2	85	85	236
3	150	155	220.5
4	145	50	159
5	130	150	230
6(新进入)	0	0	52

2.模型1: