% uza4.m
% fd, 28, 29  janvier 2006, 19 mars 2025, 08 juillet 2026 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear ; format long 

fprintf  ('uza4 \n')  
fprintf  ('opmisation sous contraintes egalites \n') 
fprintf  ('algorithme de Uzawa avec pas optimal \n  ') ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


a=0 , b=0 , al0=2 , be0=1 

hold off
dth=ellips (a,b,al0,be0);
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
title ('une ellipse d isovaleurs de la fonction cout') ;
print('04a-ellipse.jpg', '-djpg','-S891,630'); 
pause (.5) 
hold on 	     

dr=.1; rr=.1 ; 
for j=1:20;
al=rr*al0 ; be=rr*be0 ; 
dth=ellips (a,b,al,be);
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
rr=rr+dr;
end;
title ('ellipses d isovaleurs de la fonction cout') ;
print('04b-ellipses.jpg', '-djpg','-S891,630'); 
pause (.5)



fprintf  ('droite des contraintes de la forme x/p + y/q = 1   \n') ;
p=1 , q=1
  xmin = -3 ; xmax = 3 ;
ymin = q*(1-(xmin/p)) ; ymax = q*(1-(xmax/p))
zeromin = 1 - (xmin/p) - (ymin/q) 
zeromax = 1 - (xmax/p) - (ymax/q) 
  
  plot ([xmin,xmax],[ymin,ymax],'k;;')
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
title ('ellipse d isovaleurs de la fonction cout et convexe des contraintes') ;
print('04c-ellipses.jpg', '-djpg','-S891,630'); 
pause (.5)


fprintf  ('valeur initiale du multiplicateur de Lagrange \n') ;

%% lambda = -.95 
%% lambda = -.6
  lambda = .2 ;
lambda = 0. 
  
rho = .15
  
fprintf  ('point initial de l algorithme de Uzawa \n') ;
xx = a + lambda * al0**2 / p 
yy = b + lambda * be0**2 / q

plot (xx,yy ,'*k;;')
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
title ('premier point de l algorithme de Uzawa') ;
print('04d-point-initial.jpg', '-djpg','-S891,630'); 
pause   (.5)

fprintf  ('deux iterations de l algorithme de Uzawa a pas optimal\n') ;

for j = 1 : 2
x = xx ; y = yy ; 
g = x/p + y/q - 1 ;
%   recherche du rho optimal   
rho = max1d ( lambda , a, b, al0, be0 , p , q ) ;
lambda = lambda -  rho*g 
xx = a + lambda * al0**2 / p ; 
yy = b + lambda * be0**2 / q ; 
     plot (xx,yy ,'b*;;', [x xx], [y yy], 'k;;') ;
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );

title ('algorithme de Uzawa a pas optimal') ;
if (j==1)  print('04e-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end 
if (j==2)  print('04f-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end
%%%% if (j==3)  print('04g-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end
%%%% if (j==4)  print('04h-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end
%%%% if (j==5)  print('04i-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end
pause  (.5) ; 
end


A=[ 1/al0**2 , 0 , 1/p ; 0 , 1/be0**2 , 1/q ; 1/p , 1/q , 0 ] ;
bb= [ a/al0**2 ; b/be0**2 ; 1 ] ;
xi = A \ bb ; 
xmin = xi(1:1,1) 
ymin = xi(2:2,1)       
jmin = .5 * ( ((xmin-a)/al0)**2 + ((ymin-b)/be0)**2 )  ;



plot (xmin, ymin ,'r*;;')
title ('le point de minimum est en rouge') ;
print('04g-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  
pause  (.5) ; 


n=ellips1 (a,b,2*jmin*al0,2*jmin*be0 , 1000);

eps=.1 ;  pause  (.5) 
plot (xmin, ymin ,'r*;;')   
axis ( [ xmin - eps  , xmin + eps  ,  ymin - eps  , ymin + eps ] );
plot ([xmin xmin+eps*p  ],[ymin ymin-eps*q ],'b;;')
plot ([xmin xmin-eps*p  ],[ymin ymin+eps*q ],'b;;')
dth =ellips1 (a,b,sqrt(2*jmin)*al0,sqrt(2*jmin)*be0 , 10000 );
title ('premier zoom ; le point de minimum est en rouge') ;
print('04h-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  
pause  (.5)  



eps=.01 ; 
for j=1:2
axis ( [ xmin - eps  , xmin + eps  ,  ymin - eps  , ymin + eps ] );
if (j==1)
title ('second zoom ; le point de minimum est en rouge') ;
print('04i-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end 
if (j==2)
title ('troisieme zoom ; le point de minimum est en rouge') ;
print('04j-uzawa-contrainte-affine.jpg', '-djpg','-S891,630');  end 

eps = eps / 10 ; 
 pause  (.5)  
end 





hold off

