最优化算法-Matlab程序

1、CG程序代码function x,y = cg(A,b,x0)%算法 % r0 = A*x0-b;p0 = -r0;k = 0;r = r0;p = p0;x = x0;while r=0alpha = -r'*p/(p'*A*p);x = x+alpha*p;rold = r;r = rold+alpha*A*p;beta = r'*r/(rold'*rold);p = -r+beta*p; plot(k,norm(p),'.-');hold onk = k+1;end y.funcount = k;y.fval = x'*A*x/2-

2、b'*x;CG-FR程序代码function x,y = cg_FR(fun,dfun,x0)%CGJ算法 %error = 10A-5;f0 = feval(fun,x0);df0 = feval(dfun,x0);p0 = -df0;f = f0;df = df0;p = p0;x = x0;k = 0;while (norm(df)>error)&&(k<1000)f = feval(fun,x);%alpha,funcNk,exitflag = lines(fun,0.01,0.15,0.85,6,f,df'*p,x,p);% 用线搜索找下降

3、距离 if exitflag = -1disp('Break!');break;endx = x+alpha*p;dfold = df;df = feval(dfun,x);beta = df'*df/(dfold'*dfold);p = -df+beta*p;plot(k,norm(df),'.-');hold onk = k+1;endy.funcount = k;y.fval = feval(fun,x);y.error = norm(df);BFGS程序代码function x,output = bfgs(fun,dfun,x0,vara

4、rgin)%优化课程 BFGS算法程序 % epsi = 1.0e-6;k = 0; funcN = 0; rho = 0.01;l = 0.15; u = 0.85;x = x0;f = feval(fun,x,varargin:);funcN = funcN + 1;n = length(x0);H = eye(n);g = feval(dfun,x,varargin:); while norm(g)>epsi && k<=1000 itercon = true;d = -H*g;alpha_0 = 1;gd = g'*d;alpha,funcNk,ex

5、itflag = lines(fun,rho,l,u,alpha_0,f,gd,x,d,varargin:); funcN = funcN + funcNk;if exitflag = -1 itercon = false; restart = true; H = eye(n); gold = g;endif itercon s = alpha*d;x = x + s; f = feval(fun,x,varargin:); funcN = funcN + 1; gold = g;g = feval(dfun,x,varargin:);y = g - gold; hy = H*y; sy =

6、s'*y; yhy = y'*hy;if sy<0.2*yhytheta = 0.8*yhy/(yhy-sy); s = theta*s + (1 - theta)*hy; sy = 0.2*yhy;endv = sqrt(yhy)*(s/sy - hy/yhy);H = H + s*s'/sy - hy*hy'/yhy + v*v' endk = k + 1;end output.fval = f; output.iteration = k; output.funcount = funcN; output.gnorm = norm(g);Newt

7、on-CG 程序代码function x,output = newton_cg(fun,dfun,ddfun,x0)%优%化 7.1 算法程序 % n = length(x0);x = x0;k = 1;d(1,1:n) = zeros(1,n); %循环开始，直至满足条件结束%while (n orm(feval(dfu n,x)>=10A-5) &&(k<20000)epsi(k) = min(0.5,sqrt(norm(feval(dfun,x)*norm(feval(dfun,x); z(1,1:n) = zeros(1,n);B = feval(ddfun

8、,x);r(1,1:n) = feval(dfun,x)'d(1,1:n) = -r(1,1:n);p = zeros(1,n)'for j=1:k%寻找最优下降方向 %if d(j,1:n)*B*d(j,1:n)'<=0if j=1p = -r(k,1:n)' break;elsep = z(j,1:n)'break;endendalpha = r(j,1:n)*r(j,1:n)'/(d(j,1:n)*B*d(j,1:n)');z(j+1,1:n) = z(j,1:n) + alpha*d(j,1:n);r(j+1,1:n) =

9、r(j,1:n) + alpha*(B*d(j,1:n)')'if norm(r(j+1,1:n)<epsi(k)p = z(j+1,1:n)'break;endbeta = r(j+1,1:n)*r(j+1,1:n)'/(r(j,1:n)*r(j,1:n)');d(j+1,1:n) = -r(j+1,1:n) + beta*d(j,1:n);end %方向找寻结束 % %获得新的迭代点 %gd = -d(1,1:n)*d(1,1:n)'alpha1,funcNk,exitflag = lines(fun,0.01,0.15,0.85,1,

10、feval(fun,x),gd,x,d(1,1:n)');x = x + alpha1*p; k = k + 1;endoutput.fval = feval(fun,x);output.funcount = k;CG- Steihaug程序代码function x,output = cg_steihaug(fun,dfun,ddfun,x0)%优%化 7.2 算法程序 % %参数初始化 %n = length(x0);x = x0;k = 1;d(1,1:n) = zeros(1,n);%开始循环，直至满足条件%while (n orm(feval(dfu n,x)>=10A-

11、5) &&(k<2000)epsi(k) = min(0.5,sqrt(norm(feval(dfun,x)*norm(feval(dfun,x);z(1,1:n) = zeros(1,n);B = feval(ddfun,x);r(1,1:n) = feval(dfun,x)'d(1,1:n) = -r(1,1:n);delta(k) = 1;%p = 0,0'%开%始搜索最优下降方向%此%算法利用信赖域法寻找方向%if norm(r(1,1:n)<epsi(k)p = z(1,1:n)'break;endfor j=1:kif d(j,1

12、:n)*B*d(j,1:n)'<=0%delta(k)s,val,posdef,count,lambda = trustM(r(k,1:n)',B,delta(k);p = s;break;endalpha = r(j,1:n)*r(j,1:n)'/(d(j,1:n)*B*d(j,1:n)');z(j+1,1:n) = z(j,1:n) + alpha*d(j,1:n);if(norm(z(j+1,1:n)>=delta(k)%tao =(-4*z(j,1: n)*d(j,1: n)'+sqrt(2*z(j,1: n)*d(j,1: n)'F2-4*(z(j,1: n)*z(j,1: n)'-delta(k)A2)/(2*d(j,1: n)*d(j,1:n)');p = z(j,1:n)' + tao*d(j,1:n)'break;endr(j+1,1:n) = r(j,1:n) + alpha*(B*d(j,1:n)')'if norm(r(j+1,1:n)<epsi(k)p = z(j+1,1:n)'break 八endbefa H rji一5)*0+1 V-S7(r015)*A1 m)-)八dq+1 i m) nA+l iS