Procédures | Cast3M

* OPTIMISE  PROCEDUR  SP204843  26/05/07    21:15:05     12544          
'DEBP' OPTIMISE TAB1*'TABLE' ;
 
*------------------------- Lecture des entrees ------------------------*
 
* Arguments obligatoires
* ----------------------
 
* Nom procedure simulation :
simu1    = mot tab1 . simulation ;
 
* lpini1 : jeu de parametres initiaux :
lpini1   = tab1 . parametres_initiaux ;
 
* labsc1 : valeurs des abscisses :
labsc1   = tab1 . points_de_mesure ;
 
* lsref1 : vecteur objectif :
lsref1   = tab1 . valeurs_objectif ;
 
* Arguments optionnels
* --------------------
 
* Choix methode de optimisation :
si (exis tab1 methode_optimisation) ;
  ilevmar1 = ega tab1.methode_optimisation 'LEVENBERG_MARQUARDT' ;
  imoca1   = ega tab1.methode_optimisation 'MOINDRE_CARRES' ;
  si (non (ilevmar1 ou imoca1)) ;
    erre '***** OPTIMISE : methode inconnue' ;
  fins ;
sino ;
  ilevmar1 = vrai ;
  imoca1   = faux ;
fins ;
 
* liste poids :
si (exis tab1 poids) ;
  lpoid1   = tab1 . poids ;
sino ;
  lpoid1   = prog (dime lsref1) * 1. ;
fins ;
 
* Bornes parametres :
si (exis tab1 bornes_inf) ;
  lpmin1   = tab1 . bornes_inf ;
sino ;
  lpmin1   = 0.5 * lpini1 ;
fins ;
si (exis tab1 bornes_sup) ;
  lpmax1   = tab1 . bornes_sup ;
sino ;
  lpmax1   = 1.5 * lpini1 ;
fins ;
 
* crit1  : critere de precision 
si (exis tab1 precision) ;
  crit1    = tab1 . precision ;
sino ;
  crit1    = 1.e-2 ;
fins ;
 
* nmini1 : nombre iterations max. boucle minimisation :
si (exis tab1 niter_max) ;
  nmini1   = tab1 . niter_max ;
sino ;
  nmini1   = 100 ;
fins ;
 
* lambda : parametre ponderation jacobienne (vs. hessienne) :
ilambda1 = exis tab1 lambda ;
si ilambda1 ;
  lambda0  = tab1 . lambda ;
fins ;
 
* Activation dessin evolution criteres :
idess1   = faux ;
si (exis tab1 dessin_evol) ;
  idess1   = tab1 . dessin_evol ;
fins ;
 
* Activation affichages calcul :
echoff1  = vrai ;
si (exis tab1 echo_simu) ;
  echoff1  = non (tab1 . echo_simu) ;
fins ;
vimpr1   = vale impr ;
 
*-------------------------- Initialisations ---------------------------*
 
* Point dupport rigidite associee a matrice hessienne :
vdim1    = vale dime ;
si (vdim1 neg 3) ; opti dime 3 ; fins ;
prig1    = manu poi1 (0 0 0) ;
si (vdim1 > 0 ) ; opti dime vdim1 ; fins ;
 
* Noms composantes primales / duales :
npar1    = dime lpini1 ;
lprim1   = mots ;
ldual1   = mots ;
repe bpar1 npar1 ;
  moti1    = chai 'P' &bpar1 ;
  lprim1   = lprim1 et moti1 ;
  moti1    = chai 'R' &bpar1 ;
  ldual1   = ldual1 et moti1 ;
fin bpar1 ;
*list lprim1 ; list ldual1 ;
opti inco lprim1 ldual1 ;
 
* Table enregistrement listreel ecarts relatifs pour dessin :
nsol1    = dime lsref1 ;
si idess1 ;
  tdess1   = tabl ;
  repe bsol1 nsol1 ;
    tdess1 . &bsol1 = prog ;
  fin bsol1 ;
fins ;
 
*------------------------- Affichage initial --------------------------*
 
saut 1 lign ;
mess '----------------------- Debut procedure OPTIMISE -----------------------' ;
saut 1 lign ;
si ilevmar1 ;
mess ' > Methode minimisation : LEVENBERG-MARQUARDT' ;
fins ;
si imoca1 ;
mess ' > Methode minimisation : MOINDRE CARRES' ;
fins ;
saut 1 lign ;
 
*------------------------- boucle minimisation ------------------------*
 
* Initialisations :
lpi1     = lpini1 ;
Rdeno1   = (0.5 * (somm (lsref1 * lsref1 * lpoid1 * lpoid1))) ** 0.5 ;
rcprec1  = (vale prec) ** 0.5 ;
liter1   = prog ;
iconst1  = faux ;
iannul1  = faux ;
iti1     = 0 ; 
repe bmini1 (nmini1 + 1) ;
 
* Calcul vecteur solution (appel a la simulation) :
  si echoff1 ; opti impr 30 ; fins ;
  lsoli1   = simu1 labsc1 lpi1 ;
  opti impr vimpr1 ;
 
* Calcul Ecart Relatif a la solution de reference :
  lRi1     = (lsoli1 - lsref1) * lpoid1 ;
  Fi1      = 0.5 * (somm (lRi1 * lRi1)) ;
  Rni1     = (Fi1 ** 0.5) / Rdeno1 ;
 
* Impressions messages critere et param. :
  repe bpar1 npar1 ;
    motpari1 = chai '   Param. ' &bpar1*14 ;
    valpari1 = extr lpi1 &bpar1 ;
    si (&bpar1 ega 1) ;
      motpar1  = motpari1 ;
      mespar1  = chai valpari1*14 ;
    sino ;
      motpar1  = chai motpar1 motpari1*(&bpar1*14) ;
      mespar1  = chai mespar1 valpari1*(&bpar1*14) ;
    fins ;
  fin bpar1 ;
  mlambda1 = lambda1 ;
  si (&bmini1 ega 1) ;
    si ilevmar1 ;
      entete1  = chai ' It'*5 'Critere'*18  'lambda'*31 motpar1 ;
    fins ;
    si imoca1 ;
      entete1  = chai ' It'*5 'Critere'*18  motpar1 ;
    fins ;
    mess  entete1 ;
    si ilambda1 ;
      mlambda1 = lambda0 ;
    sino ;
      mlambda1 = mot '-' ;
    fins ;
  fins ;
  si (non iannul1) ;
    si ilevmar1 ;
      messi1   = chai  iti1*5      Rni1*18  mlambda1*31 mespar1 ;
    fins ;
    si imoca1 ;
      messi1   = chai  iti1*5      Rni1*18  mespar1 ;
    fins ;
    mess messi1 ;
    si idess1 ; liter1  = liter1 et iti1 ; fins ;
    iti1     = iti1 + 1 ;
  fins ;
 
* Indicateur d'arret si critere constant d'une iteration sur l'autre :
    si (&bmini1 ega 1) ;
      Rnmin1   = Rni1 ;
    sino ;
      iconst1  = (abs (Rni1 - Rnmin1)) '<' (rcprec1 * Rni1) ;
      iconst1  = iconst1 et (non iannul1) ;
*list (abs (Rni1- Rnmin1)) ;
*list iconst1 ;
    fins ;
 
* Enregistrement resultats dans table :
  si ((Rni1 &lt;EG Rnmin1) ou (&bmini1 ega 1) ou iconst1) ;
    tab1 . parametres_finaux = lpi1 ;
    tab1 . valeurs_solution  = (lRi1 / lpoid1) + lsref1 ;
    tab1 . valeur_critere    = Rni1 ;
  fins ;
 
* Critere d'arret de la boucle de minimisation :
  si ((Rni1 < crit1) ou iconst1) ;
    saut 1 lign ;
    mess '------------------------ Fin procedure OPTIMISE ------------------------' ;
    saut 1 lign ;
    quit bmini1 ;
  fins ;
 
* Pilotage de lambda si Levenberg-Marquardt :
  si ilevmar1 ;
    iannul1  = faux ;
    si (&bmini1 ega 1) ;    comm initialisation de lambda ;
      si ilambda1 ;
        lambda1  = lambda0 ;
      sino ;
        lambda1  = 1. ;
      fins ;
    sino ;
      si (Rni1 > (1.01 * Rnmin1)) ;
        lambda1  = 10. * lambda1 ;
        lpi1     = lpi1p ;
        mess ' annulation iteration : reprise valeurs precedentes' ;
        iannul1  = vrai ;
        iter bmini1 ;
      sino ;
        si (Rni1 < (0.99 * Rnmin1)) ;
          lambda1  = 0.1 * lambda1 ;
        fins ;
      fins ;
    fins ;
  fins ;
 
* Stockage valeur meilleur critere :
  si (Rni1 < Rnmin1) ; Rnmin1   = Rni1 ; fins ;
*list Rnmin1 ;
 
* Calcul de la matrice Jacobienne de la fonction cout Fi1 :
* JacoFi1 : listreel contenant la Jacobienne de la fonction
* JacoFi1 : Jkl = dFl/dpk = ((k-1)*nsol1+l)-ieme element du listreel
  JacoFi1  = prog ;
  lderiv1  = enum ;
  repe bparam1 npar1 ;
    parami1  = extr lpi1 &bparam1 ;
    izeroi1  = faux ;
    si (parami1 ega 0.) ;
      izeroi1  = vrai ;
      pmaxi1   = lpmax1 extr &bparam1 ;
      si (pmaxi1 ega 0.) ;
        parami2  = 1. ;
      sino ;
        parami2  = 0.02 * pmaxi1 ;
      fins ;
    sino ;
      parami2  = 1.02 * parami1 ;
    fins ;
    lpi2     = (lpi1 enle &bparam1) inse &bparam1 parami2 ;
    si echoff1 ; opti impr 30 ; fins ;
    lsoli2   = simu1 labsc1 lpi2 ;
    opti impr vimpr1 ;
    si izeroi1 ;
      JacoFi2  = (lsoli2 - lsoli1) / parami2 ;
    sino ;
      JacoFi2  = (lsoli2 - lsoli1) / 0.02 / parami1 ;
    fins ;
    lderiv1  = lderiv1 et JacoFi2 ;
    JacoFi1  = JacoFi1 et JacoFi2 ;
  fin bparam1 ;
 
*----- MOINDRE_CARRES :
 
  si imoca1 ;
    lpi1p    = lpi1 ;
    lpi1     = moca lpi1 lsref1 lsoli1 lderiv1 'POIDS' lpoid1 ;
  fins ;
 
*----- LEVENBERG_MARQUARDT :
 
  si ilevmar1 ;
 
*   Construction de la matrice Hessienne = H ~ tJac * Jac :
    nbhes1   = npar1 * npar1 ;
    lrigi1   = prog ;
    ldiag1   = prog ;
    repe bhes1 nbhes1 ;
      Pi1      = (&bhes1 - 1) / npar1 + 1 ;
      Qi1      = &bhes1 - ((Pi1 - 1) * npar1) ;
      Hpqi1    = 0. ;
      repe bsol1 nsol1 ;
        JacoKP1   = JacoFi1 extr ((Pi1 - 1) * nsol1 + &bsol1) ;
        JacoKQ1   = JacoFi1 extr ((Qi1 - 1) * nsol1 + &bsol1) ;
        Hpqi1    = Hpqi1 + (JacoKP1 * JacoKQ1) ;
      fin bsol1 ;
      lrigi1   = lrigi1 et Hpqi1 ;
      si (Pi1 ega Qi1) ;
        ldiag1   = ldiag1 et Hpqi1 ;
      sino ;
        ldiag1   = ldiag1 et 0. ;
      fins ;
    fin bhes1 ;
    rigi1    = manu rigi prig1 lprim1 lrigi1 ;
    diagi1   = manu rigi prig1 lprim1 ldiag1 ;
    operi1   = rigi1 et (lambda1 * diagi1) ;
 
*   Caclul du gradient de Fi1 = JacFi1 * Ri1 :
    gradFi1  = vide chpoint ;
    repe bpar1 npar1 ;
      lJacoFi1 = extr JacoFi1 (lect ((&bpar1 - 1) * nsol1 + 1) pas 1 (npar1 * nsol1)) ;
      RJacoFj1 = 0. ;
      repe bsol1 nsol1 ;
        RJacoFj1 = (lRi1 extr &bsol1) * (lJacoFi1 extr &bsol1) + RJacoFj1 ;
      fin bsol1 ;
      moti1    = ldual1 extr &bpar1 ;
      gradFij1 = manu chpo prig1 1 moti1 RJacoFj1 nature discret ;
      gradFi1  = gradFi1 et gradFij1 ;
    fin bpar1 ;
 
*   Calcul de l'increment des parametres avec RESO :
    ldpi1    = RESO operi1 gradFi1 ;
 
*   Conversion LISTREEL :
    ldpi2    = prog ;
    repe bpar1 npar1 ;
      ldpi2    = ldpi2 et (ldpi1 extr vale prig1 (lprim1 extr &bpar1)) ;
    fin bpar1 ;
    ldpi1    = ldpi2 ;
 
*   Nouvelle liste de parametres :
    lpi1p    = lpi1 ;
    lpi1     = lpi1 - ldpi1 ;
 
  fins ;
*----- fin LEVENBERG_MARQUARDT
 
* Limitation des parametres aux bornes :
  lpi2     = prog ;
  repe bparam1 npar1 ;
    pi1      = extr lpi1 &bparam1   ;
    pimin1   = extr lpmin1 &bparam1 ;
    pimax1   = extr lpmax1 &bparam1 ;
    pi1      = maxi (prog pi1 pimin1) ;
    pi1      = mini (prog pi1 pimax1) ;
    lpi2     = lpi2 et pi1 ;
  fin bparam1 ;
  lpi1     = lpi2 ;
 
* Dessin evolution critere relatif de chaque objectif :
  si idess1 ;
    evcrit1 = vide evolution ;
    repe bsol1 nsol1 ;
      tdess1 . &bsol1 = (tdess1 . &bsol1) et ((lRi1 extr &bsol1) abs) ;
      moti1   = chai 'Obj.' (&bsol1)*6 ;
      evi1    = evol manu lege moti1 marq &bsol1 'Iter' liter1 'Critere' (tdess1 . &bsol1) ;
      evi1    = evi1 coul (&bsol1 + 1) ;
      evcrit1 = evcrit1 et evi1 ;
    fin bsol1 ;
    si (iti1 &lt;EG 10) ;
      dess evcrit1 nclk lege xbor -1. (flot iti1) xgra 1. posy exce ;
    sino ;
      dess evcrit1 nclk lege xbor -1. (flot iti1) posy exce ;
    fins ;
  fins ;
 
fin bmini1 ;
 
'FINP' TAB1 ;