DFP算法及Matlab程序

作业二用DFP算法求解1212221422)(minxxxxxxf，取Tx110，10010H。一、求解：TTTgHpgxxxxxg)2,4(,)2,4(,)42,422()(00002121（1）求迭代点x1令32040)()(2000pxf，得)(0的极小值点410，所以得：.)4,3(,)5.0,1(,2,1,5.0,201001010001TTTTggyxxsgpxx于是，由DFP修正公式有4138388410010000000000001syssyHyHyyHHHTTTT下一个搜索方向为TgHp6,851111（2）求迭代点x2令5.5458)()(2111pxf，得)(1的极小值点541于是得：TTgpxx)0,0(,2,421112，所以：8,)2,4(2fxxT，因Hesse阵TGxG4222)(为正定阵，)(xf为严格凸函数，所以x为整体极小点。二、DFP算法迭代步骤如下：（1）给定初始点0x，初始矩阵0H（通常取单位阵），计算0g，令k=0，给定控制误差。（2）令kkkgHp。（3）由精确一维搜索确定步长k，)(min)(0kkkkkpxfpxf（4）令kkkkpxx1。（5）若kg，则1kxx停；否则令kkkxxs1，kkkggy1。（6）由DFP修正公式得1kH。令k=k+1，转步骤（2）三、DFP算法matlab程序实现function[best_x,best_fx,count]=DFP(x0,ess)symsx1x2t;f=x1*x1+2*x2*x2-2*x1*x2-4*x1;fx=diff(f,x1);%求表达式f对x1的一阶求导fy=diff(f,x2);%求表达式f对x2的一阶求导fi=[fxfy];%构造函数f的梯度函数%初始点的梯度和函数值g0=subs(fi,[x1x2],x0);f0=subs(f,[x1x2],x0);H0=eye(2);%输出x0,f0,g0x0f0g0xk=x0;fk=f0;gk=g0;Hk=H0;k=1;while(norm(gk)ess)%迭代终止条件||gk||=essdisp('************************************************************')disp(['第'num2str(k)'次寻优'])%确定搜索方向pk=-Hk*gk';%由步长找到下一点x(k+1)xk=xk+t*pk';f_t=subs(f,[x1x2],xk);%构造一元搜索的一元函数φ(t)%由一维搜索找到最优步长df_t=diff(f_t,t);tk=solve(df_t);iftk~=0tk=double(tk);elsebreak;end%计算下一点的函数值和梯度xk=subs(xk,t,tk)fk=subs(f,[x1x2],xk)gk0=gk;gk=subs(fi,[x1x2],xk)%DPF校正公式，找到修正矩阵yk=gk-gk0;sk=tk*pk';Hk=Hk-(Hk*yk'*yk*Hk)/(yk*Hk*yk')+sk'*sk/(yk*sk')%修正公式k=k+1;enddisp('结果如下：')best_x=xk;%最优点best_fx=fk;%最优值count=k-1;四、程序执行结果在命令窗口输入以下命令：x0=[11];ess=1e-6;[best_x,best_fx,count]=DFP(x0,ess)程序运行结果：x0=11f0=-3g0=-42************************************************************第1次寻优xk=2.00000.5000fk=-5.5000gk=-1-2Hk=0.84000.38000.38000.4100************************************************************第2次寻优xk=42fk=-8gk=00Hk=1.00000.50000.50000.5000结果如下：best_x=42best_fx=-8count=2可以看到，最优点8,)2,4(fxT，迭代次数2次，与前面结果一致。