共轭梯度法编程.docx

1、共轭梯度法编程P90 第二章 11（2）用大M法求解min w=2x1+x2-x3-x4s.t x1-x2+2x3-x4=2 2x1+x2-3x3+x4=6 x1+x2+x3+x4=7 xi0, i = 1, 2, 3, 4用matlab求解如下：f=2,1,-1,-1,200,200,200;a=1 -1 2 -1 1 0 0;2 1 -3 1 0 1 0;1 1 1 1 0 0 1;b=2 6 7;lb=zeros(7,1);x,fval,exitflag,output,lambda=linprog(f,a,b,lb)Optimization terminated.运行结果如下：x = 3

2、.0000 0.0000 1.0000 3.0000 0.0000 0.0000 0.0000fval = 2.0000exitflag = 1output = iterations: 7 algorithm: large-scale: interior point cgiterations: 0 message: Optimization terminated.lambda = ineqlin: 0x1 double eqlin: 3x1 double upper: 7x1 double lower: 7x1 double从上述运行结果可以得出：最优解为x= ，最小值约为f*=2。P151

3、第三章 26用共轭梯度算法求f(x) = (x1-1)2+5*(x2-x12)2的极小点，取初始点x0= 。用matlab求解如下：function mg=MG()%共轭梯度法求解习题三第26题%clc;clear;n=2;x=2 0;max_k=100;count_k=1;trace(1,1)=x(1);trace(2,1)=x(2);trace(3,1)=f_fun(x);k=0;g1=f_dfun(x);s=-g1;while count_k=max_k if k=n g0=f_dfun(x); s=-g0; k=0; else r_min=fminbnd(t) f_fun(x+t*s)

4、,-100,100); x=x+r_min*s; g0=g1; g1=f_dfun(x); if norm(g1)0.001 & k300) if k=0 p=-H0*g1; else vk=sk/(sk*yk)-(H0*yk)/(yk*H0*yk); w1=(yk*H0*yk)*vk*vk; H1=H0-(H0*yk*yk*H0)/(yk*H0*yk)+(sk*sk)/(sk*yk)+w1; p=-H1*g1; H0=H1; end x00=x0; result=Usearch1(f,x1,x2,df,x0,p); arf=result(1); x0=x0+arf*p; g0=g1; g1=

5、subs(df,x1,x2,x0(1,1),x0(2,1); p0=p; yk=g1-g0; sk=x0-x00; k=k+1; end; k x0 f0=subs(f,x1,x2,x0(1,1),x0(2,1)function result=Usearch1(f,x1,x2,df,x0,p)mu=0.001;sgma=0.99;a=0;b=inf;arf=1;pk=p;x3=x0;x4=x3+arf*pk;f1=subs(f,x1,x2,x3(1,1),x3(2,1);f2=subs(f,x1,x2,x4(1,1),x4(2,1);gk1=subs(df,x1,x2,x3(1,1),x3(2

6、,1);gk2=subs(df,x1,x2,x4(1,1),x4(2,1);while (f1-f20) if(gk2*pksgma*gk1*pk) a=arf;a=min(2*arf,(a+b)/2);x4=x3+arf*pk;f2=subs(f,x1,x2,x4(1,1),x4(2,1);gk2=subs(df,x1,x2,x4(1,1),x4(2,1); while (f1-f2=-mu*arf*gk1*pk) b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,x1,x2,x4(1,1),x4(2,1);gk2=subs(df,x1,x2,x4(1,1

7、),x4(2,1); end; else break; endend;result=arf;运行结果如下：k = 19x0 = 1.0000 1.0000f0 =4.4505e-013从上述运行结果可以得出：最优解为x= ，最小值约为f*=0。P229 第四章 10用SUMT内点法求解(1) min f(x) = x1+x2 s.t g1(x) = -x12+x20 g2(x) = x10用matlab求解如下：function sumt=SUMT(x0,e0,max_k0)%SUMT内点法global x s Mx=x0;e=e0;max_k=max_k0;trace(1,1)=x(1);t

8、race(2,1)=x(2);M=100;c=3;e_FR=10-10;max_FR=200;for k=0:max_k x=FR(x,e_FR,max_FR); trace(1,k+2)=x(1); trace(2,k+2)=x(2); if f_pfun(x)e break; end M=c*M;endxf=f_fun(x)ktracefunction f=FR(x0,e,max_k)global x s;count_k=1;k=0;x=x0;g1=f_dfun(x);s=-g1;while count_k=max_k if k=n g0=f_dfun(x); s=-g0; k=0; el

9、se r_min=fminbnd(t) f_fun(x+t*s),-100,100); x=x+r_min*s;% g0=g1; g1=f_dfun(x) if norm(g1)=0&b=0 g(1,1)=1; g(2,1)=1;elseif a=0 g(1,1)=1+2*M*(-x(1)2+x(2)*(-2*x(1); g(2,1)=1+2*M*(-x(1)2+x(2);elseif a=0&b=0&b=0 p=0;elseif a=0&b0 p=b2;elseif a=0 p=a2;else p=a2+b2;end%题8(1)运行结果如下：x = 1.0e-004 * -0.2058 -0

10、.2058f = -2.0576e-005k = 5trace = 100.0000 -0.0050 -0.0017 -0.0006 -0.0002 -0.0001 -0.00003.0000 -0.0050 -0.0017 -0.0006 -0.0002 -0.0001 -0.0000从上述运行结果可以得出：最优解为x= ，最小值约为f*=0。(2) min f(x) = x12+x22 S.t 2x1+x2-20 -x2+10用matlab求解如下：主程序与（1）相同；注意：以下三个函数体在运行计算习题四第8(2)题时使用%*function g=f_dfun(x)global M;a=-

11、2*x(1)-x(2)+2;b=x(2)-1;if a=0&b=0 g(1,1)=2*x(1); g(2,1)=2*x(2);elseif a=0 g(1,1)=2*x(1)+2*M*(-2*x(1)-x(2)+2)*(-2); g(2,1)=2*x(2)+2*M*(-2*x(1)-x(2)+2)*(-1);elseif a=0&b=0&b=0 p=0;elseif a=0&b0 p=b2;elseif a=0 p=a2;else p=a2+b2;end%题8(2)运行结果如下：x = -0.0000 1.0000f = 1.0000k = 6trace = 100.0000 -0.0000

12、-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.00003.0000 0.9901 0.9967 0.9989 0.9996 0.9999 1.0000 1.0000从上述运行结果可以得出：最优解为x= ，最小值约为f*=1。P232 第四章 25用梯度投影法求解下列线性约束优化问题min f(x) = x12+x1*x2+2*x22+2*x32+2*x2*x3+4*x1+6*x2+12*x3s.t x1 + x2 + x3 6x1x2 +2x3 2xi0 , i = 1, 2, 3取x1= , 用matlab求解如下：f=x12+x1*x2+2*x22+

13、2*x32+2*x2*x3+4*x1+6*x2+12*x3;a=1 1 1;1 1 -2;b=6;-2;l=zeros(3,1);x0=1 1 3;x,fval,exitflag,output,lambda,grad,hessian=fmincon(f,x0,a,b,l,)运行结果如下：x = 0 0 1fval = 14exitflag = 1output = iterations: 2 funcCount: 14 stepsize: 1 algorithm: medium-scale: SQP, Quasi-Newton, line-search firstorderopt: 0 cgit

14、erations: message: 1x144 charlambda = lower: 3x1 double upper: 3x1 double eqlin: 0x1 double eqnonlin: 0x1 double ineqlin: 2x1 double ineqnonlin: 0x1 doublegrad = 4.0000 8.0000 16.0000hessian = 1.1146 0.6771 0.6042 0.6771 3.3646 2.4792 0.6042 2.4792 3.4583从上述运行结果可以得出：最优解为x= ，最小值约为f*=14。P234 第四章 34（1）

15、用乘子法求解max f(x) = 10x1+4.4x22 +2x3s.t x1+4x2+5x332 x1+3x2+2x329 x32/2+x223 x12, x20, x30用matlab求解如下：function x,minf = ymh434(l)format long;syms x1 x2 x3f=10*(x1)+4.4*(x2)2+2*(x3);h=-x1-4*x2-5*x3+32,-x1-3*x2-2*x3+29,(x3)2)/2+(x2)2-3,x1-2,x2,x3;x0=2 2 0;v=1 0 0 0 0 0;M=2;alpha=2;gama=0.25;var=x1 x2 x3;

16、eps=1.0e-4;if nargin = 8 eps = 1.0e-4;endm1 = transpose(x0);m2 = inf;while l FE = 0; u=subs(h,x1,x2,x3,m1(1),m1(2),m1(3);for i=1:length(h) if (v(i)+M*u(i)=0 FE = FE +(v(i)+M*u(i)2-(v(i)2); else FE=FE+(v(i)2; endend SumF = f + (1/(2*M)*FE; m2,minf = minNT(SumF,transpose(m1),var,eps); Hm2 =subs(h,x1,x

17、2,x3,m2(1),m2(2),m2(3); Hm1 =subs(h,x1,x2,x3,m1(1),m1(2),m1(3); Hx2 = Funval(h,var,x2); Hx1 = Funval(h,var,x1); if norm(Hx2) = gama M = alpha*M; x1 = x2; else v = v - M*transpose(Hx2); x1 = x2; end endendminf = Funval(f,var,x);format short;运行结果如下：x = 19.75 0 2.45minf = -202.42maxf1 = 202.42从上述运行结果可以

18、得出：最优解为x= ，最大值约为f*=202.42, 最小值约为f*=-202.42。P235 第四章 35（2）用序列二次规划法求解min f(x)=-5x1-5x2-4x3-x1x3-6x4-5x5/(1+x5)-8x6/(1+x6)-10(1-2e-x7+e-2x7)s.t g1(x)=2x4+x5+0.8x6+x7-5=0 g2(x)=x22+x32+x52+x62-5=0 g3(x)=x1+x2+x3+x4+x5+x6+x710 g4(x)=x1+x2+x3+x45 g5(x)=x1+x3+x5+x62-x72-50 xi0, i=1,2,7用matlab求解如下：function

19、f=objfun35(x)f=-5*x(1)-5*x(2)-4*x(3)-x(1)*x(3)-6*x(4)-5*x(5)/(1+x(5)-8*x(6)/(1+x(6)-10*(1-2*exp(-x(7)+exp(-2*x(7);function c,ceq = confun35(x)ceq=x(2)2+x(3)2+x(5)2+x(6)2-5;c=x(1)+x(3)+x(5)+x(6)2-x(7)2-5;clcclearx0=1,1,1,1,1,1,1;aeq=0,0,0,2,1,0.8,1; beq=5; a=1,1,1,1,1,1,1;1,1,1,1,0,0,0;b=10,5; lb=0,0

20、,0,0,0,0,0; ub=;x,fval,exitflag,output,lambda,grad,hessian=fmincon(objfun35,x0,a,b,aeq,beq,lb,ub,confun35,options)c,ceq=confun35(x)运行结果如下：Iter F-count f(x) constraint Step-size derivative optimality Procedure 0 8 -31.4958 1 Infeasible start point 1 17 -41.7175 0.8642 1 -6.3 1.62 2 26 -43.0309 0.5358

21、 1 -1.09 0.924 3 35 -44.5199 0.8348 1 -1.4 0.438 4 44 -44.4543 0.04621 1 0.126 0.247 5 53 -44.4636 0.004354 1 -0.00456 0.127 6 62 -44.4697 0.005324 1 -0.0007 0.00717 7 71 -44.4687 1.631e-005 1 0.000984 0.000523 8 80 -44.4687 6.735e-009 1 3.03e-006 7.25e-006 9 89 -44.4687 4.308e-013 1 1.25e-009 4.11e

22、-007 x = 3.2418 0.0000 1.6342 0.1240 0.8896 1.2402 2.8702fval = -44.4687hessian = 0.5738 0.3772 -0.2369 0.1036 -0.1250 -0.0757 -0.1739 0.3772 0.6898 0.4769 0.0641 -0.0816 -0.0879 -0.0803 -0.2369 0.4769 1.4590 0.0628 0.0657 0.1229 -0.0055 0.1036 0.0641 0.0628 0.8811 0.0721 0.0120 0.0699 -0.1250 -0.0816 0.0657 0.0721 1.7848 -0.2129 -0.0244 -0.0757 -0.0879 0.1229 0.0120 -0.2129 1.4472 -0.1761 -0.1739 -0.0803 -0.0055 0.0699 -0.0244 -0.1761 1.0268c = -5.9342ceq = 4.3077e-013从上述运行结果可以得出：最优解为x= ，最小值约为 f*=-44.4687。