飞行管理问题的优化模型

MATLAB:

%约束条件函数

function [c,ceq]=aircon(delta)

x0=[150, 85, 150, 145, 130, 0];

y0=[140, 85, 155, 50, 150, 0];

alpha0=[243,236,220.5,159,230,52];

%第一组数据;

%x0=[60, 70, 80, 50, 40, 0];

%y0=[100, 100, 100, 100, 100, 40];

%alpha0=[270, 270, 270, 270, 270, 0];

%第二组数据;

%x0=[60, 60, 60, 60, 60, 0];

%y0=[80, 70, 60, 90, 100, 60];

%alpha0=[180, 180, 180, 180, 180, 0];

%第三组数据;

%x0=[0, 55, 90, 40, 80, 0];

%y0=[70, 5, 60, 130, 5, 60];

%alpha0=[0, 90, 180, 270, 180, 0];

alpha0=alpha0*pi/180; v=800;

co=cos(alpha0+delta);

si=sin(alpha0+delta);

for i=2:6

for j=1:i-1

b(i,j)=(x0(i)-x0(j))*(co(i)-co(j))+(y0(i)-y0(j))*(si(i)-si(j));

a(i,j)=v*((co(i)-co(j))^2+(si(i)-si(j))^2);

t(i,j)=-b(i,j)/a(i,j);

if t(i,j)<0

d(i,j)=1000;

else

d(i,j)=(x0(i)-x0(j)+v*t(i,j)*(co(i)-co(j)))^2+(y0(i)-y0(j)+v*t(i,j)*(si(i)-si(j)))^2;

end

end

end

c=64-[d(2,1), d(3,1:2), d(4,1:3), d(5,1:4), d(6,1:5)];

ceq=[ ];

/***************分割符****************************/

%目标函数

function f=airfun(delta)

f=delta*delta';

%f=sum(abs(delta));

/***************分割符****************************/

%计算程序

clear; clc

lb=-pi/18*ones(1,6); ub=pi/18*ones(1,6);

options=optimset('LargeScale','off');

n=20; fval=100;

for i=1:n

delta0=10*(rand(1,6)-0.5);

[dt,feval]=fmincon(@airfun,delta0,[ ],[ ], [ ], [ ],lb,ub,@aircon,options);

if (feval

fval=feval; dl=dt;

end

end

dl=dl*180/pi,

f=dl*dl'

%f=sum(abs(dl))

clc

x=[60, 70, 80, 50, 40, 0];

y=[100,100,100,100,100,40];

scatter(x,y,30,'r','filled');

for i=1:6

text(x(i)+2, y(i), num2str(i));

end

c=[5 4 2 1 6 3];

axis([-10, 195, -10, 170]);

grid on; hold on;

plot([0,160,160,0,0],[0,0,160,160,0],'b');

zt=[270,270,270,270,270,0];

%zt=[270, 270, 270, 270, 270, 0]+[6.04 5.56 5.08 -6.15 -5.67 4.78];

%zt=[270, 270, 270, 270, 270, 0]+ [-7.7172 5.5435 0.3477 -1.0448 6.0843 -3.4276];

zt1=zt*pi/180; vt=1;

dx=vt*cos(zt1); dy=vt*sin(zt1);

for vt=1:160

x1=x+dx; y1=y+dy;

scatter(x1,y1,11,c,'filled');

for j=1:5

for k=2:6

if k~=j

tx=x1(j)-x1(k); ty=y1(j)-y1(k);

dl=sqrt(tx^2+ty^2);

if dl<=8

fprintf('\n i=%d j=%d vt=%d', j,k,vt);

end

end

end

end

x=x1; y=y1; pause(0.1);

end

x=[60,60,60,60,60,0];

y=[80,70,60,90,100,60];

scatter(x,y,30,'r','filled');

for i=1:6

text(x(i)+2, y(i), num2str(i));

end

c=[5 4 2 1 6 3];

axis([-10, 195, -10, 170]);

grid on; hold on;

plot([0,160,160,0,0],[0,0,160,160,0],'b');

zt=[180,180,180,180,180,0];

%zt=[180,180,180,180,180,0]+[-1.29, -1.77, 9.74, -0.81 , -0.335, 5.58]; zt1=zt*pi/180; vt=1;

dx=vt*cos(zt1); dy=vt*sin(zt1);

for vt=1:160

x1=x+dx; y1=y+dy;

scatter(x1,y1,11,c,'filled');

for j=1:5

for k=2:6

if k~=j

tx=x1(j)-x1(k); ty=y1(j)-y1(k);

dl=sqrt(tx^2+ty^2);

if dl<=8

fprintf('\n i=%d j=%d vt=%d', j,k,vt);

end

end

end

end

x=x1; y=y1; pause(0.1);

end

clc

x=[0,55,90,40,80,0];

y=[70,5,60,130,5,60];

scatter(x,y,30,'r','filled');

for i=1:6

text(x(i)+2, y(i), num2str(i));

end

c=[5 4 2 1 6 3];

axis([-10, 195, -10, 170]);

grid on; hold on;

plot([0,160,160,0,0],[0,0,160,160,0],'b');

zt=[0,90,180,270,180,0];

%zt=[0,90,180,270,180,0]+[5.19, 6.13, 4.53, 0.0, 0.0, 5.67];

%zt=[0,90,180,270,180,0]+[0 -2.682928 -2.682928 0.0 0. 0.0];

zt1=zt*pi/180; vt=1;

dx=vt*cos(zt1); dy=vt*sin(zt1);

for vt=1:160

x1=x+dx; y1=y+dy;

scatter(x1,y1,11,c,'filled');

for j=1:5

for k=2:6

if k~=j

tx=x1(j)-x1(k); ty=y1(j)-y1(k);

dl=sqrt(tx^2+ty^2);

if dl<=8

fprintf('\n i=%d j=%d vt=%d', j,k,vt);

end

end

end

end

x=x1; y=y1; pause(0.1);

end

clc

x=[150, 85, 150, 145, 130, 0];

y=[140, 85, 155, 50, 150, 0];

scatter(x,y,30,'r','filled');

for i=1:6

text(x(i)+2, y(i), num2str(i));

end

c=[1 2 3 4 5 6];

%c=[5 4 2 1 6 3];

axis([-10, 195, -10, 170]);

grid on; hold on;

plot([0,160,160,0,0],[0,0,160,160,0],'b');

zt=[243,236,220.5,159,230,52];

%zt=[243,236,220.5,159,230,52]+[0, 0, 2.84, 0, 0, 0.79];

%zt=[243,236,220.5,159,230,52]+[0.0000 0.0001 0.0000 0.0000 -1.3846 -7.0638];

zt1=zt*pi/180; vt=1;

dx=vt*cos(zt1); dy=vt*sin(zt1);

for vt=1:160

x1=x+dx; y1=y+dy;

scatter(x1,y1,11,c,'filled');

for j=1:5

for k=2:6

if k~=j

tx=x1(j)-x1(k); ty=y1(j)-y1(k);

dl=sqrt(tx^2+ty^2);

if dl<=8

fprintf('\n i=%d j=%d vt=%d', j,k,vt);

%任意两架飞机距离小于8公里记录下来

end

end

end

end

x=x1; y=y1; pause(0.1);

end

%

clc; clear;

lb=-pi/18*ones(6,1); ub=pi/18*ones(6,1);

options=optimset('LargeScale','off');

n=20; fmin=100;

for i=1:n

z0=10*(rand(6,1)-0.5);

[x,fval]=fmincon(@fun1,z0,[ ], [ ], [ ], [ ], lb, ub, @fun2,options);

if fval

fmin=fval; xmin=x;

end

end

%xmin

fmin=xmin'*xmin;

xmin=xmin*180/pi

fmin=fmin*180/pi

function f=fun1(zt)

f=zt'*zt;

function [g,h]=fun2(zt)

zt1=zeros(6,1);

%x=[150, 85, 150, 145, 130, 0]';

%y=[140, 85, 155, 50, 150, 0]';

%zt0=[243,236,220.5,159,230,52]';

%第一组数据;

x=[60, 70, 80, 50, 40, 0]';

y=[100, 100, 100, 100, 100, 40]';

zt0=[270, 270, 270, 270, 270, 0]';

%第二组数据;

%x=[60, 60, 60, 60, 60, 0]';

%y=[80, 70, 60, 90, 100, 60]';

%zt0=[180, 180, 180, 180, 180, 0]';

%第三组数据;

%x=[0, 55, 90, 40, 80, 0]';

%y=[70, 5, 60, 130, 5, 60]';

%zt0=[0, 90, 180, 270, 180, 0]';

zt0=zt0*pi/180; v=800;

for i=1:6

zt1(i)=zt0(i)+zt(i);

end

for i=2:6

for j=1:i-1

dx(i,j)=x(i)-x(j); dy(i,j)=y(i)-y(j);

c(i,j)=cos(zt1(i))-cos(zt1(j));

s(i,j)=sin(zt1(i))-sin(zt1(j));

t(i,j)=-(dy(i,j)*s(i,j)+dx(i,j)*c(i,j))/(c(i,j)^2+s(i,j)^2)/v;

if (t(i,j)<0)

d(i,j)=1000;

else

d(i,j)=(dx(i,j)*s(i,j)-dy(i,j)*c(i,j))^2/(c(i,j)^2+s(i,j)^2);

end

end

end

g=64-[d(2,1), d(3,1:2), d(4,1:3), d(5,1:4), d(6,1:5)]';

h=[ ];

x=[150, 85, 150, 145, 130, 0];

y=[140, 85, 155, 50, 150, 0];

scatter(x,y,30,'r','filled');

for i=1:6

text(x(i)+2, y(i), num2str(i));

end

axis([-10, 195, -10, 170]);

grid on; hold on;

plot([0,160,160,0,0],[0,0,160,160,0],'g');

zt=[243,236,220.5,159,230,52];

zt1=zt*pi/180; a=40;

x1=x+a*cos(zt1); y1=y+a*sin(zt1);

for i=1:6

plot([x(i),x1(i)],[y(i),y1(i)],'b')

end

hold off

LEGO程序：

model:

sets:

fj/p1..p6/:zt0,dzt,zt1,x0,y0; !fj-飞机;

sj/t1..t901/:t; !sj--时间, 将时间离散化;

endsets

data:

x0=150,85,150,145,130,0;

y0=140,85,155,50,150,0;

zt0=243,236,220.5,159,230,52;

v=800; PI=3.14159265359;

enddata

min=@sum(fj:@abs(dzt)); !目标函数;

!min=@sum(fj:@abs(dzt*dzt));!目标函数;

@for(sj(k):t(k)=(k-1)/3000);

@for(fj:@bnd(-10,dzt,10)); !限制偏差角的变化范围;

!@for(fj:zt0=zt*PI/180);!把角度化为弧度;

@for(fj:zt1=zt0+dzt);

@for(sj(k):@for(fj(i)|i#lt#6:@for(fj(j )|j#gt#i:

(x0(i)-x0(j)+v*t(k)*(@cos(zt1(i)*PI/18 0)-@cos(zt1(j)*PI/180)))^2

+(y0(i)-y0(j)+v*t(k)*(@sin(zt1(i)*PI/1 80)-@sin(zt1(j)*Pi/180)))^2>=64))); End

model:

sets:

fj/p1..p6/:zt0,dzt,zt1,x0,y0; !fj-飞机;

sj/t1..t901/:t; !sj--时间, 将时间离散化;

endsets

data:

v=800; PI=3.14159265359;

!第一组数据;

x0=60,70,80,50,40,0;

y0=100,100,100,100,100,40;

zt0=270,270,270,270,270,0;

!第二组数据;

!x0=60,60,60,60,60,0;

!y0=80,70,60,90,100,60;

!zt0=180,180,180,180,180,0;

!第三组数据;

!x0=0,55,90,40,80,0;

!y0=70,5,60,130,5,60;

!zt0=0,90,180,270,180,0;

enddata

!min=@sum(fj:@abs(dzt));!目标函数;

min=@sum(fj:@abs(dzt*dzt)); !目标函数;

@for(sj(k):t(k)=(k-1)/3000);

@for(fj:@bnd(-10,dzt,10)); !限制偏差角的变化范围;

@for(fj:zt1=zt0+dzt);

@for(sj(k):@for(fj(i)|i#lt#6:@for(fj(j

)|j#gt#i:

(x0(i)-x0(j)+v*t(k)*(@cos(zt1(i)*PI/18 0)-@cos(zt1(j)*PI/180)))^2

+(y0(i)-y0(j)+v*t(k)*(@sin(zt1(i)*PI/1 80)-@sin(zt1(j)*Pi/180)))^2>=64))); end

model:

sets:

fj/p1..p6/:zt,zt0,dzt,zt1,x0,y0; !fj-飞机;

links(fj,fj):c,s,dx,dy;

endsets

data:

x0=150,85,150,145,130,0;

y0=140,85,155,50,150,0;

zt=243,236,220.5,159,230,52;

v=800;

enddata

min=@sum(fj:@abs(dzt)); !目标函数;

@for(fj:@bnd(-0.1,dxt,0.1)); !限制偏差角的变化范围;

@for(fj:zt0=zt*3.14159265359/180); !把角度化为弧度;

@for(fj:zt1=zt0+dzt);

!@for(links(i,j):c(i,j)=@cos(zt1(i))-@ cos(zt1(j)));

!@for(links(i,j):s(i,j)=@sin(zt1(i))-@ sin(zt1(j)));

!@for(links(i,j):dx(i,j)=x0(i)-x0(j)); !@for(links(i,j):dy(i,j)=y0(i)-y0(j)); @for(fj(i)|i#lt#6:@for(fj(j)|j#gt#i:(( x0(i)-x0(j))*(@cos(zt1(i))-@cos(zt1(j) ))

+(y0(i)-y0(j))*(@sin(zt1(i))-@sin(zt1( j))))/((@cos(zt1(i))-@cos(zt1(j)))^2 +(@sin(zt1(i))-@sin(zt1(j)))^2)/v<=0)) ;

@for(fj(i)|i#lt#6:@for(fj(j)|j#gt#i:(( x0(i)-x0(j))*(@sin(zt1(i))-@sin(zt1(j)

))

-(y0(i)-y0(j))*(@cos(zt1(i))-@cos(zt1( j))))^2/((@cos(zt1(i))-@cos(zt1(j)))^2 +(@sin(zt1(i))-@sin(zt1(j)))^2)>=64)); !@for(links(i,j)|i#lt#6#and#j#gt#i:-(c (i,j)*(x0(i)-x0(j))+s(i,j)*(y0(i)-y0(j )))/(c(i,j)^2+s(i,j)^2)/v>=0);

!@for(links(i,j)|i#lt#6#and#j#gt#i:(c( i,j)*(y0(i)-y0(j))-s(i,j)*(x0(i)-x0(j) ))^2/(c(i,j)^2+s(i,j)^2)>=64);

end