$$$$ MOCA_MMA * 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 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 ;