机械优化设计MATLAB程序.doc

完整word版)机械优化设计MATLAB程序

机械优化设计作业
1. 用二次插值法求函数极小值，精度e=0.01。
在MATLAB的M文件编辑器中编写的M文件，如下：
f=inline('(t+1)*(t-2)^2','t')
a=0;b=3;epsilon=0.01;
t1=a;f1=f(t1);
t3=b;f3=f(t3);
t2=0.5*(t1+t3);f2=f(t2);
c1=(f3-f1)/(t3-t1);
c2=((f2-f1)/(t2-t1)-c1)/(t2-t3);
t4=0.5*(t1+t3-c1/c2);f4=f(t4);
k=0;
while(abs(t4-t2)>=epsilon)
if t2f4
f1=f2;t1=t2;
t2=t4;f2=f4;
else
f3=f4;t3=t4;
end
else
if f2>f4
f3=f2;t3=t2;
t2=t4;f2=f4;
else
f1=f4;t2=t4;
end
end
c1=(f3-f1)/(t3-t1);
c2=((f2-f1)/(t2-t1)-c1)/(t2-t3);
t4=0.5*(t1+t3-c1/c2);f4=f(t4);
k=k+1;
end
%输出最优解
if f2>f4
t=t4;f=f(t4);
else
t=t2;f=f(t2);
end
fprintf(1,'迭代计算k=%3.0f\n',k)
fprintf(1,'极小点坐标t=%3.0f\n',t)
fprintf(1,'函数值f=%3.4f\n',f)
运行结果如下：
迭代计算k= 7
极小点坐标t= 2
函数值f=0.0001
2. 用黄金分割法求函数的极小值，精度e=0.01。
在MATLAB的M文件编辑器中编写的M文件，如下：
f=inline('t^(2/3)-(t^2+1)^(1/3)','t');
a=0;b=3;epsilon=0.01;
t1=b-0.618*(b-a);f1=f(t1);
t2=a+0.618*(b-a);f2=f(t2);
k=1;
while abs(b-a)>=epsilon
if f1t1);
else
a=t1;t1=t2;f1=f2;
t2=a+0.618*(b-a);f2=f(t2);
end
t=0.5*(b+a);
k=k+1;
f0=f(t);
end
fprintf(1,'迭代次数k=% 3.0f\n',k)
fprintf(1,'迭代区间-左端a=%3.4f\n',a)
fprintf(1,'试点1坐标值t1=%3.4f\n',t1)
fprintf(1,'函数值f1=%3.4f\n',f(t1))
fprintf(1,'迭代区间-右端b=%3.4f\n',b)
fprintf(1,'试点2坐标值t2=%3.4f\n',t2)
fprintf(1,'函数值f2=%3.4f\n',f(t2))
fprintf(1,'区间中点t=%3.4f\n',t)
fprintf(1,'函数值f0=%3.4f\n',f(t))
运行结果如下：
迭代次数k= 13
迭代区间-左端a=0.0000
试点1坐标值t1=0.0036
函数值f1=-0.9767
迭代区间-右端b=0.0093
试点2坐标值t2=0.0058
函数值f2=-0.9679
区间中点t=0.0047
函数值f0=-0.9721
由黄金分割法在初始区间[0，3]求得的极小值点为t=0.0047，极小值为-0.9721。
3. 用牛顿法、阻尼牛顿法及变尺度法求函数的极小点。
(1)在用牛顿法在MATLAB的M文件编辑器中编写的M文件，如下：
function [x,fx,k]=niudunfa(x0)
syms x1 x2
f=(x1-2)^4+(x1-2*x2)^2;
fx=0;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
G=jacobian(df,v);
epson=1e-12;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1x2},{x0(1,1),x0(2,1)});
k=0;
p=-G1\g1;
x0=x0+p;
while(norm(g1)>epson)
p=-G1\g1;
x0=x0+p;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
k=k+1;
end
x=x0;
fx=subs(f,{x1,x2},{x(1,1),x(2,1)});
运行结果如下：
>> [x,fx,k]=niudunfa([1;1])

x =1.9999554476059523381489991377897
0.99997772380297616907449956889483
fx =0.0000000000000000039398907941382470301534502947647
k =23
(2)用阻尼牛顿法在MATLAB的M文件编辑器中编写的M文件，如下：
function [x,fx,k]=zuniniudunfa(x0)%阻尼牛顿法
syms x1 x2
f=(x1-2)^4+(x1-2*x2)^2;
fx=0;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
G=jacobian(df,v);
epson=1e-12;%停机原则
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
k=0;%迭代次数
p=-G1\g1;
a0=-p'*g1/(p'*G1*p);
x0=x0+a0*p;
while(norm(a0*p)>epson)
p=-G1\g1;
a0=-p'*g1/(p'*G1*p);
x0=x0+a0*p;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
k=k+1;
end
x=x0;
fx=subs(f,{x1,x2},{x0(1,1),x0(2,1)});

运行结果如下：
>>[x,fx,k]=zuniniudunfa([1;1])

x=1.9999554476059523381489991377897
0.99997772380297616907449956889483
fx=0.0000000000000000039398907941382470301534502947647
k=23
(3)用变尺度法在MATLAB的M文件编辑器中编写的M文件，如下：
4. 用共轭梯度法求函数的极小点
（1）用共轭梯度法在MATLAB的M文件编辑器中编写的M文件，如下：

function[y,x,k]=CG(A,b,c,x0)
%共轭梯度法解minf（x）=0.5*X'*A*X+b'x+c
eps=1e-6;%迭代停机原则
%fx=0.5*x0'.*A.*x0+b'.*x0+c;
r0=A*x0+b;
if norm(r0)<=eps
x=x0;
y=0.5*x'*A*x+b'*x+c;
k=0;
end
p0=-r0;
a=-r0'*p0/(p0'*A*p0);
x1=x0+a*p0;
r1=A*x1+b;
k=0;
while norm(r1)>eps
beta=(r1'*r1)/(r0'*r0);
p1=-r1+beta*p0;
alpha=-(r1'*p1)/(p1'*A*p1);
x1=x1+alpha*p1;
r2=A*x1+b;
p0=p1;
r0=r1;
r1=r2;
k=k+1;
end
x=x1;
y=0.5*x'*A*x+b'*x+c;
运行结果如下：
[y,x,k]=CG([3 -1;-1 1],[-2;0],0,[21])
y = -1
x = 1.0000
1.0000
k = 1
(2) 用变尺度法在MATLAB的M文件编辑器中编写的M文件，如下：
function [x,fx,k]=bianchidufa(A,b,c,x0)
%用变尺度法求fx=0.5*x'*A*x+b'*x+c;
epson=1e-12;
g0=A*x0+b;
G0=A;
H0=eye(2);
k=0;
d0=-H0*g0;
a0=-d0'*g0/(d0'*G0*d0);
s0=a0*d0; %x(k+1)-x(k);
y0=A*a0*d0; %g(k+1)-g(k);
x1=x0+a0*d0;
while (norm(s0)>=epson)
switch k
case{10}
x0=x1;
g0=A*x0+b;
H0=eye(2);
k=0;
d0=-H0*g0;
a0=-d0'*g0/(d0'*G0*d0);
s0=a0*d0;
x1=x0+a0*d0;
break
otherwise
g1=A*x1+b;
y0=A*a0*d0;
s0=a0*d0;
% H1=H0+s0*s0'/(s0'*y0)-H0*y0*y0'*H0/(y0'*H0*y0);
H1=H0+((1+y0'*H0*y0/(s0'*y0))*s0*s0'-H0*y0*s0'-s0*y0'*H0)/(s0'*y0);
k=k+1;
d1=-H1*g1;
a1=-d1'*g1/(d1'*G0*d1);
a0=a1;
d0=d1;
H0=H1;
s0=a0*d0;
x1=x1+a0*d0;
break
end
end
x=x1;
fx=0.5*x1'*A*x1+b'*x1+c;
运行结果如下：
》 [x,fx,k]=bianchidufa([3 -1;-1 1],[-2;0],0,[2;1])
H1 =
　0.4031 0.2578
0.2578 0.8945
fx = -1
x =
1.0000
1.0000
fx = -1
k = 1
　故函数极小点是点（1,1）
5. 用鲍威尔法求函数的极小点。

用鲍威尔法在MATLAB的M文件编辑器中编写的M文件，如下：
function [x,fx,k]=bowell(A,b,c,x0)%鲍威尔法
d01=[1;0];
d02=[0;1];
x02=[0;0];
esp=1e-12;%停机原则
k=0;%迭代次数
while norm(x0-x02)>=esp
k=k+1;
g01=A*x0+b;
a01=-d01'*g01/(d01'*A*d01);
x01=x0+a01*d01;
g02=A*x01+b;
a02=-d02'*g02/(d02'*A*d02);
x02=x01+a02*d02;
d10=x02-x0;
g10=A*x02+b;
a10=-d10'*g10/(d10'*A*d10);
x10=x0+a01*d01;
d01=d02;
d02=d10;
x0=x10;
end
x=x0;
fx=0.5*x'*A*x+b'*x+c;
运行结果如下：
[x,fx,k]=bowell([2 -2;-2 4],[-4;0],0,[2;1])
fx =
-8
x =
4
2
fx =
-8
k =
3

6. 用单纯形法求线性规划问题

用单纯形法在MATLAB的M文件编辑器中编写的M文件，如下：
%单纯形法matlab程序-danchunxingfa
% 求解标准型线性规划:max c*x; s.t. A*x=b; x>=0
% 本函数中的A是单纯初始表，包括:最后一行是初始的检验数，最后一列是资源向量b
% N是初始的基变量的下标
% 输出变量sol是最优解, 其中松弛变量（或剩余变量）可能不为0
% 输出变量val是最优目标值，kk是迭代次数
function [sol,val,kk]=danchunxingfa(A,N)
[mA,nA]=size(A);
kk=0; % 迭代次数
flag=1;
while flag
kk=kk+1;
if A(mA,:)<=0 % 已找到最优解
flag=0;
sol=zeros(1,nA-1);
for i=1:mA-1
sol(N(i))=A(i,nA);
end
val=-A(mA,nA);
else
for i=1:nA-1
if A(mA,i)>0&A(1:mA-1,i)<=0 % 问题有无界解
disp('have infinite solution!');
flag=0;
break;
end
end
if flag % 还不是最优表，进行转轴运算
temp=0;
for i=1:nA-1
if A(mA,i)>temp
temp=A(mA,i);
inb=i; % 进基变量的下标
end
end
sita=zeros(1,mA-1);
for i=1:mA-1
if A(i,inb)>0
sita(i)=A(i,nA)/A(i,inb);
end
end
temp=inf;
for i=1:mA-1
if sita(i)>0&sita(i)utb
N(i)=inb;
end
end
% 以下进行转轴运算
A(outb,:)=A(outb,:)/A(outb,inb);
for i=1:mA
if i~=outb
A(i,:)=A(i,:)-A(outb,:)*A(i,inb);
end
end
end
end
end ;

运行结果如下：

>> A=[1 1 1 0 4;
1 2 2.5 3 5;
1.1 2.2 -3.3 4.4 0];N=[3;4];[sol,val,kk]=danchunxingfa(A,N)
sol =
0 0 4.0000 1.6667
val =
7.3333
kk =
2
所以，
7. 求解线性规划问题

用单纯形法在MATLAB的M文件编辑器中编写的M文件，如下：
%单纯形法matlab程序-danchunxingfa
% 求解标准型线性规划:max c*x; s.t. A*x=b; x>=0
% 本函数中的A是单纯初始表，包括:最后一行是初始的检验数，最后一列是资源向量b
% N是初始的基变量的下标
% 输出变量sol是最优解, 其中松弛变量（或剩余变量）可能不为0
% 输出变量val是最优目标值，kk是迭代次数
function [sol,val,kk]=danchunxingfa(A,N)
[mA,nA]=size(A);
kk=0; % 迭代次数
flag=1;
while flag
kk=kk+1;
if A(mA,:)<=0 % 已找到最优解
flag=0;
sol=zeros(1,nA-1);
for i=1:mA-1
sol(N(i))=A(i,nA);
end
val=-A(mA,nA);
else
for i=1:nA-1
if A(mA,i)>0&A(1:mA-1,i)<=0 % 问题有无界解
disp('have infinite solution!');
flag=0;
break;
end
end
if flag % 还不是最优表，进行转轴运算
temp=0;
for i=1:nA-1
if A(mA,i)>temp
temp=A(mA,i);
inb=i; % 进基变量的下标
end
end
sita=zeros(1,mA-1);
for i=1:mA-1
if A(i,inb)>0
sita(i)=A(i,nA)/A(i,inb);
end
end
temp=inf;
for i=1:mA-1
if sita(i)>0&sita(i) N(i)=inb;
end
end
% 以下进行转轴运算
A(outb,:)=A(outb,:)/A(outb,inb);
for i=1:mA
if i~=outb
A(i,:)=A(i,:)-A(outb,:)*A(i,inb);
end
end
end
end
end ;
运行结果如下：

>> A=[9 4 1 0 0 360;
4 5 0 1 0 200;
3 10 0 0 1 300;
7 12 0 0 0 0]; N=[3;4;5]; [sol,val,kk]=danchunxingfa(A,N)
sol =
20 24 84 0 0
val =
420
kk =
3
所以，经3次转轴运算，得到的最优解为