Télécharger moca_mma.procedur

Retour à la liste

* MOCA_MMA  PROCEDUR  FD218221  26/04/29    21:15:02     12530          
DEBP MOCA_MMA ;
 
* Description :
* -------------
* Ajustement des parametres d'une fonction pour coller au mieux a
* une serie de points
* Methode des moindres carres ponderee, sous contraintes,
* basee sur la Methode des Asymptotes Mobiles (MMA)
*
* Notations :
* -----------
* x = [x1,x2,...,xn] vecteur des parametres (les inconnues du probleme)
* n                  dimension de x
* xmin  xmax         vecteurs des bornes sur les parametres
* P = [P1,P2,...,Pm] vecteur des points connus
* H = [H1,H2,...,Hm] vecteur des valeurs connues
* w = [w1,w2,...,wm] vecteur des ponderations associees aux points P
* m                  dimension de P, H et w
* h(Pi,x)            fonction a ajuster (decrite par la procedure "nomproc")
* gj(x) j=1,...,r    ensemble de fonctions contraintes sur les parametres
* r                  nombre de relations gj imposees (optionnelles)
*
* Objectif :
* ----------
*                                   m
*                                  ----
*                              1   \
*       Minimiser :    S(x) = ---  /     wi * [h(Pi,x) - Hi]^2
*       sur x                  2   ----
*                                  i = 1
*
*       tel que   :    gj(x) < 0   j = 1,...,r
*                      xmin < xi < xmax
 
 
 
************************************************************************
*          A C Q U I S I T I O N   D E S   P A R A M E T R E S         *
*                            D' E N T R E E                            *
************************************************************************
 
* Lecture des arguments obligatoires par couple ( 'MOT_CLEF' ,  ARG1 )
nargobl  = 4 ;
REPE b1 nargobl ;
* Lecture du mot clef
  ARGU mot1*'MOT' ;
* Variables de depart : x0
  SI (EGA mot1 'X0') ;
    ARGU x0*'LISTREEL' ;
  FINSI ;
* Points connus (Pi,Hi) : lp,lh
  SI (EGA mot1 'ABSC') ;
    ARGU lp*'LISTREEL' ;
  FINSI ;
  SI (EGA mot1 'ORDO') ;
    ARGU lh*'LISTREEL' ;
  FINSI ;
* Procedure decrivant la fonction h(lp,x) a caler
  SI (EGA mot1 'PROC') ;
    ARGU mot2*'MOT' ;
    nomproc  = MOT (TEXT (CHAI mot2)) ;
  FINSI ;
FIN b1 ;
 
* Dimension des donnees
n        = DIME x0 ;
m        = DIME lp ;
 
* Quelques verifications
SI (NEG (DIME lh) m) ;
  ERRE 'Les listes ''ABSC'' et ''ORDO'' n''ont pas la meme dimension' ;
FINSI ;
 
* Valeurs par defaut des arguments optionnels
xmin     = PROG n*0. ;
xmax     = PROG n*1000. ;
lw       = PROG m*1. ;
icont    = FAUX ;
move     = 0.1 ;
nitmax   = 50 ;
crit1    = 1.E-4 ;
iecho    = 0 ;
bvisu    = FAUX ;
 
* Lecture des arguments optionnels par couple ( 'MOT_CLEF' ,  ARG1 )
nargopt  = 9 ;
REPE b1 nargopt ;
* Lecture du mot clef
  ARGU mot1/'MOT' ;
  SI (EGA (TYPE mot1) 'ANNULE') ;
    QUIT b1 ;
  FINSI ;
* Bornes min/max sur les variables : xmin,xmax
  SI (EGA mot1 'XMIN') ;
    ARGU xmin*'LISTREEL' ;
  FINSI ;
  SI (EGA mot1 'XMAX') ;
    ARGU xmax*'LISTREEL' ;
  FINSI ;
* Ponderations wi : lw
  SI (EGA mot1 'POID') ;
    ARGU lw*'LISTREEL' ;
  FINSI ;
* Procedure decrivant les fonctions contraintes g(x)
  SI (EGA mot1 'CONT') ;
    ARGU mot2*'MOT' ;
    nomcont  = MOT (TEXT (CHAI mot2)) ;
    icont    = VRAI ;
  FINSI ;
* Critere de convergence sur le residu kkt : crit1
  SI (EGA mot1 'CRIT') ;
    ARGU crit1*'FLOTTANT' ;
  FINSI ;
* Nombre d'iterations max. pour la MMA : nitmax
  SI (EGA mot1 'ITER') ;
    ARGU nitmax*'ENTIER' ;
  FINSI ;
* Limite d'increment pour MMA : move
  SI (EGA mot1 'INCR') ;
    ARGU move*'FLOTTANT' ;
  FINSI ;
* Niveau d'affichage
  SI (EGA mot1 'ECHO') ;
    ARGU iecho*'ENTIER' ;
  FINSI ;
* Visualisation ?
  SI (EGA mot1 'VISU') ;
    bvisu    = VRAI ;
    ARGU l1plot*'LISTREEL' ;
  FINSI ;
FIN b1 ;
 
* Quelques verifications
SI (NEG (DIME xmin) n) ;
  ERRE 'Les listes ''X0'' et ''XMIN'' n''ont pas la meme dimension' ;
FINSI ;
SI (NEG (DIME xmax) n) ;
  ERRE 'Les listes ''X0'' et ''XMAX'' n''ont pas la meme dimension' ;
FINSI ;
SI (NEG (DIME lw) m) ;
  ERRE 'Les listes ''ABSC'' et ''POID'' n''ont pas la meme dimension' ;
FINSI ;
SI (move &lt;EG 0.) ;
  ERRE 'La limite d''increment ''INCR'' doit etre > 0' ;
FINSI ;
SI (crit1 < 0.) ;
  ERRE 'Le critere de convergence doit etre > 0.' ;
FINSI ;
SI (nitmax < 1) ;
  ERRE 'Le nombre d''iterations doit etre > 0' ;
FINSI ;
SI ((iecho < 0) OU (iecho > 2)) ;
  ERRE 'L''entier ''ECHO'' doit etre compris entre 0 et 2' ;
FINSI ;
 
* Pour la visualisation
SI bvisu ;
  evmes    = EVOL 'ROUG' 'MANU' 'STYL' 'NOLI' 'MARQ' 'ROND' lp lh ;
FINSI ;
 
 
 
************************************************************************
*          I N I T I A L I S A T I O N   D E S   D O N N E E S         *
*                P O U R   L' O P E R A T E U R   M M A                *
************************************************************************
 
* On applique l'operateur MMA (voir notice) avec le parametrage suivant :
*    f0(x)      = 0
*    fi(x)      =  h(Pi,x) - Hi   pour i dans [1,m]
*    f{m+i}(x)  = -h(Pi,x) + Hi   pour i dans [1,m]
*    f{2m+j}(x) =  gj             pour j dans [1,r]
*            a0 = 0
*            ai = 0               pour i dans [1,2m+r]
*            ci = 0               pour i dans [1,2m]
*       c{2m+j} = "grand"         pour j dans [1,r]
*            di = wi              pour i dans [1,2m]
*       d{2m+j} = 0               pour j dans [1,r]
* avec
* m : nombre de points connus (Pi,Hi)
* r : nombre de relations supplementaires (gj) sur les parametres x
* q : nombre total de contraintes imposees pour l'operateur MMA
 
* Calcul des valeurs h(Pi,x0) et des dh(Pi,x0)/dx1 ... dh(Pi,x0)/dxn
f0       = 0. ;
df0      = PROG n*0. ;
h0 dh0   = (TEXT nomproc) lp x0 ;
f        = (h0 - lh) ET ((-1. * h0) + lh) ;
df       = TABL ;
REPE b1 m ;
  df . &b1 = EXTR dh0 &b1 ;
FIN b1 ;
REPE b1 m ;
  df . (&b1 + m) = -1. * (df . &b1) ;
FIN b1 ;
 
* Calcul des contraintes gj(x0) et des dgj(x0)/dx1 ... dgj(x0)/dxn
r        = 0 ;
SI icont ;
  g0 dg0   = (TEXT nomcont) x0 ;
  r        = DIME g0 ;
  f        = f ET g0 ;
  REPE b1 r ;
    df . (&b1 + (2 * m)) = EXTR dg0 &b1 ;
  FIN b1 ;
FINSI ;
q        = (2 * m) + r ;
 
* Calcul de l'ecart quadratique S(x0)
mse0     = 0.5 * (SOMM (lw * ((h0 - lh) ** 2))) ;
* Initialisation de la table pour l'operateur MMA :
t        = TABL ;
t . 'X'     = x0 ;
t . 'XMIN'  = xmin ;
t . 'XMAX'  = xmax ;
t . 'F0VAL' = f0 ;
t . 'DF0DX' = df0 ;
t . 'FVAL'  = f ;
t . 'DFDX'  = df ;
t . 'A0'    = 0. ;
t . 'A'     = PROG q*0. ;
t . 'C'     = PROG (2 * m)*0. ;
t . 'D'     = lw ET lw ;
SI icont ;
  t . 'C'     = (t . 'C') ET (PROG r*1.E5) ;
  t . 'D'     = (t . 'D') ET (PROG r*0.) ;
FINSI ;
t . 'MOVE'  = move ;
 
* Pour l'affichage
SI (iecho > 0) ;
  txt0     = CHAI 'It' *5 ;
  txt1     = CHAI  '0' *5 ;
  REPE b1 n ;
    txt0     = CHAI txt0 (CHAI 'x' &b1) <11 ;
    txt1     = CHAI 'FORMAT' '(F10.5)' txt1 (EXTR x0 &b1) <11 ;
  FIN b1 ;
  txt0     = CHAI txt0 'mse' <11 'kkt' <11 ;
  txt1     = CHAI 'FORMAT' '(F10.5)' txt1 mse0 <11 ;
  MESS txt0 ;
  MESS txt1 ;
* Historique des valeurs de x selon les iterations
  SI (iecho > 1) ;
    lmse     = PROG mse0 ;
    lx       = ENUM x0 ;
  FINSI ;
FINSI ;
* Pour la visualisation
SI bvisu ;
  l2plot   = (TEXT nomproc) l1plot x0 'VISU' ;
  ev0      = EVOL 'GRIS' 'MANU' l1plot l2plot ;
  DESS (evmes ET ev0) 'TITR' (CHAI 'Iteration :' ' ' 0) ;
FINSI ;
 
* Boucle d'optimisation par MMA
REPE loop nitmax ;
* Calcul d'un nouveau jeu de variables xnew par MMA
  MMA t ;
  xnew     = t . 'X' ;
* Mise a jour des valeurs des fonctions
  h1 dh1   = (TEXT nomproc) lp xnew ;
  f        = (h1 - lh) ET ((-1. * h1) + lh) ;
  df       = TABL ;
  REPE b1 m ;
    df . &b1 = EXTR dh1 &b1 ;
  FIN b1 ;
  REPE b1 m ;
    df . (&b1 + m) = -1. * (df . &b1) ;
  FIN b1 ;
  SI icont ;
    g1 dg1   = (TEXT nomcont) xnew ;
    f        = f ET g1 ;
    REPE b1 r ;
      df . (&b1 + (2 * m)) = EXTR dg1 &b1 ;
    FIN b1 ;
  FINSI ;
  t . 'FVAL'  = f ;
  t . 'DFDX'  = df ;
* Calcul de l'ecart quadratique S(x)
  mse1     = 0.5 * (SOMM (lw * ((h1 - lh) ** 2))) ;
* Calcul du residu pour les conditions KKT
  res kkt2 kktinf = KKT_MMA t ;
* Pour l'affichage
  SI (iecho > 0) ;
    txt1     = CHAI &loop *5 ;
    REPE b1 n ;
      txt1     = CHAI 'FORMAT' '(F10.5)' txt1 (EXTR xnew &b1) <11 ;
    FIN b1 ;
    txt1     = CHAI 'FORMAT' '(F10.5)' txt1 mse1 <11 kkt2 <11 ;
    MESS txt1 ;
  FINSI ;
* Historique des valeurs de x selon les iterations
  SI (iecho > 1) ;
    lmse     = lmse ET mse1 ;
    lx       = lx   ET xnew ;
  FINSI ;
* Pour la visualisation
  SI bvisu ;
    l2plot   = (TEXT nomproc) l1plot xnew 'VISU' ;
    ev1      = EVOL 'GRIS' 'MANU' l1plot l2plot ;
    DESS (evmes ET ev1) 'TITR' (CHAI 'Iteration :' ' ' &loop) ;
  FINSI ;
* Critere d'arret
  SI (kkt2 < crit1) ;
    QUIT loop ;
  FINSI ;
FIN loop ;
 
* Ecriture des arguments de sortie
RESP xnew ;
SI (iecho > 1) ;
  RESP lmse lx ;
FINSI ;
 
FINP ; 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales