基于MATLAB的可行方向法求极值问题讲解

基于MATLAB 可行方向法求极值的实现姓名：xxx学号：xxx（北京理工大学机械与车辆学院车辆工程，北京 100081）摘要：在工程实际的优化设计中，随着设计变量数和约束条件数的增加，随机方向搜索法和复合形法等直接优化解法的求解效率会偏低。

可行方向法，顾名思义，一种始终在可行域内寻找下降方向的搜索法，以其收敛速度快、效果好的优点已成为求解约束非线性问题的一种有代表性的直接解法，同时也是求解大型约束优化问题的主要方法之一。

本文将简单介绍可行方向法的数学思想，采用线性规划法和约束最优步长法编写MATLAB 程序，最后通过算例完成对优化问题的求解。

关键字：可行方向法；MATLAB ；优化方法。

1. 可行方向法的基本数学思想1.1可行方向法的搜索策略可行方向法迭代的第一步都是从可行域的某一初始点(0)X 出发，沿负梯度(0)()f X -∇方向移至某一个或J 个起作用约束面的交集()k X 上。

以后的搜索路线和迭代计算可根据约束函数和目标函数的不同性状，分别采用以下三种不同策略继续搜索。

1）由点()k X 出发，沿可行方向作一维最优化搜索，若所得新点(1)k X +在可行域内，则再沿(1)()k f X +-∇方向作一维最优化搜索；若所得的新点不在可行域内，则将它移至约束面上再反复重复上述步骤，若(1)()k f X ε+∇≤，则停止迭代，如图1.1所示。

2）由点()k X 出发，沿可行方向作一维最优化搜索，若所得新点(1)k X +在可行域外，则沿可行方向以最大步长到达另一个约束面上一点，将该点作为迭代点(1)k X +进行反复搜索，直至满足给出的K-T 条件，如图1.2所示。

3） 沿着约束面进行搜索。

对于只具有线性约束的非线性规划问题，如图1.3所示，从点()k X 出发，沿约束面移动，在有限的几步内即可搜索到约束最优点；对于非线性约束函数，如图1.4所示情况，就是沿着约束面的切线移动，但这样将会进入非可行域，使问题变得复杂。

此时，应将进入非可行域的新点X 设法调整回约束边界。

调整方法为：先规定约束面容差0δ>，建立新的约束边界，然后将已经离开约束面的X 点沿起作用约束的负梯度方向()g X -∇返回约束面上，计算公式为(1)()k t X X g X α+=-∇式中，t α为调整步长，可以用试探法确定，也可以用下式估计，使其回到约束边界上。

[]()()()t Tg X g X g X α=∇∇ (1-1)图1.1图1.2图1.3图1.4不管采用哪种搜索方法都要进行两条决策：一是产生一个适用的可行方向()k d；二是沿()k d方向确定一个不会越出可行域的适当的步长因子。

1.2产生可行方向和步长的数学原理可行方向要保证沿该方向作微小移动后所得的新点是可行点，且目标函数值有所下降。

显然可行方向应满足可行和下降两个条件。

这里从数学原理的角度分析可行性和下降性。

1.2.1可行性设()k X 为可行域D 中的一个点，即()k X D ∈。

对于某一方向()k d来说，若存在数00A >，使其对于任意的A 0(0)A A ≤≤，式1-2均成立，则称方向()k d 对()k X 点满足可行性条件。

()()k k X AdD +∈ （1-2）若点()k X 在可行域内，即()k X 为内点，则任何方向都满足可行性条件。

当点()k X 在某一约束边界上时，如图1.5所示，该约束面就是起作用约束。

当点()k X 位于一个约束面上时，如图1.6所示，作出()k X 初起作用约束的梯度()g X ∇和切线，可见只要()k d与起作用约束的梯度成直角或钝角，则可指向可行域内，满足可行性条件。

当()k X 同时处于J 个约束面时，就要求()k d与这J 个约束面的梯度均垂直或交成钝角，如图1.6所示，写成数学表达式如下：()()()0(1,2,,)Tk k j g X dj J ⎡⎤∇≤=⎣⎦式中，J 为起作用约束的个数。

上式称为可行性条件。

图1.5图1.61.2.2下降性沿一方向搜索时，要求下降愈快愈好。

对某一点来说，负梯度方向为最速下降方向。

如果负梯度方向是可行方向，那么负梯度方向是最有利的方向；否则为了保证目标函数值有所下降，至少要使搜索方向()k d 与目标函数的负梯度方向成锐角或与梯度方向成钝角，写成数学表达式为()()()0Tk k f X d⎡⎤∇<⎣⎦该式称为下降性条件。

如图1.7所示，可行下降方向显然位于点()k X 约束面的切线与目标函数等值线的切线所围成的扇形区域内，推广到一般的情况就是可行方向在目标函数超等值面的超切面和J 个起作用约束的超切面所围成的超锥体内。

该区域称为可行下降方向区。

图1.7综上所述，当点()k X 位于J 个起作用的约束面上时，满足()()()()()0(1,2,,)()0T k k j T k k g X d j J f X d ⎧⎡⎤∇≤=⎪⎣⎦⎨⎡⎤⎪∇<⎣⎦⎩ （1-3）式中，()k d即为可行方向。

1.3可行方向的产生本文采用线性规划法来产生可行方向，用约束最优步长确定最优步长，为后面编制MATLAB 程序做理论基础。

本节先阐述线性规划法，下节再对约束最优步长法进行描述。

线性规划法对包含线性和非线性的不等式约束的最优化问题都适用，但不允许有等式约束。

其基本原理是将具有一阶连续偏导数的目标函数和约束条件在()k X 点用Taylor 展开式展成线性近似函数（一次项），并用这些线性近似函数代替目标函数和它的约束条件，使得问题线性化。

这样min ()f X 成为()()()min{()[()]()}()k k T k n f X f X X X X R +∇-∈约束条件()0u g X ≤变为()()()()[()]()0(1,2,,)k k T k u u g X g X X X u m +∇-≤=用(1)k X +代替上式中的X ，得到()()()min{()[()]}k k k T k f X f X d α+∇s.t. ()()()()[()]0(1,2,,)k k k T u k u g X g X du m α+∇≤=式中，()()k f X 、()()k u g X 为常数，只有搜索方向()k d和k α是未知量。

当迭代点()k X 位于J 个起作用约束边界的交点上时，问题就变成求解线性规划问题：()()min[()]n k k TX Rf X d ∈∇s.t. ()()[()]0(1,2,,J)k k T u g X du ∇≤=式中，()k d只起到方向作用，故规定其向量的模是有界的，也就是规定()k d 每个分量的绝对值不大于1，即()j||1(1,2,,)k d j n ≤=由于可行方向还得满足式1-3的可行性和下降性条件，故确定可行方向的线性规划数学模型如下：()()()min ()[()]nk k k T X RZ d f X d∈=∇s.t. ()()[()]0k k T f X d ∇<()()[()]0(1,2,,J)k k T u g X du ∇≤=()j11(1,2,,)k d j n -≤≤= (1-4)求解后，可得到迭代方向()k d，作为下一步的搜索方向。

1.4最优步长的确定由线性规划法求的可行方向()k d后，由()k X 出发，找到该方向上的可行新点(1)k X +。

新点必须满足两个条件，一是新点(1)k X +必须还是可行点，二是目标函数值具有最大的下降量。

约束最优步长即可满足这两点，可用以下一维搜索法求解：()()min ()k k f X dαα+s.t. max 0αα≤≤式中，()()max sup{|()0(1,2,,)}k k u g X d u m αλα=+≤=。

2. 可行方向法的迭代过程和算法流程图2.1迭代过程1)输入收敛精度ε；2)形成可行初始点(0)X ，令k=1;3)确定起作用约束下标集合，()()(){|()0(1)}k k i E X i g X i m ==≤≤； 4)判断E 是否为空集，若是空集则说明()k X 在可行域内，转步骤5）；否则说明()k X 在约束边界上，转步骤6)；5)判断是否满足终止条件()()k f X ε∇≤，若满足，转步骤10)；若不满足，则取搜索方向为()()()k k df X =-∇，转步骤8)；6)采用线性规划法，用式1-4求解后，得到可行方向()k d ；7)判断0k Z =是否成立，若成立，则转步骤10)输出，停机结束；若不成立，则继续下一步；8)采用约束最优步长法进行一维搜索，得到新点(1)k X +; 9)令k=k+1，转步骤3)；10)输出*()k X X =，**()f f X =，停机结束。

2.2算法流程图算法流程图如图2-1。

()()()()()()()()()jmin ()[()]..[()]0[()]0(1,2,,J)11(1,2,,)k k k k Tk k T k k T u k d Z d f X ds t f X d g X d u d j n =∇∇<∇≤=-≤≤=求解线性规划，确定可行方向()()()k k df X =-∇k ()()max ()()max min ()..0sup{|()0(1,2,,)}k k k k u f X d s t g X d u m ααααααλα+≤≤=+≤=求解一维搜索问题，得到约束步长其中()(1)()1k k k k X X d k k α+⇐+⇐+图 2-1 可行方向法的算法流程图3. 可行方向法的MATLAB 程序function kxfxf(x0) eps=1.0e-6;%刚开始给的x0为行向量 x0=transpose(x0); funcsz=length(x0); G1 = G; for k=1:1:50%'f(x)梯度：'sf=diff_val(x0);sf=eval('sf');%'f(x)梯度范数：'norm_s=norm(diff_val(x0));GL=length(G);G1 = G;G_copy=G;for i=1:1:GLg=subs(G1(i,:),x,x0);G1(i,:)=g;endG_zero=eval('G1');for i=GL:-1:1if abs(G_zero(i,:))>0.1G_zero(i,:)=[];G_copy(i,:)=[];endendI=size(G_zero,1);add=-ones(I,1);%根据E是否为空集(I==0)确定不同情况下的搜索方向if I==0if norm_s<=3'最优解：''自变量值:'x0'目标函数值：'f=fval(x0)'迭代次数：'k'======================================================' returnelsed0=-diff_val(x0);%线性搜索d0=transpose(d0);vmax = gmin(x0,d0);h=fmin(x0,d0,vmax);x0=x0+h*d0;endelse%线性规划grad=jacobian(G_copy,x);G_zero=subs(grad,x,x0);Ac=[sf;G_zero];lb=-1*ones(sz,1);ub=ones(sz,1);p=zeros(I+1,1);dz = size(1,sz);c=sf*dz';[dz,mn,m1,m2,m3]=linprog(c,Ac,p,[],[],lb,ub);z0 = sf*dz;z0=abs(z0);if z0<0.01'最优解：''自变量值:'x0'目标函数值：'f=fval(x0)'迭代次数：'k'======================================================'returnelsed0=dz;%线性搜索vmax = gmin(x0,d0);h=fmin(x0,d0,vmax);x0=x0+h*d0;endendend程序中用到一些回调函数，简单列举如下，.m文件详细内容见附录一。