Accueil | Cast3M

* fichier :  excel4.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 de la géométrie d'un treillis de 2 barres de section et *
* d'écartement variables                                               *
*                                                                      *
*                                                                      *
*                                                                      *
*                          F                                           *
*                        /                                             *
*                       /                                              *
*                      /                                               *
*     ʌ             3 *                                                *
*     |             ⟋  ⟍                                              *
*     |           ⟋      ⟍  x1 : section des barres (cm2)             *
*  1m |         ⟋          ⟍                                          *
*     |       ⟋              ⟍                                        *
*     v   1 ⟋                  ⟍ 2                                    *
*      ----*---------------------*-----                                *
*           <-------->                                                 *
*               x2 : 1/2 distance entre les barres (m)                 *
*                                                                      *
* Les 2 barres sont rotulées en bas et soumises à une force au noeud 3 *
* Fx = 24.8kN   Fy = 198.4kN                                           *
*                                                                      *
* Le problème d'optimisation est le suivant :                          *
* - Trouver (x1 x2) pour minimiser la masse du treillis                *
* - tels que la contrainte de traction dans chaque barre doit être     *
*   inférieure à 100 N/mm2                                             *
*                                                                      *
* Au final, le problème à résoudre se ramène à :                       *
* Minimiser :   w(x1,x2) = c1 * x1 + (1 + x2**2)**0.5                  *
* avec :       s1(x1,x2) = c2 * (1 + x2**2)**0.5 * (8/x1 + 1/x1x2) < 1 *
*              s2(x1,x2) = c2 * (1 + x2**2)**0.5 * (8/x1 - 1/x1x2) < 1 *
*              0.2 cm2 < x1 < 4 cm2                                    *
*              0.1 m   < x2 < 1.6 m                                    *
*                                                                      *
* si on choisit c1 = 1 et c2 = 0.124 alors la solution est :           *
* x1 = 1.41 cm2   x2 = 0.38 m                                          *
************************************************************************
 
* Options d'affichage
OPTI 'ECHO' 0 ;
idess = FAUX ;
 
* Choix des paramètres du problème mécanique c1 et c2
c1 = 1. ;
c2 = 0.124 ;
 
 
 
************************************************************************
*                  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 w(x) et de ses
* dérivées partielles
DEBP FF tx*'TABLE' ;
  x1 = tx . 1 ;
  x2 = tx . 2 ;
  tf = TABL 'VECTEUR' ;
  ax2 = (1. + (x2**2))**0.5 ;
  tf . 0 = c1 * x1 * ax2 ;
  tf . 1 = c1 * ax2 ;
  tf . 2 = c1 * x1 * x2 / ax2 ;
FINP tf ;
 
* Procédure FS1 : calcul de la fonction contrainte s1(x1,x2) et ses
* dérivées partielles
DEBP FS1 tx*'TABLE' ;
  x1 = tx . 1 ;
  x2 = tx . 2 ;
  tc = TABL 'VECTEUR' ;
  ax2 = (1. + (x2**2))**0.5 ;
  tc . 0 = c2 * ax2 * ((8./x1) + (1./(x1*x2))) ;
  tc . 1 = c2 * ax2 * ((-8./(x1*x1)) - (1./(x1*x1*x2))) ;
  tc . 2 = c2 * ((x2 * ((8./x1) + (1./(x1*x2))) / ax2) - (ax2/(x1*x2*x2))) ;
FINP tc ;
* Procédure FS2 : calcul de la fonction contrainte s2(x1,x2) et ses
* dérivées partielles
DEBP FS2 tx*'TABLE' ;
  x1 = tx . 1 ;
  x2 = tx . 2 ;
  tc = TABL 'VECTEUR' ;
  ax2 = (1. + (x2**2))**0.5 ;
  tc . 0 = c2 * ax2 * ((8./x1) - (1./(x1*x2))) ;
  tc . 1 = c2 * ax2 * ((-8./(x1*x1)) + (1./(x1*x1*x2))) ;
  tc . 2 = c2 * ((x2 * ((8./x1) + (1./(x1*x2))) / ax2) + (ax2/(x1*x2*x2))) ;
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 x1 et x2
to . 'VX0' = TABL 'VECTEUR' ;
to . 'VX0' . 1 = 0.25 ;
to . 'VX0' . 2 = 1.5 ;
* -- valeur initiale de la fonction f(x1,x2)
to . 'VF' = FF (to . 'VX0') ;
* -- valeurs initiales des contraintes s1(x1,x2) et s2(x1,x2)
to . 'MC' = TABL ;
to . 'MC' . 1 = FS1 (to . 'VX0') ;
to . 'MC' . 2 = FS2 (to . 'VX0') ;
* -- bornes pour x1 et x2
to . 'VXMIN' = TABL 'VECTEUR' ;
to . 'VXMIN' . 1 = 0.2 ;
to . 'VXMIN' . 2 = 0.1 ;
to . 'VXMAX' = TABL 'VECTEUR' ;
to . 'VXMAX' . 1 = 4. ;
to . 'VXMAX' . 2 = 1.6 ;
* -- bornes pour les fonctions contraintes
to . 'VCMAX' = TABL 'VECTEUR' ;
to . 'VCMAX' . 1 = 1. ;
to . 'VCMAX' . 2 = 1. ;
* -- choix de la methode ('STA' ou 'MOV')
to . 'METHODE' = 'MOV' ;
*to . 'T0' = 2./3. ;
*to . 'S0' = 0.5 ;
 
* Listes pour stocker les resultats intermediaires :
* fonction objectif, infaisabilité et variables de conception
lf  = PROG (to . 'VF'  . 0) ;
g1 = to . 'MC' . 1 . 0 ;
g2 = to . 'MC' . 2 . 0 ;
li = PROG (MAXI (PROG 0. ((g1 - 1.) / 1.) ((g2 - 1.) / 1.))) ;
lx1 = PROG (to . 'VX0' . 1) ;
lx2 = PROG (to . 'VX0' . 2) ;
 
* Boucle d'optimisation
REPE bop 20 ;
* appel a l'operateur EXCE --> nouveau couple (x1,x2) optimisé
  to_new = EXCE to ;
* récuperation des nouvelles valeurs des variables x1 et x2
  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
  tg1_new = FS1 tx_new ;
  tg2_new = FS2 tx_new ;
  g1 = tg1_new . 0 ;
  g2 = tg2_new . 0 ;
* mise à jour de la table d'optimisation
  to . 'VX0'    = tx_new ;
  to . 'VF'     = tf_new ;
  to . 'MC' . 1 = tg1_new ;
  to . 'MC' . 2 = tg2_new ;
* calcul de l'infaisabilité (écart sur les fonctions contraintes)
  infi = MAXI (PROG 0. ((g1 - 1.) / 1.) ((g2 - 1.) / 1.)) ;
* 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) ;
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 x1,x2 (article de Svanberg)
fref = 1.51 ;
x1ref = 1.41 ;
x2ref = 0.38 ;
* Comparaison avec les valeurs calculées par Cast3M
tol1 = 1.E-2 ;
ierr = FAUX ;
MESS ;
MESS ' Variable de | Valeur    | Valeur    | Erreur' ;
MESS ' conception  | calculee  | attendue  | relative' ;
MESS '-----------------------------------------------' ;
x1 = tx_new . 1 ;
x2 = tx_new . 2 ;
err1 = ABS ((x1 - x1ref) / x1ref) ;
err2 = ABS ((x2 - x2ref) / x2ref) ;
mot1 = CHAI 'FORMAT' '(F6.4)' 'x1' *5 '|' *14 x1 /16 '|' *26 x1ref /28 '|' *38 err1 /40 ;
mot2 = CHAI 'FORMAT' '(F6.4)' 'x2' *5 '|' *14 x2 /16 '|' *26 x2ref /28 '|' *38 err2 /40 ;
SI (err1 > tol1) ;
  ierr = VRAI ;
  mot1 = CHAI mot1 ' <-- NOK' ;
FINSI ;
MESS mot1 ;
SI (err2 > tol1) ;
  ierr = VRAI ;
  mot2 = CHAI mot2 ' <-- NOK' ;
FINSI ;
MESS mot2 ;
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 'x1,x2' lx1 ;
  x1r   = EXTR lxref 1 ;
  evx1r = EVOL 'ROUG' 'MANU' 'Iteration' (PROG 0. nit) 'F' (PROG x1ref x1ref) ;
  evx2  = EVOL 'ORAN' 'MANU' 'Iteration' lit 'x1,x2' lx2 ;
  x2r   = EXTR lxref 2 ;
  evx2r = EVOL 'ORAN' 'MANU' 'Iteration' (PROG 0. nit) 'F' (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 x1 et x2 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 ;