- 积分
- 83
- 注册时间
- 2003-11-14
- 仿真币
-
- 最后登录
- 1970-1-1
|
本帖最后由 bainhome 于 2013-3-23 22:52 编辑
从前写过坐标轮换法多维无约束优化程序,当时用6,没匿名函数,加之水平眼界有限,符号函数实现,又笨又慢。最近因接触现代设计方法这本书,看到这个算法流程,来了兴致重写一遍。
言归正传,坐标轮换法是利用函数值信息直接寻优最基础的一个算法,迭代步数较多,效能取决于目标函数的形式,目标函数等值线椭圆形长短轴与坐标轴不平行时迭代计算曲折,收敛速度慢。对很多问题搜索不到合理解,算法本身没太多值得提倡的东西,不过匿名函数的用法倒还不错,另外算法实现方面,rocwoods这货今天给出重要提示,因此拿来和大家共享。
- function [ outX,outY,count ] = CoodinateOpt( fx,x0,epsx,epsf )
- % ================================
- % 用于解决多维函数无约束优化问题
- % 算法:坐标轮换法
- % 用户测试代码:
- % fx=@(x)100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
- % [ outX,outY,count ] = CoodinateOpt( fx,[1,1],1e-2,1e-3 )
- % ================================
- clc
- tic
- lenX=length(x0);
- x00=x0(:);
- xStart=x0(:);
- xn1=100*eye(lenX,1);
- dir1=eye(lenX);
- fx0=fx(x0);
- count=0;
- fn=1000;
- while norm(xn1-x00)>epsx&&abs(fn-fx0)>epsf&&count<=5000
- xIter=xStart;
- x00=xStart;
- fx0=fx(xStart);
- xmin=0;
- for i=1:lenX
- f=@(a) fx(xIter+a*dir1(:,i));
- xmin=Xquad(f,xmin,.01,.001);
- xIter=xIter+xmin*dir1(:,i);
- end
- xStart=xIter;
- xn1=xIter;
- fn=fx(xn1);
- count=count+1;
- end
- outX=xn1;
- outY=fx(xn1);
- toc
- function xmin=Xquad(f,x1,h,t)
- % ***********************************************
- % 二次插值法计算最优步长子程序
- % ***********************************************
- clc
- [ax,bx,countjintui]=jintui(f,x1,h);
- count=0;
- a1=(ax+bx)/2;
- xx=[ax,a1,bx];
- for i=1:3
- Y(i)=f(xx(i));
- end
- c1=(Y(3)-Y(1))/(bx-ax);
- c2=((Y(2)-Y(1))/(a1-ax)-c1)/(a1-bx);
- xp=.5*(ax+bx-c1/c2);
- if abs(f(xp)-f(a1))<=t
- xmin=min([xp,a1]);
- fmin=f(xmin);
- step_opt=xmin-x1;
- else
- while abs(f(xp)-f(a1))>=t
- [ax,a1,xp,bx]=quadCompute(f,ax,a1,xp,bx);
- end
- xmin=min([xp,a1]);
- fmin=f(xmin);
- step_opt=xmin-x1;
- count=count+1;
- end
- count;
- function [ax,a1,xp,bx]=quadCompute(f,ax,a1,xp,bx)
- if f(a1)<f(xp)&&xp>a1
- bx=xp;
- elseif f(a1)>f(xp)&&xp>a1
- ax=a1;
- a1=xp;
- elseif f(a1)<f(xp)&&xp<a1
- ax=xp;
- elseif f(a1)>f(xp)&&xp<a1
- bx=a1;
- a1=xp;
- end
- xx=[ax,a1,bx];
- for i=1:3
- Y(i)=f(xx(i));
- end
- c1=(Y(3)-Y(1))/(bx-ax);
- c2=((Y(2)-Y(1))/(a1-ax)-c1)/(a1-bx);
- xp=.5*(ax+bx-c1/c2);
- function [a,b,countjintui]=jintui(f,x1,h)
- % ================================
- % 进退法确定搜索区间
- % ================================
- f1=f(x1);
- x2=x1+h;
- f2=f(x2);
- countjintui=0;
- if f1<f2
- h=-h;
- x2=x1;
- f2=f1;
- x3=x2+h;
- f3=f(x3);
- else
- h=2*h;
- x3=x2+h;
- f3=f(x3);
- end
- while f2>f3
- if countjintui<5000
- x1=x2;
- x2=x3;
- x3=x2+h;
- f1=f2;
- f2=f3;
- f3=f(x3);
- countjintui=countjintui+1;
- else
- disp('该搜索区间不存在“高-低-高”的正确搜索形态,请修改后再试!')
- return
- end
- end
- a=min(x1,x3);
- b=max(x1,x3);
复制代码 基本思想是将一个无约束n维优化问题转化为依次沿着相应的n个坐标轴方向的一维优化问题求解。
以二维目标函数f(X)=f(x1,x2)为例,从某个初始点X^(0)作为第一轮迭代起始点X0^(1)。
1.以x1坐标轴的单位向量$e_1=[1,0]^T作为搜索方向,用一维优化方法确定其最优步长k1^(1),获得第一轮第一个迭代点X1^(1)=X0^(1)+k1^(1)e1;
2.以X1^(1)为新起点,改以x2坐标轴的单位向量e2=[0,1]^T为搜索方向,用一维优化方法确定其最优步长k2^(1),获得第一轮第二个迭代点X2^(1)=X1^(1)+k2^(1)e2
3.对二维问题,经过依次沿着e1和e2的两次搜索,完成第一轮迭代,得到第一轮迭代极小点X^(1)=X2^(1);
4.如经过第一轮迭代计算后仍然没有得到收敛到目标函数f(X)的极小点X*,则以第一轮迭代的极小点X^{1}作为第二轮迭代的起始点X2^(0)=X^(1),依次再沿着e1和e2方向搜索,按上述方法做第3,4,...搜索,直至迭代逼近目标函数f(X)的极小点X*。
难点:例如如下二维无约束问题目标函数应如何使其被降为一维搜索问题,以确定每一次的最优步长因子。
- fx=@(x) 100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
复制代码 解决办法是构造如下匿名函数:
- fx=@(x)100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
- d1=[1:5;1:5]';
- d2=.5*d1;
- x=arrayfun(@(i) fminsearch(@(a) fx(d1(i,:)+a*d2(i,:)),0),1:length(d1))
复制代码 计算结果:
- x =
- 1.5826 -0.20869 -0.80581 -1.1044 -1.2835
复制代码 匿名函数的确是MATLAB中比较有特色的编程方法。
用该算法计算上述二维无约束极值问题,结果如下:
- Elapsed time is 0.429963 seconds.
- outX =
- 1.7913
- 1.7914
- outY =
- 9.6461e-08
- count =
- 274
复制代码
fminunc命令的计算结果:
x =
1.7913 1.7913
fval =
7.3758e-12
exitflag =
1
output =
iterations: 17
funcCount: 69
stepsize: 1
firstorderopt: 5.3476e-08
algorithm: 'medium-scale: Quasi-Newton line search'
message: [1x438 char]
Elapsed time is 0.251491 seconds.
fminunc采用拟牛顿法17次迭代到1e-8容差,自编程序迭代274次,再次证明梯度法收敛速度的确不是盖的。
ps1:从前倒看过很多优化问题书籍,很多优化算法仍在使用inline函数定义目标函数,这个命令的效率比较低。
ps2:很久没来,冒个泡再急速下潜...
|
评分
-
1
查看全部评分
-
|