Sources Esope | Cast3M

gcmmasub
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 &lt;= 0,  i = 1,...,m
C                xmin_j &lt;= x_j &lt;= 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 &lt;= b_i,  i = 1,...,m
C                  alfa_j &lt;= x_j &lt;= 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