matlab最优化实验

6.1知识要点与背景

6.1.1 由简入繁： 最佳水槽断面问题的推广

6.1.2 微分法求最大和最小

◆问题4      (1)求受检点： ，

zxy6_1.m

【    syms x1 x2                                          %定义符号变量。

f=x1^3-x2^3+3*x1^2+3*x2^2-9*x1;   % 函数z。

v=[x1 x2];df=jacobian(f,v)                %计算雅可比 。

[X,Y]=solve(df(1),df(2));[X,Y]         % 用指令solve求驻点。             】

zxy6_2.m

【    clf,xmin=-5;xmax=3.5;ymin=-3;ymax=5;

x1=linspace(xmin,xmax,30);x2=linspace(ymin,ymax,30);

[X1,X2]=meshgrid(x1,x2);Z= X1.^3.-X2.^3+3*X1.^2+3*X2.^2-9*X1;

contour(X1,X2,Z,60);,hold on,  xp=[-3,1,-3,1];yp=[0 0 2 2];

plot(xp,yp,'ro'),axis([xmin xmax ymin ymax]),colorbar

xlabel('x_1'),ylabel('x_2'),

for i=1:length(xp)

text(xp(i),yp(i),['\leftarrow (',num2str(xp(i)),',',num2str(yp(i)),')'] )

end              】

6.2  实验与观察(Ⅰ)：模拟盲人下山的迭代寻优法

zxy6_3.m(盲人下山的模拟)

【      clf, a=-2;b=4;

xmin=a;xmax=b;ymin=a;ymax=b;    %设置变量范围和坐标轴显示范围。

x1=linspace(xmin,xmax,100);x2=linspace(ymin,ymax,100);

[X1,X2]=meshgrid(x1,x2);

[Z,DZ1,DZ2]=zxy6_3f(X1,X2); %计算函数和梯度向量。

contour(X1,X2,Z,30),               %画等值线图。

axis([xmin xmax ymin ymax]),hold on,

axis equal,                              %该命令将使横轴、纵轴具有相同比例，避免失真。

plot([1.46808510638298],[1.148936170212776],'o'),  %标注最优点。

axis([xmin xmax ymin ymax])

x=[];y=[];            %开始用鼠标选点，按左键选点，按右键中止选点过程。

disp('Select a point by put on mouse left-key')

%disp指令,在命令窗口显示文字。

disp('Stop selecting point by put on mouse right-key')

button=1;         %button和ginput命令结合使用可用鼠标选点， 按左键时button＝1。

x=[];y=[];

while button==1

[xi,yi,button]=ginput(1);

%ginput(n)用鼠标选n个点，xi,yi分别为点的横坐标和纵坐标。

plot([xi],[yi],'r.','MarkerSize',10)，hold on,             %画所选的点。

[zi,dz1,dz2]=zxy6_3f(xi,yi);                       %计算函数值和梯度方向。

v=zi;

contour(X1,X2,Z,[v v],'-'),             %在点所在的高度画一条等高线。

axis([xmin xmax ymin ymax]),

x=[x,xi];y=[y,yi];

H_line2=plot(x,y);                                  %画已走的路径连线。

set(H_line2,'color','red','linewidth',2);     %设置颜色和线宽。

xt=xi-dz1;yt=yi-dz2;

H_line=plot([xi xt],[yi yt],'k:','linewidth',1);    %画最速下降方向路径。

end           %若按左键button=1，继续循环。若按右键，button~=1,循环终止 。 】

zxy6_3f.m(模拟山谷的二次函数程序)

【     function [f,df1,df2]=zxy6_4f(x1,x2)

f=8*x1.*x1+9*x2.*x2-10*x1.*x2-12*x1-6*x2;       %计算函数值。

if nargout > 1

df1=2*8*x1-10*x2-12*ones(size(x1));                %计算梯度向量。

df2=2*9*x2-10*x1-6*ones(size(x2));

end               】

6.3 .实验与观察(Ⅱ)：Matlab优化工具箱简介

6.3.1多元函数无约束优化指令fminunc和fminsearch

1. 观察：运行香蕉函数的优化程序bandemo.m

2. 使用fminunc和fminsearch指令

◆ 观察：用inline生成函数。

【      f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2'),

x=[2,2],y=f(x),    %代入一个点计算看看效果 。             】

3. bandemo.m的简化和剖析

zxy6_4.m

【 clf;  clear

%以下程序段是画香蕉函数图形。

xx = [-2:0.125:2]'; yy = [-1:0.125:3]'; [x,y]=meshgrid(xx',yy') ;

meshd = 100.*(y-x.*x).^2 + (1-x).^2; conts = exp(3:20);

xlabel('x1'),ylabel('x2'),title('Minimization of the Banana function')

contour(xx,yy,meshd,conts), hold on

plot(-1.9,2,'ro'),text(-1.9,2,'Start Point')

plot(1,1,'ro'),text(1,1,'Solution')

%优化程序段开始。

x0=[-1.9,2];        %赋初值。

l=1;

while l           %while 语句是可以重复运行下面的程序段，直至l=0退出循环。

clc                 %清除命令窗口的全体内容 。

%以下程序段是在命令窗口显示相应的文字内容 。

disp(' ')

disp('   Choose any of the following methods to minimize the …

banana   function')

disp('')

disp('    UNCONSTRAINED:    1) BFG direction ')

disp('                                     2) DFP direction')

disp('                                     3) Steepest Descent direction')

disp('                                     4) Simplex Search')

disp('                                     0) Quit')

method=input('Select method : ');  % input 从键盘输入控制变量method数据。

switch method         %Switch体开始。

case 0             %当method＝0，终止程序。

hold off

disp('End of demo')

break                         %break指令：中止程序。

case 1              %当method＝1，采用BFGS法。

clf,hold on    %每一个case中重新画等值线图,下面的程序段是重新画图。

xlabel('x1'),ylabel('x2'),

title('Minimization of the    Banana function')

contour(xx,yy,meshd,conts)

plot(-1.9,2,'ro'), text(-1.9,2,'Start Point')

plot(1,1,'ro'), text(1,1,'Solution')

% 这里是学习的重点： OPTIONS是控制fminunc和fminsearch指令的重要参数，

%用optimset('属性'，'属性值',…)指令改变设置，可以容易地控制算法。

OPTIONS=optimset('LargeScale','off');

%fminunc默认的大规模算法是“信赖域方法”，这是一种有效的算法；

%将LargeScale的属性设置为off时,fminunc的默认中等规模的算法就是BFGS方法。

OPTIONS = optimset(OPTIONS,'gradobj','on'); %使用解析梯度。

%定义梯度函数和画图函数banplot6_4。

GRAD=inline('[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;…

100*(2*x(2)-2*x(1)^2); banplot6_4(x)]');

f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2'); %定义目标函数。

disp('[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);');

%(调用fminunc指令，输出x,fval分别为最优点和最优函数值，exitflag和output

% 提供算法的一些信息,读者可在程序结束后，键入output或exitflag查看这些信息)

[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);

hold off

disp(' ')

disp('Strike any key for menu')

pause

case 2      %当method＝2，采用DFP法。

clf,  xlabel('x1'),ylabel('x2'),

title('Minimization of the Banana function')

contour(xx,yy,meshd,conts),  hold on

plot(-1.9,2,'ro'),   text(-1.9,2,'Start Point')

plot(1,1,'ro'),     text(1,1,'Solution')

OPTIONS=optimset('LargeScale','off');

OPTIONS = optimset(OPTIONS,'gradobj','on');

OPTIONS=optimset(OPTIONS,'HessUpdate','dfp');

% 将HessUpdate属性设置为dfp就使fminunc指令采用DFP法。

GRAD=inline('[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;…

100*(2*x(2)-2*x(1)^2); banplot6_4(x)]');

f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');

disp('[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);');

[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);

hold off

disp(' ')

disp('Strike any key for menu')

pause

case 3   %当method＝3，采用最速下降法。

clf,   xlabel('x1'),ylabel('x2'),

title('Minimization of the Banana function')

contour(xx,yy,meshd,conts)

hold on

plot(-1.9,2,'ro'),  text(-1.9,2,'Start Point')

plot(1,1,'ro'),     text(1,1,'Solution')

OPTIONS=optimset('LargeScale','off');

OPTIONS = optimset(OPTIONS,'gradobj','on');

OPTIONS=optimset(OPTIONS,'HessUpdate','steepdesc');

%将HessUpdate属性设置为steepdesc就使fminunc指令采用最速下降法。

GRAD=inline('[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;…

100*(2*x(2)-2*x(1)^2); banplot6_4(x)]');

f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');

disp('[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);');

[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);

hold off

disp(' ')

disp('Strike any key for menu')

pause

case 4          %当method＝4，采用单纯形方法。

clf,hold on,  xlabel('x1'),ylabel('x2'),

title('Minimization of the Banana function')

contour(xx,yy,meshd,conts),

plot(-1.9,2,'ro'),  text(-1.9,2,'Start Point')

plot(1,1,'ro'),      text(1,1,'Solution')

OPTIONS=optimset('LargeScale','off');

OPTIONS = optimset(OPTIONS,'gradobj','off');

%该方法不使用导数，所以要设置gradobj属性为off。

f=inline('[100*(x(2)-x(1)^2)^2+(1-x(1))^2; banplot6_4(x)]');

%如果要画迭代过程的中间图，就要编制一个画图程序 banplot6_4，

% 套用本程序的格式定义目标函数。

disp('[x,fval,exitflag,output] = fminsearch(f,x0,OPTIONS);');

[x,fval,exitflag,output] = fminsearch(f,x0,OPTIONS);

%fminsearch 是多变量函数寻优的单纯形法指令，用法和fminunc是类似的。

hold off

disp(' ')

disp('Strike any key for menu')

pause

end

end                                       】

banplot6_4.m

【function out =  banplot6_4(x)

plot(x(1),x(2),'O','Erasemode','none')

drawnow;         % Draws current graph now

out = [];                              】

6.3.2其它的优化算法指令

1.多变量约束优化指令fmincon

2. 线性规划linprog指令

【    f = [-5; -4; -6]; A =  [1 -1  1;3  2  4;3 2 0];b = [20; 42; 30];

lb = zeros(3,1); x0=[1,1,0];

options=optimset('Diagnostics','on', 'largescale','off');

%查看诊断信息并采用单纯形算法

[x,fval,exitflag,output,lambda]= …linprog(f,A,b,[],[],lb,[],x0,options)

%没有等式约束且变量无上界,故需置Aeq=[ ];Aeb=[ ]；ub=[ ];               】

3. 二次规划quadprog指令

4. 一元函数寻优fminbnd指令

5.非线性最小二乘指令lsqnonlin和非线性数据拟合指令lsqcurvefit

程序如下(zxy6_5.m)

【    clf;x=1:10; y=2+2*x;                                    %选直线上的10个点。

a0=[0.1,0.4]; y1=exp(x*a0(1))+exp(x*a0(2));  %计算一条曲线。

[a,resnorm,residual] = lsqnonlin('zxy6_5f',a0); % 求最优解初始点a0。

disp('a='), disp(a);

disp('resnorm='), disp(resnorm)

y2=exp(x*a(1))+exp(x*a(1));

plot( －x,y,x,y1,'r:',x,y2,'o-',x,residual,'.-'),grid on

legend('直线','猜测的曲线','解曲线','残差')                         】

函数子程序为(zxy6_5f.m)

【    function F = zxy6_5f(a)

x = 1:10;

F = 2 + 2*x-exp(a(1)*x)-exp(a(2)*x);                    】

{     Optimization terminated successfully:

Norm of the current step is less than OPTIONS.TolX

a=    0.2578    0.2578

6.4 应用、思考与练习

6.4.1 .计算最佳水槽断面面积

zxy6_6S.m

【      function s=zxy6_6S(x)

l=24;a(1)=x(3);a(2)=x(4);   xs0=(0.5*l-x(1)-x(2));

xs1=xs0+x(1)*cos(a(1));     xs2=xs1+x(2)*cos(a(2));

h1=x(1)*cos(a(1));             h2=x(2)*cos(a(2));

s=(xs0+xs1)*h1+(xs1+xs2)*h2;s=-s;                    】

zxy6_6.m

【    clf,A=[1,1,0,0;-1,-1,0.0,0.0];b=[12,0]';

lb=[0,0,0,0]';ub=[100,100,pi/2,pi/2]';x0=[4,4,pi/3,pi/3]';

[s,fval] = fmincon('zxy6_6S',x0,A,b,[],[],lb,ub)       %最优化计算

%以下是绘制最优断面的图形，首先计算坐标点。将底边放在x轴上，并让断面关于

%y轴对称。逆时针计算坐标点，使之成为一个封闭的图形)

x(1)=(24-2*s(1)-2*s(2))/2; y(1)=0;

x(2)=x(1)+s(1)*cos(s(3));  y(2)=s(1)*sin(s(3));

x(3)=x(2)+s(2)*cos(s(4));  y(3)=y(2)+s(2)*sin(s(4));

x(4)=-x(3); y(4)=y(3);x(5)=-x(2); y(5)=y(2);x(6)=-x(1); y(6)=y(1);

x(7)=x(1);y(7)=y(1);        %首尾相接。

%plot(x,y),axis equal       %用这命令可画出封闭图形。

patch(x,y,'y'); axis equal,  %用patch命令画块对象并填充颜色。     】

{      s =    4.8000    4.8000    0.6283    1.2566

fval =  -88.6373                                            }

6.4.2  对约束优化的讨论

6.4.3.工程优化问题的计算

1. 啤酒配方问题: 线性规划

2.  储能飞轮的设计

3.  齿轮减速器设计

200多个MATLAB经典教程和MATLAB论文请查看：matlab教程