基于Matlab的机械优化设计课后题

基于Matlab的机械优化设计课后题1.一维搜索法

说明：采用0.618法进行编程

代码如下：syms t

f=（需要计算的函数）

x1=0;h=2;

f1=subs(f,x1);x2=x1+h;

f2=subs(f,x2);

f3=f2-1;t=1;

if (f1-f2)>0

while f3

f3=subs(f,x2+t*h);

t=t+1;

end

x3=x2+t*h;

else

f3=f2;f2=f1;t=1;f1=f2-1;x3=x2;x2=x1;

while f1

f1=subs(f,x2-t*h);

t=t+1;

end

x1=x2-t*h;

end

a=x1;b=x3;

e=1e-05;k=0.618;

a1=b-k*(b-a);

a2=a+k*(b-a);

f1=subs(f,a1);

f2=subs(f,a2);

c=(b-a)/2;

while c>e

if f1>f2

a=a1;a1=a2;

a2=a+k*(b-a);

f1=f2;

f2=subs(f,a2);

else

b=a2;a2=a1;

a1=b-k*(b-a);

f2=f1;

f1=subs(f,a1);

end

c=(b-a)/2;

end

实际运行结果如下：

1. t=a1=5；f=11;

2. t=a1=

3.2796;f=22.6590;

3. t=a1=2;f=1.1122e-011;

2.无约束的优化设计

说明：一共采用了三种算法的编程

首先要建立步长函数，步长函数代码如下：

function az=buchang(x,d,f)

b=symvar(f);

syms aa

xn=x-d*aa;

pa=subs(f,b,xn);

pd=diff(pa,aa);

aa=solve(pd);

aa=real(aa);

aa=subs(aa);

n=size(aa,1);z=zeros(1,n);

for i=1:n

xn(:,i)=x-aa(i)*d;

end

for i=1:n

z(i)=subs(f,b,xn(:,i));

end

s=z(1);c=1;

for i=1:n-1

if s>z(i+1)

s=z(i+1);c=i+1;

end

end

az=aa(c);

end

1.DFP：·¨

syms （函数中的变量）

f=（所需要优化的函数）

x0=?????;b=symvar(f);

n=size(b,2);g=cell(n,1);G=cell(n,n);H0=eye(n,n);

for i=1:n

g{i,1}=diff(f,b(i));

end

for i=1:n %2úéúo￡è????ó

for j=1:n

G{j,i}=diff(g{i},b(j));

end

end

g0=subs(g,b,x0);

d0=H0*g0;

a=buchang(x0,d0,f);

x1=x0-a*d0;

xcha=x1-x0;

dis=mo(xcha);k=1

while dis>1e-5

g1=subs(g,b,x1);

if g1==0

Gy=subs(G,b,x1)

break

end

y0=g1-g0;

s0=x1-x0;

H1=H0+(s0*s0')/(s0'*y0)-(H0*y0*y0'*H0)/(y0'*H0*y0); d1=H1*g1;

a=buchang(x1,d1,f);

x2=x1-a*d1;

g0=g1;x0=x1;x1=x2;H0=H1;

xcha=x1-x0;dis=mo(xcha);

x1

subs(f,b,x1)

k=k+1

end

2.坐标轮换法：

syms t1t2

f=(t1^2+t2-11)^2+(t1+t2^2-7)^2;

b=symvar(f);n=size(b,2);

A=eye(n,n);x0=[1;1];

dis=1;k=0;

while dis>1e-5

if k==30

break

end

e=A(:,1);

a=buchang(x0,e,f);

x1=x0-e*a;i=2;

while i

e=A(:,i);

a=buchang(x1,e,f);

x2=x1-e*a;

x1=x2;i=i+1;

end

xcha=x1-x0;dis=mo(xcha);x0=x1;k=k+1

x0

end

3.鲍威尔法：

syms t1t2t3t4

f=(t1+10*t2)^2+5*(t3-t4)^2+(t2-2*t3)^4+10*(t1-t4)^4;

x0=[3;-1;0;1];b=symvar(f);n=size(b,2);

A=eye(n,n);dis=1;

while dis>1e-6

F0=subs(f,b,x0);

e=A(:,1);

a=buchang(x0,e,f);

x1=x0-e*a;i=2;fz=zeros(1,n);

while i

e=A(:,i);

fd=subs(f,b,x1);

fz(1,i-1)=fd;

a=buchang(x1,e,f);

x2=x1-a*e;x1=x2;

i=i+1;

end

xcha=x1-x0;dis=mo(xcha);

d=x1-x0;x=x1+d;

F2=subs(f,b,x1);F3=subs(f,b,x);

fz=[F0,fz];cz=zeros(1,n);

for i=1:n

cz(i)=fz(i)-fz(i+1);

end

czm=max(cz);m=find(czm);

P1=F3-F0;P2=(F0-2*F2+F3)*(F0-F2-czm)^2-0.5*czm*(F0-F3)^2;

if P1<0&&P2<0

a=buchang(x1,d,f);

x0=x1-a*d;

A(:,m)=[];

A=[A,d];

else

if F2

x0=x1;

else

x0=x;

end

end

x0

end

实际运行结果如下：

1.x=[5;6];f=0

2.x=[3;2];f=0

3.x=[3;0.5];f=0.52978;

4.x=[0.1239;0.2844];f=

5.8969;

5.x=[0;0;0;0];f=0;

3.约束优化方法程序

说明：本题采用matlab自带的fmincon函数来解决非线性优化问题具体代码如下：

1. function y=fun1(x)

y=(x(1)-2).^2+(x(2)-1).^2;

end

function [c,ceq]=gt1(x)

c=x(1).^2-x(2);

ceq=[];

end

A=[1,1];b=2;

x0=[1;2];lb=[];ub=[];

[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun1,x0,A,b,[], [],lb,ub,@gt1);

运行结果如下：

x=[1;1];fval=1

2. function y=fun2(x)

y=x(2).^3*((x(1)-3).^2-9)/(27*sqrt(3));

end

A=[-1/sqrt(3),1;-1,1/sqrt(3);1,1/sqrt(3)];b=[0;0;6];

x0=[2;3];lb=[0;0];ub=[];

[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun2,x0,A,b,[], [],lb,ub)

运行结果如下:

X=[4.5000;2.5981];f=-2.5313;

备注：引用函数并未能够达到最优点，原因不明；

3. function y=fun3(x)

y=1000-x(1).^2-2*x(2).^2-x(3)^2-x(1).*x(2)-x(1).*x(3);

end

function [c,cev]=gt2(x)

c=[];

cev(1)=x(1).^2+x(2).^2+x(3).^2-25;

cev(2)=8*x(1)+14*x(2)+7*x(3)-56;

end

A=[];b=[];Aq=[];bq=[];lb=[0;0;0];ub=[];

x0=[2;2;2];

[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun3,x0,A,b,Aq, bq,lb,ub,'gt2')

运行结果如下：

x= [3.5121;0.2170;3.5522];f= 961.7152;exitflag =4

备注：引用函数由于算法在迭代过程中产生了NaN，迭代被迫终止4. function y=fun4(x)

y=100*(x(2)-x(1).^2).^2+(1-x(1)).^2+90*(x(4)-x(3).^2).^2+(1-x(3))^2+1 0*((x(2)-1).^2+(x(4)-...

1).^2)+19.8*(x(2)-1)*(x(4)-1);

end

A=[];b=[];Aq=[];bq=[];lb=[-10;-10;-10;-10];ub=[10;10;10;10];x0=[-3;-1 ;-3;-1];

[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun4,x0,A,b,Aq, bq,lb,ub)

运行结果如下：

X=[1.0001;1.0002;0.9999;0.9997];f=0;exitflag=5 备注：原始对偶问题没有可行解