西安交大结构优化设计实验报告

结构优化设计实验报告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。