Exemples de jeux de données Cast3M

* fichier :  mma_05.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*      Test de l'operateur MMA : Methode des Asymptotes Mobiles        *
*       Application a l'optimisation d'un treillis a 3 barres          *
*                                                                      *
* Reference :                                                          *
* 2008 M. Bruggi                                                       *
* On an alternative approach to stress constraints relaxation in       *
* topology optimization                                                *
* Struct Multidisc Optim, 36:125–141                                   *
* https://link.springer.com/article/10.1007/s00158-007-0203-6          *
*                                                                      *
* On considere un treillis de 3 barres reparties entre 4 noeuds        *
* Les noeuds p2, p3 et p4 sont encastres, le noeud p1 est soumis a     *
* 2 cas de chargement f1 et f2                                         *
*                                                                      *
*                              Ʌ f1 = [0 1.5]^T                        *
*                              ‖                                       *
*                              ‖                                       *
*                              ‖                                       *
*                           p1 *====> f2 = [1 0]^T                     *
*                             /|\                                      *
*                            / | \                                     *
*                           /  |  \                                    *
*                       x1 /   |   \ x3                                *
*                         /  x2|    \                                  *
*                        *     *     *                                 *
*                       p2     p3    p4                                *
*                                                                      *
* Les variables d'optimisation sont les densites des barres x1, x2, x3 *
* Le problème d'optimisation consiste a minimiser le volume de la      *
* structure tel que les contraintes ne depassent pas une valeur seuil  *
* En choisissant :                                                     *
* - des caracteristiques geometriques/mecaniques unitaires             *
* - des longueur de barre l1 = l3 = 1. et l2 = sqrt(2)                 *
* - une limite en contrainte sigy = 3. pour chaque barre               *
* - une loi puissance reliant le module d'Young a la densite (SIMP)    *
* - d'eliminer la variable x3 (avec x3=x1), car la contrainte dans la  *
*   barre 3 est superflue pour ce chargement                           *
* alors, le probleme d'optimisation s'ecrit :                          *
*                                                                      *
*  Minimiser : f0(x1,x2) = 2sqrt(2) x1 + x2                            *
*  sur x1,x2                                                           *
*                                0.75 x1^p                             *
*       avec : f1(x1,x2) = ----------------------  -  3 x1^p  <= 0     *
*                           x1^p / sqrt(2) + x2^p                      *
*                                                                      *
*                                 1.5 x2^p                             *
*              f2(x1,x2) = ----------------------  -  3 x2^p  <= 0     *
*                           x1^p / sqrt(2) + x2^p                      *
*                                                                      *
*              f3(x1,x2) = 1/sqrt(2) - 3 x1^p  <= 0                    *
*                                                                      *
*              0 < x1 <= 1        0 <= x2 <= 1                         *
*                                                                      *
* Le probleme admet un optimum global en   x1 = 0.708  x2 = 0.         *
*                   un optimum local  en   x1 = 0.617  x2 = 0.7        *
************************************************************************
 
* Options
OPTI 'ECHO' 0 ;
itrac = FAUX ;
 
* Parametres du probleme
p      = 3. ;
q      = 3. ;
 
* Procedures pour le calcul de la fonction objectif (f0) des fonctions limitations (f1 f2 f3)
* et de leurs derivees partielles
* Ici, l'exposant q est utilise pour penaliser les contraintes relaxer le probleme original
* Si q < p --> probleme relaxe
* Si q = p --> probleme initial (singulier, avec espaces degeneres)
 
* Juste les valeurs des fonctions
DEBP F0123 x1 x2 p q ;
  rac2   = 2. ** 0.5 ;
  f0     = (2. * rac2 * x1) + x2 ;
  xdeno  = ((x1 ** p) / rac2) + (x2 ** p) ;
  f1     = (0.75 * (x1 ** p) / xdeno) - (3. * (x1 ** q)) ;
  f2     = (1.5  * (x2 ** p) / xdeno) - (3. * (x2 ** q)) ;
  f3     = (1. / rac2) - (3. * (x1 ** q)) ;
FINP f0 f1 f2 f3 ;
 
* Les fonctions et leurs derivees
DEBP FONC lx*'LISTREEL' p*'FLOTTANT' q*'FLOTTANT' ;
  x1     = EXTR lx 1 ;
  x2     = EXTR lx 2 ;
  rac2   = 2. ** 0.5 ;
  f0 f1 f2 f3 = F0123 x1 x2 p q ;
  f      = PROG f1 f2 f3 ;
  df0dx  = PROG (2. * rac2) 1. ;
  xdeno  = ((x1 ** p) / rac2) + (x2 ** p) ;
  dfdx   = TABL ;
  dfdx . 1 = PROG ((0.75 * p * (x1 ** (p - 1.)) * (x2 ** p) / (xdeno ** 2.)) - (3. * q * (x1 ** (q - 1.))))
                  (-0.75 * p * (x2 ** (p - 1.)) * (x1 ** p) / (xdeno ** 2.)) ;
  dfdx . 2 = PROG ((-1.5 / rac2) * p * (x1 ** (p - 1.)) * (x2 ** p) / (xdeno ** 2.))
                  ((( 1.5 / rac2) * p * (x2 ** (p - 1.)) * (x1 ** p) / (xdeno ** 2.)) - (3. * q * (x2 ** (q - 1.)))) ;
  dfdx . 3 = PROG (-1. * 3. * q * (x1 ** (p - 1.)))
                  0. ;
FINP f0 df0dx f dfdx ;
 
* Solution de reference
x1ref = 0.708 ;
x2ref = 0. ;
xref  = PROG x1ref x2ref ;
f0ref df0dxref fval0 dfdxref = FONC xref p q ;
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 0.95 0.45 ;
 
* Calcul de f0(x0), df0dx(x0), fval(x0), dfdx(x0)
f0 df0dx fval dfdx = FONC x0 p q ;
 
* Initialisation de la table pour la MMA
t = TABL ;
t . 'X'     = x0 ;
t . 'XMIN'  = PROG 1.E-6 0. ;
t . 'XMAX'  = 1. ;
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 p q ;
  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)) ;
MESS 'Erreur 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' ;
  xmin  = MAXI (t . 'XMIN') ;
  xmax  = t . 'XMAX' ;
  msh   = (DROI 100 (xmin xmin) (xmax xmin)) TRAN 100 (0. (xmax - xmin)) ;
  cmsh  = CONT msh ;
  x y   = COOR msh ;
  f0msh f1msh f2msh f3msh = F0123 x y p q ;
  lf1 toto = @ISOSURF msh (PROG 0.) f1msh ;
  lf2 toto = @ISOSURF msh (PROG 0.) f2msh ;
  lf3 toto = @ISOSURF msh (PROG 0.) f3msh ;
  path  = QUEL 'SEG2' lx1 lx2 ;
  pdep  = x01 x02 ;
  pfin  = x1 x2 ;
  cdep  = (CERC 10 'ROTA' 360. (pdep PLUS (0.01 0.)) pdep) COUL 'VIOL' ;
  cfin  = (CERC 10 'ROTA' 360. (pfin PLUS (0.01 0.)) pfin) COUL 'ROUG' ;
  annd  = ANNO 'ETIQ' pdep 'VIOL' 'NE' 0.1 VRAI 'Depart' ; 
  annf  = ANNO 'ETIQ' pfin 'ROUG' 'NE' 0.1 VRAI 'Arrivee' ; 
  TRAC f0msh msh (cmsh ET path ET lf1 ET lf2 ET lf3 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-3) ;
  ERRE 'Erreur dans le calcul d''optimisation' ;
SINON ;
  MESS 'Cas test passe avec succes !' ;
FINSI ;
 
 
FIN ;