基于MATLAB的鲍威尔法求极值问题

基于MATLAB的鲍威尔法求极值问题姓名：xxx 学号：xxx（北京理工大学机械与车辆学院车辆工程，北京 100081）摘要：无约束优化方法主要有七种，按照求导与否把这些方法分为间接法和直接法。

牛顿法的成败与初始点选择有极大关系，其可靠性最差；坐标轮换法、单纯形法和最速下降法对于高维优化问题计算效率很低，有效性差；由于编制变尺度法程序复杂，其简便性不足。

综合考虑后，鲍威尔法、共轭梯度法具有较好的综合性能。

本文首先对鲍威尔法的原理进行阐述，根据其迭代过程给出流程图，并编写MATLAB程序。

最后用此MATLAB程序求解实际的极值问题，并对求解结果进行简要分析。

1.鲍威尔法的基本思想1.1其他优化方法对鲍威尔法形成的影响通过对鲍威尔法的学习，可以很明显看出来其迭代思想中汲取了其他几种优化方法的核心思想。

为了更全面、更深入的学习鲍威尔法，很有必要对其他有影响的优化思想进行学习和梳理。

由最基本的数学基础知识可知，梯度方向是函数增加最快的方向，负梯度方向是函数下降最快的方向，于是，利用这个下降最快方向产生了最速下降法。

每次迭代都沿着负梯度方向进行一维搜索，直到满足精度要求为止。

其特点是相邻两个搜索方向互相正交，所以很明显的一个现象就是刚开始搜索步长比较大，愈靠近极值点其步长愈小，收敛速度愈慢，特别当二维二次目标函数的等值线是较扁的椭圆时，迭代速度更慢。

这时，倘若目标函数是等值线长、短轴都平行于坐标轴的椭圆形，则通过坐标轮换法可以很高效的解决问题。

通过两次分别沿坐标轴进行一维搜索，便可达到极值点。

但对于目标函数的等值线椭圆的长、短轴倾斜于坐标轴时，坐标轮换法的搜索效率也显得极低。

抛开这两种特殊情况，对于一般形态的目标函数，如果在某些明显可以直达最优点的情况下（一般为靠近极值点区域），迭代过程完全可以不沿负梯度方向搜索，取而代之的是找到直达最优点的方向，一步到位。

但这样的直达方向应该如何去找呢？共轭梯度法由此产生。

其基本原理是：任意形式的目标函数在极值点附近的特性都近似一个二次函数，其等值线在极值点附近为近似的同心椭圆簇，而同心椭圆簇有一个特性便是任意两条平行线与椭圆簇切点的连线必通过椭圆的中心。

而这个连线方向便是所寻找的直达方向。

通过对其迭代过程的分析，很明显可以看出需大量的求目标函数的一阶和二阶偏导数。

对于一些实际的机械优化问题，目标函数可能复杂到难以求取其偏导数或者根本就不存在，求取极值问题就显得无从下手。

所以，一个效率高的、适应性强的优化方法急需出现，而鲍威尔法便是这么一种综合的方法。

1964年，鲍威尔提出了对共轭方向法的改进方法——鲍威尔共轭方向法。

一维搜索法、共轭方向和坐标轮换法的思想在鲍威尔法中体现的淋漓尽致。

下面就对鲍威尔法的基本原理进行阐述。

1.2鲍威尔法的数学原理通过前文可知，鲍威尔法也算一种共轭方向法，但与共轭梯度法相比，不需要对函数求导，而是在迭代过程中逐次构造出用于搜索的共轭方法。

1). 对于二维无约束优化问题，采用鲍威尔法求解的迭代过程如图1-1所示。

任选一初始点0X ,令00X X =，按照坐标轮换法，选择两个单位向量1110d e ⎡⎤==⎢⎥⎣⎦和2201d e ⎡⎤==⎢⎥⎣⎦，以此作为搜索方向进行第一轮搜索得到02X 点。

2). 用00X 和02X 的连线方向构成新的搜索方向1d 。

从02X 点出发，沿1d 方向一维搜索得到10X 点，作为下一轮搜索的初始点。

3). 从10X 出发，依次沿1201d e ⎡⎤==⎢⎥⎣⎦和12d d =方向进行一维搜索，得到点11X 和点12X 。

4). 用点10X 和点12X 的连线方向21120d X X =-构成新的搜索方向。

10X 和12X 是从两个不同点出发沿相同方向1d 搜索得到的，所以1d 与10X 和12X 的连线方向2d 互为共轭方向。

从12X 点出发，沿2d 方向一维搜索得到20X 点。

因20X 是从02X 点出发依次沿两个互为共轭的方向1d 和2d 进行两次一维搜索得到的，所以20X 就是该二维二次函数的极小点*X 。

图1-1 二维情况下鲍威尔法的迭代过程将上述二维优化问题扩展到n 维的情况，得到鲍威尔法的基本迭代过程：从初始点出发依次沿n 个线性无关的方向组进行一维搜索得到一个终点，沿初始点和终点的连线方向一维搜索得到下一轮迭代的初始点，并以这个方向作为下一轮迭代方向组中的最后一个方向，同时去掉第一个方向，组成的新方向组进行第一轮迭代。

若目标函数是个n 维的正定二次函数，则经过这样的n 轮迭代以后，就可以收敛到最优点。

但是这种方法有一个缺陷，通过这种方法产生的n 个新方向，有可能是线性相关或近似线性相关的。

因为新方向01nkk k i ni i d X X d α==-=∑，如果其中出现了0i α=，则k d 就可以表示为2d ，3d ，…，n d 的线性组合，这样用新方向kd 替换1d 后，新的坐标轮换方向组就成为线性相关的一组向量，以后的各轮迭代计算将在维数下降了的空间内进行，这将导致算法收敛不到真正的最优点。

针对此现象，通过适当改进，产生了新的鲍威尔法。

1.3改进的鲍威尔法由上一小节讨论可知，如果对新方向组的形成加以适当的选择，防止其线性相关，则可以避免鲍威尔法的“退化”。

在这里直接给出是否用kd 方向替代m d 来组成新的搜索方向组的判别条件，供后面编程所用。

设()10k f f X =，()2k n f f X =，()302k k n f f X X =-。

其中最大者11max k k k m i m m i nf f -<<∆=∆=-。

若31f f <且221231213(2)()0.5()k k m m f f f f f f f -+--∆<∆-成立，则用kd 方向替代m d 方向，否则仍用原来的n 个搜索方向。

这样，改进后的鲍威尔法保证了对于非线性函数计算时收敛的可靠性。

2. 迭代过程和算法流程图2.1迭代过程1). 给定初值点0X ，n 个线性无关的初始向量组0100[,,,]n D d d d =及精度0ε>，置0k =。

2). 从0k k X X =点出发沿01[,,,]n k k D d d d =中的方向进行一维搜索：111111110()min ()k kj k kj j j j j k k j j j t X X X X d f X d f X d ααα--------≥⎧=⎪⎪=+⎨⎪+=+⎪⎩，其中1,2,,;j n =3). 判断是否满足收敛准则，若0k k n X X ε-<成立，则迭代停止，输出*k n X X =，*X 即为最优点；否则跳至第4步。

4). 计算第k 轮迭代中相邻两点目标函数值的下降量，并找出下降量最大者。

1k k k i i i f f -∆=-11max k k km i m m i nf f -<<∆=∆=-相应的方向：1kk km m m d X X -=-5). 沿共轭方向10kk kn n d X X +=-计算反射点。

102k k kn n X X X +=-6). 检验判别条件。

令()10kf f X =，()2k n f f X =，()31k n f f X +=，若同时满足：31f f <221231213(2)()0.5()k k m m f f f f f f f -+--∆<∆-则转第6步；否则转第7步。

7). 由k n X 出发沿1kn d +方向进行一维搜索，求出此方向极小点k X ，并令下一轮迭代初始点10k kX X +=。

同时令1111212111[,,,][,,,,,,]k k k kkkkkn m m n d d d d d d d d +++-++=，同时置1k k =+。

并转第2步。

8). 若32f f <，则101k k n X X ++=；否则10k kn X X +=。

置1k k =+，转第2步。

2.2算法流程图鲍威尔法的算法流程图如图2-1。

kki i d α沿方向一维搜索求最优步长1kk k k ii i i X X d α-=+1kk km m m d X X -=-n+1kd 沿方向一维搜索求最优步长1k n d X +=+110n 1kk k kn n X X dα++=+11(,1,,)k kii d d i m m n ++⇐=+1(1,2,,1)k k iidd i m +⇐=-图 2-1 鲍威尔法的算法流程图3. 鲍威尔法的MATLAB 程序function [x,minf,k,allX,allY] = Powell_Graphic(f,x0,D,var,eps) format long;if nargin == 4 %默认精度为1.0e-2 eps = 1.0e-2;endn = length(var)+1;k = 0;syms l;while 1tx = zeros(size(D));tx(:,1) = x0;for i=1:n-1 %在每个方向上开始一维搜索xv = tx(:,i) + l*D(:,i);fx = Funval(f, var,xv);[a,b] = minJT(fx,0,0.1); %一维搜索进退法确定搜索区间tl = minHJ(fx,a,b); %黄金分割法求最优点tx(:,i+1) = tx(:,i) + tl*D(:,i);endD(:,n) = tx(:,n) - tx(:,1); %判断收敛性是否满足精度要求if norm(D(:,n)) <= epsx = tx(:,n);allX(:,k*n+1:(k+1)*n) = tx(:,1:n);for j=1:nallY(:,k*n+j) = Funval(f, var,tx(:,j));endbreak;elsefor j=1:nFX(j) = Funval(f, var,tx(:,j));endmaxDF = -inf;m = 0;for j=1:n-1df = FX(j) - FX(j+1);if df > maxDFmaxDF = df;m = j+1;endendtmpF = Funval(f, var,2*tx(:,n)-tx(:,1));fl=(FX(1)-2*FX(n-1)+tmpF)*(FX(1)-FX(n-1)-maxDF)^2;if fl<0.5*maxDF*(FX(1)-tmpF)^2 && tmpF<FX(1)%检验判断条件决定换方向 xv = tx(:,n) + l*D(:,n);fx = Funval(f, var,xv);[a,b] = minJT(fx,0,0.1);tl = minHJ(fx,a,b);x0 = tx(:,n) + tl*D(:,n); D(:,m:(n-1)) = D(:,(m+1):n); elseif tmpF < FX(n-1)x0 = 2*tx(:,n)-tx(:,1); elsex0 = tx(:,n); end end k=k+1;allX(:,(k-1)*n+1:k*n) = tx(:,1:n); for j=1:nallY(:,(k-1)*n+j) = Funval(f, var,tx(:,j)); end end endminf = Funval(f,var,x); format short;程序中用到的函数minJT 、minHJ 和Funval 分别为进退法确定搜索区间函数、黄金分割法求极值函数和计算函数值函数，直接调用即可，由于篇幅限制，故不在此列出其代码。