$$$$ HASOFER * HASOFER PROCEDUR AM 09/12/07 21:15:37 6578 DEBPROC HASOFER TAB1*TABLE; *----------------------------------------------------------------- * Calcul de l'indice de fiabilite d'Hasofer-Lind * Entree : * - nombre de variables : tab1.nbre_variables * - parametres des variables NATAF : tab1.param_va * - matrice de covariance : tab1.matcov; * - coefficient de majoration 1 : tab1.major1 * - coefficient de majoration 2 : tab1.major2 * - precision relative du BETA : tab1.prec * - nombre max d'iterations : tab1.itmax * - point de depart P0 : tab1.depart * - parametres modele MECA : tab1.parametres * Sortie : * - point d'arrivee P* : tab1.arrivee * - indices de fiabilite : tab1.beta * - nombres d'appels à la FEL : tab1.appels * - points P*(iter) : tab1.beta_point * - valeurs de G(iter) : tab1.glim * - ecarts |UP*i(iter+1)-UP*i(iter)| : tab1.ecu * - increments DELTA_X(iter) : tab1.delx * - indice d'erreur : err1 *----------------------------------------------------------------- * Recuperation des entrees *------------------------- nva1=tab1.nbre_variables; para1=tab1.param_va; matc1=tab1.matcov; vecu1=tab1.depart; SI (EXISTE tab1.parametres); para2=tab1.parametres; FINSI; * Controle des entrees *--------------------- err1=0; SI ((DIME vecu1) NEG nva1); MESS 'HASOFER : incoherence DIME depart / nva'; err1=1; QUIT HASOFER ; FINSI; SI ((DIME para1) NEG nva1); MESS 'HASOFER : incoherence DIME param_va / nva'; err1=1; QUIT HASOFER ; FINSI; naux1 = 0; REPE bouc1 nva1; naux1=naux1 + &bouc1; FIN bouc1; SI ((DIME matc1) NEG naux1); MESS 'HASOFER : incoherence DIME matcov / nva'; err1=1; QUIT HASOFER ; FINSI; * Declarations et initialisations *-------------------------------- * Coefficient d'amplification c1 maj1=10.; * Coefficient d'amplification c2 maj2=10.; * Precision sur beta (inutile en pratique) epsilon_beta prec1=0.005; * Precision sur la distance entre deux P* consecutifs epsilon_p prec2=0.05; * Nombre iteration maximal itmax1=100; tab1.beta=TABLE; tab1.beta_pointu=TABLE; tab1.beta_pointx=TABLE; tab1.glim=TABLE; tab1.ecu=TABLE; tab1.delx=TABLE; tab1.distance=TABLE; tab1.gradu=TABLE; tab1.gradx=TABLE; tab1.cal_grad=TABLE; tabx1=TABLE; tabx1.nbre_variables=nva1; SI (EXISTE para2); tabx1.parametres=para2; FINSI; tabg1=TABLE; tabg1.nbre_variables=nva1; SI (EXISTE para2); tabg1.parametres=para2; FINSI; tab1.appels=0; ch1=MOTS; id_cal1=LECT; REPE bouc1 nva1; id_cal1=id_cal1 ET (LECT 1); ch1=ch1 ET (MOTS (CHAIN 'V' &bouc1)); FIN bouc1; * Calcul de BETA beta1 = NORV1 vecu1; *------------------------------------------------------------ * Boucle principale *------------------ REPE boucp1; SI (&boucp1 > itmax1); mess 'HASOFER : depassement de itmax'; err1=1; QUIT boucp1; FINSI; * Point courant en espace physique *--------------------------------- tabn1=TABLE; tabn1.transformation_directe=FAUX; tabn1.points_espace_reference=vecu1; tabn1.noms_des_variables=ch1; tabn1.matcov=matc1; tabn1.param_va=para1; NATAF tabn1; vecx1=tabn1.points_espace_physique; * Valeur de la fonction d'EL au point courant *-------------------------------------------- tabx1.variables=vecx1; err1=FLIM1 tabx1; SI (EGA err1 0); g1=tabx1.valeur_G; precr1=tabx1.precision_R; precs1=tabx1.precision_S; r1=tabx1.valeur_R; s1=tabx1.valeur_S; SINON; QUIT boucp1; FINSI; tab1.appels=tab1.appels+1; * Valeurs de la fonction d'EL aux bornes de la perturbation delta_X *------------------------------------------------------------------ * Increments sur les axes en espace physique *------------------------------------------- SI (EGA &boucp1 1); delx1=maj1*(precr1+precs1)*vecx1; SINON; delx1=delx2; FINSI; tabg1.vecx=vecx1; tabg1.delx=delx1; tabg1.calcul=id_cal1; err1=GDFLIM1 tabg1; SI (EGA err1 0); delg1=tabg1.valeurs_G; appel1=tabg1.appels; SINON; QUIT boucp1; FINSI; tab1.appels=tab1.appels + appel1; * Gradient en espace physique X *------------------------------ gradx1=PROG; REPE bouc2 nva1; dg1=EXTR &bouc2 delg1; dx1=EXTR &bouc2 delx1; SI (EGA (EXTR &bouc2 id_cal1) 1); gradx1=gradx1 ET (PROG ((dg1-g1)/dx1)); SINON; gradx1=gradx1 ET (PROG (EXTR &bouc2 gradx2)); FINSI; FIN bouc2; * Gradient en espace standard U *------------------------------ * Increment delta_Ui pour calculer le vecteur gradient de X par rapport a Ui delu1=1.E-7; gradu1=PROG; REPE bouc3 nva1; * Gradient de X par rapport a Ui *------------------------------- ui1=EXTR &bouc3 vecu1; v_delu1=INSERER vecu1 &bouc3 (ui1+delu1); v_delu1=ENLEVER v_delu1 (&bouc3+1); tabn1=TABLE; tabn1.transformation_directe=FAUX; tabn1.points_espace_reference=v_delu1; tabn1.noms_des_variables=ch1; tabn1.matcov=matc1; tabn1.param_va=para1; NATAF tabn1; v_delx1=tabn1.points_espace_physique; gradxu1=(v_delx1 - vecx1)/delu1; * Gradient de G par rapport a Ui *------------------------------- gradui1=LTL gradx1 gradxu1; gradu1=gradu1 ET (PROG gradui1); FIN bouc3; * Nouveau point en espace standard (Rackwitz-Fiessler) *----------------------------------------------------- alpu1=(-1.)*gradu1/(NORV1 gradu1); vecu2=((LTL vecu1 alpu1)+(g1/(NORV1 gradu1)))*alpu1; beta2=NORV1 vecu2; * Nouveau point en espace physique *--------------------------------- tabn1=TABLE; tabn1.transformation_directe=FAUX; tabn1.points_espace_reference=vecu2; tabn1.noms_des_variables=ch1; tabn1.matcov=matc1; tabn1.param_va=para1; NATAF tabn1; vecx2=tabn1.points_espace_physique; * Stockage des resultats *----------------------- tab1.beta.&boucp1=beta2; tab1.beta_pointu.&boucp1=vecu2; tab1.beta_pointx.&boucp1=vecx2; tab1.glim.&boucp1=g1; tab1.ecu.&boucp1=ABS (vecu2-vecu1); tab1.distance.&boucp1=NORV1 (vecu2-vecu1); tab1.gradu.&boucp1=gradu1; tab1.gradx.&boucp1=gradx1; tab1.iterations=&boucp1; tab1.delx.&boucp1=delx1; tab1.cal_grad.&boucp1 = id_cal1; * Verification du critere de convergence entre 2 P* consecutifs * (crtitere plus severe que sur 2 BETA consecutifs) *-------------------------------------------------------------- indc=1; REPETER boucc1 nva1; SI ((EXTR &boucc1 tab1.ecu.&boucp1) < prec2); indc=indc*1; SINON; indc=indc*0; FINSI; FIN boucc1; SI (EGA indc 1); * Convergence : fin ******************** tab1.arrivee = vecu2; QUIT boucp1; SINON; * Non convergence ***************** * Increments optimaux sur les axes en espace physique *---------------------------------------------------- * Distance entre 2 points P* consecutifs *--------------------------------------- d1_x12=PROG; REPE bouc4 nva1; x1=(EXTR &bouc4 vecx1); x2=(EXTR &bouc4 vecx2); SI (EGA (EXTR &bouc4 id_cal1) 1); d1_x12=d1_x12 ET (PROG (x2 - x1)); SINON; d1_x12=d1_x12 ET (PROG ((x2 - x1) + (EXTR &bouc4 d2_x12))); FINSI; FIN bouc4; * Indice de calcul des derivees partielles: 1 calcul, 0 pas de calcul *------------------------------------------------------------------------- id_cal1=LECT; REPE bouc5 nva1; dx1=(EXTR &bouc5 delx1); SI (OU ((ABS (EXTR &bouc5 d1_x12)) > dx1) (EGA &boucp1 1)); id_cal1=id_cal1 ET (LECT 1); SINON; id_cal1=id_cal1 ET (LECT 0); FINSI; FIN bouc5; * Calcul des gradients et increments *----------------------------------- gradx1=PROG; delx2=PROG; REPE bouc7 nva1; dx1=(EXTR &bouc7 delx1); SI (EGA (EXTR &bouc7 id_cal1) 1); dg1=(EXTR &bouc7 delg1); gradx1=gradx1 ET (PROG ((dg1-g1)/dx1)); SINON; gradx1=gradx1 ET (PROG (EXTR &bouc7 gradx2)); FINSI; dx1=maj2*((precr1*r1)+(precs1*s1))/(EXTR &bouc7 gradx1); delx2=delx2 ET (ABS (PROG dx1)); FIN bouc7; * Mises à jour pour l'iteration suivante *--------------------------------------- gradx2=gradx1; beta1=beta2; vecu1=vecu2; d2_x12=d1_x12; FINSI; FIN boucp1; *------------------ * Fin Boucle principale *------------------------------------------------------------ FINPROC err1;