% uza2.m
% fd, 28  janvier 2006, 19 mars 2025.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear ; format long 

fprintf  ('uza2 \n')  
fprintf  ('opmisation sous contraintes egalites \n') 
fprintf  ('gradient projete a pas optimal \n') ;
fprintf  ('fonction cout quadratique \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('02a-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);
rr=rr+dr;
end;
title ('ellipses d isovaleurs de la fonction cout') ;
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
%%%  print('02b-isovaleurs-fonction-cout.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('02c-contrainte-affine.jpg', '-djpg','-S891,630'); 
pause (.5)

fprintf  ('point initial de l algorithme du gradient   \n') ;
     x0=-2.5 ; y0=2;
%     x0=-4; y0=-2;
%      x0=.7  ; y0=.3;
plot (x0,y0 ,'*k;;')
%%%  axis ( [ a-(al0 +.5) , a+(al0+.5) , b-(be0+.5) , b+(be0+.5)] , 'square' );
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
title ('point initial de l algorithme du gradient projete a pas optimal') ;
print('02a-point-initial.jpg', '-djpg','-S891,630'); 
pause   (.5)

[xx,yy]=proj1(x0,y0,p,q) ;
zerop = (xx/p) + (yy/q) - 1
plot (xx,yy ,'*r;;', [x0 xx], [y0 yy], 'r;;') ;
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
title ('projection sur la droite des contraintes') ;
print('02b-projection.jpg', '-djpg','-S891,630'); 
pause   (.5) 


for j=1:30 ; 
gx=(xx-a)/al0**2 ; gy=(yy-b)/be0**2;

rho = ( gx**2 + gy**2) / ( (gx/al0)**2 +  (gy/be0)**2 ) ; 
			  
     x=xx-rho*gx; 
     y=yy-rho*gy; 
     plot ([x xx], [y yy], 'b;;') ; 
     [xx,yy]=proj1(x,y,p,q) ; 
     plot (xx,yy ,'b*;;', [x xx], [y yy], 'r;;') ; 
%% axis ( [ a-(al0 +.5) , a+(al0+.5) , b-(be0+.5) , b+(be0+.5)] , 'square' );
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );
end 
title ('algorithme du gradient projete a pas optimal') ;
print('02c-gradient-projete-contrainte-affine.jpg', '-djpg','-S891,630'); 
pause  (.5) 


fprintf  ('resolution analytique du probleme   \n') ;
x0=x ; y0=y ; 
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*;;')
axis ( [ -3.1, 3.1, -3.1, 3.1, ] , 'square' );


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;;')
axis ( [ xmin - eps  , xmin + eps  ,  ymin - eps  , ymin + eps ] );
title ('zoom') ;
%%% print('02g-gradient-projete-contrainte-affine.jpg', '-djpg','-S891,630'); 

dth =ellips1 (a,b,sqrt(2*jmin)*al0,sqrt(2*jmin)*be0 , 10000 );

eps=.01 ;  pause (.5)
axis ( [ xmin - eps  , xmin + eps  ,  ymin - eps  , ymin + eps ] );
%%% print('02h-gradient-projete-contrainte-affine.jpg', '-djpg','-S891,630'); 

eps=.001 ;  pause (.5)  
axis ( [ xmin - eps  , xmin + eps  ,  ymin - eps  , ymin + eps ] );
%%% print('02i-gradient-projete-contrainte-affine.jpg', '-djpg','-S891,630'); 

eps=.0001 ;  pause (.5) 
axis ( [ xmin - eps  , xmin + eps  ,  ymin - eps  , ymin + eps ] );
title ('l algorithme propose un point  maintenant tres proche du point optimal') ;
%%% print('02j-gradient-projete-contrainte-affine.jpg', '-djpg','-S891,630'); 

hold off


