lingo

昨天做完了1996年的数学建模题最优捕鱼策略，感觉基理题可以说难也可以说不难。

开始建模的时候模型的制定我们过于简单，经过老师的提醒我们在第一天下午建立了微分方程模型，思路挺清晰的，但是在求解方面停留了很久，lingo求解比较方便但是有些错误会看不出，建模的思路要有条理的展现在代码上。主要lingo软件能力其实还是比不上MATLAB，遇到大数据就会不稳定，这次题目上有十几个数量级的，计算时容易出现一些错误。下次计算时最好通过降级处理得到结果。

下面展示第一题的代码

max=17.86*k3*x03/(k3+r)*(1-@exp(-2/3*(k3+r)))+22.99*k4*x04/(k4+r)*(1-@exp(-2/3*(k4+r)));
n4=1.109e+5;
k3=0.42*k4;
r=0.8;
n=0.5*n4*x03*@exp(-2/3*(r+k3))+n4*x04*@exp(-2/3*(r+k4));
c=1.22e+11/(1.22e+11+n);
c*n=x01;
x01*@exp(-r)=x02;
x02*@exp(-r)=x03;
x03*@exp(-2/3*k3-r)+x04*@exp(-4/3*k4-r)=x04;
第二题的结果不是标准答案，感觉思路也是对的，结果少了一点，捕捞系数过小

sets:
yu/1,2,3,4/:thega;
shiqi/1,2/:;
nian/1,2,3,4,5,6/:;
link(nian,shiqi):y,r,san,si;
endsets

max=w3*@sum(nian(i)|i#lt#6:k3*san(i,1)/(k3+s)*(1-@exp(-2/3*(k3+s))))+w4*@sum(nian(j)|j#lt#6:k4*si(j,1)/(k4+s)*(1-@exp(-2/3*(k4+s))));
y(1,1)=122e+1;r(1,1)=29.7e+1;san(1,1)=10.1e+1;si(1,1)=3.29e+1;s=0.8;k3=0.42*k4;a=0.5*n4*@exp(-2/3*(s+k3));b=n4*@exp(-2/3*(s+k4));
n4=1.109e+5;w3=17.86;w4=22.99;
thega(1)=@exp(-s);thega(2)=@exp(-s);thega(3)=@exp(-s-2/3*k3);thega(4)=@exp(-s-2/3*k4);
@for(nian(i)|i#lt#6:y(i,2)=y(i,1)*thega(1));
@for(nian(i)|i#lt#6:r(i,2)=r(i,1)*thega(2));
@for(nian(i)|i#lt#6:san(i,2)=san(i,1)*thega(3));
@for(nian(i)|i#lt#6:si(i,2)=si(i,1)*thega(4));
@for(nian(i)|i#gt#1:y(i,1)=(san(i-1,1)*a+si(i-1,1)*b)*1.22e+3/(1.22e+3+san(i-1,1)*a+si(i-1,1)*b));
@for(nian(i)|i#gt#1:r(i,1)=y(i-1,2));
@for(nian(i)|i#gt#1:san(i,1)=r(i-1,2));
@for(nian(i)|i#gt#1:si(i,1)=san(i-1,2)+si(i-1,2));
!y(6,1)>0.75*y(1,1);!r(6,1)>0.75*r(1,1);san(6,1)>0.75*san(1,1);si(6,1)>0.75*si(1,1);
y1=y(6,1)/y(1,1);r1=r(6,1)/r(1,1);san1=san(6,1)/san(1,1);si1=si(6,1)/si(1,1);
对于这种非线性优化问题MATLAB可以求解，但是第二题就比较麻烦；下午学习了MATLAB的非线性编程，

function feixingguangli
clear;clc
x0=[0 0 0 0 0 0 ];
sl=-pi/6*ones(1,6);
su=pi/6*ones(1,6);
[x,fval]=fmincon(@feixfun,x0,[],[],[],[],sl,su,@feixfun1)
end
function f=feixfun(x)
f=sum(x.^2);
end
function [f,g]=feixfun1(x)
a=[];
b=[];
c=[];
h=[];
v=800;
t=0.283;
X=[150,85,150,145,130,0];
Y=[140,85,155,50,150,0];
A=[243 236 220.5 159 230 52];
A=A*pi/180;
for i=1:5
    for j=i+1:6
        t1=v.^2*(cos(A(i)+x(i))-cos(A(j)+x(j))).^2+v.^2*(sin(A(i)+x(i))-sin(A(j)+x(j))).^2;
        a=[a t1];
        t2=2*v*(cos(A(i)+x(i))-cos(A(j)+x(j)))*(X(i)-X(j))+2*v*(sin(A(i)+x(i))-sin(A(j)+x(j)))*(Y(i)-Y(j));
        b=[b t2];
        t3=(X(i)-X(j)).^2+(Y(i)-Y(j)).^2;
        c=[c t3];
    end
end
for k=1:15
    if  -b(k)/(2*a(k))<=0
        h=[h;64-c(k)];
    else if -b(k)/(2*a(k))>0.283
        h=[h;64-a(k)*t^2-b(k)*t-c(k)];
    else h=[h;64-(4*a(k).*c(k)-b(k).^2)/(4*a(k)+eps)];
        end
    end
end
f=h;g=[];
end

这是1995年建模题飞行管理问题。