C GCMMASUB SOURCE FD218221 26/10/07 21:15:06 12662 SUBROUTINE GCMMASUB(M,N,XVAL,XMIN,XMAX,LOW,UPP, & F0VAL,DF0DX,FVAL,DFDX, & RAA0,RAA,ALBEFA, & ALFA,BETA,P0,Q0,P,Q,B,R0) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C ---------------------------------------------------------------------- C Ce code est une adaptation en FORTRAN/ESOPE de l'algorithme de la MMA C (Methode des Asymptotes Mobiles) propose initialement en langage C Matlab par K. Svarberg (sources : https://www.smoptit.se/) C ---------------------------------------------------------------------- C Cette subroutine definie un sous-probleme d'optimisation convexe C selon la methode des asymptotes mobiles globalement convergente (GC-MMA) C C Le probleme d'optimisation initial est le suivant : C Minimiser : f_0(x) + a_0*z + SUM[ c_i*y_i + 0.5*d_i*(y_i)^2 ] C soumis a : f_i(x) - a_i*z - y_i <= 0, i = 1,...,m C xmin_j <= x_j <= xmax_j, j = 1,...,n C z >= 0, y_i >= 0, i = 1,...,m C C Le probleme est transforme en une suite de sous-problemes convexes : C Minimiser : SUM[ p0_j/(upp_j-x_j) + q0_j/(x_j-low_j) ] + a_0*z + C + SUM[ c_i*y_i + 0.5*d_i*(y_i)^2 ] C Soumis a : SUM[ p_ij/(upp_j-x_j) + q_ij/(x_j-low_j) ] - a_i*z - y_i <= b_i, i = 1,...,m C alfa_j <= x_j <= beta_j, j = 1,...,n C z >= 0, y_i >= 0, i = 1,...,m C C Cette subroutine calcule les matrices/vecteurs decrivant un sous-probleme C Elle necessite d'etre iteree pour approcher la solution du probleme initial C C Entrees : C --------- C M = Nombre de contraintes C N = Nombre de variables x_j C XVAL = Vecteur (N) des variables x_j C XMIN/XMAX = Vecteur (N) des bornes inferieures/superieures pour les variables x_j C LOW/UPP = Vecteurs (N) des asymptotes inferieures/superieures du sous-probleme C DF0DX = Vecteur (N) des valeurs du gradient de la fonction objectif f_0, C par rapport aux variables x_j, calcules en XVAL C FVAL = Vecteur (M) des valeurs des fonctions contraintes f_i, C calculees en XVAL C DFDX = Matrice (M x N) des valeurs du gradient des fonctions contraintes f_i, C par rapport aux variables x_j, calcules en XVAL C DFDX(i,j) = gradient de f_i par rapport a x_j C ALBEFA = Pourcentage d'ecart par rapport aux asymptotes alfa et beta C RAA0 = Parametre de conservatisme / courbure de l'objectif C RAA = Vecteur (M) des coefficients de conservatisme / courbure des contraintes C C Sorties : C --------- C ALFA/BETA = Vecteurs (N) des bornes inferieures/superieures mises a jour C P0/Q0 = Vecteurs (N) des coefficients des asymptotes superieure/inferieure C dans la fonction objectif C P/Q = Matrices (M x N) des coefficients des asymptotes superieure/inferieure C des contraintes C B = Vecteur (M) des seconds membres des contraintes C R0 = Parametre d'ajustement de la fonction objectif approximee C ---------------------------------------------------------------------- C Arguments INTEGER M,N REAL*8 RAA0,ALBEFA,F0VAL,R0 REAL*8 XVAL(N),XMIN(N),XMAX(N),LOW(N),UPP(N),RAA(M) REAL*8 DF0DX(N),FVAL(M),DFDX(M,N) REAL*8 ALFA(N),BETA(N),P0(N),Q0(N),P(M,N),Q(M,N),B(M) C Variables locales INTEGER I,J C Mise a jour des asymptotes alfa et beta DO 10 I = 1, N ALFA(I) = MAX(LOW(I)+ALBEFA*(XVAL(I)-LOW(I)),XMIN(I)) BETA(I) = MIN(UPP(I)-ALBEFA*(UPP(I)-XVAL(I)),XMAX(I)) 10 CONTINUE C Calcul de p0, q0 et r0 R0 = F0VAL DO 20 I = 1, N DF0 = DF0DX(I) XMAMI = MAX(XMAX(I)-XMIN(I), 0.00001D0) P0(I) = MAX( DF0,0.D0) + 0.001D0*ABS(DF0) + RAA0/XMAMI Q0(I) = MAX(-DF0,0.D0) + 0.001D0*ABS(DF0) + RAA0/XMAMI P0(I) = P0(I) * ((UPP(I) - XVAL(I))**2) Q0(I) = Q0(I) * ((XVAL(I) - LOW(I))**2) R0 = R0 - P0(I)/(UPP(I)-XVAL(I)) - Q0(I)/(XVAL(I)-LOW(I)) 20 CONTINUE C Calcul de p, q et b DO 30 J = 1, M B(J) = -FVAL(J) DO 40 I = 1, N DFJ = DFDX(J,I) P(J,I) = MAX(0.D0, DFJ) + 0.001D0*ABS(DFJ) & + (RAA(J)/(XMAX(I)-XMIN(I))) Q(J,I) = MAX(0.D0,-DFJ) + 0.001D0*ABS(DFJ) & + (RAA(J)/(XMAX(I)-XMIN(I))) P(J,I) = P(J,I) * ((UPP(I)-XVAL(I))**2) Q(J,I) = Q(J,I) * ((XVAL(I)-LOW(I))**2) B(J) = B(J) + P(J,I) / (UPP(I) - XVAL(I)) & + Q(J,I) / (XVAL(I) - LOW(I)) 40 CONTINUE 30 CONTINUE RETURN END