(数学建模)非线性规划

定义


如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问 题。
非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。
一般形式:






线性规划与非线性规划的区别


如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。


极值约束问题



二次规划


若某非线性规划的目标函数为自变量x的二次函数，约束条件又全是线性的，就称 这种规划为二次规划。
Matlab 中二次规划的数学模型可表述如下：



Matlab 中求解二次规划的命令是:


[x,fval]=quadprog(H,f,A,b,Aeq,Beq,lb,ub,x0,options)



例: 求解二次规划




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




罚函数法


利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题， 因而也称这种方法为序列无约束小化技术，简记为 SUMT .
罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函 数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法，下面介绍外罚函数法。



例: 求下列非线性规划



(1)定义增广目标函数,编写 M 文件 test.m


function g=test1(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); 



或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下：


function g=test2(x); 
M=50000; 
f=x(1)^2+x(2)^2+8; 
g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-x(2),0)+M*abs(-x(1)-x(2)^2+2); 



也可以修改罚函数的定义，编写test.m如下


function g=test3(x); 
M=50000; 
f=x(1)^2+x(2)^2+8; 
g=f-M*min(min(x),0)-M*min(x(1)^2-x(2),0)+M*(-x(1)-x(2)^2+2)^2; 



(2)求增广目标函数的极小值,在 Matlab 命令窗口输入


[x,y]=fminunc('test3',rand(2,1)) 



由于是非线性问题,很难求得问题的全局最优解,因此只能求得局部最优解,并且每次运行结果都不一样



  注:
  
  (1)如果非线性规划问题要求实时算法,则可以使用罚函数方法,但计算精度较低;
  (2)如果非线性规划问题不要求实时算法,但要求精度高,则可以使用Lingo软件编程求解或使用Matlab的fmincon命令



Matlab优化工具箱中的optimtool命令提供了优化问题的用户图形界面解法。 optimtool可应用到所有优化问题的求解，计算结果可以输出到Matlab工作空间中。


应用实例


在约10000m高空的某边长160km的正方形区域内，常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均又计算机记录其数据，以便进行飞行管理。当一家欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机相撞。如果会相撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，为避免碰撞。先假定条件如下：



  1.不碰撞的标准位任意两架飞机的距离大于8km;
   2.飞机飞行方向角调整的幅度不应超过30度
   3.所有飞机飞行速度均为每小时800km
   4.进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60km以上
   5.最多考虑六架飞机
   6.不必考虑飞机离开此区域后的情况



请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过 0.01 度），要求飞机飞行方向角调整的幅度尽量小。
 设该区域 4 个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)。
 记录数据见下表。


截屏2020-02-27上午8.12.11




  注：方向角指飞行方向与x轴正向的夹角。为方便之后的讨论，我们引进如下记号。
  
  $D$为飞行管理区域的边长。
   $Omega$为飞行管理区域，取直角坐标系使其为$[0, D] 	imes[0, D]]$为飞机飞行速度，a=800 mathrm{km} / mathrm{h}
  $left(x_{i}^{0}, y_{i}^{0}
ight)$为第$i$架飞机的初始位置。
   $(x_i(t),y_i(t))$为第$i$架飞机在$t$时刻的位置
   $	heta_{i}^{0}$为第$i$架飞机的原飞行方向角，即飞行方向与x轴夹角：$0 leq 	heta_{i}^{0}<2 pi$为第$i$架飞机的方向角调整，$-frac{pi}{6} leq Delta 	heta_{i} leq frac{pi}{6}$
   $theta_{i}=	heta_{i}^{0}+Delta 	heta_{i}$为第i架飞机调整后的飞行航向角




模型构建过程


根据相对运动的观点在考察两架飞机i和j的飞行时，可以将飞机i视为不动而飞机j以相对速度

截屏2020-02-27上午8.23.05

相对于飞机i运动，对上式进行和差化积,不妨设	heta_{j} geq 	heta_{i}，此时相对飞行方向角为eta_{i j}=frac{pi}{2}+frac{	heta_{i}+	heta_{j}}{2}

由于两架飞机的初始距离是：
r_{i j}(https://math.jianshu.com/math?formula=r_{i j}(0)%3Dsqrt{left(x_{i}^{0}-x_{j}^{0}
ight)^{2}%2Bleft(y_{i}^{0}-y_{j}^{0}
ight)^{2}})=sqrt{left(x_{i}^{0}-x_{j}^{0}
ight)^{2}+left(y_{i}^{0}-y_{j}^{0}
ight)^{2}}

截屏2020-02-27上午8.24.19

假设eta_{i j}^{0}为调整前j架飞机相对于第i架飞机的相对速度（矢量）与这两架飞机连线（从j指向i的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负）。得两架飞机不碰撞的条件为：
 left|eta_{i j}^{0}+frac{1}{2}left(Delta 	heta_{i}+Delta 	heta_{j}
ight)
ight|>alpha_{i j}^{0}
 其中：
 eta_{m n}^{0}=相对速度v_{mn}的幅角-从n指向m的连线矢量的幅角
 =arg frac{e^{i 	heta_{n}}-e^{i 	heta_{n}}}{left(x_{m}+i y_{m}
ight)-left(x_{n}+i y_{n}
ight)}
 （注意eta_{m n}^{0}表达式中的i表示虚数单位）这里我们利用复数的幅角，可以很方便地计算角度eta_{m n}^{0}(m, n=1,2, cdots, 6)

本问题中的优化目标函数可以有不同的形式：如使所有飞机的最大调整量最小；所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标函数，可以得到数学规划模型。


Matlab求解


利用matlab程序计算alpha_{i j}^{0}和eta_{i j}^{0}：


clc,clear
x0=[150 85 150 145 130 0];
y0=[140 85 155 50 150 0];
q=[243 236 220.5 159 230 52];
xy0=[x0; y0];
d0=dist(xy0); %求矩阵各个列向量之间的距离
d0(find(d0==0))=inf;
a0=asind(8./d0) %以度为单位的反函数
xy1=x0+i*y0
xy2=exp(i*q*pi/180)
for m=1:6
    for n=1:6
        if n~=m
            b0(m,n)=angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n)));
        end
    end
end
b0=b0*180/pi;
dlmwrite('txt1.txt',a0,'delimiter', '	','newline','PC');
fid=fopen('txt1.txt','a');
fwrite(fid,'~','char'); %往纯文本文件中写 LINGO 数据的分割符
dlmwrite('txt1.txt',b0,'delimiter', '	','newline','PC','-append','roffset', 1)



matlab会输出一个txt文件，将这个txt文件输入Lingo软件。


model:
sets:
plane/1..6/:delta;
link(plane,plane):alpha,beta;
endsets
data:
alpha=@file('txt1.txt'); !需要在alpha的数据后面加上分隔符"~";
beta=@file('txt1.txt');
enddata
min=@sum(plane:@abs(delta));
@for(plane:@bnd(-30,delta,30));
@for(link(i,j)|i#ne#j:@abs(beta(i,j)+0.5*delta(i)+0.5*delta(j))>a
lpha(i,j));
end



求解的最优解Delta 	heta_{3}=2.83858^{circ}, Delta 	heta_{6}=0.7908^{circ}，其他调整角度为0.


参考文章


https://www.jianshu.com/p/06ca582e8abe

https://blog.csdn.net/sunyueqinghit/article/details/81505409?depth_1-utm_source=distribute.pc_relevant.none-task&utm_source=distribute.pc_relevant.none-task

查看原文：https://upcwsh.top/ahhhh/190/
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/pteromyini/p/12370345.html