Capítulo 4: OPTIMIZACIÓN Y
RESOLUCIÓN DE SISTEMAS
DE ECUACIONES NO LINEALES

Por defecto se asume la resolución del sistema (vectores en negrita):   f(x)=0.
Índice:

  1. Ajuste de una nube de puntos por mínimos cuadrados
  2. Método iterativo multivariable
  3. Método de Newton multivariable
  4. Método de Newton modificado
  5. Método de Newton multivariable con diferencias finitas
  6. Método cuasi-Newton de Broyden
  7. Estudio de convergencia de los métodos iterativos
  8. Método de minimización del máximo descenso
  9. Método de minimización del máximo descenso con backtracking
  10. Método de minimización de gradientes conjugados
  11. Método de minimización de Newton
  12. Método de minimización de Newton con diferencias finitas
  13. Método de minimización cuasi-Newton DFP
  14. Método de minimización cuasi-Newton BFGS
  15. Estudio de convergencia de los métodos de minimización
  16. Función fminunc

Volver al índice del libro: ERRORES, OPTIMIZACIÓN Y RESOLUCIÓN NUMÉRICA DE SISTEMAS



Ajuste por mínimos cuadrados de una nube de puntos:

Se ajusta por una parábola (polinomio de grado 2):
x=[0;0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9;1.0];
y=[-0.14;1.98;3.28;6.16;7.08;7.34;7.66;9.56;9.48;9.30;11.2];
b=y; A=[x.^0,x.^1,x.^2]; X=(A'*A)\(A'*b)
plot(x,y,'r*'), hold on
xx=(0:0.02:1)'; yy=[xx.^0,xx.^1,xx.^2]*X; plot(xx,yy,'b')

Ir al inicio

Método iterativo multivariable:

Se resuelve  f(x)=0  reescribiéndolo como  x=g(x) con control de convergencia:
xx=[0.1;0.1;-0.1];JJ=[];syms x y z;var=[x;y;z];itermax=7;
g=[1/3*cos(y*z)+1/6; 1/9*sqrt(x^2+sin(z)+1.06)-0.1; ...
-1/20*exp(-x*y)-(10*pi-3)/60]; Jg=jacobian(g,var);
for i=1:itermax
jacb=norm(eval(subs(Jg,var,xx(:,i))),1); JJ=[JJ,jacb];
x1=eval(subs(g,var,xx(:,i))); xx=[xx,x1];
end; errit=sum(abs(diff(xx')')), xit=xx, JJ

Ir al inicio

Método de Newton multivariable:

xx=[0.1;0.1;-0.1];syms x y z;var=[x;y;z];itermax=7;
f=[3*x-cos(y*z)-1/2; x^2-81*(y+0.1)^2+sin(z)+1.06; ...
exp(x*y)+20*z+(10*pi-3)/3]; Jf=jacobian(f,var);
for i=1:itermax
fxant=eval(subs(f,var,xx(:,i))); jacb=eval(subs(Jf,var,xx(:,i)));
x1=xx(:,i)-jacb\fxant; xx=[xx,x1];
end; errne=sum(abs(diff(xx')')), xx

Ir al inicio

Método de Newton modificado:

Se actualiza la matriz jacobiana una vex cada cuatro iteraciones.
xx=[0.1;0.1;-0.1];syms x y z;var=[x;y;z];itermax=7;
f=[3*x-cos(y*z)-1/2; x^2-81*(y+0.1)^2+sin(z)+1.06; ...
exp(x*y)+20*z+(10*pi-3)/3]; Jf=jacobian(f,var);
for i=1:itermax
if mod(i-1,4)==0 jacb=eval(subs(Jf,var,xx(:,i)));end;
fxant=eval(subs(f,var,xx(:,i)));
x1=xx(:,i)-jacb\fxant;xx=[xx,x1];
end; errnm=sum(abs(diff(xx')')),xx

Ir al inicio


Método de Newton multivariable con aproximación de la matriz jacobiana mediante diferencias finitas:

xx=[0.1;0.1;-0.1];syms x y z;var=[x;y;z];itermax=7;
f=[3*x-cos(y*z)-1/2; x^2-81*(y+0.1)^2+sin(z)+1.06; ...
exp(x*y)+20*z+(10*pi-3)/3];
n=length(xx);E=eye(n); h=1e-7;jacb=zeros(n,n);
for i=1:itermax
fxant=eval(subs(f,var,xx(:,i)));
for j=1:n;jacb(:,j)=(eval(subs(f,var,xx(:,i)+h*E(:,j))-fxant))/h;end;
x1=xx(:,i)-jacb\fxant; xx=[xx,x1];
end; errnedf=sum(abs(diff(xx')')),xx

Ir al inicio

Método cuasi-Newton de Broyden:

xx=[0.1;0.1;-0.1];syms x y z;var=[x;y;z];itermax=7;
f=[3*x-cos(y*z)-1/2; x^2-81*(y+0.1)^2+sin(z)+1.06;...
exp(x*y)+20*z+(10*pi-3)/3]; Jf=jacobian(f,var);
jacb1=(eval(subs(Jf,var,xx(:,1))))^(-1);
fxant=eval(subs(f,var,xx(:,1)));
for i=1:itermax
x1=xx(:,i)-jacb1*fxant; xx=[xx,x1]; pp=x1-xx(:,i);
fx1=eval(subs(f,var,x1)); yy=fx1-fxant; fxant=fx1;
jacb1=jacb1+((pp-jacb1*yy)*pp'*jacb1)/(pp'*jacb1*yy);
end; errbr=sum(abs(diff(xx')')), xx

Ir al inicio

Estudio de convergencia de los métodos iterativos:

semilogy((1:7)',errit,'r'),hold on,grid,semilogy((1:7)',errit,'sr')
semilogy((1:7)',errnm,'r'),semilogy((1:7)',errnm,'or')
semilogy((1:7)',errne,'b'),semilogy((1:7)',errne,'xb')
semilogy((1:7)',errbr,'b'),semilogy((1:7)',errbr,'vb')  

Ir al inicio

Método de minimización del máximo descenso:

xx=[-1;1]; syms x y a;var=[x;y];itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);Jf=jacobian(f,var);
fmd=eval(subs(f,var,xx(:,1)));
jacb=-eval(subs(Jf,var,xx(:,1)))';vmd=-jacb;
for i=1:itermax
ff=eval(subs(f,var,(xx(:,i)+a*jacb)));
Jff=diff(ff);sol=real(double(solve(Jff)));fsol=eval(subs(ff,a,sol));
[ffun,indx]=min(fsol); aa=sol(indx);
x1=xx(:,i)+aa*jacb;xx=[xx,x1];fmd=[fmd,ffun];
jacb=-eval(subs(Jf,var,x1))';vmd=[vmd,-jacb];
end; xmd=xx,sum(abs(diff(xmd')')),fmd,sum(abs(vmd))

Para la representación de la minimización unidimensional, insertar las siguientes líneas antes del end de la línea anterior:
t=0:0.02:0.5;ffunt=double(subs(ff,a,t));
figure,plot(t,ffunt,'b'),hold on,plot(aa,ffun,'r*');

Representación gráfica 3D:
figure,ezsurfc(f,[min(xx(1,:))-0.2,max(xx(1,:))+0.2,...
min(xx(2,:))-0.2,max(xx(2,:))+0.2]),hold on,xmdt=xmd';axis equal
plot3(xmdt(:,1),xmdt(:,2),fmd,'r')
plot3(xmdt(:,1),xmdt(:,2),fmd,'r*')

 Ir al inicio

Método de minimización del máximo descenso con backtracking:

xx=[-1;1]; syms x y a; var=[x;y]; itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);Jf=jacobian(f,var);
fmdx=eval(subs(f,var,xx(:,1)));
rho=0.5; beta=0.2; jacb=-eval(subs(Jf,var,xx(:,1)))';vmdx=jacb;
for i=1:itermax
ff=eval(subs(f,var,(xx(:,i)+a*jacb)));
aa=0.5;iter=0;fant=eval(subs(ff,a,aa));gradv=jacb'*vmdx(:,i);
while(fant>fmdx(i)+aa*rho*gradv&&iter<itermax)
aa=beta*aa;fant=eval(subs(ff,a,aa)); iter=iter+1; end;iter
ffun=eval(subs(ff,a,aa));x1=xx(:,i)+aa*jacb;xx=[xx,x1];
fmdx=[fmdx,ffun];
jacb=-eval(subs(Jf,var,x1))'; vmdx=[vmdx,jacb];
end; sum(abs(diff(xx')')),xmd=xx,fmdx,sum(abs(vmdx)) 

Ir al inicio

Método de minimización de gradientes conjugados:

xx=[-1;1]; syms x y a; var=[x;y]; itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);Jf=jacobian(f,var);
fgc=eval(subs(f,var,xx(:,1)));jacb0=(eval(subs(Jf,var,xx(:,1))))';
jacb1=jacb0;vdir=-jacb0;vgc=jacb1;
for i=1:itermax
ff=eval(subs(f,var,(xx(:,i)+a*vdir)));
Jff=diff(ff);sol=real(double(solve(Jff)));fsol=eval(subs(ff,a,sol));
[ffun,indx]=min(fsol);aa=sol(indx);
x1=xx(:,i)+aa*vdir;xx=[xx,x1];fgc=[fgc,ffun];
jacb1=(eval(subs(Jf,var,x1)))';vgc=[vgc,jacb1];
vdir=-jacb1+vdir*(jacb1'*jacb1)/(jacb0'*jacb0);
end; xgc=xx,sum(abs(diff(xgc')')),fgc,sum(abs(vgc)) 

Ir al inicio

Método de minimización de Newton:

xx=[-1;1]; syms x y a; var=[x;y]; itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);Jf=(jacobian(f,var));
Hf=(jacobian(Jf,var)); fne=eval(subs(f,var,xx(:,1)));
jacb=(eval(subs(Jf,var,xx(:,1))))';vne=jacb;
for i=1:itermax
Hess=eval(subs(Hf,var,xx(:,i))); vdir=-Hess\jacb;
ff=eval(subs(f,var,eval(xx(:,i)+a*vdir)));
Jff=diff(ff);sol=real(double(solve(Jff,a)));fsol=eval(subs(ff,a,sol));
[ffun,indx]=min(fsol); aa=sol(indx);
x1=xx(:,i)+aa*vdir; xx=[xx,x1];fne=[fne,ffun];
jacb=(eval(subs(Jf,var,x1)))';vne=[vne,jacb];
end; xne=xx,sum(abs(diff(xne')')),fne,sum(abs(vne))

Ir al inicio

Método de minimización de Newton con diferencias finitas:

xx=[-1;1]; syms x y a; var=[x;y]; itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);
n=length(xx);E=eye(n);h=1e-7;jacb=zeros(n,1);Hess=zeros(n,n);
fne=eval(subs(f,var,xx(:,1)));ffun=fne;hh=h*h;fh=zeros(n,1);
for j=1:n;fh(j)=eval(subs(f,var,xx+h*E(:,j)));jacb(j)=(fh(j)-ffun)/h;
end; vne=jacb;%valor del gradiente por diferencias finitas
for i=1:itermax
for k=1:n,for j=1:n;%Hessiana por diferencias finitas:
Hess(k,j)=(eval(subs(f,var,xx(:,i)+h*E(:,k)+h*E(:,j)))...
-fh(k)-fh(j)+ffun)/hh; end;end;
vdir=-Hess\jacb; ff=eval(subs(f,var,eval(xx(:,i)+a*vdir)));
Jff=diff(ff);sol=real(double(solve(Jff,a)));
fsol=eval(subs(ff,a,sol));[ffun,indx]=min(fsol); aa=sol(indx);
x1=xx(:,i)+aa*vdir; xx=[xx,x1];fne=[fne,ffun];
for j=1:n;fh(j)=eval(subs(f,var,x1+h*E(:,j)));jacb(j)=(fh(j)-ffun)/h;
end; vne=[vne,jacb];%valor del gradiente por diferencias finitas
end; xne=xx,sum(abs(diff(xne')')),fne,sum(abs(vne))

Ir al inicio

Método de minimización cuasi-Newton DFP:

xx=[-1;1];syms x y a;var=[x;y];itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);Jf=jacobian(f,var);
Hf=jacobian(Jf,var); fcn=eval(subs(f,var,xx(:,1)));
jacb=(eval(subs(Jf,var,xx(:,1))))';vcn=[jacb];
Hess=eval(subs(Hf,var,xx(:,1)));Hess=inv(Hess);
for i=1:itermax
vdir=-Hess*jacb; vdir=vdir/norm(vdir);%se normaliza
ff=eval(subs(f,var,(xx(:,i)+a*vdir)));Jff=diff(ff);
sol=real(double(solve(Jff)));fsol=eval(subs(ff,a,sol));
[ffun,indx]=min(fsol); aa=sol(indx);
x1=xx(:,i)+aa*vdir; xx=[xx,x1]; fcn=[fcn,ffun];
jacb0=jacb;jacb=(eval(subs(Jf,var,x1)))';
vcn=[vcn,jacb];p=x1-xx(:,i);q=jacb-jacb0;
Hess=Hess-(Hess*q*q'*Hess)/(q'*Hess*q)+(p*p')/(q'*p); %DFP
end; xcn=xx,sum(abs(diff(xcn')')),fcn,sum(abs(vcn))

Ir al inicio

Método de minimización cuasi-Newton BFGS:

xx=[-1;1];syms x y a;var=[x;y];itermax=10;
f=exp(x+3*y-0.1)+exp(x-3*y-0.1)+exp(-x-0.1);Jf=jacobian(f,var);
Hf=jacobian(Jf,var); fcn=eval(subs(f,var,xx(:,1)));
jacb=(eval(subs(Jf,var,xx(:,1))))';vcn=[jacb];
Hess=eval(subs(Hf,var,xx(:,1)));Hess=inv(Hess);
for i=1:itermax
vdir=-Hess*jacb; vdir=vdir/norm(vdir);%se normaliza
ff=eval(subs(f,var,(xx(:,i)+a*vdir)));Jff=diff(ff);
sol=real(double(solve(Jff)));fsol=eval(subs(ff,a,sol));
[ffun,indx]=min(fsol); aa=sol(indx);
x1=xx(:,i)+aa*vdir; xx=[xx,x1]; fcn=[fcn,ffun];
jacb0=jacb;jacb=(eval(subs(Jf,var,x1)))';
vcn=[vcn,jacb];p=x1-xx(:,i);q=jacb-jacb0;
Hess=Hess+(1+q'*Hess*q)/(q'*p)*(p*p')/(q'*p)... %BFGS
-(p*q'*Hess+Hess*q*p')/(q'*p); %BFGS
end; xcn=xx,sum(abs(diff(xcn')')),fcn,sum(abs(vcn))

Ir al inicio

Estudio de convergencia de los métodos de minimización:

gvmd=sum(abs(vmd))'; %modulo del gradiente
gvmdx=sum(abs(vmdx))';gvgc=sum(abs(vgc))';
gvne=sum(abs(vne))';gvcn=sum(abs(vcn))';
ermd=sum(abs(diff(xmd')')); %error en coordenadas
ermdx=sum(abs(diff(xmdx')'));ergc=sum(abs(diff(xgc')'));
erne=sum(abs(diff(xne')'));ercn=sum(abs(diff(xcn')'));
efmd=abs(diff(fmd')');%error o variacion en la funcion
efmdx=abs(diff(fmdx')');efgc=abs(diff(fgc')');
efne=abs(diff(fne')');efcn=abs(diff(fcn')');
figure,semilogy(1:length(efmd),efmd,'r*')%%Graf.de.converg.en.|f|:
hold on,semilogy(1:length(efmd),efmd,'r'),grid
semilogy(1:length(efmdx),efmdx,'r+'),semilogy(1:length(efmdx),efmdx,'r')
semilogy(1:length(efgc),efgc,'mv'),semilogy(1:length(efgc),efgc,'m')
semilogy(1:length(efne),efne,'bs'),semilogy(1:length(efne),efne,'b')
semilogy(1:length(efcn),efcn,'bo'),semilogy(1:length(efcn),efcn,'b')
figure,semilogy(0:length(gvmd)-1,gvmd,'r*')%Graf.converg.en.|Grad.f|:
hold on,semilogy(0:length(gvmd)-1,gvmd,'r'),grid
semilogy(0:length(gvmdx)-1,gvmdx,'r+')
semilogy(0:length(gvmdx)-1,gvmdx,'r')
semilogy(0:length(gvgc)-1,gvgc,'mv'),semilogy(0:length(gvgc)-1,gvgc,'m')
semilogy(0:length(gvne)-1,gvne,'bs'),semilogy(0:length(gvne)-1,gvne,'b')
semilogy(0:length(gvcn)-1,gvcn,'bo'),semilogy(0:length(gvcn)-1,gvcn,'b')
figure,semilogy(0:length(ermd)-1,ermd,'r*')%Graf.converg.en.|x(k)-x(k-1)|:
hold on,semilogy(0:length(ermd)-1,ermd,'r'),grid
semilogy(0:length(ermdx)-1,ermdx,'r+')
semilogy(0:length(ermdx)-1,ermdx,'r')
semilogy(0:length(ergc)-1,ergc,'mv'),semilogy(0:length(ergc)-1,ergc,'m')
semilogy(0:length(erne)-1,erne,'bs'),semilogy(0:length(erne)-1,erne,'b')
semilogy(0:length(ercn)-1,ercn,'bo'),semilogy(0:length(ercn)-1,ercn,'b')

Ir al inicio

Función fminunc:

f=@(x) exp(x(1)+3*x(2)-0.1)+exp(x(1)-3*x(2)-0.1)+exp(-x(1)-0.1)
x=fminunc(f,[-1;1]), f(x)

Ir al inicio