* fichier :  mma_00.dgibi
************************************************************************
* Section : Mathematiques Fonctions
************************************************************************

************************************************************************
*      Test de l'opérateur MMA : Méthode des Asymptotes Mobiles        *
*                       Globalement Convergente                        *
*    Application à une fonction à 1 variable et sans contraintes       *
*                                                                      *
* Le problème s'écrit :                                                *
*    Minimiser :  f0(x) = (x-50)² + 25                                 *
*      sur x                                                           *
*         avec :  1 <= x <= 100                                        *
*                                                                      *
* La solution (évidente) est : x = 50                                  *
************************************************************************

* Options
OPTI 'ECHO' 0 ;
itrac = FAUX ;

* Procédure pour calculer f0 ainsi que la derivée partielle df0/dx
DEBP FONC a*'LISTREEL' ;
  x = EXTR a 1 ;
  f = ((x - 50.) ** 2) + 25. ;
  df = 2. * (a - 50.) ;
  g = PROG 0. ;
  dgda = TABL ;
  dgda . 1 = PROG 0. ;
FINP f df g dgda ;

* Solution de référence
x1ref = 50. ;
xref  = PROG x1ref ;
f0ref df0dxref fval0 dfdxref = FONC xref ;
MESS ;
MESS 'Solution de reference' ;
MESS '        x1          f0' ;
MESS (CHAI 'FORMAT' '(F10.5)' x1ref /7 f0ref >3) ;

* Choix d'un point initial pour les inconnues
x0 = PROG 1. ;

* Calcul de f0(x0), df0dx(x0), fval(x0), dfdx(x0)
f0 df0dx fval dfdx = FONC x0 ;

* Initialisation de la table pour la MMA
t = TABL ;
t . 'X'     = x0 ;
t . 'XMIN'  = 1. ;
t . 'XMAX'  = 100. ;
t . 'F0VAL' = f0 ;
t . 'DF0DX' = df0dx ;
t . 'FVAL'  = fval ;
t . 'DFDX'  = dfdx ;
t . 'C'     = PROG 1.E7 ;
t . 'GC'    = VRAI ;

* Itérations de la MMA
x01 = EXTR x0 1 ;
lx1 = PROG x01 ;
lit = PROG 0. ;
lf0 = PROG f0 ;
f1  = EXTR fval 1 ;
li  = PROG (MAXI 0. f1) ;
MESS ;
MESS 'Optimisation par GC-MMA' ;
MESS 'It sIt  x1          f0         kktnorm' ;
MESS (CHAI 'FORMAT' '(F10.5)' 0 0 /6 x01 >1 f0 >3) ;
* Boucle d'optimisation
REPE loop 20 ;
* Appel à MMA
  lit = lit ET &loop ;
* Boucle des sous itérations de la GCMMA
  REPE sloop 15 ;
    MMA t ;
*   Si le sous problème est conservatif
*   mise à jour des valeurs des fonctions et arret des sous itérations
    SI (t . 'CONSERV') ;
      xmma = t . 'X' ;
      f0 df0dx fval dfdx = FONC xmma ;
      t . 'F0VAL'  = f0 ;
      t . 'DF0DX'  = df0dx ;
      t . 'FVAL'   = fval ;
      t . 'DFDX'   = dfdx ;
      QUIT sloop ;
*   Sinon, mise à jour des valeurs des fonctions (pas besoin des gradients)
*   dans WTABLE et nouvelle sous itération
    SINON ;
      xmma = t . 'WTABLE' . 'X' ;
      f0 df0dx fval dfdx = FONC xmma ;
      t . 'WTABLE' . 'F0VALNEW'  = f0 ;
      t . 'WTABLE' . 'FVALNEW'   = fval ;
    FINSI ;
  FIN sloop ;
* Calcul du résidu pour les conditions KKT
  res kkt2 kktinf = KKT_MMA t ;
* Bilan de l'itération
  sit = t . 'WTABLE' . 'ITER' ;
  x1 = EXTR xmma 1 ;
  lx1 = lx1 ET x1 ;
  lf0  = lf0 ET f0 ;
  f1   = EXTR fval 1 ;
  li   = li ET (MAXI 0. f1) ;
  MESS (CHAI 'FORMAT' '(F10.5)' &loop sit /6 x1 >1 f0 >3 kkt2 >2) ;
FIN loop ;
nit = (DIME lf0) - 1 ;
MESS ;

* Vérification du résultat
errmax = MAXI 'ABS' ((xmma - xref) / xref) ;
MESS 'Erreur relative max (sur x)' ;
MESS errmax ;
MESS ;

* Évolutions temporelles de f, de l'infaisabilité et des variables
* en fonction des itérations d'optimisation
SI itrac ;
  labs  = PROG (t . 'XMIN') 'PAS' 1. (t . 'XMAX') ;
  lord  = ((labs - 50.) ** 2) + 25. ;
  evf0  = EVOL        'MANU' 'x' labs 'F0' lord ;
  path  = EVOL 'ROUG' 'MANU' 'STYL' 'POIN' 'MARQ' 'ROND' 'x' lx1  'F0' lf0 ;
  DESS (evf0 ET path) 'TITR' 'Fonction objectif F0 et chemin au cours de l''optimisation' ;
  lit   = PROG 0. 'PAS' 1. nit ;
  evf0  = EVOL 'ROUG' 'MANU' 'Iteration' lit 'F0'  lf0 ;
  evfr  = EVOL 'ROUG' 'MANU' 'Iteration' (PROG 0. nit) 'F0' (PROG f0ref f0ref) ;
  tl    = TABL ;
  tl . 2 = 'TIRR' ;
  tl . 'TITRE' = TABL ;
  tl . 'TITRE' . 1  = 'F0 calculee' ;
  tl . 'TITRE' . 2  = 'F0 ref.' ;
  DESS (evf0 ET evfr) 'TITR' 'Fonction objectif F0(x) VS Iterations' 'LEGE' tl ;
  evi   = EVOL 'VERT' 'MANU' 'Iteration' lit 'Inf' li ;
  tl . 'TITRE' . 1  = 'Infaisabilite' ;
  DESS evi 'TITR' 'Infaisabilite VS Iterations' 'LEGE' tl ;
  evx1  = EVOL 'ROUG' 'MANU' 'Iteration' lit 'x' lx1 ;
  evx1r = EVOL 'ROUG' 'MANU' 'Iteration' (PROG 0. nit) 'x' (PROG x1ref x1ref) ;
  tl . 'TITRE' . 1  = 'x1 calcule' ;
  tl . 'TITRE' . 2  = 'x1 ref.' ;
  DESS (evx1 ET evx1r) 'TITR' 'Valeurs x VS Iterations' 'LEGE' tl ;
FINSI ;

* Sortie en erreur si l'écart aux valeurs de références est trop important
SI (errmax > 1.E-6) ;
  ERRE 'Erreur dans le calcul d''optimisation' ;
SINON ;
  MESS 'Cas test passe avec succes !' ;
FINSI ;


FIN ;