我用线性规划的时候约束函数不起作用,算出来的频率不在限定的范围之内,导致求出来的刚度也不在约束的范围之内。
这是我的求解过程,但是算不出来,不知道哪里需要修改
这是求解的:
x0=[570,108,330,570,108,330,319.5,15,15];
lb=[300,50,100,300,50,200,200,10,10];
ub=[600,300,500,600,300,500,500,100,100];
option=optimset;
option.MaxFunEvals=50000000;
option.MaxIter=20000000;
[x,fval,exitflag,output]=fmincon(@opt_jieou_myfun,x0,[],[],[],[],lb,ub,@opt_jieou_mycon,option)
目标函数:
function f=opt_jieou_myfun(x)
m=323.98;%kg
Jx=25.706;Jy=13.607;Jz=23.014;Jxy=-3.029;Jyz=3.359;Jzx=-0.876;%kg*m*m
Xz=72.44;Yz=-431.92;Zz=148.82;
Xy=4.34;Yy=452.98;Zy=250.22;
Xh=199.64;Yh=-7.72;Zh=-245.93;
%mm在下面将位移坐标变成m
% Kzx;Kzy;Kzz;
% Kyx;Kyy;Kyz;
% Khx;Khy;Khz;%N/m
ALFAzu=0;BETAzu=90*pi/180;GAMAzu=90*pi/180;
ALFAzv=90*pi/180;BETAzv=0;GAMAzv=90*pi/180;
ALFAzw=90*pi/180;BETAzw=90*pi/180;GAMAzw=0;
ALFAyu=0;BETAyu=90*pi/180;GAMAyu=90*pi/180;
ALFAyv=90*pi/180;BETAyv=0;GAMAyv=90*pi/180;
ALFAyw=90*pi/180;BETAyw=90*pi/180;GAMAyw=0;
ALFAhu=8*pi/180;BETAhu=90*pi/180;GAMAhu=82*pi/180;
ALFAhv=90*pi/180;BETAhv=0;GAMAhv=90*pi/180;
ALFAhw=98*pi/180;BETAhw=90*pi/180;GAMAhw=8*pi/180;
Th=[cos(ALFAhu),cos(BETAhu),cos(GAMAhu);cos(ALFAhv),cos(BETAhv),cos(GAMAhv);cos(ALFAhw),cos(BETAhw),cos(GAMAhw)];
Tz=[cos(ALFAzu),cos(BETAzu),cos(GAMAzu);cos(ALFAzv),cos(BETAzv),cos(GAMAzv);cos(ALFAzw),cos(BETAzw),cos(GAMAzw)];
Ty=[cos(ALFAyu),cos(BETAyu),cos(GAMAyu);cos(ALFAyv),cos(BETAyv),cos(GAMAyv);cos(ALFAyw),cos(BETAyw),cos(GAMAyw)];
Fh=[1e3,0,0,0,Zh,-Yh;0,1e3,0,-Zh,0,Xh;0,0,1e3,Yh,-Xh,0]*0.001;
Fz=[1e3,0,0,0,Zz,-Yz;0,1e3,0,-Zz,0,Xz;0,0,1e3,Yz,-Xz,0]*0.001;
Fy=[1e3,0,0,0,Zy,-Yy;0,1e3,0,-Zy,0,Xy;0,0,1e3,Yy,-Xy,0]*0.001;
%%此处将位移坐标转变成m
kz=[x(1),0,0;0,x(2),0;0,0,x(3)];
ky=[x(4),0,0;0,x(5),0;0,0,x(6)];
kh=[x(7),0,0;0,x(8),0;0,0,x(9)];
M=[m,0,0,0,0,0;0,m,0,0,0,0;0,0,m,0,0,0;0,0,0,Jx,-Jxy,-Jzx;0,0,0,-Jxy,Jy,-Jyz;0,0,0,-Jzx,-Jyz,Jz];
% D=[1200,0,0,0,0,0;0,1600,0,0,0,0;0,0,2000,0,0,0;0,0,0,60,0,0;0,0,0,0,80,0;0,0,0,0,0,100];
K=Fh'*Th'*kh*Th*Fh+Fz'*Tz'*kz*Tz*Fz+Fy'*Ty'*ky*Ty*Fy;
[v,d]=eig(K,M);
g=diag(sqrt(d)/2/pi)'
for j=1:6
for k=1:6
for l=1:6
ENERGY(k,l)=M(k,l)*v(k,j)*v(l,j);
end
end
qq=sum(ENERGY);
qqt=sum(qq);
dig=[qq/qqt];
EGEN(:,j)=dig;
ep1=max(dig);
ep2=sum(dig);
eper(j)=ep1/ep2*100;
end
eper;
EGEN=EGEN*100;
f=-(1.3*max(EGEN(1)))-(3*max(EGEN(2)))-(1.5*max(EGEN(3)))-(2*max(EGEN(4)))-(3*max(EGEN(5)))-(0.5*max(EGEN(6)));
约束函数:
function [c,ceq]=opt_jieou_mycon(x)
% m=323.98;%kg
% Jx=25.706;Jy=13.607;Jz=23.014;Jxy=-3.029;Jyz=3.359;Jzx=-0.876;%kg*m*m
% Xz=72.44;Yz=-431.92;Zz=148.82;
% Xy=4.34;Yy=452.98;Zy=250.22;
% Xh=199.64;Yh=-7.72;Zh=-245.93;
% %mm在下面将位移坐标变成m
% ALFAzu=0;BETAzu=90*pi/180;GAMAzu=90*pi/180;
% ALFAzv=90*pi/180;BETAzv=0;GAMAzv=90*pi/180;
% ALFAzw=90*pi/180;BETAzw=90*pi/180;GAMAzw=0;
% ALFAyu=0;BETAyu=90*pi/180;GAMAyu=90*pi/180;
% ALFAyv=90*pi/180;BETAyv=0;GAMAyv=90*pi/180;
% ALFAyw=90*pi/180;BETAyw=90*pi/180;GAMAyw=0;
% ALFAhu=8*pi/180;BETAhu=90*pi/180;GAMAhu=82*pi/180;
% ALFAhv=90*pi/180;BETAhv=0;GAMAhv=90*pi/180;
% ALFAhw=98*pi/180;BETAhw=90*pi/180;GAMAhw=8*pi/180;
% Th=[cos(ALFAhu),cos(BETAhu),cos(GAMAhu);cos(ALFAhv),cos(BETAhv),cos(GAMAhv);cos(ALFAhw),cos(BETAhw),cos(GAMAhw)];
% Tz=[cos(ALFAzu),cos(BETAzu),cos(GAMAzu);cos(ALFAzv),cos(BETAzv),cos(GAMAzv);cos(ALFAzw),cos(BETAzw),cos(GAMAzw)];
% Ty=[cos(ALFAyu),cos(BETAyu),cos(GAMAyu);cos(ALFAyv),cos(BETAyv),cos(GAMAyv);cos(ALFAyw),cos(BETAyw),cos(GAMAyw)];
% Fh=[1e3,0,0,0,Zh,-Yh;0,1e3,0,-Zh,0,Xh;0,0,1e3,Yh,-Xh,0]*0.001;
% Fz=[1e3,0,0,0,Zz,-Yz;0,1e3,0,-Zz,0,Xz;0,0,1e3,Yz,-Xz,0]*0.001;
% Fy=[1e3,0,0,0,Zy,-Yy;0,1e3,0,-Zy,0,Xy;0,0,1e3,Yy,-Xy,0]*0.001;
% %%此处将位移坐标转变成m
% % x(1);x(2);x(3);
% % x(4);x(5);x(6);
% % x(7);x(8);x(9);%N/m
% kz=[x(1),0,0;0,x(2),0;0,0,x(3)];
% ky=[x(4),0,0;0,x(5),0;0,0,x(6)];
% kh=[x(7),0,0;0,x(8),0;0,0,x(9)];
% M=[m,0,0,0,0,0;0,m,0,0,0,0;0,0,m,0,0,0;0,0,0,Jx,-Jxy,-Jzx;0,0,0,-Jxy,Jy,-Jyz;0,0,0,-Jzx,-Jyz,Jz];
% % D=[1200,0,0,0,0,0;0,1600,0,0,0,0;0,0,2000,0,0,0;0,0,0,60,0,0;0,0,0,0,80,0;0,0,0,0,0,100];
% K=Fh'*Th'*kh*Th*Fh+Fz'*Tz'*kz*Tz*Fz+Fy'*Ty'*ky*Ty*Fy;
% [v,d]=eig(K,M);
% g=diag(sqrt(d)/2/pi)'
% for j=1:6
% for k=1:6
% for l=1:6
% ENERGY(k,l)=M(k,l)*v(k,j)*v(l,j);
% end
% end
% qq=sum(ENERGY);
% qqt=sum(qq);
% dig=[qq/qqt];
% EGEN(:,j)=dig;
% ep1=max(dig);
% ep2=sum(dig);
% eper(j)=ep1/ep2*100;
% end
% eper;
% EGEN=EGEN*100;
c(1)=g(1)-8;
c(2)=g(2)-15;
c(3)=g(3)-11;
c(4)=g(4)-18;
c(5)=g(5)-12;
c(6)=g(6)-17;
c(2)=7-g(1);
c(8)=7-g(2);
c(9)=8-g(3);
c(10)=8-g(4);
c(11)=8-g(5);
c(12)=8-g(6);
% ceq=[];
每次算的时候求出来的频率,(这里我改成了g)不在约束的范围之内,并且经常出现评价次数不够或者迭代次数不够,实在没办法。
约束函数不会写,目标函数我改了一下,设置了几个加权值。但是还是不对,算出来的g都不在范围之内。 |