可行方向法

可行方向法摘要可行方向法是求解最优化问题的重要方法，在可行方向法求解过程中，一般需要构造一个求解可行下降方向的子问题，而可行方向法的不同取决于所采用的求解可行下降方向的子问题，它具有如下特点：迭代过程中所采用的搜索方向为可行方向，所产生的迭代点列是中在可行域内，目标函数值单调下降，由此可见，很多方法都可以归入可行方向法一类，本文主要介绍Frank-Wolf 方法。

一、问题形式min ().. 0f x Ax b s t x ≥⎧⎨≥⎩ （11.1）其中A 为m n ⨯矩阵，m b R ∈，n x R ∈。

记{},0,n D x Ax b x x R =≥≥∈并设()f x 一阶连续可微。

二、算法基本思想D 是一个凸多面体，任取0x D ∈，将()f x 在0x 处线性展开 000()()()()()T L f x f x f x x x f x ≈+∇-= 用min ().. 0L f x Ax b s t x ≥⎧⎨≥⎩ 或 0min ().. 0T f x xAx b s t x ∇≥⎧⎨≥⎩ （11.2）逼近原问题，这是一个线性规划问题，设0y D ∈是其最优解。

1）若000()()0T f x y x ∇-=，则0x 也是线性规划问题（11.2）的最优解，此时可证0x 为原问题的K-T 点。

2）若000()()0T f x y x ∇-≠，则由0y 是（11.2）的最优解，故必有0000()()T T f x y f x x ∇<∇从而 000()()0T f x y x ∇-<即00y x -为()f x 在0x 处的下降方向，沿此方向作有约束的一维搜索 00001min (())f x y x λλ≤≤+-设最佳步长因子为0λ，令100000000()(1)()x x y x y x D λλλ=+-=+-∈当λ充分小时100000001()min (())(())f x f x y x f x y x λλλ≤≤=+-≤+-00000()()()()()T f x f x y x o f x λλ=+∇-+< 用1x 取代0x ，重复以上计算过程。

三、算法迭代步骤1）给定初始点0x D ∈，允许误差0ε>，:0k =。

2）求解线性规划问题min ()k T x Df x x ∈∇，得最优解k y D ∈。

3）若()()k T k k f x y x ε∇-≤，Stop,*k x x ≈；否则go to 4)。

4）进行一维搜索01min (())k k k f x y x λλ≤≤+-，得最优步长因子k λ；令1()k k k k k x x y x λ+≈+-，:1k k =+，go to 2) 四、算法收敛性定理定理11.1设非线性规划问题（11.1）的最优解存在，且对算法产生的点列{}k x ，线性规划问题min ().. 0k T f x xAx b s t x ∇≥⎧⎨≥⎩ （11.3）的最优解总存在。

则1） 若迭代到某步0k ，有000()()0k k k T f x y x ∇-=，则0k x 为问题（11.1）的K-T 点； 2） 若情形1）总不发生，则算法产生一有界无穷点列{}k x ，其任意极限点都是原问题（11.1）的K-T 点。

证明：若情形1）出现，则0k x 也是问题（11.3）的最优解，故0k x 满足（11.3）的K-T 条件：00000000000()0()00,0,0,k T kT k T k k f x A u v u Ax b v x u v x Ax b ⎧∇--=⎪-=⎪⎨=⎪⎪≥≥≥≥⎩（11.4） 而(11.1)的K-T 条件：()0()00,0,0,T TTf x A u v u Ax b v x u v x Ax b ⎧∇--=⎪-=⎪⎨=⎪⎪≥≥≥≥⎩（11.5） （11.4）表明0k x ，0u ，0v 一起满足K-T 条件（11.5），故0k x 是原问题的K-T 点。

2）由点列{}k x 包含在D 的极点的凸组合中，而{}k y 均为D 的极点，故{}k x 、{}k y 均为有界点列。

设x ∞为{}k x 的任一极限点，即存在子列{}i k x ，使得： lim i i k k x x ∞→∞=注意到点列{}k x 满足：1()k k k k k x x y x λ+=+- 考虑点列{}i k y 、{}i k λ、1{}i k x +，不失一般性，设lim i i k k y y ∞→∞=， lim i i k k λλ∞→∞=， 11lim i i k k x x +∞+→∞=否则，可以通过反复抽取子序列，使上式对某个子序列成立。

由 i k y 是min ()i k T x Df x x ∈∇的最优解，故x D ∀∈，有()()i i i k k k T T f x x f x y ∇≥∇且 ()()0i i i k k k T f x y x ∇-<再由 (())(())i i i i i i i k k k k k k k f x y x f x y x λλ+-≥+- []0,1λ∀∈ 及 1()i i i i i k k k k k x x y x λ+=+-取极限，有 1()(),()()0(())(())()T T T f x x f x y f x y x f x y x f x y x x x y x λλλ∞∞∞∞∞∞∞∞∞∞∞∞∞∞+∞∞∞∞⎧∇≥∇∇-≤⎪+-≥+-⎨⎪=+-⎩（11.6）不等式组（11.6）等价于：011min ()()()()0min (())(())()T T x DT f x x f x y f x y x f x y x f x y x x x y x λλλλ∞∞∞∈∞∞∞∞∞∞∞∞∞∞≤≤∞+∞∞∞∞⎧∇=∇⎪⎪∇-≤⎪⎨+-=+-⎪⎪=+-⎪⎩ （11.7） 若能证明 ()()0T f x y x ∞∞∞∇-= 即x ∞为问题 min ()T x Df x x ∞∈∇的最优解，由本定理的第一部分可知，x ∞为原问题K-T 点。

下证：()()0T f x y x ∞∞∞∇-=，若不然，由（11.7）即知必有()()0T f x y x ∞∞∞∇-<故y x ∞∞-为x ∞处的下降方向，因而当λ充分小时，有 (())()f x y x f x λ∞∞∞∞+-<进而有： 1()()f x f x ∞+∞< （11.8） 但{()}k f x 为单调下降的有界序列，故lim ()k k f x →∞存在而且 lim ()lim ()()i i k k k k f x f x f x ∞→∞→∞==11lim ()lim ()()i i k k k k f x f x f x +∞+→∞→∞==即 1()()f x f x ∞∞+=与（11.8）矛盾，故必有()()0T f x y x ∞∞∞∇-=。

五、流程图min ()()11,k Tk T j i f X g X d ηη⎧⎪≤⎪⎨-⎪⎪-≤≤⎩η≤}n六、部分算法说明及代码编写目标函数文件myf.m如下：function f=myf(x)%目标函数f=4*x(1)^2+x(2)^2-32*x(1)-34x(2);编写约束函数文件mycf.m如下：function[c,ceq]=confun(x)%约束函数%非线性不等式约束c=[x(1);x(2);2-x(1);6-x(1)-2*x(2)];%非线性等式约束ceq=[];%初始条件x0=[0,0]；%目标函数的梯度文件dmyf.mfunction g=dmyf(x)g=[8*x(1)-32;2*x(2)-34];%约束函数的梯度文件dmycf.mfunction [dc,dceq]=dmycf(x)dc=[-1;-2];dceq=[];然后在命令窗口输入如下命令：x0=[0;0];function Zoutendiijkclear;x0=[0; 0];A=[2.0 -1.0;1.0 1.0;-1.0 0.0;0.0 -1.0]; b=[1.0;2.0;0.0;0.0];c=0;m=0;while c<5c=c+1;k=0;j=0;m=m+1;for i=1:4c=A(i,:)*x0;if c>=b(i)-1e-3k=k+1;A1(k,:)=A(i,:);b1(k,1)=b(i);endif c<b(i)-1e-3j=j+1;A2(j,:)=A(i,:);b2(j,1)=b(i);endendif (length(A1(:,:))==0)breakend%pausef=dfx(x0);b1=[-1 -1];bu=[1 1];b0=zeros(size(b1));d=linprog(f,A1,b0,[],[],b1,bu)if(abs(d)<1e-5)breakenddd=A2*d;bb=b2-A2*x0;lamdmax=max(bb./dd);options=optimset('display','off');lamda=fminbnd(@(lamda)fun(lamda,d,x0),0,lamdmax,options) if (length(lamda(:))==0)breakendx0=x0+lamda*d;endx0=[0,1]';[x,fval]=fmincon(@ f1,x0,A,b,[],[],[],[],[],options); x,fvalfunction f=f1(x)f=x(1)^2+x(2)^2-2*x(1)-4*x(2)+6;endfunction f=fun(lamda,d,x)xx=x+lamda*d;f=xx(1)^2+xx(2)^2-2*xx(1)-4*xx(2)+6;function dfx=dfx(x)dfx=[2*x(1)-2,2*x(2)-4];end七、数值试验及结果分析Optimization terminated.d =1.00001.0000lamda =0.9999Optimization terminated.d =-1.00001.0000lamda =0.5000Optimization terminated.d =1.0e-009 *-0.7803-0.5941Warning: Trust-region-reflective algorithm does not solve this type of problem, using active-setalgorithm. You could also try the interior-point or sqp algorithms: set the Algorithm option to'interior-point' or 'sqp' and rerun. For more help, see Choosing the Algorithm in thedocumentation.> In fmincon at 472In Zoutendiijk at 49x =0.50001.5000 fval =1.5000。