Accueil | Cast3M

* fichier :  mma_05.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*      Test de l'opérateur MMA : Méthode des Asymptotes Mobiles        *
*       Application à l'optimisation d'un treillis à 3 barres          *
*                                                                      *
* Référence :                                                          *
* 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 considère un treillis de 3 barres réparties entre 4 noeuds        *
* Les noeuds p2, p3 et p4 sont encastrés, le noeud p1 est soumis à     *
* 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 densités des barres x1, x2, x3 *
* Le problème d'optimisation consiste à minimiser le volume de la      *
* structure tel que les contraintes ne depassent pas une valeur seuil  *
* En choisissant :                                                     *
* - des caractéristiques géometriques/mécaniques unitaires             *
* - des longueurs 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 à la densité (SIMP)    *
* - d'éliminer la variable x3 (avec x3=x1), car la contrainte dans la  *
*   barre 3 est superflue pour ce chargement                           *
* alors, le problème d'optimisation s'écrit :                          *
*                                                                      *
*  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 problème 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 ;
 
* Paramètres du problème
p      = 3. ;
q      = 3. ;
 
* Procédures pour le calcul de la fonction objectif (f0) des fonctions
* limitations (f1 f2 f3) et de leurs dérivées partielles
* Ici, l'exposant q est utilisé pour pénaliser les contraintes et
* relaxer le problème original
* Si q < p --> problème relaxé
* Si q = p --> problème initial (singulier, avec espaces dégénérés)
 
* 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 dérivées
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 référence
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 à 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 à 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 résidu pour les conditions KKT
  res kkt2 kktinf = KKT_MMA t ;
* Bilan de l'itération
  MESS (CHAI 'FORMAT' '(F10.5)' &loop x1 /4 x2 /16 f0 /28 kkt2 /40) ;
FIN loop ;
nit = (DIME lf0) - 1 ;
MESS ;
 
* Vérification du résultat
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 itérations 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 ;