ntontrust5duowei08.m

解无约束优化问题的锥模型拟牛顿信赖域算法。不同与传统的二次模型方法。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  1.出现矩阵奇异现象,需要限制s'y>=1.0e-8;
%  2.出现一致性系数大于1对吗;应该也是对的.
%  3.出现趋于渐进线方向怎么办?.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2008年5月修改时候对四维的问题进行测试 


function x=newtontrust2   %   (fun,x0,e)
% this a algorithm for a quasi-newton type trust region method based on the conic model.
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%  10.Brown and Dennis function, n=4,m=20，n=5 m=10 可以%%%%%
x0=[25;5;-5;-1];
disp('x0=[25;5;-5;-1]');
%x0=10*x0;
%x0=10*x0;
n=length(x0);
syms x1 x2 x3 x4;
x=[x1;x2;x3;x4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  m=5,10,20 %
fun=((x1+0.2*x2-exp(0.2))^2+(x3+x4*sin(0.2)-cos(0.2))^2)^2 ...
+((x1+0.4*x2-exp(0.4))^2+(x3+x4*sin(0.4)-cos(0.4))^2)^2 ...
+((x1+0.6*x2-exp(0.6))^2+(x3+x4*sin(0.6)-cos(0.6))^2)^2 ...
+((x1+0.8*x2-exp(0.8))^2+(x3+x4*sin(0.8)-cos(0.8))^2)^2 ...
+((x1+x2-exp(1))^2+(x3+x4*sin(1)-cos(1))^2)^2 ...
+((x1+1.2*x2-exp(1.2))^2+(x3+x4*sin(1.2)-cos(1.2))^2)^2 ...   
+((x1+1.4*x2-exp(1.4))^2+(x3+x4*sin(1.4)-cos(1.4))^2)^2 ...   
+((x1+1.6*x2-exp(1.6))^2+(x3+x4*sin(1.6)-cos(1.6))^2)^2 ...   
+((x1+1.8*x2-exp(1.8))^2+(x3+x4*sin(1.8)-cos(1.8))^2)^2 ...   
+((x1+2*x2-exp(2))^2+(x3+x4*sin(2)-cos(2))^2)^2 ...
+((x1+2.2*x2-exp(2.2))^2+(x3+x4*sin(2.2)-cos(2.2))^2)^2 ...  
+((x1+2.4*x2-exp(2.4))^2+(x3+x4*sin(2.4)-cos(2.4))^2)^2 ...  
+((x1+2.6*x2-exp(2.6))^2+(x3+x4*sin(2.6)-cos(2.6))^2)^2 ...  
+((x1+2.8*x2-exp(2.8))^2+(x3+x4*sin(2.8)-cos(2.8))^2)^2 ...   
+((x1+3*x2-exp(3))^2+(x3+x4*sin(3)-cos(3))^2)^2 ...  
+((x1+3.2*x2-exp(3.2))^2+(x3+x4*sin(3.2)-cos(3.2))^2)^2 ...  
+((x1+3.4*x2-exp(3.4))^2+(x3+x4*sin(3.4)-cos(3.4))^2)^2 ...   
+((x1+3.6*x2-exp(3.6))^2+(x3+x4*sin(3.6)-cos(3.6))^2)^2 ...   
+((x1+3.8*x2-exp(3.8))^2+(x3+x4*sin(3.8)-cos(3.8))^2)^2 ...   
+((x1+4*x2-exp(4))^2+(x3+x4*sin(4)-cos(4))^2)^2 ;     
disp('fun=((x1+0.2*x2-exp(0.2))^2+(x3+x4*sin(0.2)-cos(0.2))^2)^2 ...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% 9.Extended Powell singular function m=n=4 其中 n 需要是4的整数倍 
%%%%%%%%%%%%%%%%%%%%%%%%%%%可以了
%x0=[3;-1;0;1];
%disp('x0=[3;-1;0;1]');
%x0=10*x0;
%x0=10*x0;
%n=length(x0);
%syms x1 x2 x3 x4;
%x=[x1;x2;x3;x4];
%fun=(x1+10*x2)^2 +5*(x3-x4)^2 +(x2-2*x3)*(x2-2*x3)*(x2-2*x3)^2+10*(x1-x4)*(x1-x4)*(x1-x4)^2;
%disp('(x1+10*x2)^2 +5*(x3-x4)^2 +(x2-2*x3)*(x2-2*x3)*(x2-2*x3)^2+10*(x1-x4)*(x1-x4)*(x1-x4)^2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%  4.Penalty-I function%%%%%

%x0=[1;2;3;4];
%disp('x0=[1;2;3;4]');
%x0=10*x0;
%x0=10*x0;

%n=length(x0);
%syms x1 x2 x3 x4;
%x=[x1;x2;x3;x4];
%arf=1.0e-005;
%fun=arf*(x1-1)^2+arf*(x2-1)^2+arf*(x3-1)^2+arf*(x4-1)^2+(x1^2+x2^2+x3^2+x4^2-0.25)^2;
%disp('fun=arf*(x1-1)^2+arf*(x2-1)^2+arf*(x3-1)^2+arf*(x4-1)^2+(x1^2+x2^2+x3^2+x4^2-0.25)^2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  5.Penalty-II function%%%%%
%x0=[1/2;1/2;1/2;1/2];

%disp('x0=[1/2;1/2;1/2;1/2]');
%x0=10*x0;
%x0=10*x0;
%n=length(x0);
%syms x1 x2 x3 x4;
%x=[x1;x2;x3;x4];
%arf=1.0e-005;

%fun1=(x1-0.2)*(x1-0.2);
%fun2=arf*(exp(x2/10.0)+exp(x1/10.0)-exp(2/10.0)-exp(1/10.0))*(exp(x2/10.0)+exp(x1/10.0)-exp(2/10.0)-exp(1/10.0));
%fun3=arf*(exp(x3/10.0)+exp(x2/10.0)-exp(3/10.0)-exp(2/10.0))*(exp(x3/10.0)+exp(x2/10.0)-exp(3/10.0)-exp(2/10.0));
%fun4=arf*(exp(x4/10.0)+exp(x3/10.0)-exp(4/10.0)-exp(3/10.0))*(exp(x4/10.0)+exp(x3/10.0)-exp(4/10.0)-exp(3/10.0));
%fun5=arf*(exp(x2/10.0)-exp(-1/10.0))*(exp(x2/10.0)-exp(-1/10.0));
%fun6=arf*(exp(x3/10.0)-exp(-1/10.0))*(exp(x3/10.0)-exp(-1/10.0));
%fun7=arf*(exp(x4/10.0)-exp(-1/10.0))*(exp(x4/10.0)-exp(-1/10.0));
%fun8=(4*x1*x1+3*x2*x2+2*x3*x3+x4*x4-1)*(4*x1*x1+3*x2*x2+2*x3*x3+x4*x4-1);

%fun=fun1+fun2+fun3+fun4+fun5+fun6+fun7+fun8;
%disp('fun=fun1+fun2+fun3+fun4+fun5+fun6+fun7+fun8;');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% 2.Wood function

%x0=[-3;-1;-3;-1];
%x0=10*x0;
%x0=10*x0;
%disp('x0=[-3;-1;-3;-1]');

%n=length(x0);
%syms x1 x2 x3 x4;
%x=[x1;x2;x3;x4];
%fun=10*(x2-x1^2)^2 +(1-x1)^2 +90*(x4-x3^2)^2+(1-x3)^2+10*(x2+x4-2)^2+0.1*(x2-x4)^2;
%disp('10*(x2-x1^2)^2 +(1-x1)^2 +90*(x4-x3^2)^2+(1-x3)^2+10*(x2+x4-2)^2+0.1*(x2-x4)^2^2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
gg=jacobian(fun,x);
g1=gg(1);
g2=gg(2);
g3=gg(3);
g4=gg(4);

%%%%%%%%%%%%%%%%  set a serial of parameters.
eps=1.000e-004;
yita1=0.25;
yita2=0.75;
gama1=0.25;
gama2=2.00;
rmax=10.00;
r0=1.5;
A0=eye(n,n);
a0=[0.0;0.0;0.0;0.0];
f0=funt(fun,x0);
g0=gunt(g1,g2,g3,g4,x0);

xc=x0;
rc=r0;
Ac=A0;
ac=a0;
fc=f0;
gc=g0;
k=0;
sc=0;
gd=g0;
yc=1.0;
ared=1.0;nf=0;
while norm(gd)>=eps
     if k==0
          sc=trustvale0(Ac,ac,gd,rc);%% calculate a solution of the subproblem.
          %disp(sc);
          %fc=funt(fun,xc);
          fd=funt(fun,xc+sc);nf=nf+1;
          %disp(fc);disp(fd);
          ared=fc-fd;
          %disp(ared);
          pred=-(gd'*sc)/(1-ac'*sc)-(sc'*Ac*sc)/(2*(1-ac'*sc)^2);
          %disp(pred);
          seta=ared/pred;
          %disp(seta);
          while (seta<yita1)
          rc=gama1*norm(sc);
          sc=trustvale0(Ac,ac,gd,rc);%% calculate a solution of the subproblem.
          %disp(sc);
          %fc=funt(fun,xc);
          fd=funt(fun,xc+sc);nf=nf+1;
          %disp(fc);disp(fd);
          ared=fc-fd;
          %disp(ared);
          pred=-(gd'*sc)/(1-ac'*sc)-(sc'*Ac*sc)/(2*(1-ac'*sc)^2);
          %disp(pred);
          seta=ared/pred;
          disp(seta);
          end
          if (seta>=yita2)&(norm(sc)==rc)
             rc=min(rc*gama2,rmax);
          else
             rc=rc;
          end
          %gc=gunt(g1,g2,g3,g4,xc);
          xc=xc+sc;
          gd=gunt(g1,g2,g3,g4,xc); 
          k=k+1;
       end
    %%%%%%%%%%%%% update
     if (k>=1)    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面这部分注释掉就是二次模型的算法 begin
           rop=(fc-fd)^2-(gc'*sc)*(gd'*sc);
           if rop>0
             beta=-(fc-fd + sqrt(rop))/(gc'*sc);
           else 
             beta=1.0;
           end
           ac=(beta-1.0)/(gc'*sc)*gc;
           if norm(ac)==0.0
              yc=gd-gc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 除了这一行这部分就变为二次模型的
           else
              yc=beta*gd-beta^3*gc;
           end         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%上面这部分注释掉的就是二次模型的 end
           if (sc'*yc > 1.0e-8)   
               Ac=Ac+(yc*yc')/(sc'*yc)-(Ac*sc*sc'*Ac)/(sc'*Ac*sc); 
           end
           %%%%%%%%%%%% After updating, calculate the snbproblem. 
            sc=trustvale0(Ac,ac,gd,rc);%% calculate a solution of the subproblem.
               %disp(sc);
               fc=fd;%fc=funt(fun,xc);
               fd=funt(fun,xc+sc);nf=nf+1;
               %disp(fc);disp(fd);
               ared=fc-fd;
               %disp(ared);
               pred=-(gd'*sc)/(1-ac'*sc)-(sc'*Ac*sc)/(2*(1-ac'*sc)^2);
               %disp(pred);
               seta=ared/pred     
           while (seta<yita1)%%如果一致性太差，就缩小半径再计算子问题。直到一致性满足最小即可。
               rc=gama1*norm(sc);
               sc=trustvale0(Ac,ac,gd,rc)%% calculate a solution of the subproblem.
               %disp(sc);
               %fc=funt(fun,xc);
               fd=funt(fun,xc+sc);nf=nf+1;
               %disp(fc);disp(fd);
               ared=fc-fd;
               %disp(ared);
               pred=-(gd'*sc)/(1-ac'*sc)-(sc'*Ac*sc)/(2*(1-ac'*sc)^2);
               %disp(pred);
               seta=ared/pred;
               %disp(seta);
           end
           if (seta>=yita2)&(norm(sc)==rc)%%%%如果一致性比较好且到达边界则扩大信赖域半径。
              rc=min(rc*gama2,rmax);
           else%55一致性一般则信赖域不变
              rc=rc;
           end
           gc=gd;%gc=gunt(g1,g2,g3,g4,xc);%%%记录迭代前该点的梯度
           xc=xc+sc;%%%此处一致性都满足最小要求所以迭代一步。
           gd=gunt(g1,g2,g3,g4,xc); %%%记录迭代后该点的梯度
           k=k+1;
      end 
end
iters=k;
disp('iters  n=ng');
disp(k);
disp('nf');
disp(nf);
x=xc;

disp('||g||');
disp(norm(gd));
disp('f');
disp(fd);

disp('x');
%%%%%%%%%%%%%% subfuntion

%%%%%%%%%%%% value of function.
function f1=funt(fun,x)
x1=x(1);x2=x(2);x3=x(3);x4=x(4);
f1=eval(fun);

%%%%%%%%%%%% value of gradient.
function g=gunt(g1,g2,g3,g4,x)
x1=x(1);x2=x(2);x3=x(3);x4=x(4);
g=[eval(g1);eval(g2);eval(g3);eval(g4)]; 

%%%%%%%%%%%% subproblem.
function sd=trustvale0(A,a,g,r)
% to solve a trust region subproblem for conic model
% a:  the gauge vector.
% A:  a positive definte matrix.
% g:  gradient.
% r:  trust region radius.
% bar_e0:   a small positive number. 
% e0:   a small positive number between 0 and 1.
% sk:
% s_*:
%  mnsn:
HN=1-(a'/A)*g;
if HN==0.0%55555   If the conic model is not succeed,
%error('shi bai');%555   we use the quadratic model. 
a=[0.0;0.0;0.0;0.0];
%return
end
%if norm(g)==0
 %   disp('0 is !');return
 %end
sn=-(A\g)/(HN);
dsn=1/abs(1-(a'/A)*g);
HC=g'*A*g-a'*g*g'*g;
mnsn=min(1,dsn);
if HC>0
    scp=-g'*g*g/HC;
    dscp=(g'*A*g)/HC;
   bar_e0=min(dscp,mnsn);
else
   bar_e0=mnsn; 
end
e0=0.9*bar_e0;
sd=zihan(A,e0,a,g,r);
%preds=-g'*y/(1-a'*y)-y'*A*y/(2*(1-a'*y)^2);

%% 111111111111111111111111111111111111111111111111111111111111111111111111
function sk=zihan(A,e0,a,g,r)
HN=1-(a'/A)*g;
sn=-(A\g)/HN;
if norm(sn)<=r
    sk=sn;  %disp('在牛顿点得到');
    return
end
HC=g'*A*g-a'*g*g'*g;

if HC>0
    scp=-g'*g*g/HC;
end
tokc=((norm(g))^3)/(r*HC);
if HC>0
    tok=min(1,tokc);
else
    tok=1;
end
sc=-tok*r*g/norm(g);
if tok==1
    sk=sc;return %disp('在限制最速下降点信赖域边界得到');
end
p=[(norm(sn-scp))^2,2*(sn-scp)'*scp,norm(scp)^2-r^2];
h=roots(p);%解出牛顿点和柯西点连线与信赖域的交点t1和其反向延长线与信赖域交点t2。
t1=h(1);t2=h(2);
if t1<=t2%其根是二维向量正的分量为t1，负的分量为t2。
   kt=t1;t1=t2; t2=kt;
end
st1=t1*sn+(1-t1)*scp;
st2=t2*sn+(1-t2)*scp;
%%%%%%%%%%%  if abs(a)=0,that is it is a quadratic model. 
if norm(a)==0
  sk=st1;%disp('一般信赖域')
  return
end
%%%%%%%%%%%  if a is not a zero vector,that is it is a conic model. 
if HN>0
    t=t1;
   %disp('两者同侧，在柯西方点和牛顿点连线与信赖域交点上得到');
else
    te02=(e0*a'*g*g'*g-(e0+1)*g'*A*g*HN)/((a'/A)*g*g'*A*g-a'*g*g'*g);%%workkkkkkkkkkkkkkkkkkkkkk柯西方点和牛顿点连线与负超平面交点
    vst1=g'*st1/(1-a'*st1)+st1'*A*st1/(2*(1-a'*st1)^2);
    vst2=g'*st2/(1-a'*st2)+st2'*A*st2/(2*(1-a'*st2)^2);
    if (te02<=t1)&(vst1<=vst2);%牛点与原点异侧时候在异侧的与信赖域的交点是可行点，且其函数值比在同侧与信赖域交点所得可行点函数值小，
                                    %//即st1和st2都可行，且st1点比st2点好，在e0限制后的算法一定有te02<1。
        t=t1;%disp('两者异侧，在柯西方点和牛顿点连线与信赖域交点上得到，在异侧时候的交点');
    else
        t=t2;%disp('两者异侧，在柯西方点和牛顿点反向延长线与信赖域交点上得到');
    end    
end
sk=t*sn+(1-t)*scp;