您所在位置：主题：讨论
zongmei
[专家分：0] 发布于 2008-04-08 20:50:00
我编了一个优化程序
这个是目标函数（保存名为proportionalfun.m文件）
function f=proportionalfun(x,a,b,h,G,k1,k2,L)
%alpha Define the limit of variable x1.
%Ps Define the limit of variable x2.
%Ps0 Define the limit of variable x3.
%beta Define the limit of variable x4.
%z Define the limit of variable x5.
G=52528;
a=2358;
b=0942;
L=3300;
h=1400;
k1=2228;
k2=1805;
Fs1=k1*x(2);
Fs2=k2*x(2);
betaF=atan(tan(x(1))*tan(x(4)));
q=Fs1-Fs2*tan(atan(tan(x(1))*tan(x(4))));
q0=Fs2-Fs1*tan(atan(tan(x(1))*tan(x(4))));
zA=(Fs1+Fs2)/G;
zB=(G*(a-b*tan(betaF))+sqrt(G^2*(a-b*tan(betaF))^2-4*h*G*L*q0*(1+tan(betaF))))/2*h*G*(1+tan(betaF));
zC=(a-b*tan(x(1)))/h*(1+tan(x(1)));
zD=(G*(a-b*tan(betaF))-sqrt(G^2*(a-b*tan(betaF))^2-4*h*G*L*q*(1+tan(betaF))))/2*h*G*(1+tan(betaF));
%满载时的相关数据
if zA<=zC
phiA=zA*L/(1+tan(x(1)))*(b+zA*h);
elseif zA>zC
phiA=zA*L*tan(x(1))/(1+tan(x(1)))*(a-zA*h);
end
phiB=zB;
phiC=zC;
phiD=zD;
phif1=x(5)*L/((1+tan(x(1)))*(b+x(5)*h));
phir1=x(5)*L*tan(x(1))/((1+tan(x(1)))*(a-x(5)*h));
phif2=(L*(G*x(5)+tan(betaF)*Fs1-Fs2))/G*(b+x(5)*h)*(1+tan(betaF));
phir2=(L*(tan(betaF)*G*x(5)-tan(betaF)*Fs1+Fs2))/G*(a-x(5)*h)*(1+tan(betaF));
A1=1/tan(x(1));
Ef1=b/(A1*L-phif1*h);
Er1=a/((1-A1)*L+phir1*h);
Ef2=(phif2*G(1+tan(betaF))*b-(Fs1*tan(betaF)-Fs2)*L)/(L*G-G*(1+tan(betaF))*phif2*h)*phif2;
Er2=(phir2*G(1+tan(betaF))*a+(Fs1*tan(betaF)-Fs2)*L)/(L*G*tan(betaF)+G*(1+tan(betaF))*phir2*h)*phir2;
%满载时的目标函数
syms phi;
if zA<zC
if zB<=0.8
f=0.6-(int(Ef1,'phi',phiA,0.2)+int(Ef2,'phi',phiB,phiA)+int(Er2,'phi',0.8,0.2));
elseif zB>0.8
f=0.6-(int(Ef1,'phi',phiA,0.2)+int(Ef2,'phi',0.8,phiA));
end
elseif zA>=zC
if (zD<0.8)&(zB<0.8)
if zC<=0.2
f=0.6-(int(Er1,'phi',phiA,0.2)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',phiB,phiD)+int(Er2,'phi',0.8,phiB));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',phiA,phiC)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',phiB,phiD)+int(Er2,'phi',0.8,phiB));
end
elseif (zD<0.8)|(zB>=0.8)
if zC<=0.2
f=0.6-(int(Er1,'phi',phiA,0.2)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',0.8,phiD));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',phiA,phiC)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',0.8,phiD));
end
elseif (zD>=0.8)&(zA<0.8)
if zC<=0.2
f=0.6-(int(Er1,'phi',phiA,0.2)+int(Er1,'phi',0.8,phiA));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',phiA,phiC)+int(Er2,'phi',0.8,phiA));
end
else zA>=0.8
if zC<=0.2
f=0.6-(int(Er1,'phi',0.8,0.2));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',0.8,phiC));
end
end
end
这个是约束条件（保存名为proportionalfun.m文件）
function [c,ceq]=proportionalcon(x,a,b,h,G,k1,k2,L)
%alpha Define the limit of variable x1.
%Ps Define the limit of variable x2.
%Ps0 Define the limit of variable x3.
%beta Define the limit of variable x4.
%z Define the limit of variable x5.
%满载时的数据
G=52528;
a=2358;
b=942;
L=3300;
h=1400;
k1=2228;
k2=1805;
Fs1=k1*x(2);
Fs2=k2*x(2);
Fp11=0.15*G*(b+0.15*h)/L;
Fp12=0.15*G*(a-0.15*h)/L;
Fp21=0.3*G*(b+0.3*h)/L;
Fp22=0.3*G*(a-0.3*h)/L;
Fp31=0.45*G*(b+0.45*h)/L;
Fp32=0.45*G*(a-0.45*h)/L;
Fp41=0.8*G*(b+0.8*h)/L;
Fp42=0.8*G*(a-0.8*h)/L;
Fp51=G*(b+0.61*h)*(0.61+0.07)/0.85*L;
Fp52=0.61*G-Fp51;
Fp61=G*(b+0.15*h)*(0.15+0.07)/0.85*L;
Fp62=0.15*G-Fp61;
gamma1=atan(Fp12/Fp11);
gamma2=atan(Fp22/Fp21);
gamma3=atan(Fp32/Fp31);
gamma4=atan(Fp42/Fp41);
gamma5=atan(Fp52/Fp51);
gamma6=atan(Fp62/Fp61);
for i=1,4,5,6
if Fpi1>Fs1
deltai=atan((Fpi2-Fs2)/(Fpi1-Fs1));
elseif Fpi1<=Fs1
deltai=atan((Fpi2-Fs2)/(Fpi1-Fs1))+pi;
end
end
A1=G*((b+0.07*h)^2-4*0.07*b*h)+4*0.85*L*h*Fs1;
A2=G*(2*(b+0.07*h)*(b+0.07*h-0.85*L)-8*0.07*b*h)+4*0.85*L*h*(Fs1-Fs2);
A3=G*((b+0.07*h-0.85*L)^2-4*0.07*b*h)-4*0.85*L*h*Fs2;
A41=G*(1+tan(delta81))/0.85*L;
A42=G*(1+tan(delta82))/0.85*L;
delta81=atan((-A2+sqrt((A2)^2-4*A1*A3))/2*A1);
delta82=atan((-A2-sqrt((A2)^2-4*A1*A3))/2*A1);
z81=-(A41*b+0.07*A41*h-G)/2*A41*h;
z82=-(A42*b+0.07*A42*h-G)/2*A42*h;
if z81>=z82
delta8=delta81;
elseif z81<z82
delta8=delta82;
end
if z81>=z82
z8=z81;
elseif z81<z82
z8=z82;
end
A71=G*((b+0.07*h)^2-4*0.07*b*h);
A72=G*(2*(b+0.07*h)*(b+0.07*h-0.85*L)-8*0.07*b*h);
A73=G*((b+0.07*h-0.85*L)^2-4*0.07*b*h);
A741=G*(1+tan(delta71))/0.85*L;
A742=G*(1+tan(delta72))/0.85*L;
delta71=atan((-A72+sqrt((A72)^2-4*A71*A73))/2*A71);
delta72=atan((-A72-sqrt((A72)^2-4*A71*A73))/2*A71);
z71=-(A741*b+0.07*A741*h-G)/2*A741*h;
z72=-(A742*b+0.07*A742*h-G)/2*A742*h;
if z71>=z72
delta7=delta71;
elseif z71<z72
delta7=delta72;
end
if z71>=z72
z7=z71;
elseif z71<z72
z7=z72;
end
if x(1)>gamma1
Fe1=Fp11;
Fe2=Fp12;
elseif (x(1)>=gamma2)&(x(1)<=gamma1)|(x(1)>=gamma4)&(x(1)<=gamma3)
ze=(a-b*tan(x(1)))/h*(1+tan(x(1)));
Fe1=ze*G*(b+ze*h)/L;
Fe2=ze*G*(a-ze*h)/L;
else (x(1)>=gamma3)&(x(1)<=gamma2)
q0=Fp32-tan(tau)*Fp31;
tan(tau)=(Fp42-Fp32)/(Fp41-Fp31);
Fe1=q0/(tan(x(1))-tan(tau));
Fe2=q0*tan(x(1))/(tan(x(1))-tan(tau));
else x(1)<=gamma4
Fe1=+\infty;
Fe2=-\infty;
end
%c(i)(i=[1:13])非线性不等式约束
c(1)=x(1)/gamma1-1.0;
if (x(5)>=0.15)&(x(5)<=0.61)
c(2)=gamma7/x(1)-1.0;
elseif (z7<0.15)||(z7>0.61)
c(2)=max(gamma5,gamma6)/x(1)-1.0;
end
c(3)=Fs2/Fe2-1.0;
c(4)=-Fs2;
c(5)=betaF/delta4-1.0;
if (z8>=0.15)&(z8<=0.61)
c(6)=delta8/betaF-1.0;
elseif (z8<0.15)|(z8>0.61)
c(6)=max(delta5,delta6)/\betaF-1.0;
end
c(7)=betaF/x(1)-1.0;
c(8)=-betaF;
c(9)=-x(5);
c(10)=x(5)-1;
然后我调用x0=[atan(k2/k1),1.3835e6,8.2e6,0.552,0.15]；
x=fmincon(@proportionalfun1,x0,[],[],[],[],[],[],@proportionalfun1)
出现这种情况
??? Error using ==> fmincon
FMINCON cannot continue because user supplied objective function failed with the following error:
Attempted to access G(1.49894); index must be a positive integer or logical.
这是怎么回事啊！
[em18][em18][em18]
这个是目标函数（保存名为proportionalfun.m文件）
function f=proportionalfun(x,a,b,h,G,k1,k2,L)
%alpha Define the limit of variable x1.
%Ps Define the limit of variable x2.
%Ps0 Define the limit of variable x3.
%beta Define the limit of variable x4.
%z Define the limit of variable x5.
G=52528;
a=2358;
b=0942;
L=3300;
h=1400;
k1=2228;
k2=1805;
Fs1=k1*x(2);
Fs2=k2*x(2);
betaF=atan(tan(x(1))*tan(x(4)));
q=Fs1-Fs2*tan(atan(tan(x(1))*tan(x(4))));
q0=Fs2-Fs1*tan(atan(tan(x(1))*tan(x(4))));
zA=(Fs1+Fs2)/G;
zB=(G*(a-b*tan(betaF))+sqrt(G^2*(a-b*tan(betaF))^2-4*h*G*L*q0*(1+tan(betaF))))/2*h*G*(1+tan(betaF));
zC=(a-b*tan(x(1)))/h*(1+tan(x(1)));
zD=(G*(a-b*tan(betaF))-sqrt(G^2*(a-b*tan(betaF))^2-4*h*G*L*q*(1+tan(betaF))))/2*h*G*(1+tan(betaF));
%满载时的相关数据
if zA<=zC
phiA=zA*L/(1+tan(x(1)))*(b+zA*h);
elseif zA>zC
phiA=zA*L*tan(x(1))/(1+tan(x(1)))*(a-zA*h);
end
phiB=zB;
phiC=zC;
phiD=zD;
phif1=x(5)*L/((1+tan(x(1)))*(b+x(5)*h));
phir1=x(5)*L*tan(x(1))/((1+tan(x(1)))*(a-x(5)*h));
phif2=(L*(G*x(5)+tan(betaF)*Fs1-Fs2))/G*(b+x(5)*h)*(1+tan(betaF));
phir2=(L*(tan(betaF)*G*x(5)-tan(betaF)*Fs1+Fs2))/G*(a-x(5)*h)*(1+tan(betaF));
A1=1/tan(x(1));
Ef1=b/(A1*L-phif1*h);
Er1=a/((1-A1)*L+phir1*h);
Ef2=(phif2*G(1+tan(betaF))*b-(Fs1*tan(betaF)-Fs2)*L)/(L*G-G*(1+tan(betaF))*phif2*h)*phif2;
Er2=(phir2*G(1+tan(betaF))*a+(Fs1*tan(betaF)-Fs2)*L)/(L*G*tan(betaF)+G*(1+tan(betaF))*phir2*h)*phir2;
%满载时的目标函数
syms phi;
if zA<zC
if zB<=0.8
f=0.6-(int(Ef1,'phi',phiA,0.2)+int(Ef2,'phi',phiB,phiA)+int(Er2,'phi',0.8,0.2));
elseif zB>0.8
f=0.6-(int(Ef1,'phi',phiA,0.2)+int(Ef2,'phi',0.8,phiA));
end
elseif zA>=zC
if (zD<0.8)&(zB<0.8)
if zC<=0.2
f=0.6-(int(Er1,'phi',phiA,0.2)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',phiB,phiD)+int(Er2,'phi',0.8,phiB));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',phiA,phiC)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',phiB,phiD)+int(Er2,'phi',0.8,phiB));
end
elseif (zD<0.8)|(zB>=0.8)
if zC<=0.2
f=0.6-(int(Er1,'phi',phiA,0.2)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',0.8,phiD));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',phiA,phiC)+int(Er2,'phi',phiD,phiA)+int(Ef2,'phi',0.8,phiD));
end
elseif (zD>=0.8)&(zA<0.8)
if zC<=0.2
f=0.6-(int(Er1,'phi',phiA,0.2)+int(Er1,'phi',0.8,phiA));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',phiA,phiC)+int(Er2,'phi',0.8,phiA));
end
else zA>=0.8
if zC<=0.2
f=0.6-(int(Er1,'phi',0.8,0.2));
elseif zC>0.2
f=0.6-(int(Ef1,'phi',phiC,0.2)+int(Er1,'phi',0.8,phiC));
end
end
end
这个是约束条件（保存名为proportionalfun.m文件）
function [c,ceq]=proportionalcon(x,a,b,h,G,k1,k2,L)
%alpha Define the limit of variable x1.
%Ps Define the limit of variable x2.
%Ps0 Define the limit of variable x3.
%beta Define the limit of variable x4.
%z Define the limit of variable x5.
%满载时的数据
G=52528;
a=2358;
b=942;
L=3300;
h=1400;
k1=2228;
k2=1805;
Fs1=k1*x(2);
Fs2=k2*x(2);
Fp11=0.15*G*(b+0.15*h)/L;
Fp12=0.15*G*(a-0.15*h)/L;
Fp21=0.3*G*(b+0.3*h)/L;
Fp22=0.3*G*(a-0.3*h)/L;
Fp31=0.45*G*(b+0.45*h)/L;
Fp32=0.45*G*(a-0.45*h)/L;
Fp41=0.8*G*(b+0.8*h)/L;
Fp42=0.8*G*(a-0.8*h)/L;
Fp51=G*(b+0.61*h)*(0.61+0.07)/0.85*L;
Fp52=0.61*G-Fp51;
Fp61=G*(b+0.15*h)*(0.15+0.07)/0.85*L;
Fp62=0.15*G-Fp61;
gamma1=atan(Fp12/Fp11);
gamma2=atan(Fp22/Fp21);
gamma3=atan(Fp32/Fp31);
gamma4=atan(Fp42/Fp41);
gamma5=atan(Fp52/Fp51);
gamma6=atan(Fp62/Fp61);
for i=1,4,5,6
if Fpi1>Fs1
deltai=atan((Fpi2-Fs2)/(Fpi1-Fs1));
elseif Fpi1<=Fs1
deltai=atan((Fpi2-Fs2)/(Fpi1-Fs1))+pi;
end
end
A1=G*((b+0.07*h)^2-4*0.07*b*h)+4*0.85*L*h*Fs1;
A2=G*(2*(b+0.07*h)*(b+0.07*h-0.85*L)-8*0.07*b*h)+4*0.85*L*h*(Fs1-Fs2);
A3=G*((b+0.07*h-0.85*L)^2-4*0.07*b*h)-4*0.85*L*h*Fs2;
A41=G*(1+tan(delta81))/0.85*L;
A42=G*(1+tan(delta82))/0.85*L;
delta81=atan((-A2+sqrt((A2)^2-4*A1*A3))/2*A1);
delta82=atan((-A2-sqrt((A2)^2-4*A1*A3))/2*A1);
z81=-(A41*b+0.07*A41*h-G)/2*A41*h;
z82=-(A42*b+0.07*A42*h-G)/2*A42*h;
if z81>=z82
delta8=delta81;
elseif z81<z82
delta8=delta82;
end
if z81>=z82
z8=z81;
elseif z81<z82
z8=z82;
end
A71=G*((b+0.07*h)^2-4*0.07*b*h);
A72=G*(2*(b+0.07*h)*(b+0.07*h-0.85*L)-8*0.07*b*h);
A73=G*((b+0.07*h-0.85*L)^2-4*0.07*b*h);
A741=G*(1+tan(delta71))/0.85*L;
A742=G*(1+tan(delta72))/0.85*L;
delta71=atan((-A72+sqrt((A72)^2-4*A71*A73))/2*A71);
delta72=atan((-A72-sqrt((A72)^2-4*A71*A73))/2*A71);
z71=-(A741*b+0.07*A741*h-G)/2*A741*h;
z72=-(A742*b+0.07*A742*h-G)/2*A742*h;
if z71>=z72
delta7=delta71;
elseif z71<z72
delta7=delta72;
end
if z71>=z72
z7=z71;
elseif z71<z72
z7=z72;
end
if x(1)>gamma1
Fe1=Fp11;
Fe2=Fp12;
elseif (x(1)>=gamma2)&(x(1)<=gamma1)|(x(1)>=gamma4)&(x(1)<=gamma3)
ze=(a-b*tan(x(1)))/h*(1+tan(x(1)));
Fe1=ze*G*(b+ze*h)/L;
Fe2=ze*G*(a-ze*h)/L;
else (x(1)>=gamma3)&(x(1)<=gamma2)
q0=Fp32-tan(tau)*Fp31;
tan(tau)=(Fp42-Fp32)/(Fp41-Fp31);
Fe1=q0/(tan(x(1))-tan(tau));
Fe2=q0*tan(x(1))/(tan(x(1))-tan(tau));
else x(1)<=gamma4
Fe1=+\infty;
Fe2=-\infty;
end
%c(i)(i=[1:13])非线性不等式约束
c(1)=x(1)/gamma1-1.0;
if (x(5)>=0.15)&(x(5)<=0.61)
c(2)=gamma7/x(1)-1.0;
elseif (z7<0.15)||(z7>0.61)
c(2)=max(gamma5,gamma6)/x(1)-1.0;
end
c(3)=Fs2/Fe2-1.0;
c(4)=-Fs2;
c(5)=betaF/delta4-1.0;
if (z8>=0.15)&(z8<=0.61)
c(6)=delta8/betaF-1.0;
elseif (z8<0.15)|(z8>0.61)
c(6)=max(delta5,delta6)/\betaF-1.0;
end
c(7)=betaF/x(1)-1.0;
c(8)=-betaF;
c(9)=-x(5);
c(10)=x(5)-1;
然后我调用x0=[atan(k2/k1),1.3835e6,8.2e6,0.552,0.15]；
x=fmincon(@proportionalfun1,x0,[],[],[],[],[],[],@proportionalfun1)
出现这种情况
??? Error using ==> fmincon
FMINCON cannot continue because user supplied objective function failed with the following error:
Attempted to access G(1.49894); index must be a positive integer or logical.
这是怎么回事啊！
[em18][em18][em18]