1=length(mu); d21=d2la(x,mu);dh=dh1(x); N=[d21,-dh';-dh,zeros(1,1)]; %%%%目标函数f(x)%%%%%%%% function f=f1(x) s=x(1)*x(2)*x(3)*x(4)*x(5); f=exp(s)-0.5*(x(1)^3+x(2)3+1)2; %%%%约束函数h(x)%%%%%%% function h=h1(x) h=[x(1)2+x(2)^2+x(3)2+x(4)^2+x(5)2-10;x(2)*x(3)-5*x(4)*x(5);x(1)3+x(2)^3+1]; %%%%%目标函数f(x)的梯度%%%%%% function df=df1(x) s=x(1)*x(2)*x(3)*x(4)*x(5); df(1)=s/(x(1)*exp(s)-3*(x(1)3+x(2)3+1)*x(1)2; df(2)=s/(x(2)*exp(s)-3*(x(1)^3+x(2)*3+1)*x(2)2; df(3)=s/(x(3)*exp(s);df1(4)=s/(x(4))*exp(s); df(5)=s/(x(5)*exp(s); df=df (:) %%%%约束函数h(x)的Jacobi矩阵A(x)%%%%%% function dh=dh1(x) Back Close
11/84 JJ II J I Back Close l=length(mu); d2l=d2la(x,mu); dh=dh1(x); N=[d2l, -dh’; -dh, zeros(l,l)]; %%%%%%%%% 8IºÍ f(x) %%%%%%%%%%%%% function f=f1(x) s=x(1)*x(2)*x(3)*x(4)*x(5); f=exp(s)-0.5*(x(1)^3+x(2)^3+1)^2; %%%%%%%%% ÂºÍ h(x) %%%%%%%%%%%%% function h=h1(x) h=[x(1)^2+x(2)^2+x(3)^2+x(4)^2+x(5)^2-10;x(2)*x(3)-5*x(4)*x(5);x(1)^3+x(2)^3+1]; %%%%%%%%% 8IºÍ f(x) F›%%%%%%%%%%%%% function df=df1(x) s=x(1)*x(2)*x(3)*x(4)*x(5); df(1)=s/(x(1))*exp(s)-3*(x(1)^3+x(2)^3+1)*x(1)^2; df(2)=s/(x(2))*exp(s)-3*(x(1)^3+x(2)^3+1)*x(2)^2; df(3)=s/(x(3))*exp(s); df1(4)=s/(x(4))*exp(s); df(5)=s/(x(5))*exp(s); df=df(:); %%%%%%%%% ÂºÍ h(x) Jacobi› A(x)%%%%%%%%%%%%% function dh=dh1(x)
dh=[2*x(1),2*x(2),2*x(3),2*x(4),2*x(5);0,x(3),x(2),-5*x(5),-5*x(4);. 3*x(1)2,3*x(2)2,0,0,0]; %%%%%目标函数f(x)的Hesse阵%%%%%%% function d2f=d2f1(x) s=x(1)*x(2)*x(3)*x(4)*x(5); d2f=[(s/(x(1))~2*exp(s)-6*x(1)*(x(1)^3+x(2)^3+1)-9*x(1)^4,. (1+s)*x(3)*x(4)*x(5)*exp(s)-9*x(1)^2*x(2)^2,(1+s)*x(2)*x(4)*x(5)*exp(s),.· (1+s)*x(2)*x(3)*x(5)*exp(s),(1+s)*x(2)*x(3)*x(4)*exp(s);. (1+s)*x(3)*x(4)*x(5)*exp(s)-9*x(1)^2*x(2)~2,. (s/(x(2))2*exp(s)-6*x(2)*(x(1)^3+x(2)^3+1)-9*x(2)4,·.· (1+s)*x(1)*x(4)*x(5)*exp(s),(1+s)*x(1)*x(3)*x(5)*exp(s),. (1+s)*x(1)*x(3)*x(4)*exp(s);. (1+s)*x(2)*x(4)*x(5)*exp(s),(1+s)*x(1)*x(4)*x(5)*exp(s),s^2/(x(3)*exp(s),.: (1+s)*x(1)*x(2)*x(5)*exp(s),(1+s)*x(1)*x(2)*x(4)*exp(s);. (1+s)*x(2)*x(3)*x(5)*exp(s),(1+s)*x(1)*x(3)*x(5)*exp(s),. (1+s)*x(1)*x(2)*x(5)*exp(s),s^2/(x(4)*ep(s),(1+s)*x(1)*x(2)*x(3)*ep(s);.· (1+s)*x(2)*x(3)*x(4)*exp(s),(1+s)*x(1)*x(3)*x(4)*exp(s),. (1+s)*x(1)*x(2)*x(4)*exp(s),(1+s)*x(1)*x(2)*x(3)*exp(s),s^2/(x(5)*exp(s)]'; %%%%%约束函数h(x)的Hesse阵%%%%%% function [d2h1,d2h2,d2h3]=d2h(x) Back Close
12/84 JJ II J I Back Close dh=[2*x(1),2*x(2),2*x(3),2*x(4),2*x(5);0,x(3),x(2),-5*x(5),-5*x(4);. 3*x(1)^2,3*x(2)^2,0,0,0]; %%%%%%%%% 8IºÍ f(x) Hesse %%%%%%%%%%%%% function d2f=d2f1(x) s=x(1)*x(2)*x(3)*x(4)*x(5); d2f=[(s/(x(1)))^2*exp(s)-6*x(1)*(x(1)^3+x(2)^3+1)-9*x(1)^4,. (1+s)*x(3)*x(4)*x(5)*exp(s)-9*x(1)^2*x(2)^2,(1+s)*x(2)*x(4)*x(5)*exp(s),. (1+s)*x(2)*x(3)*x(5)*exp(s), (1+s)*x(2)*x(3)*x(4)*exp(s); . (1+s)*x(3)*x(4)*x(5)*exp(s)-9*x(1)^2*x(2)^2, . (s/(x(2)))^2*exp(s)-6*x(2)*(x(1)^3+x(2)^3+1)-9*x(2)^4,. (1+s)*x(1)*x(4)*x(5)*exp(s),(1+s)*x(1)*x(3)*x(5)*exp(s),. (1+s)*x(1)*x(3)*x(4)*exp(s);. (1+s)*x(2)*x(4)*x(5)*exp(s),(1+s)*x(1)*x(4)*x(5)*exp(s),s^2/(x(3))*exp(s),. (1+s)*x(1)*x(2)*x(5)*exp(s),(1+s)*x(1)*x(2)*x(4)*exp(s); . (1+s)*x(2)*x(3)*x(5)*exp(s),(1+s)*x(1)*x(3)*x(5)*exp(s), . (1+s)*x(1)*x(2)*x(5)*exp(s),s^2/(x(4))*exp(s),(1+s)*x(1)*x(2)*x(3)*exp(s);. (1+s)*x(2)*x(3)*x(4)*exp(s),(1+s)*x(1)*x(3)*x(4)*exp(s), . (1+s)*x(1)*x(2)*x(4)*exp(s),(1+s)*x(1)*x(2)*x(3)*exp(s),s^2/(x(5))*exp(s)]’; %%%%%%%%% ÂºÍ h(x) Hesse %%%%%%%%%%%%% function [d2h1,d2h2,d2h3]=d2h(x)
d2h1=[20000;02000;00200;00020;00002]'; d2h2=[00000;00100;01000;0000-5;000-50]'; d2h3=[6*x(1)0000;06*x(2)000;00000;00000;00000]’; 利用上面的程序，取乘子向量的初值为0=(0,0,0)，终止准则值 取为川VL(x,)川2≤10-12，对于不同的初始点得到计算结果如下表 所示 初始点(o) 迭代次数()f()的值 Ih(x)的值 (-1.71.51.8-0.6-0.6)7 11 0.0539 3.2894e-011 （-1.7,1.6,1.8.-0.7,-0.7)T 11 0.0539 5.9272e-012 (-1.8,1.7,1.9,-0.8.-0.8)T 9 0.0539 3.0606e-011 (-2,1.5,2,-1,-1)7 14 0.0539 1.4766e-011 (-3,2,3,-2,-2)T 16 0.0539 3.6545e-011 说明Matlab调用方式为：在命令窗口依次输入如下命令并回车 即得计算结果， Back x0=[-1.7,1.6,1.8,-0.7,-0.7]’; Close
13/84 JJ II J I Back Close d2h1=[2 0 0 0 0; 0 2 0 0 0; 0 0 2 0 0; 0 0 0 2 0; 0 0 0 0 2]’; d2h2=[0 0 0 0 0; 0 0 1 0 0; 0 1 0 0 0; 0 0 0 0 -5; 0 0 0 -5 0]’; d2h3=[6*x(1) 0 0 0 0;0 6*x(2) 0 0 0;0 0 0 0 0;0 0 0 0 0;0 0 0 0 0]’; |^˛°ßS, ¶fï˛–äè µ0 = (0, 0, 0), ™éOKä è k∇L(xk, µk)k 2 ≤ 10−12 , Èuÿ”–©:Oé(JXeL §´. –©: (x0) SìgÍ (k) f(xk) ä kh(xk)k ä (−1.71.51.8 − 0.6 − 0.6)T 11 0.0539 3.2894e-011 (−1.7, 1.6, 1.8, −0.7, −0.7)T 11 0.0539 5.9272e-012 (−1.8, 1.7, 1.9, −0.8, −0.8)T 9 0.0539 3.0606e-011 (−2, 1.5, 2, −1, −1)T 14 0.0539 1.4766e-011 (−3, 2, 3, −2, −2)T 16 0.0539 3.6545e-011 `² Matlab N^ê™è: 3·-Iùùg—\Xe·-ø£ê =Oé(J. x0=[-1.7,1.6,1.8,-0.7,-0.7]’;
mu0=[0.10.10.1]’; [x,mu,val,mh,k]=newtlagr(x0,mu0) §12.2SQP方法的算法模型 §12.2.1基于拉格朗日函数Hesse阵的SQP方法 前一节介绍的牛顿-拉格朗日法，由于每一迭代步求解方程组(12.6) 数值上不是很稳定，因此这一方法并不实用.但它有一个重要的贡献， 就是以它为基础发展了序列二次规划方法(SQP方法).鉴于方程组 (12.6)的求解数值不稳定，故考虑将它转化为一个严格凸二次规划问 题.转化的条件是(12.2)的解点x*处最优性二阶充分条件成立，即对 满足A(x*)Td=0的任一向量d≠0，成立 dw(x*,ud >0. Back Close
14/84 JJ II J I Back Close mu0=[0.1 0.1 0.1]’; [x,mu,val,mh,k]=newtlagr(x0,mu0) §12.2 SQP ê{é{ . §12.2.1 ƒu.ÇKFºÍ Hesse SQP ê{ cò!0 ⁄Ó-.ÇKF{, duzòSì⁄¶)êß| (12.6) Íä˛ÿ¥È½, œd˘òê{øÿ¢^. ßkòááz, “¥±ßèƒ:u– Sg5yê{ (SQP ê{). Åuêß| (12.6) ¶)Íäÿ½, ƒÚß=zèòáÓLJg5yØ K. =z^ᥠ(12.2) ): x ∗ ?Å`5 ø©^᧷, =È ˜v A(x ∗ ) T d = 0 ?òï˛ d 6= 0, §· d TW(x ∗ , µ∗ )d > 0.
这时，由引理93知，当T>0充分小时，有 W(c,u）+27AxYA 正定.考虑(12.6)中的W(xk,)用一个正定矩阵来代替，记 B)=W,)+2万APA2, 则当(xk,4)→(x*,μ)时，矩阵B(x*,u*)正定.注意到(12.6)的第一 个展开式为 W(Zk:pk)dk -A(k)Tvk=-Vf(k)+A(k)T pk; 将上式变形为 W()AG)AGd-A(- 令 4=以+4+27Ad, Back Close
15/84 JJ II J I Back Close ˘û, d⁄n 9.3 , τ > 0 ø©û, k W(x ∗ , µ∗ ) + 1 2τ A(x ∗ ) TA(x ∗ ) ½. ƒ (12.6) • W(xk, µk) ^òá½› 5ìO, P B(xk, µk) = W(xk, µk) + 1 2τ A(xk) TA(xk), K (xk, µk) → (x ∗ , µ∗ ) û, › B(x ∗ , µ∗ ) ½. 5ø (12.6) 1ò á–m™è W(xk, µk)dk − A(xk) T νk = −∇f(xk) + A(xk) Tµk, Ú˛™C/è h W(xk, µk)+ 1 2τ A(xk) TA(xk) i dk−A(xk) T h µk+νk+ 1 2τ A(xk)dk i = −∇f(xk). - µ¯k := µk + νk + 1 2τ A(xk)dk,