maoving
2012-07-04, 12:19
哪位大神帮忙看看这段课程设计的优化程序,运行不出来,貌似就一个小问题。。。急啊!:)
%主程序
clear all
clc
%膜片弹簧结构参数值
E=210000;%材料弹性模量
miu=0.3;%泊松比
ds=2;%磨损极限
dt=3;%推力行程
D=300;d=175;
r0=39;rf=40.3;%结构参数膜片弹簧小端内半径和分离作用半径
x0=[5.8 2.93 145.7 116.8 143.66 116.1 4.8];%分别是H h R r R1 r1 lamda
x7=0:0.1:9;
aa=pi*E.*x0(2).*x7/(6*(1-miu.^2));
bb=log(x0(3)./x0(4))./((x0(5)-x0(6)).^2);
cc=x0(1)-x7.*(x0(3)-x0(4))./(x0(5)-x0(6));
dd=x0(1)-0.5.*x7.*(x0(3)-x0(4))./(x0(5)-x0(6));
ee=x0(2).^2;
F=aa.*bb.*(cc.*dd+ee);
plot(x7,F,'b')%绘制原始膜片弹簧弹性特性曲线
hold on
%设计变量的上下限
Lb=[4 2 140 115 135 115 4];%设计变量下限
Ub=[7 4 150 125 145 125 6];%设计变量上限
%线性不等式约束系数矩阵和常数
% H h R r R1 r1 lamda
A=[1,-2.2,0,0,0,0,0;
-1,1.7,0,0,0,0,0;
1,0,-pi/15,pi/15,0,0,0,;
-1,0,pi/20,-pi/20,0,0,0;
0,0,1,-1.35,0,0,0;
0,0,-1,1.2,0,0,0;
0,-50,1,0,0,0,0;
0,35,-1,0,0,0,0;
0,0,0,0,1,0,0;
0,0,0,0,-1,0,0;
0,0,1,0,-1,0,0;
0,0,-1,0,1,0,0;
0,0,0,-1,0,1,0;
0,0,0,1,0,-1,0;];
b=[0 0 0 0 0 0 0 0 D./2-(D+d)./4 7 -1 6 0];
%线性等式约束
Aeq=[];
beq=[];
options=optimset('largescale','off','display','iter');
[x,fval,exitflag,out]=fmincon(@ objfun,x0,A,b,Aeq,beq,Lb,Ub,@ confun,options);
[c]=confun(x);
x7=0:0.1:9;
aa=pi*E.*x(2).*x7/(6*(1-miu.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc=x(1)-x7.*(x(3)-x(4))./(x(5)-x(6));
dd=x(1)-0.5.*x7.*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
F=aa.*bb.*(cc.*dd+ee);
plot(x7,F,'r')%绘制优化后膜片弹簧弹性特性曲线
%%目标函数的确定
FUNCTION f=obufun(x)
E=210000;%材料弹簧模量
miu=0.3;%泊松比
rf=40.3;%分离轴承推力半径
%%将膜片弹簧特性公式分成aa、bb、cc、dd、ee五部分表示
aa=pi*E.*x(2).*x(7)/(6*(1-miu.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc=x(1)-x(7).*(x(3)-x(4))./(x(5)-x(6));
dd=x(1)-0.5.*x(7).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
ff=x((5)-x(6))./(x(6)-rf);%ff=(R1-r1)/(r1-rf)
%%%%磨损后的公式参数变化
ds=2;%磨损极限在1.6-2.2之间,取2mm
aa1=pi*E.*x(2).*(x(7)-ds)/(6*(1-miu.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc1=x(1)-(x(7)-ds).*(x(3)-x(4))./(x(5)-x(6));
dd1=x(1)-0.5.*(x(7)-ds).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
%%%%分离过程公式参数变化
dt=3;%分离行程取值为三毫米
aa2=pi*E.*x(2).*(x(7)+dt)/(6*(1-miu.^2));
bb2=log(x(3)./x(4))./((x(5)-x(6)).*(x(5)-rf));
cc2=x(1)-(x(7)+dt).*(x(3)-x(4))./(x(5)-x(6));
dd2=x(1)-0.5.*(x(7)+dt).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
%%%%双目标函数表达式
f1=abs(aa.*bb.*(cc.*dd+ee)-aa1.*bb.*(cc1.*dd1+ee));%第一个目标函数:磨损极限内正压力的变化值
f2=aa2.*bb2.*(cc2.*dd2+ee);%第二个目标函数:膜片弹簧在分离位置时的弹力
fac=0.6;%加权因子
f=fac.*f1+(1-fac)*f2;%总体目标函数
%%%建立非线性约束条件
FUNCTION [c,ceq]=confun(x)
aa=pi*210000.*x(2).*x(7)/(6*(1-0.3.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc=x(1)-x(7).*(x(3)-x(4))./(x(5)-x(6));
dd=x(1)-0.5.*x(7).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
miu=3;
kk=210000/((1-miu.^2).*x(4));
e=(x(3)-x(4))./(log(x(3)./x(4)));%中性点半径
tt=0.5.*(e-x(4));
alfa=atan(x(1)./(x(3)-x(4)));%膜片弹簧锥角的计算
fa=alfa+0.5.*x(2)./(e-x(4));%切向压应力达到最大值时的膜片转角
thegatb=abs(kk.*(tt.*fa.^2-(2.*tt.*alfa+x(2)./2).*fa));%膜片弹簧危险部位切向压应力计算
%%%%%%%%%%%%%%%%%%
rf=40.3;
ff=(x(5)-x(6))./(x(6)-rf);%ff=(R1-r1)/(r1-rf) F2/F1力的比值
fff=aa.*bb.*(cc.*dd+ee)*ff;%F2
%%%%%%%%%%%%%%%%%%%
n=18;%分离指个数
b=10;%分离指根部的宽度
thegarb=abs(6*(x(4)-rf)*fff./(x(2).^2.*n.*b));%膜片弹簧危险部位弯曲应力计算
%%%%%%%%%%%%%%%%%%%%%%%%%
T=700*1000;%离合器所要求传递最大扭矩,单位转化为N*mm
a=300/2;b=175/2;z=4;fz=0.3;
c(1)=sqrt(thegarb.^2+thegatb.^2)-1500;%膜片弹簧危险点最大当量应力约束,非线性不等式1
c(2)=T/(z.*fz.*(2./3).*(a.^3-b.^3)/(a.^2-b.^2))-(aa.*bb.*(cc.*dd+ee));%膜片弹簧产生压紧力的约束,非线性不等式2
ceq=[];
%主程序
clear all
clc
%膜片弹簧结构参数值
E=210000;%材料弹性模量
miu=0.3;%泊松比
ds=2;%磨损极限
dt=3;%推力行程
D=300;d=175;
r0=39;rf=40.3;%结构参数膜片弹簧小端内半径和分离作用半径
x0=[5.8 2.93 145.7 116.8 143.66 116.1 4.8];%分别是H h R r R1 r1 lamda
x7=0:0.1:9;
aa=pi*E.*x0(2).*x7/(6*(1-miu.^2));
bb=log(x0(3)./x0(4))./((x0(5)-x0(6)).^2);
cc=x0(1)-x7.*(x0(3)-x0(4))./(x0(5)-x0(6));
dd=x0(1)-0.5.*x7.*(x0(3)-x0(4))./(x0(5)-x0(6));
ee=x0(2).^2;
F=aa.*bb.*(cc.*dd+ee);
plot(x7,F,'b')%绘制原始膜片弹簧弹性特性曲线
hold on
%设计变量的上下限
Lb=[4 2 140 115 135 115 4];%设计变量下限
Ub=[7 4 150 125 145 125 6];%设计变量上限
%线性不等式约束系数矩阵和常数
% H h R r R1 r1 lamda
A=[1,-2.2,0,0,0,0,0;
-1,1.7,0,0,0,0,0;
1,0,-pi/15,pi/15,0,0,0,;
-1,0,pi/20,-pi/20,0,0,0;
0,0,1,-1.35,0,0,0;
0,0,-1,1.2,0,0,0;
0,-50,1,0,0,0,0;
0,35,-1,0,0,0,0;
0,0,0,0,1,0,0;
0,0,0,0,-1,0,0;
0,0,1,0,-1,0,0;
0,0,-1,0,1,0,0;
0,0,0,-1,0,1,0;
0,0,0,1,0,-1,0;];
b=[0 0 0 0 0 0 0 0 D./2-(D+d)./4 7 -1 6 0];
%线性等式约束
Aeq=[];
beq=[];
options=optimset('largescale','off','display','iter');
[x,fval,exitflag,out]=fmincon(@ objfun,x0,A,b,Aeq,beq,Lb,Ub,@ confun,options);
[c]=confun(x);
x7=0:0.1:9;
aa=pi*E.*x(2).*x7/(6*(1-miu.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc=x(1)-x7.*(x(3)-x(4))./(x(5)-x(6));
dd=x(1)-0.5.*x7.*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
F=aa.*bb.*(cc.*dd+ee);
plot(x7,F,'r')%绘制优化后膜片弹簧弹性特性曲线
%%目标函数的确定
FUNCTION f=obufun(x)
E=210000;%材料弹簧模量
miu=0.3;%泊松比
rf=40.3;%分离轴承推力半径
%%将膜片弹簧特性公式分成aa、bb、cc、dd、ee五部分表示
aa=pi*E.*x(2).*x(7)/(6*(1-miu.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc=x(1)-x(7).*(x(3)-x(4))./(x(5)-x(6));
dd=x(1)-0.5.*x(7).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
ff=x((5)-x(6))./(x(6)-rf);%ff=(R1-r1)/(r1-rf)
%%%%磨损后的公式参数变化
ds=2;%磨损极限在1.6-2.2之间,取2mm
aa1=pi*E.*x(2).*(x(7)-ds)/(6*(1-miu.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc1=x(1)-(x(7)-ds).*(x(3)-x(4))./(x(5)-x(6));
dd1=x(1)-0.5.*(x(7)-ds).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
%%%%分离过程公式参数变化
dt=3;%分离行程取值为三毫米
aa2=pi*E.*x(2).*(x(7)+dt)/(6*(1-miu.^2));
bb2=log(x(3)./x(4))./((x(5)-x(6)).*(x(5)-rf));
cc2=x(1)-(x(7)+dt).*(x(3)-x(4))./(x(5)-x(6));
dd2=x(1)-0.5.*(x(7)+dt).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
%%%%双目标函数表达式
f1=abs(aa.*bb.*(cc.*dd+ee)-aa1.*bb.*(cc1.*dd1+ee));%第一个目标函数:磨损极限内正压力的变化值
f2=aa2.*bb2.*(cc2.*dd2+ee);%第二个目标函数:膜片弹簧在分离位置时的弹力
fac=0.6;%加权因子
f=fac.*f1+(1-fac)*f2;%总体目标函数
%%%建立非线性约束条件
FUNCTION [c,ceq]=confun(x)
aa=pi*210000.*x(2).*x(7)/(6*(1-0.3.^2));
bb=log(x(3)./x(4))./((x(5)-x(6)).^2);
cc=x(1)-x(7).*(x(3)-x(4))./(x(5)-x(6));
dd=x(1)-0.5.*x(7).*(x(3)-x(4))./(x(5)-x(6));
ee=x(2).^2;
miu=3;
kk=210000/((1-miu.^2).*x(4));
e=(x(3)-x(4))./(log(x(3)./x(4)));%中性点半径
tt=0.5.*(e-x(4));
alfa=atan(x(1)./(x(3)-x(4)));%膜片弹簧锥角的计算
fa=alfa+0.5.*x(2)./(e-x(4));%切向压应力达到最大值时的膜片转角
thegatb=abs(kk.*(tt.*fa.^2-(2.*tt.*alfa+x(2)./2).*fa));%膜片弹簧危险部位切向压应力计算
%%%%%%%%%%%%%%%%%%
rf=40.3;
ff=(x(5)-x(6))./(x(6)-rf);%ff=(R1-r1)/(r1-rf) F2/F1力的比值
fff=aa.*bb.*(cc.*dd+ee)*ff;%F2
%%%%%%%%%%%%%%%%%%%
n=18;%分离指个数
b=10;%分离指根部的宽度
thegarb=abs(6*(x(4)-rf)*fff./(x(2).^2.*n.*b));%膜片弹簧危险部位弯曲应力计算
%%%%%%%%%%%%%%%%%%%%%%%%%
T=700*1000;%离合器所要求传递最大扭矩,单位转化为N*mm
a=300/2;b=175/2;z=4;fz=0.3;
c(1)=sqrt(thegarb.^2+thegatb.^2)-1500;%膜片弹簧危险点最大当量应力约束,非线性不等式1
c(2)=T/(z.*fz.*(2./3).*(a.^3-b.^3)/(a.^2-b.^2))-(aa.*bb.*(cc.*dd+ee));%膜片弹簧产生压紧力的约束,非线性不等式2
ceq=[];