结构优化设计实验报告1.实验背景结构优化能在保证安全使用的前提下保证工程结构减重,提高工程的经济效益,这也是课程练习的有效补充。
2.实验课题问题1:考察最速下降法、拟牛顿法(DFP,BFGS)、单纯形法的性能,使用matlab中的fminunc 和fminsearch 函数。
●目标函数1: 目标函数,多元二次函数其中,,,,初值运行结果:最小值迭代次数最速下降法-1.864 0.455 -0.046 5.750 20牛顿DFP法-1.864 0.455 -0.046 5.750 9牛顿BFGS法-1.864 0.455 -0.046 5.750 7单纯性法-2.010 0.503 0.003 5.750 57●目标函数2运行结果:最小值迭代次数2 最速下降法-3.596 -3.800 1.4011.3牛顿DFP法-0.329 -0.380 1.869 0.446 7牛顿BFGS法-0.329 -0.380 1.869 0.446 7单纯性法-55.85 -137.8 218.2 0 43结果分析:从上述结果可以看出牛顿法具有较好的稳定性,最速下降法和单纯形法在求解超越函数时稳定性不佳,最速下降法迭代次数最少,单纯形法迭代次数最多。
问题2:使用matlab中的linprog和quadprog函数验证作业的正确性。
用单纯形法求解线性规划问题的最优解●目标函数16,运行结果:单纯形法的解析解用两相法求解线性规划问题的最优解●目标函数2,运行结果:单纯形法的解析解求解二次规划问题的最优解●目标函数2,,运行结果:问题3:用Matlab命令函数fmincon求解非线性约束规划问题●目标函数1运行结果:迭代次数:8●目标函数2运行结果:迭代次数:16问题4:用Matlab命令函数fmincon求解人字形钢管架优化问题。
已知:2F = 600kN,2B = 6 m,T=5 mm,钢管材料E = 210 GPa,密度=,许用应力[ ]=160MPa,根据工艺要求2m ≤ h≤6m ,20mm ≤ D≤300mm 。
求h , D 使总重量W为最小。
求目标函数1运行结果:迭代次数:8问题5:修改满应力程序opt4_1.m和齿形法程序opt4_2.m,自行设计一个超静定桁架结构,并对其进行优化。
要求:(1)设计变量数目不小于2;(2)给出应力的解析表达式;(3)建立以重量最小为目标函数、应力为约束的优化模型。
分别用满应立法和齿轮法求解图2超静定结构,已知材料完全相同,,σ,=σ,1500=,2000满应力法和齿轮法运行结果:附录目标函数1的代码:●function f=objfun1_2(x)●f=exp(5*(x(1))+5*(x(2)))*(3*(x(1))^2+4*(x(2))^2+3*(x(3))^2+...● 4*x(1)*x(2)+2*x(2)*x(3)+5);●●x0=[-1 1 2];●options1=optimset('HessUpdate','bfgs','LineSearchType','quadcubic');●[x,fval,exitflag,output]=fminunc('objfun1_1',x0 ,options1)●●options2=optimset('HessUpdate','dfp','LineSearchType','quadcubic');●[x,fval,exitflag,output]=fminunc('objfun1_1', x0,options2)●●options3=optimset('HessUpdate','steepdesc');●[x,fval,exitflag,output]=fminunc('objfun1_1', x0,options3)●●[x,fval,exitflag,output]=fminsearch('objfun1_1', x0)目标函数2的代码●function f=objfun1_2(x)●f=exp(5*(x(1))+5*(x(2)))*(3*(x(1))^2+4*(x(2))^2+3*(x(3))^2+...● 4*x(1)*x(2)+2*x(2)*x(3)+5);●●x0=[1 1 2];●options3=optimset('HessUpdate','steepdesc');●[x,fval,exitflag,output]=fminunc('objfun1_2', x0,options3)●●options2=optimset('HessUpdate','dfp','LineSearchType','quadcubic');●[x,fval,exitflag,output]=fminunc('objfun1_2', x0,options2)●●options1=optimset('HessUpdate','bfgs','LineSearchType','quadcubic');●[x,fval,exitflag,output]=fminunc('objfun1_2',x0 ,options1)●[x,fval,exitflag,output]=fminsearch('objfun1_2', x0)问题2目标函数1的代码●x0=[0 0]';● c=[-2 -6]';● A=[-1 1; 2 1];● b=[1 2]';● lb=[0 0]';● options = optimset('LargeScale', 'off', 'Simplex', 'on')● [x,fval,exitflag,output] =linprog(c,A,b,[],[],lb ,[],x0,options)目标函数2的代码●x0=[0 0]';●c=[3 4]';●A=[-5 -4; 3 5;5 4];●b=[-100 150 200]';●lb=[0 0]';●options = optimset('LargeScale', 'off', 'Simplex', 'on')●[x,fval,exitflag,output] =linprog(c,A,b,[],[],lb ,[],x0,options)目标函数3的代码●H=[6 2; 2 2];●c= [-30;-14];●A=[1 1;2 -1; 0 1];●b=[3; 4; 2];●lb=[0; 0];●%有效集法(active-set)●[x,fval,exitflag,output ,lambda] =quadprog (H,c,A,b,[],[],lb)●问题3目标函数1●●function opt3_1●x0=[-1 1]●options=optimset('LargeScale', 'off');●[x,fval,exitflag,output]=fmincon('objfun3_1',x0,[],[],[],[],[],[],'confun3_1',options)●%目标函数●function f=objfun3_1(x)●f=(x(1))^4+(x(2))^4-14*(x(1))^2-38*(x(2))^2-24*x(1)+120*x(2);●%约束函数●function [c,ceq]=confun3_1(x)●%非线性不等式约束函数●c=[ x(1)^2+x(2)^2-18 ; x(1)-x(2)-5 ];●%非线性等式约束函数●ceq=[];目标函数2●function opt3_2●clc●x0=[1 1 1 1 1 1 1];●lb=[0;0;0;0;0;0;0];●options=optimset('LargeScale', 'off');●[x,fval,exitflag,output]=fmincon('objfun3_2',x0,[],[],[],[],lb,[],'confun3_2',options)●%目标函数●function f=objfun3_2(x)●e=2.7183;●f=-5*x(1)-5*x(2)-4*x(3)-6*x(4)-x(1)*x(3)-5*x(5)/(1+x(5))-8*x(6)/(1+x(6))...● -10*(1-2*e^(-x(7)))+2*e^(-x(7));●%约束函数●function [c,ceq]=confun3_2(x)●%非线性不等式约束函数●c=[ x(1)+x(2)+x(3)+x(4)-5 ;●x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)-10;● x(1)+x(3)+x(5)+(5)^2-(x(7))^2-5];●非线性等式约束函数●ceq=[2*x(4)+x(5)+0.8*x(6)+x(7)-5;(x(2))^2+(x(3))^2+(x(5))^2+(x(6))^2-5];问题4●function opt4_1●clc●x0=[4 0.2];●lb=[2;0.02];●ub=[6;0.3];●options=optimset('LargeScale', 'off');●[x,fval,exitflag,output]=fmincon('objfun4_1',x0,[],[],[],[],lb,ub,'confun4_1',options)●●%目标函数●function f=objfun4_1(x)●f=245*x(2)*sqrt(9+(x(1))^2);●●%约束函数●function [c,ceq]=confun4_1(x)●%非线性不等式约束●c=[300*sqrt(9+(x(1))^2)/(15.71*x(1)*x(2))-160;● 300*sqrt(9+(x(1))^2)/(15.71*x(1)*x(2))-...● (259077.12*((x(2))^2+0.000025))/(9+(x(1))^2)];●%线性等式约束●ceq=[];问题5●function [sr]=stressratio(x,pb)●% 三杆超静定问题——应力计算●% 计算应力比●switch pb●case 1●% 两个工况● stress=[3000/x(1);● 3000/x(3);●2000*sqrt(2)*x(1)*x(3)/(x(1)*((x(1)*x(2)+x(2)*x(3)+sqrt(2)*x(3)*x(1))));●2000*(x(1)+x(3))/(sqrt(2)*(x(1)*x(3)+x(2)*x(3)+x(1)*x(2)));●2000*sqrt(2)*x(1)*x(3)/(x(3)*((x(1)*x(2)+x(2)*x(3)+sqrt(2)*x(3)*x(1))));]● astress=[ 2000;1500;2000;2000;2000];● sr=stress./astress;●End●function [w]=weight(x)●%目标函数●w=10*(sqrt(2)*x(1)+x(2)+sqrt(2)*x(3));●function [ksi]=getksi(sr,pb)●% 三杆超静定问题●%计算迭代算子●switch pb●case 1●% 两个工况● ksi=[max([sr(1) sr(3)]);sr(4);max([sr(2) sr(5)])];●end●function opt4_1●% 满应力法●clc●beta=1.0;●x0=[1.0;1.0;1.0];●miter=300;●pb=1●iter=0●while iter<miter● [sr]=stressratio(x0,pb)● [ksi]=getksi(sr,pb)● w=weight(x0)● x1=x0.*ksi.^beta● xn=norm(x1-x0);●if iter>1 & xn<1.0e-5●break;●end● x0=x1;● iter=iter+1●end●sr=stressratio(x1,pb);●%齿轮形法●function opt4_2●clc●beta=1.0; %松弛因子●alpha=0.8; %满应力步长因子●x0=[1.0;1.0;1.0]; %初值●w0=1.0e10;●miter=100; %对大迭代次数●pb=1 %问题1●iter=0●while iter<miter● [sr]=stressratio(x0,pb)● [ksi]=getksi(sr,pb)● eta=max(ksi)●%射步法● x0=x0*eta;● sr=sr/eta;● [ksi1]=getksi(sr,pb)● w=weight(x0)● wtol=abs(w-w0);●if iter>1 & (w>w0 | wtol<=1.0e-5)●break;●end●%忙应力步● x1=alpha*x0.*ksi1.^beta+(1.0-alpha)*x0● x0=x1;● w0=w;● iter=iter+1●end●sr=stressratio(x1,pb);11。