工程结构优化设计编程作业.docx
《工程结构优化设计编程作业.docx》由会员分享，可在线阅读，更多相关《工程结构优化设计编程作业.docx（14页珍藏版）》请在冰豆网上搜索。
工程结构优化设计编程作业
工程结构优化设计—编程作业
1.用黄金分割法求方程
在区间[-1,1]上的解。
在MATLAB中没有专门的函数实现黄金分割法求解线性方程，可通过编写glodf.m函数实现黄金分割法求解，其代码如下:
functionx=glodf（f,a,b,eps）
if（nargin==3）
eps=1.0e-4;
end
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
x=a;
end
if（f2==0）
x=b;
end
if（f1*f2>0）
disp（’两端点函数值乘积大于0!
’）;
return;
else
t1=a+（b-a）*0.382;
t2=a+（b-a）*0.618;
f_1=subs（sym（f）,findsym（sym（f））,t1）;
f_2=subs（sym（f）,findsym（sym（f））,t2）;
tol=abs（t1-t2）;
while（tol>eps）
if（f_1*f_2<0）
a=t1;
b=t2;
else
fa=subs（sym（f）,findsym（sym（f））,a）;
if（f_1*fa>0）
a=t2;
else
b=t1;
end
end
t1=a+（b-a）*0.382;
t2=a+（b-a）*0.618;
f_1=subs（sym（f）,findsym（sym（f））,t1）;
f_2=subs（sym（f）,findsym（sym（f））,t2）;
tol=abs（t2-t1）;
end
x=（t1+t2）/2;
end
其实现的MATLAB代码如下:
>>clearall;
>>tic
>>x=glodf（'x^3-3*x+1',-1,1,1e-6）
x=
0.3473
2.用抛物线法求解
的极值点，取值范围[0.5，1]，允许误差为1e-6。
在MATLAB中没有专门的函数实现抛物线法求解线性方程，可通过编写paowuxian.m函数实现抛物线法求解，其代码如下:
functionx=paowuxian（f,x0,x1,x2,e）
ifnargin<4,e=1e-4;
end
x=x2;y=x1;z=x0;
whileabs（x-y）>e
h1=y-z;
h2=x-y;
c1=（feval（f,y）-feval（f,z））/h1;
c2=（feval（f,x）-feval（f,y））/h2;
d=（c1-c2）/（h2+h1）;
w=c2+h2*d;
xi=x-（2*feval（f,x））/（w+（w/abs（w））*sqrt（w^2-4*feval（f,x）*d））;
z=y;
y=x;
x=xi;
end
其实现的MATLAB代码如下:
>>fun=inline（'x^3+2*x-0.2'）;
>>paowuxian（fun,0,0.5,1,1e-6）,formatshort
ans=
0.0995
3.用齿形法对下图的结构进行杆件截面优化设计：
节点处作用两种工况：
（1）P1=2000kN，P2=0kN；
（2）P1=0kN，P2=2000kN。
各杆采用同一材料制成，弹性模量E为常数，容许应力分别为：
[σ+]=2╳107kPa，[σ-]=1.5╳107kPa。
此结构为对称结构，工况也对称，所以A1=A3，计算时只取两个变量A1和A2，一种工况P1=2000kN，P2=0kN。
目标函数取为结构总重量：
在MATLAB中没有专门的函数实现齿形法求解线性方程，可通过编写Zigzag.m函数实现齿形法求解，其代码如下:
function[newA1,newA2,newA3,W,exitflag]=Zigzag（A1,A2,maxstep）
ifnargin<3
maxstep=96;
end
A=[A1;A2];
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
w=sqrt
（2）*2*A1+A2;
w1=w;
w2=w;
newA1=A1;
newA2=A2;
newA3=A1;
symspl;
W=vpa（w*p*l）;
A11=A1;
A13=A1;
A22=A2;
A24=A2;
exitflag=-1;
fork=1:
maxstep
ifw2>w1
newA1=A11;
newA2=A22;
newA3=A11;
W=vpa（w1*p*l,4）;
exitflag=1;
return;
end
u1=max（U）;
A=A.*u1;
A11=A
（1）;
A22=A
（2）;
U=U./u1;
U=roundn（U,-4）;
w=sqrt
（2）*2*A11+A22;
w1=w;
ifw1>w2
newA1=A13;
newA2=A24;
newA3=A13;
W=vpa（w2*p*l,4）;
exitflag=0;
return;
end
A=A.*U;
A1=A
（1）;
A2=A
（2）;
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
u1=max（U）;
A=A.*u1;
A13=A
（1）;
A24=A
（2）;
U=U./u1;
U=roundn（U,-4）;
w=sqrt
（2）*2*A13+A24;
w2=w;
A=A.*U;
A1=A
（1）;
A2=A
（2）;
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
end
end
其实现的MATLAB代码如下:
>>clear;
A1=1e-4;
A2=1e-4;
[Anew1,Anew2,Anew3,W,exitflag]=Zigzag（A1,A2）
Anew1=
7.7346e-05
Anew2=
4.5309e-05
Anew3=
7.7346e-05
W=
0.0002641*l*p
exitflag=
0
4.用修改齿形法对结构进行杆件截面优化设计。
在MATLAB中没有专门的函数实现修改齿形法求解线性方程，可通过编写IZigzag.m函数实现修改齿形法求解，其代码如下:
function[newA1,newA2,newA3,W,exitflag]=IZigzag（A1,A2,v,maxstep）
ifnargin<5
maxstep=96;
end
A=[A1;A2];
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
w=sqrt
（2）*2*A1+A2;
w1=w;
w2=w;
newA1=A1;
newA2=A2;
newA3=A1;
symspl;
W=vpa（w*p*l）;
A11=A1;
A13=A1;
A22=A2;
A24=A2;
exitflag=-1;
fork=1:
maxstep
ifw2>w1
newA1=A11;
newA2=A22;
newA3=A11;
W=vpa（w1*p*l,4）;
exitflag=1;
return;
end
u1=max（U）;
A=A.*u1;
A11=A
（1）;
A22=A
（2）;
U=U./u1;
x1=U
（1）;
x2=U
（2）;
y1=1-v*（1-x1）;
y2=1-v*（1-x2）;
U=[y1;y2];
U=roundn（U,-4）;
w=sqrt
（2）*2*A11+A22;
w1=w;
ifw1>w2
newA1=A13;
newA2=A24;
newA3=A13;
W=vpa（w2*p*l,4）;
exitflag=0;
return;
end
A=A.*U;
A1=A
（1）;
A2=A
（2）;
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
u1=max（U）;
A=A.*u1;
A13=A
（1）;
A24=A
（2）;
U=U./u1;
x1=U
（1）;
x2=U
（2）;
y1=1-v*（1-x1）;
y2=1-v*（1-x2）;
U=[y1;y2];
U=roundn（U,-4）;
w=sqrt
（2）*2*A13+A24;
w2=w;
A=A.*U;
A1=A
（1）;
A2=A
（2）;
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
end
end
其实现的MATLAB代码如下:
>>clear;
A1=1e-4;
A2=1e-4;
v=0.5;
[Anew1,Anew2,Anew3,W,exitflag]=IZigzag（A1,A2,v）
Anew1=
7.8168e-05
Anew2=
4.2839e-05
Anew3=
7.8168e-05
W=
0.0002639*l*p
exitflag=
0
5.用满应力法对结构进行杆件截面优化设计
在MATLAB中没有专门的函数实现满应力法求解线性方程，可通过编写Full_Stressed.m函数实现满应力法求解，其代码如下:
function[newA1,newA2,newA3,exitflag]=Full_Stressed（A1,A2,wucha,maxstep）
ifnargin<4
maxstep=96;
end
ifnargin<3
wucha=1e-6;
end
A=[A1;A2];
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
newA1=A1;
newA2=A2;
newA3=A1;
exitflag=-1;
fork=1:
maxstep
m=U
（1）;
h=1-m;
n=U
（2）;
j=1-n;
ifabs（h）newA1=A1;
newA2=A2;
newA3=A1;
exitflag=0;
return;
end
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
U=roundn（U,-4）;
A=A.*U;
A1=A
（1）;
A2=A
（2）;
A=[A1;A2];
end
newA1=A1;
newA2=A2;
newA3=A1;
end
其实现的MATLAB代码如下:
>>clearall;closeall;
A1=1e-4;
A2=1e-4;
[newA1,newA2,newA3,exitflag]=Full_Stressed（A1,A2）
newA1=
9.8988e-05
newA2=
1.4373e-06
newA3=
9.8988e-05
exitflag=
-1
6.用加松弛因子满应力法对结构进行杆件截面优化设计
在MATLAB中没有专门的函数实现加松弛因子满应力法求解线性方程，可通过编写SRF_Full_Stressed.m函数实现满应力法求解，其代码如下:
function[newA1,newA2,newA3,exitflag]=SRF_Full_Stressed（A1,A2,B,wucha,maxstep）
ifnargin<5
maxstep=96;
end
ifnargin<3
wucha=1e-6;
end
A=[A1;A2];
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
newA1=A1;
newA2=A2;
newA3=A1;
exitflag=-1;
fork=1:
maxstep
m=U
（1）;
h=1-m;
n=U
（2）;
j=1-n;
ifabs（h）newA1=A1;
newA2=A2;
newA3=A1;
exitflag=0;
return;
end
U=[（1e-4）*（sqrt
（2）*A1+A2）/（sqrt
（2）*A1^2+2*A1*A2）;（1e-4）*（sqrt
（2）*A1）/（sqrt
（2）*A1^2+2*A1*A2）];
U=roundn（U,-4）;
U=U.^B;
A=A.*U;
A1=A
（1）;
A2=A
（2）;
A=[A1;A2];
end
newA1=A1;
newA2=A2;
newA3=A1;
end
其实现的MATLAB代码如下:
>>clearall;closeall;
A1=1e-4;
A2=1e-4;
B=1.1;
wucha=0.05;
[newA1,newA2,newA3,exitflag]=SRF_Full_Stressed（A1,A2,B,wucha）
newA1=
9.5527e-05
newA2=
6.3262e-06
newA3=
9.5527e-05
exitflag=
0
升级会员 