Accueil | Cast3M

* fichier :  excel3.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*      Test de l'opérateur EXCE : Méthode des Asymptotes Mobiles       *
*     Reproduction d'un cas test de l'article originel de Svanberg     *
*                                                                      *
*  Krister Svanberg: The method of moving asymptotes - a new method    *
*  for structural optimization                                         *
*  International Journal for Numerical Methods in Engineering,         *
*  vol. 24, 359-373 (1987)                                             *
*  https://doi.org/10.1002/nme.1620240207                              *
*                                                                      *
*                                                                      *
*                                                                      *
* Optimisation d'une poutre en flexion de section variable             *
*                                                                      *
* |‾‾‾‾‾‾‾‾‾‾‾‾|____________                                           *
* |            |            |____________                              *
* |            |            |            |____________                 *
* |            |            |            |            |____________    *
* |            |            |            |            |            |   *
* |     1      |     2      |     3      |     4      |     5      |   *
* |            |            |            |            |‾‾‾‾‾‾‾‾‾‾‾‾|   *
* |            |            |            |‾‾‾‾‾‾‾‾‾‾‾‾             |   *
* |            |            |‾‾‾‾‾‾‾‾‾‾‾‾                          |   *
* |            |‾‾‾‾‾‾‾‾‾‾‾‾                                      \|/  *
*  ‾‾‾‾‾‾‾‾‾‾‾‾                                                    v   *
*                                                                      *
* La poutre est divisée en 5 troncons de même longueur et de sections  *
* x1 x2 x3 x4 et x5 qui seront les variables de conception.            *
* Elle est encastrée à gauche et soumise à une force verticale à       *
* droite.                                                              *
*                                                                      *
* Le problème d'optimisation est le suivant :                          *
* - Trouver x=(x1 x2 x3 x4 x5) pour minimiser la masse de la poutre    *
* - tels que la flèche soit inférieure à une limite donnée             *
*                                                                      *
* Au final, le problème à résoudre se ramène à :                       *
* Minimiser :  f(x) = c1 * (x1 + x2 + x3 + x4 + x5)                    *
* avec :       g(x) = 61/x1**3 + 37/x2**3 + 19/x3**3                   *
*                              +  7/x4**3 +  1/x5**3  < c2             *
*              xi > 0                                                  *
*                                                                      *
* si on choisit c1 = 0.0624 et c2 = 1 alors la solution est :          *
* x1 = 6.016   x2 = 5.309   x3 = 4.494   x4 = 3.502   x5 = 2.153       *
************************************************************************
 
* Options d'affichage
OPTI 'ECHO' 0 ;
idess = FAUX ;
 
* Choix des paramètres du problème mécanique c1 et c2
c1 = 0.0624 ;
c2 = 1. ;
 
 
 
************************************************************************
*                  F O N C T I O N   O B J E C T I F                   *
*                                 E T                                  *
*              F O N C T I O N S   C O N T R A I N T E S               *
************************************************************************
 
* Procédure FF : calcul de la fonction objectif f(x) et de ses dérivées
* partielles
DEBP FF tx*'TABLE' ;
  x1 = tx . 1 ;
  x2 = tx . 2 ;
  x3 = tx . 3 ;
  x4 = tx . 4 ;
  x5 = tx . 5 ;
  tf = TABL 'VECTEUR' ;
  tf . 0 = c1 * (x1 + x2 + x3 + x4 + x5) ;
  tf . 1 = c1 ;
  tf . 2 = c1 ;
  tf . 3 = c1 ;
  tf . 4 = c1 ;
  tf . 5 = c1 ;
FINP tf ;
 
* Procedure FG : calcul de la fonction contrainte g(x) et ses dérivées
* partielles
DEBP FG tx*'TABLE' ;
  x1 = tx . 1 ;
  x2 = tx . 2 ;
  x3 = tx . 3 ;
  x4 = tx . 4 ;
  x5 = tx . 5 ;
  tc = TABL 'VECTEUR' ;
  tc . 0 = (61. / (x1**3)) + (37. / (x2**3)) + (19. / (x3**3)) + (7. / (x4**3)) + (1. / (x5**3)) ;
  tc . 1 = -3. * 61. / (x1**4) ;
  tc . 2 = -3. * 37. / (x2**4) ;
  tc . 3 = -3. * 19. / (x3**4) ;
  tc . 4 = -3. *  7. / (x4**4) ;
  tc . 5 = -3. *  1. / (x5**4) ;
FINP tc ;
 
 
 
 
 
************************************************************************
*                       O P T I M I S A T I O N                        *
************************************************************************
 
* Initialisation de la table d'optimisation
to = TABL ;
* -- valeurs initales de x
to . 'VX0' = TABL 'VECTEUR' ;
to . 'VX0' . 1 = 5. ;
to . 'VX0' . 2 = 5. ;
to . 'VX0' . 3 = 5. ;
to . 'VX0' . 4 = 5. ;
to . 'VX0' . 5 = 5. ;
* -- valeur initiale de la fonction f(x)
to . 'VF' = FF (to . 'VX0') ;
* -- valeurs initiales des contraintes g(x)
to . 'MC' = TABL ;
to . 'MC' . 1 = FG (to . 'VX0') ;
* -- bornes pour les xi
to . 'VXMIN' = TABL 'VECTEUR' ;
to . 'VXMIN' . 1 = 0. ;
to . 'VXMIN' . 2 = 0. ;
to . 'VXMIN' . 3 = 0. ;
to . 'VXMIN' . 4 = 0. ;
to . 'VXMIN' . 5 = 0. ;
to . 'VXMAX' = TABL 'VECTEUR' ;
to . 'VXMAX' . 1 = 50. ;
to . 'VXMAX' . 2 = 50. ;
to . 'VXMAX' . 3 = 50. ;
to . 'VXMAX' . 4 = 50. ;
to . 'VXMAX' . 5 = 50. ;
* -- bornes pour la fonction contrainte (c2)
to . 'VCMAX' = TABL 'VECTEUR' ;
to . 'VCMAX' . 1 = c2 ;
* -- choix de la methode ('STA' ou 'MOV')
to . 'METHODE' = 'MOV' ;
to . 'T0' = 2./3. ;
 
* Listes pour stocker les resultats intermediaires :
* fonction objectif, infaisabilité et variables de conception
lf  = PROG (to . 'VF'  . 0) ;
g1 = to . 'MC' . 1 . 0 ;
li = PROG (MAXI ((PROG 0.) ET ((g1 - c2) / c2))) ;
lx1 = PROG (to . 'VX0' . 1) ;
lx2 = PROG (to . 'VX0' . 2) ;
lx3 = PROG (to . 'VX0' . 3) ;
lx4 = PROG (to . 'VX0' . 4) ;
lx5 = PROG (to . 'VX0' . 5) ;
 
* Boucle d'optimisation
REPE bop 20 ;
* appel a l'operateur EXCE --> nouveaux xi optimisés
  to_new = EXCE to ;
* récuperation des nouvelles valeurs des variables xi
  tx_new = to_new . 'VX0' ;
* calcul de la nouvelle valeur de la fonction objectif et de ses dérivées
  tf_new = FF tx_new ;
  fi = tf_new . 0 ;
* calcul des nouvelles valeurs des fonctions contraintes et de leurs dérivées
  tg_new = FG tx_new ;
  gi = tg_new . 0 ;
* mise à jour de la table d'optimisation
  to . 'VX0'    = tx_new ;
  to . 'VF'     = tf_new ;
  to . 'MC' . 1 = tg_new ;
* calcul de l'infaisabilité (écart sur les fonctions contraintes)
  infi = MAXI (PROG 0. ((gi - c2) / c2)) ;
* remplissage des listes pour visualition des itérations
  lf  = lf  ET fi ;
  li  = li  ET infi ;
  lx1 = lx1 ET (tx_new . 1) ;
  lx2 = lx2 ET (tx_new . 2) ;
  lx3 = lx3 ET (tx_new . 3) ;
  lx4 = lx4 ET (tx_new . 4) ;
  lx5 = lx5 ET (tx_new . 5) ;
FIN bop ;
nit = (DIME lf) - 1 ;
 
 
 
 
 
************************************************************************
*                       V É R I F I C A T I O N                        *
************************************************************************
 
* Valeurs de référence attendues pour la fonction objectif et les
* variables de conception xi (article de Svanberg)
fref = 1.340 ;
lxref = PROG 6.016 5.309 4.494 3.502 2.153 ;
* Comparaison avec les valeurs calculée par Cast3M
tol1 = 1.E-3 ;
ierr = FAUX ;
MESS ;
MESS ' Variable de | Valeur    | Valeur    | Erreur' ;
MESS ' conception  | calculee  | attendue  | relative' ;
MESS '-----------------------------------------------' ;
REPE b1 5 ;
  xi = tx_new . &b1 ;
  xref = EXTR lxref &b1 ;
  erri = ABS ((xi - xref) / xref) ;
  moti = CHAI 'FORMAT' '(F6.4)' 'x' *5 &b1 *6 '|' *14 xi /16 '|' *26 xref /28 '|' *38 erri /40 ;
  SI (erri > tol1) ;
    ierr = VRAI ;
    moti = CHAI moti ' <-- NOK' ;
  FINSI ;
  MESS moti ;
FIN b1 ;
MESS ;
MESS ' Fonction    | Valeur    | Valeur    | Erreur' ;
MESS ' objectif    | calculee  | attendue  | relative' ;
MESS '-----------------------------------------------' ;
erri = ABS ((fi - fref) / fref) ;
moti = CHAI 'FORMAT' '(F6.4)' 'f' *5 '|' *14 fi /16 '|' *26 fref /28 '|' *38 erri /40 ;
SI (erri > tol1) ;
  ierr = VRAI ;
  moti = CHAI moti ' <-- NOK' ;
FINSI ;
MESS moti ;
 
 
 
 
 
************************************************************************
*                    P O S T   T R A I T E M E N T                     *
************************************************************************
 
* Évolutions temporelles de f, de l'infaisabilité et des variables
* de conception en fonction des iterations d'optimisation
SI idess ;
  lit   = PROG 0. 'PAS' 1. nit ;
  evf   = EVOL 'ROUG' 'MANU' 'Iteration' lit 'F'  lf ;
  evfr  = EVOL 'ROUG' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG fref fref) ;
  tl    = TABL ;
  tl . 2 = 'TIRR' ;
  tl . 'TITRE' = TABL ;
  tl . 'TITRE' . 1  = 'F calculee' ;
  tl . 'TITRE' . 2  = 'F ref.' ;
  DESS (evf ET evfr) 'TITR' 'Fonction objectif F(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 ;
  x1r   = EXTR lxref 1 ;
  evx1r = EVOL 'ROUG' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG x1r x1r) ;
  evx2  = EVOL 'ORAN' 'MANU' 'Iteration' lit 'x' lx2 ;
  x2r   = EXTR lxref 2 ;
  evx2r = EVOL 'ORAN' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG x2r x2r) ;
  evx3  = EVOL 'JAUN' 'MANU' 'Iteration' lit 'x' lx3 ;
  x3r   = EXTR lxref 3 ;
  evx3r = EVOL 'JAUN' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG x3r x3r) ;
  evx4  = EVOL 'VERT' 'MANU' 'Iteration' lit 'x' lx4 ;
  x4r   = EXTR lxref 4 ;
  evx4r = EVOL 'VERT' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG x4r x4r) ;
  evx5  = EVOL 'BLEU' 'MANU' 'Iteration' lit 'x' lx5 ;
  x5r   = EXTR lxref 5 ;
  evx5r = EVOL 'BLEU' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG x5r x5r) ;
  tl . 4  = 'TIRR' ;
  tl . 6  = 'TIRR' ;
  tl . 8  = 'TIRR' ;
  tl . 10 = 'TIRR' ;
  tl . 'TITRE' . 1  = 'x1 calcule' ;
  tl . 'TITRE' . 2  = 'x1 ref.' ;
  tl . 'TITRE' . 3  = 'x2 calcule' ;
  tl . 'TITRE' . 4  = 'x2 ref.' ;
  tl . 'TITRE' . 5  = 'x3 calcule' ;
  tl . 'TITRE' . 6  = 'x3 ref.' ;
  tl . 'TITRE' . 7  = 'x4 calcule' ;
  tl . 'TITRE' . 8  = 'x4 ref.' ;
  tl . 'TITRE' . 9  = 'x5 calcule' ;
  tl . 'TITRE' . 10 = 'x5 ref.' ;
  DESS (evx1 ET evx1r ET evx2 ET evx2r ET evx3 ET evx3r ET evx4 ET evx4r ET evx5 ET evx5r)
    'TITR' 'Valeurs x VS Iterations' 'LEGE' tl ;
FINSI ;
 
* Sortie en erreur si l'écart aux valeurs de références est trop important
SI ierr ;
  ERRE 'Erreur dans le calcul d''optimisation' ;
FINSI ;
 
FIN ;