《数值最优化方法》课程教学课件（讲稿，打印版）序列二次规划法

最优化方法及其atlab程序设计 第十二章序列二次规划法 Back Close

1/84 JJ II J I Back Close Å`zê{9Ÿ Matlab ßSO 1õŸ Sg5y{

本章考虑求解一般非线性优化问题 min f(x), s.t.h(x)=0,i∈E={1,.,}, (12.1) g(x)≥0，i∈I={1,.,m} 的序列二次规划法(SQP,Sequential Quadratic Programming).SQP方 法是求解约束优化问题最有效的算法之一，其基本思想是：在每一迭 代步通过求解一个二次规划子问题来确立一个下降方向，以减少价值 函数来取得步长，重复这些步骤直到求得原问题的解.这也是之所以 称为序列二次规划法的由来.本章主要介绍SQP方法的基本思想、迭 代步骤和收敛性分析.首先介绍求解等式约束优化问题的牛顿-拉格朗 日法 Back Close

2/84 JJ II J I Back Close Ÿƒ¶)òÑöÇ5`zØK    min f(x), s.t. hi(x) = 0, i ∈ E = {1, · · · , l}, gi(x) ≥ 0, i ∈ I = {1, · · · , m} (12.1) Sg5y{ (SQP, Sequential Quadratic Programming). SQP ê {¥¶)Â`zØKÅké{Éò, Ÿƒgé¥: 3zòS ì⁄œL¶)òág5yfØK5(·òáe¸êï, ±~dä ºÍ5⁄, ­E˘ ⁄½Ü¶ØK). ˘è¥É§± °èSg5y{d5. ŸÃá0 SQP ê{ƒgé!S ì⁄½⁄¬Ò5©¤. ƒk0¶)™Â`zØK⁄Ó-.ÇK F{

§12.1牛顿-拉格朗日法 §12.1.1牛顿-拉格朗日法的基本理论 考虑纯等式约束的优化问题 min f(x), (12.2) s.t.h(x)=0,i∈E={1,.,l} 其中f:RnR,h,:Rn→R(i∈E)都是二阶连续可微的实函数. 记h(x)=(h1(x),·,hu(x)T,则不难写出该问题的拉格朗日函数为 L(z,L)=f(a)->uh:(a)=f(a)-jTh(a), i=1 7 其中μ=（山1，·，4）T为拉格朗日乘子向量.约束函数h(x)的梯度 矩阵为 Vh(x)=Vhi(x),.,Vhi(a), Back Close

3/84 JJ II J I Back Close §12.1 ⁄Ó-.ÇKF{ §12.1.1 ⁄Ó-.ÇKF{ƒnÿ ƒX™Â`zØK    min f(x), s.t. hi(x) = 0, i ∈ E = {1, · · · , l}, (12.2) Ÿ• f : R n → R, hi : R n → R (i ∈ E) —¥ÎYåᢺÍ. P h(x) = (h1(x), · · · , hl(x))T , KÿJ—TØK.ÇKFºÍè L(x, µ) = f(x) − X l i=1 µihi(x) = f(x) − µ T h(x), Ÿ• µ = (µ1, · · · , µl) T è.ÇKF¶fï˛. ÂºÍ h(x) F› › è ∇h(x) = ∇h1(x), · · · , ∇hl(x) ,

则h()的Jacobi矩阵为A(x)=Vh(x)T.根据问题(12.2)的KT条 件，可以得到如下的方程组 Vf(x)-A(a)Tu =0. (12.3) -h(c) 现在考虑用牛顿法求解上述的非线性方程组(12.3).记函数7L(x,4 的Jacobi矩阵为 N(c,四= -A() (12.4 其中 W(,四=VL(x,）=V2f)-Vh(c) 是拉格朗日函数L(x,)关于x的Hesse阵.(12.4)所定义的矩阵 N(a,)亦称之为KT矩阵.对于给定的点=(xk,),牛顿法的迭 Back Close

4/84 JJ II J I Back Close K h(x)  Jacobi › è A(x) = ∇h(x) T . ä‚ØK (12.2)  KT ^ á, å±Xeêß| ∇L(x, µ) =   ∇xL(x, µ) ∇µL(x, µ)   =   ∇f(x) − A(x) Tµ −h(x)   = 0. (12.3) y3ƒ^⁄Ó{¶)˛„öÇ5êß| (12.3). PºÍ ∇L(x, µ)  Jacobi › è N(x, µ) =   W(x, µ) −A(x) T −A(x) 0   , (12.4) Ÿ• W(x, µ) = ∇2 xxL(x, µ) = ∇2 f(x) − X l i=1 µi∇2hi(x) ¥.ÇKFºÍ L(x, µ) 'u x  Hesse . (12.4) §½¬› N(x, µ) ½°Éè KT › . Èuâ½: zk = (xk, µk), ⁄Ó{S

代格式为 2k+1=2k+Pk, (12.5) 其中pk=(d,V)满足下面的线性方程组： N(Tk;Bk)Pk =-VL(Tk:Lk), 即 W(k; -A(Zk) h(Zk) (12.6) 不难发现，只要矩阵A(xk)行满秩且W(xk,)是正定的，那么方 程组(12.6)的系数矩阵是非奇异的，且该方程有唯一解.由于KT条 件(12.3)是拉格朗日函数稳定点的条件，所以人们通常把基于求解方 程(12.3)的优化方法称为拉格朗日方法.特别地，如果用牛顿法求解 Back Close

5/84 JJ II J I Back Close ìÇ™è zk+1 = zk + pk, (12.5) Ÿ• pk = (dk, νk) ˜ve°Ç5êß|: N(xk, µk)pk = −∇L(xk, µk), =   W(xk, µk) −A(xk) T −A(xk) 0     dk νk   =   −∇f(xk) + A(xk) Tµk h(xk)   . (12.6) ÿJuy, êá› A(xk) 1˜ùÖ W(xk, µk) ¥½, @oê ß| (12.6) XÍ› ¥ö¤., ÖTêßkçò). du KT ^ á (12.3) ¥.ÇKFºÍ­½:^á, §±<Çœ~rƒu¶)ê ß (12.3) `zê{°è.ÇKFê{. AO/, XJ^⁄Ó{¶)

该方程组，则称之为牛顿拉格朗日方法.因此，根据牛顿法的性质，该 方法具有局部二次收敛性质 下面写出牛顿-拉格朗日方法的详细算法步骤 算法12.1（牛顿-拉格朗日方法） 步0选取x0∈Rm,0∈R,p,Y∈(0,1)，0≤e≤1.令k=0. 步1计算7L(x,)川的值.若VL(xk,)川≤E,停算.否则， 转步2. 步2解方程组(12.6)得pk=（(dk,V) 步3若 I川VL(ck+dk,k+)川2≤(1-Y)川7L(ck,)l川2， (12.7) 则置=1，转步5；否则，转步4. Back Close

6/84 JJ II J I Back Close Têß|, K°Éè⁄Ó-.ÇKFê{. œd, ä‚⁄Ó{5ü, T ê{‰k¤‹g¬Ò5ü. e°—⁄Ó-.ÇKFê{ç[é{⁄½. é{ 12.1 (⁄Ó-.ÇKFê{) ⁄ 0 ¿ x0 ∈ R n , µ0 ∈ R l , ρ, γ ∈ (0, 1), 0 ≤ ε  1. - k := 0. ⁄ 1 Oé k∇L(xk, µk)k ä. e k∇L(xk, µk)k ≤ ε, é. ƒK, =⁄ 2. ⁄ 2 )êß| (12.6)  pk = (dk, νk). ⁄ 3 e k∇L(xk + dk, µk + νk)k 2 ≤ (1 − γ)k∇L(xk, µk)k 2 , (12.7) Kò αk := 1, =⁄ 5; ƒK, =⁄ 4

步4令mk是使下面的不等式成立的最小非负整数m: VL(ck +pmdk,uk+pmvk2<(1-Ypm)L(xk,uk)2,(12.8) 置ak=pmk. 步5令xk+1=xk十Qkdk,k+1=k十QVk.置k=k十1，转 步1. S12.1.2牛顿-拉格朗日法的Matlab程序 本小节通过一个具体的例子来介绍算法12.1（牛顿-拉格朗日方 法)的Matlab实现. Back Close

7/84 JJ II J I Back Close ⁄ 4 - mk ¥¶e°ÿ™§·ÅöKÍ m : k∇L(xk + ρ mdk, µk + ρ mνk)k 2 ≤ (1 − γρm)k∇L(xk, µk)k 2 , (12.8) ò αk = ρ mk . ⁄ 5 - xk+1 = xk + αkdk, µk+1 = µk + αkνk. ò k := k + 1, = ⁄ 1. §12.1.2 ⁄Ó-.ÇKF{ Matlab ßS !œLòá‰N~f50é{ 12.1 (⁄Ó-.ÇKFê {)  Matlab ¢y.

例12.1用算法12.1编程计算下列最优化问题的极小点 min f(c)=ei2243425- ++1只 1 s.t. 晚+吃+喝+x+喝-10=0， T2C3-5x4x5=0, x3+x号+1=0. 该问题有最优解x*=(-1.71,1.59,1.82,-0.763,-0.763)T,最优值f(x*)= 0.0539. 解编制Matlab程序如下 function [x,mu,val,mh,k]=newtlagr(x0,mu0) %功能：用牛顿-拉格朗日法求解约束优化问题： %minf(x),s.t.h_i(x)=0,=1,.,1. %输入：x0是初始点，mu0是乘子向量的初始值 %输出：x,m山分别是近似最优点及相应的乘子， %val是最优值，mh是约束函数的模，k是迭代次数. Back maxk=500;%最大迭代次数 Close

8/84 JJ II J I Back Close ~ 12.1 ^é{ 12.1 ?ßOéeÅ`zØK4:    min f(x) = ex1x2x3x4x5 − 1 2 (x 3 1 + x 3 2 + 1)2 , s.t. x2 1 + x 2 2 + x 2 3 + x 2 4 + x 2 5 − 10 = 0, x2x3 − 5x4x5 = 0, x 3 1 + x 3 2 + 1 = 0. TØKkÅ`) x ∗ = (−1.71, 1.59, 1.82, −0.763, −0.763)T , Å`ä f(x ∗ ) = 0.0539. ) ?õ Matlab ßSXe function [x,mu,val,mh,k]=newtlagr(x0,mu0) % ıU:^⁄Ó-.ÇKF{¶)Â`zØK: % min f(x), s.t. h_i(x)=0, i=1,., l. %—\:x0¥–©:, mu0¥¶fï˛–©ä %——: x,mu©O¥CqÅ`:9ÉA¶f, %val¥Å`ä,mh¥ºÍ,k¥SìgÍ. maxk=500; %ÅåSìgÍ

n=length(x0);1=length(mu0); rho=0.5;gamma=0.4; x=x0;mu=mu0; k=0;epsilon=1e-8; while(k<maxk) dl=d1a(x,mu);%计算乘子函数的梯度 if(norm(dl)<epsilon),break;end%检验终止准则 N=N1(x,mu);%计算拉格朗日矩阵 dz=-N\dl; %解方程组得搜索方向 dx=dz(1:n);du=dz(n+1:n+l) m=0;mk=0; while(m<20) %Armijo搜索 if (norm(dla(x+rho m*dx,mu+rho m*du))2<=(1-gamma*rho"m)*norm(dl)2) mk=m;break; end m=m+1; end x=x+rho'mk*dx;mu=mu+rho mk*du; k=k+1; end Back Close

9/84 JJ II J I Back Close n=length(x0); l=length(mu0); rho=0.5; gamma=0.4; x=x0; mu=mu0; k=0; epsilon=1e-8; while(k<maxk) dl=dla(x,mu); %Oé¶fºÍF› if(norm(dl)<epsilon), break; end %u™éOK N=N1(x,mu); % Oé.ÇKF› dz=-N\dl; %)êß||¢êï dx=dz(1:n); du=dz(n+1:n+l); m=0; mk=0; while(m<20) % Armijo|¢ if(norm(dla(x+rho^m*dx,mu+rho^m*du))^2<=(1-gamma*rho^m)*norm(dl)^2) mk=m; break; end m=m+1; end x=x+rho^mk*dx; mu=mu+rho^mk*du; k=k+1; end

val=f1(x); mh=norm(h1(x),inf); %%%%%拉格朗日函数L(x,mu)%%%%%%%% function l=la(x,mu) f=f1(x);%调用目标函数文件 h=h1(x);%调用约束函数文件 1=f-mu'*h;%计算乘子函数 %%%%%拉格朗日函数的梯度%%%%%% function dl=dla(x,mu) df=df1(x);%调用目标函数梯度文件 h=h1(x);%调用约束函数文件 dh=dh1(x);%调用约束函数Jacobi:矩阵文件 dl=[df-dh'*mu;-h];%计算乘子函数梯度文件 %%%%%拉格朗日函数的Hesse阵%%%%%% function d21=d2la(x,mu) d2f=d2f1(x);%调用目标函数Hesse阵文件 [d2h1,d2h2,d2h3]=d2h(x);%调用约束函数二阶导数文件 d21=d2f-mu(1)*d2h1-mu(2)*d2h2-mu(3)*d2h3;%计算乘子函数的Hesse阵 %%%%%系数矩阵N(x,mu)%%%%%% function N=N1(x,mu)%计算拉格朗日矩阵 Back Close

10/84 JJ II J I Back Close val=f1(x); mh=norm(h1(x),inf); %%%%%%%%% .ÇKFºÍ L(x,mu) %%%%%%%%%%%%% function l=la(x,mu) f=f1(x); %N^8IºÍ©á h=h1(x); %N^ºͩá l=f-mu’*h; % Oé¶fºÍ %%%%%%%%% .ÇKFºÍF› %%%%%%%%%%%%% function dl=dla(x,mu) df=df1(x); %N^8IºÍF›©á h=h1(x); %N^ºͩá dh=dh1(x); %N^ºÍJacobi› ©á dl=[df-dh’*mu; -h]; %Oé¶fºÍF›©á %%%%%%%%% .ÇKFºÍHesse %%%%%%%%%%% function d2l=d2la(x,mu) d2f=d2f1(x); %N^8IºÍHesse ©á [d2h1,d2h2,d2h3]=d2h(x); %N^ºÍÍ©á d2l=d2f-mu(1)*d2h1-mu(2)*d2h2-mu(3)*d2h3; %Oé¶fºÍHesse %%%%%%%%% XÍ› N(x,mu) %%%%%%%%%%% function N=N1(x,mu) %Oé.ÇKF›