a03共軛梯度法編程

1、MK7FCCK282Z3SVMK6CZBW26Q6LZ9 MKHMNXTJJQ4S83MKMDE2RESG33HSP90 第二章 11 ( 2)用大M法求解min w=2X1+X2-X3-X4 s.t x1-x2+2x3-x4=22x1+x2-3x3+x4=6 x1+x2+x3+x4=7Xi > 0, 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

2、,lambda=linprog(f,a,b,lb) Optimization terminated.運(yùn)行結(jié)果如下：3.00000.00001.00003.00000.00000.00000.0000 fval =2.0000 exitflag =1 output =iterations: 7algorithm: 'large-scale: interior point' cgiterations: 0message: 'Optimization terminated.' lambda =ineqlin: 0X1 double eqlin: 3X1 double

3、upper: 7X1 double lower: 7X1 double從上述運(yùn)行結(jié)果可以得出：最優(yōu)解為 X= ，最小值約為 f*=2 。P151第三章 26用共軛梯度算法求 f(x) = (x1-1)A2+5*(x2-x1A2)A2的極小點(diǎn)，取初始點(diǎn)xO=。26 題用 matlab 求解如下： function mg=MG() % %共軛梯度法求解習(xí)題三第 %clc; clear;n=2;x=2 O'max_k=1OO; count_k=1;trace(1,1)=x(1);trace(2,1)=x(2); trace(3,1)=f_fun(x); k=O;g1=f_dfun(x);s=

4、-g1;while count_k<=max_k if k=ng0=f_dfun(x);s=-g0;k=0;elser_min=fminbnd(t) f_fun(x+t*s),-100,100);x=x+r_min*s;g0=g1;g1=f_dfun(x);if norm(g1)<10A(-6)break;endm=( norm(g1)A2)/( norm(gO)A2);s=-g1+m*s;count_k=count_k+1;trace(1,count_k)=x(1);trace(2,count_k)=x(2);trace(3,count_k)=f_fun(x);k=k+1;end

5、 end count_k f=f_fun(x) function g=f_dfun(x)g(1,1)=20*x(1)A3-20*x(1)*x (2)+2*x(1)-2;g(2,1)=10*x (2)-10*x(1)A2;function f=f_fun(x)f=5*x(1)A4-10*x(1)A2*x (2)+x(1)A2-2*x(1)+1+5*x (2)人2; 運(yùn)行結(jié)果如下：k =59x0 =1.00001.0000f0 =0從上述運(yùn)行結(jié)果可以得出：最優(yōu)解為 x= ，最小值約為 f*=0 。P151 第三章 26用 BFGS算法求 f(x) = (x1-1)2+5*(x2-x12)2的極小點(diǎn)

6、，取初始點(diǎn)x0=。用 matlab 求解如下：syms x1 x2;f=(x1-1)A2+5*(x2-x1A2)A2;v=x1,x2; df=jacobian(f,v);df=df.' x0=2,0'g1=subs(df,x1,x2,x0(1,1),x0(2,1);k=0;H0=1,0;0,1; while(norm(g1)>0.001 & k<300)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*

7、yk*yk'*H0)/(yk'*H0*yk)+(sk*sk')/(sk'*yk)+w1;p=-H1*g1;H0=H1;endx00=x0;result=Usearch1(f,x1,x2,df,x0,p);arf=result(1);x0=x0+arf*p;g0=g1;g1=subs(df,x1,x2,x0(1,1),x0(2,1); p0=p;yk=g1-g0; sk=x0-x00;k=k+1;end;kx0 f0=subs(f,x1,x2,x0(1,1),x0(2,1)function result=Usearch1(f,x1,x2,df,x0,p)mu=0.

8、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,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(d

9、f,x1,x2,x4(1,1),x4(2,1);end;while (1>0)if(gk2'*pk<sgma*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(d f

10、,x1,x2,x4(1,1),x4(2,1);end;elsebreak;end end; result=arf; 運(yùn)行結(jié)果如下： k =19 x0 =1.00001.0000f0 = 4.4505e-013從上述運(yùn)行結(jié)果可以得出：最優(yōu)解為 x= ，最小值約為 f*=0 。P229 第四章 10用SUM-內(nèi)點(diǎn)法求解(1) min f(x) = x1+x2s.t g1(x) = -x1A2+x2> 0g2(x) = x1 > 0用 matlab 求解如下：function sumt=SUMT(x0,e0,max_k0)%SUM內(nèi)點(diǎn)法 global x s M x=x0;e=e0; m

11、ax_k=max_k0;trace(1,1)=x(1); trace(2,1)=x(2); M=100;c=3;e_FR=10A-10;max_FR=200;for k=0:max_kx=FR(x,e_FR,max_FR);trace(1,k+2)=x(1);trace(2,k+2)=x(2);if f_pfun(x)<ebreak;endM=c*M;endx f=f_fun(x) k tracefunction f=FR(xO,e,max_k)global x s;count_k=1;k=O;x=xO;g1=f_dfun(x);s=-g1;while count_k<=max_k

12、if k=ng0=f_dfun(x); s=-g0;k=0;elser_min=fminbnd(t) f_fun(x+t*s),-100,100); x=x+r_min*s;%g0=g1; g1=f_dfun(x) if norm(g1)<e break;endm=( norm(g1)A2)/( norm(gO)A2); s=-g1+m*s;count_k=count_k+1; k=k+1;endend count_k=count_k-1 f=x注意：以下三個(gè)函數(shù)體在運(yùn)行計(jì)算習(xí)題四第 8(1) 題時(shí)使用 %* function f=f_fun(x)global M;f=x(1)+x(2)

13、+M*f_pfun(x);% 題 8(1)function g=f_dfun(x)global M;a=-x(1)A2+x(2);b=x(1);if a>=O&b>=Og(1,1)=1;g(2,1)=1; elseif a<0&b>=0g(1,1)=1+2*M*(-x(1)A2+x(2)*(-2*x(1);g(2,1)=1+2*M*(-x(1)A2+x (2);elseif a>=0&b<0g(1,1)=1+2*M*x(1);g(2,1)=1;elseg(1,1)=1+2*M*(-x(1)A2+x(2)*(-2*x(1)+2*M*x(

14、1); g(2,1)=1+2*M*(-x(1)A2+x(2);end懣 8(1) function p=f_pfun(x)%罰函數(shù)a=-x(1)A2+x(2); b=x(1);if a>=0&b>=0p=0;elseif a>=0&b<0p=b2;elseif a<0&b>=0p=a2;elsep=a2+b2;end懣 8(1) 運(yùn)行結(jié)果如下：1.0e-004 *-0.2058-0.2058f = -2.0576e-005 trace =100.0000 -0.0050 -0.0017 -0.0006 -0.0002 -0.0001

15、-0.0000 3.0000 -0.0050 -0.0017 -0.0006 -0.0002 -0.0001 -0.0000從上述運(yùn)行結(jié)果可以得出：最優(yōu)解為 x= ，最小值約為 f*=0 。(2) min f(x) = x12+x22S.t 2x1+x2-2 w 0-x2+1 w 0用 matlab 求解如下： 主程序與( 1 )相同； 注意：以下三個(gè)函數(shù)體在運(yùn)行計(jì)算習(xí)題四第 8(2) 題時(shí)使用 %*function g=f_dfun(x) global M;a=-2*x(1)-x(2)+2; b=x(2)-1; if a>=0&b>=0 g(1,1)=2*x(1); g(

16、2,1)=2*x(2);elseif a<0&b>=0g(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<0g(1,1)=2*x(1);g(2,1)=2*x(2)+2*M*(x(2)-1);elseg(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)+2*M*(x(2)-1);end懣 8(2) %*functio

17、n f=f_fun(x)global M;f=x(1)A2+x(2)A2+M*1_pfun(x);% 題 8(2)function p=f_pfun(x)a=-2*x(1)-x(2)+2;b=x(2)-1;if a>=0&b>=0p=0;elseif a>=0&b<0p=b2;elseif a<0&b>=0p=a2;elsep=a2+b2;end懣 8(2) 運(yùn)行結(jié)果如下：x =-0.00001.00001.0000trace =100.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.00000.9

18、999 1.0000 1.0000-0.0000 3.0000 0.9901 0.9967 0.9989 0.9996從上述運(yùn)行結(jié)果可以得出：最優(yōu)解為 x=，最小值約為 f*=1 。P232第四章 25用梯度投影法求解下列線性約束優(yōu)化問(wèn)題min f(x) = x1A2+x1*x2+2*x2A2+2*x3A2+2*x2*x3+4*x1+6*x2+12*x3s.t x1 + x2 + x3 w 6-x1 - x2 +2x3> 2 xi > 0 , i = 1,2, 3取 x1= , 用 matlab 求解如下：f= x12+x1*x2+2*x22+2*x32+2*x2*x3+4*x1+

19、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,)運(yùn)行結(jié)果如下：fval =14 exitflag =1 output =iterations: 2funcCount: 14stepsize: 1algorithm: 'medium-scale: SQP, Quasi-Newton, line-search' firstorderopt: 0cgiterations: message: 1x14

20、4 char lambda =lower: 3x1 doubleupper: 3x1 double eqlin: 0x1 double eqnonlin: 0x1 double ineqlin: 2x1 double ineqnonlin: 0x1 double grad =4.00008.000016.0000hessian =1.1146 0.67710.60420.6771 3.3646 2.47920.6042 2.4792 3.4583從上述運(yùn)行結(jié)果可以得出：最優(yōu)解為 x= ，最小值約為 f*=14 。P234 第四章 34( 1)用乘子法求解max f(x) = 10x1+4.4x

21、22 +2x3s.t x1+4x2+5x3x1+3x2+2x3x32/2+x22x1 > 2, x2用 matlab 求解如下： function x,minf = ymh434(l) format long;syms x1 x2 x3f=10*(x1)+4.4*(x2)A2+2*(x3);h=-x1-4*x2-5*x3+32,-x1-3*x2-2*x3+29,(x3)A2)/2+(x2)A2-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; eps=1.0e-4; if nargin

22、= 8 eps = 1.0e-4;end m1 = transpose(x0); m2 = inf;while lFE = 0;u=subs(h,x1,x2,x3,m1(1),m1(2),m1(3);for i=1:length(h)if (v(i)+M*u(i)<=0FE = FE +(v(i)+M*u(i)A2-(v(i)A2);elseFE=FE+(v(i)A2;endendSumF = f + (1/(2*M)*FE;m2,minf = minNT(SumF,transpose(m1),var,eps);Hm2 =subs(h,x1,x2,x3,m2(1),m2(2),m2(3)

23、;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) < epsx = x2;break;else if Hx2/Hx1 >= gamaM = alpha*M; x1 = x2;elsev = v - M*transpose(Hx2); x1 = x2;endend endminf = Funval(f,var,x); format short;運(yùn)行結(jié)果如下：x =19.75 0 2.45 minf =-202.42maxf1 =202.4

24、2從上述運(yùn)行結(jié)果可以得出： 最優(yōu)解為 x= ，最大值約為 f*=202.42, 最小值約為 f*=-202.42 。P235 第四章 35( 2)用序列二次規(guī)劃法求解min f(x)=-5x1-5x2-4x3-x1x3-6x4-5x5/(1+x5)-8x6心+x6)-10(1-2eA-x7+eA-2x7) s.t g1(x)=2x4+x5+0.8x6+x7-5=0g2(x)=x2A2+x3A2+x5A2+x6A2-5=0g3(x)=x1+x2+x3+x4+x5+x6+x7g4(x)=x1+x2+x3+x4 w 5g5(x)=x1+x3+x5+x62-x72-5xi > 0, i=1,2,

25、7用 matlab 求解如下： function 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)A2-x(7)A2-5;clc clear x0=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,0,0,0,0,0;ub=;x,fval,exitflag,output,lambda,grad,hessian=fmincon('objfun35',x0,a,b,aeq,beq,l b,ub,'confun35',options) c,ceq=confun35(x)運(yùn)行結(jié)果如下：derivat