Accueil | Cast3M

* fichier :  mma_04.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*      Test de l'opérateur MMA : Méthode des Asymptotes Mobiles        *
*    Application a une fonction a 2 variables avec 5 contraintes       *
*                                                                      *
* Le problème s'écrit :                                                *
*    Minimiser :  f0(x) = (x1-8)² + (x2-9)²                            *
*      sur x                                                           *
*         avec :  f1(x) = ((x1-7)/5)² + ((x2-3)/2)² >= 1  (ellipse)    *
*                 f2(x) = (  x1/4  )² + ((x2-8)/5)² >= 1  (ellipse)    *
*                 f3(x) =   (x1-6)²   +   (x2-7)²   >= 1  (cercle)     *
*                 f4(x) =   (x1-7)²   +  (x2-10)²   >= 1  (cercle)     *
*                 f5(x) =  (x1-10)²   +   (x2-7)²   >= 1  (cercle)     *
*                 0 <= xi <= 10   pour i=1,2                           *
*                                                                      *
* La solution (evidente) est : x1 = 8   x2 = 9                         *
************************************************************************
 
* Options
OPTI 'ECHO' 0 ;
itrac = FAUX ;
 
* Procedures pour calculer f0 ainsi que les derivees partielles df0/dx
* Juste les valeurs des fonctions
DEBP F012345 x1 x2 ;
  f0     = ((x1 - 8.) ** 2) + ((x2 - 9.) ** 2) ;
  f1     = (-1. * ((((x1 -  7.) / 5.) ** 2) + (((x2 -  3.) / 2.) ** 2))) + 1. ;
  f2     = (-1. * ((((x1 -  0.) / 4.) ** 2) + (((x2 -  8.) / 5.) ** 2))) + 1. ;
  f3     = (-1. * ((((x1 -  6.) / 1.) ** 2) + (((x2 -  7.) / 1.) ** 2))) + 1. ;
  f4     = (-1. * ((((x1 -  7.) / 1.) ** 2) + (((x2 - 10.) / 1.) ** 2))) + 1. ;
  f5     = (-1. * ((((x1 - 10.) / 1.) ** 2) + (((x2 -  7.) / 1.) ** 2))) + 1. ;
FINP f0 f1 f2 f3 f4 f5 ;
 
* Les fonctions et leurs derivees
DEBP FONC lx*'LISTREEL' ;
  x1 = EXTR lx 1 ;
  x2 = EXTR lx 2 ;
  f0 f1 f2 f3 f4 f5 = F012345 x1 x2 ;
  f = PROG f1 f2 f3 f4 f5 ;
  df0dx = 2. * (lx - (PROG 8. 9.)) ;
  dfdx = TABL ;
  dfdx . 1 = -2. * (lx - (PROG  7.  3.)) / ((PROG 5. 2.) ** 2) ;
  dfdx . 2 = -2. * (lx - (PROG  0.  8.)) / ((PROG 4. 5.) ** 2) ;
  dfdx . 3 = -2. * (lx - (PROG  6.  7.)) / ((PROG 1. 1.) ** 2) ;
  dfdx . 4 = -2. * (lx - (PROG  7. 10.)) / ((PROG 1. 1.) ** 2) ;
  dfdx . 5 = -2. * (lx - (PROG 10.  7.)) / ((PROG 1. 1.) ** 2) ;
FINP f0 df0dx f dfdx ;
 
* Solution de reference
x1ref = 8. ;
x2ref = 9. ;
xref  = PROG x1ref x2ref ;
f0ref df0dxref fval0 dfdxref = FONC xref ;
MESS ;
MESS 'Solution de reference' ;
MESS '     x1          x2          f0' ;
MESS (CHAI 'FORMAT' '(F10.5)' x1ref /4 x2ref /16 f0ref /28) ;
 
* Choix d'un point initial pour les inconnues
x0 = PROG 1. 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'  = 0. ;
t . 'XMAX'  = 10. ;
t . 'F0VAL' = f0 ;
t . 'DF0DX' = df0dx ;
t . 'FVAL'  = fval ;
t . 'DFDX'  = dfdx ;
t . 'A0'    = 1. ;
t . 'A'     = PROG (DIME dfdx)*0. ;
t . 'C'     = PROG (DIME dfdx)*1.E5 ;
t . 'D'     = PROG (DIME dfdx)*1. ;
t . 'MOVE'  = 0.01 ;
 
* Iterations de la MMA
x01 = EXTR x0 1 ;
x02 = EXTR x0 2 ;
lx1 = PROG x01 ;
lx2 = PROG x02 ;
lit = PROG 0. ;
lf0 = PROG f0 ;
li  = PROG (MAXI (fval ET 0.)) ;
MESS ;
MESS 'Optimisation par MMA' ;
MESS 'It   x1          x2          f0         kktnorm' ;
MESS (CHAI 'FORMAT' '(F10.5)' 0 x01 /4 x02 /16 f0 /28) ;
* Boucle d'optimisation
REPE loop 150 ;
* Appel a MMA
  lit = lit ET &loop ;
  MMA t ;
  xmma = t . 'X' ;
  x1 = EXTR xmma 1 ;
  x2 = EXTR xmma 2 ;
  lx1 = lx1 ET x1 ;
  lx2 = lx2 ET x2 ;
* Mise a jour des valeurs des fonctions
  f0 df0dx fval dfdx = FONC xmma ;
  lf0  = lf0 ET f0 ;
  li   = li ET (MAXI (fval ET 0.)) ;
  t . 'F0VAL'  = f0 ;
  t . 'DF0DX'  = df0dx ;
  t . 'FVAL'   = fval ;
  t . 'DFDX'   = dfdx ;
* Calcul du residu pour les conditions KKT
  res kkt2 kktinf = KKT_MMA t ;
* Bilan de l'iteration
  MESS (CHAI 'FORMAT' '(F10.5)' &loop x1 /4 x2 /16 f0 /28 kkt2 /40) ;
FIN loop ;
nit = (DIME lf0) - 1 ;
MESS ;
 
* Verification du resultat
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 iterations d'optimisation
SI itrac ;
  OPTI 'DIME' 2 'ELEM' 'QUA8' ;
  msh   = (DROI 50 ((t . 'XMIN') (t . 'XMIN')) ((t . 'XMAX') (t . 'XMIN'))) TRAN 50 (0. ((t . 'XMAX') - (t . 'XMIN'))) ;
  cmsh  = CONT msh ;
  x y   = COOR msh ;
  f0msh f1msh f2msh f3msh f4msh f5msh = F012345 x y ;
  lf1 toto = @ISOSURF msh (PROG 0.) f1msh ;
  lf2 toto = @ISOSURF msh (PROG 0.) f2msh ;
  lf3 toto = @ISOSURF msh (PROG 0.) f3msh ;
  lf4 toto = @ISOSURF msh (PROG 0.) f4msh ;
  lf5 toto = @ISOSURF msh (PROG 0.) f5msh ;
  path  = QUEL 'SEG2' lx1 lx2 ;
  pdep  = x01 x02 ;
  pfin  = x1 x2 ;
  cdep  = (CERC 10 'ROTA' 360. (pdep PLUS (0.1 0.)) pdep) COUL 'VIOL' ;
  cfin  = (CERC 10 'ROTA' 360. (pfin PLUS (0.1 0.)) pfin) COUL 'ROUG' ;
  annd  = ANNO 'ETIQ' pdep 'VIOL' 'NO' 1. VRAI 'Depart' ; 
  annf  = ANNO 'ETIQ' pfin 'ROUG' 'NE' 1. VRAI 'Arrivee' ; 
  TRAC f0msh msh (cmsh ET path ET lf1 ET lf2 ET lf3 ET lf4 ET lf5 ET cdep ET cfin) 25
    (annd ET annf) 'TITR' 'Isovaleurs de la 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) ;
  evx2  = EVOL 'ORAN' 'MANU' 'Iteration' lit 'x' lx2 ;
  evx2r = EVOL 'ORAN' 'MANU' 'Iteration' (PROG 0. nit) 'x' (PROG x2ref x2ref) ;
  tl . 4  = 'TIRR' ;
  tl . 'TITRE' . 1  = 'x1 calcule' ;
  tl . 'TITRE' . 2  = 'x1 ref.' ;
  tl . 'TITRE' . 3  = 'x2 calcule' ;
  tl . 'TITRE' . 4  = 'x2 ref.' ;
  DESS (evx1 ET evx1r ET evx2 ET evx2r) '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-2) ;
  ERRE 'Erreur dans le calcul d''optimisation' ;
SINON ;
  MESS 'Cas test passe avec succes !' ;
FINSI ;
 
 
FIN ;