Sources Esope | Cast3M

mmasub
C MMASUB    SOURCE    FD218221  25/09/03    21:15:04     12351          
      SUBROUTINE MMASUB(M,N,ITER,XVAL,XMIN,XMAX,XOLD1,XOLD2,
     &                  F0VAL,DF0DX,FVAL,DFDX,A0,A,C,D,MOVE,
     &                  XMMA,YMMA,ZMMA,LAM,XSI,ETA,MU,ZET,S,LOW,UPP)
 
C     Typages implicites habituels
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
C ------------------------------------------------------------------------------------------------
C Ce code est une adaptation en FORTRAN de l'algorithme de la MMA (Methode des Asymptotes Mobiles)
C propose initialement en langage Matlab par K. Svarberg
C Les sources sont disponibles ici : https://www.smoptit.se/
C ------------------------------------------------------------------------------------------------
C    Cette subroutine effectue une iteration de la MMA, afin de resoudre un probleme
C    d'optimisation non lineaire suivant :
C
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  Entrees :
C  ---------
C   m    = Nombre de contraintes
C   n    = Nombre de variables x_j
C  iter  = Numero de l'iteration courante ( =1 la premiere fois que mmasub est appellee)
C  xval  = Vecteur colonne (n x 1) des variables x_j
C  xmin  = Vecteur colonne (n x 1) des bornes inferieures pour les variables x_j
C  xmax  = Vecteur colonne (n x 1) des bornes superieures pour les variables x_j
C  xold1 = xval, a 1 iteration precedente (a condition que iter>1)
C  xold2 = xval, a 2 tterations precedentes (a condition que iter>2)
C  f0val = Valeur de la fonction objectif f_0 en xval
C  df0dx = Vecteur colonne (n x 1) des valeurs du gradient de la fonction objectif f_0,
C          par rapport aux variables x_j, calcules en xval
C  fval  = Vecteur colonne (m x 1) des valeurs des fonctions contraintes f_i,
C          calculees en at 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  low   = Vecteur colonne (n x 1) des valeurs des asymptotes inferieures
C          a l'iteration precedente (a condition que iter>1)
C  upp   = Vecteur colonne (n x 1) des valeurs des asymptotes superieures
C          a l'iteration precedente (a condition que iter>1)
C  a0    = Constante a_0
C  a     = Vecteur colonne (m x 1) des constantes a_i
C  c     = Vecteur colonne (m x 1) des constantes c_i
C  d     = Vecteur colonne (m x 1) des constantes d_i
C  move  = Limite de variation pour les variables
C
C  Sorties :
C  ---------
C  xmma  = Vecteur colonne (n x 1) des variables x_j optimisees
C          pour le sous probleme courant de la MMA
C  ymma  = Vecteur colonne (m x 1) des variables y_j optimisees
C          pour le sous probleme courant de la MMA
C  zmma  = Valeur de la variable z optimisee
C          pour le sous probleme courant de la MMA
C  lam   = Vecteur colonne (m x 1) des multiplicateurs de Lagrange
C          des contraintes generales f_i de la MMA
C  xsi   = Vecteur colonne (n x 1) des multiplicateurs de Lagrange
C          des contraintes sur les bornes inferieures alfa_j - x_j &lt;= 0
C  eta   = Vecteur colonne (n x 1) des multiplicateurs de Lagrange
C          des contraintes sur les bornes superieures x_j - beta_j &lt;= 0
C   mu   = Vecteur colonne (m x 1) des multiplicateurs de Lagrange
C          des contraintes de non negativite -y_i &lt;= 0
C  zet   = Multiplicateurs de Lagrange de la contrainte de non negativite -z &lt;= 0
C   s    = Vecteur colonne (m x 1) des variables d'ecart
C          des contraintes generales de la MMA f_i
C  low   = Vecteur colonne (n x 1) des valeurs des asymptotes inferieures
C          mises a jour
C  upp   = Vecteur colonne (n x 1) des valeurs des asymptotes superieures
C          mises a jour
C ------------------------------------------------------------------------------------------------
 
C     Arguments de sortie
      REAL*8 XMMA(N,1),YMMA(M,1),ZMMA
      REAL*8 LAM(M,1),XSI(N,1),ETA(N,1),MU(M,1),ZET,S(M,1)
C     Arguments d'entree
      INTEGER M,N,ITER
      REAL*8 A0,F0VAL,FVAL(M,1),XVAL(N,1),XMIN(N,1),XMAX(N,1)
      REAL*8 XOLD1(N,1),XOLD2(N,1),DF0DX(N,1),DFDX(M,N)
      REAL*8 LOW(N,1),UPP(N,1),MOVE
      REAL*8 A(M,1),C(M,1),D(M,1)
C     Arguments internes
      REAL*8 RAA0,ALBEFA,ASYINIT,ASYINCR,ASYDECR
      REAL*8 B(M,1),EEEN(N,1),EEEM(M,1),ZERON(N,1)
      REAL*8 ZZZ(N,1),ZZZ1(N,1),ZZZ2(N,1),FACTOR(N,1)
      REAL*8 LOWMIN(N,1),LOWMAX(N,1),UPPMIN(N,1),UPPMAX(N,1)
      REAL*8 ALFA(N,1),BETA(N,1),XMAMI(N,1)
      REAL*8 XMAMIEPS(N,1),XMAMIINV(N,1),UX1(N,1),UX2(N,1)
      REAL*8 XL1(N,1),XL2(N,1),UXINV(N,1),XLINV(N,1)
      REAL*8 P0(N,1),Q0(N,1),PQ0(N,1),P(M,N),Q(M,N),PQ(M,N)
      REAL*8 TEMP(M,N)
C     Valeurs des parametres
      EPSIMIN = 1.D-7
      RAA0 = 1.D-5
      ALBEFA = 0.1D0        
      ASYINIT = 0.5D0      
      ASYINCR = 1.2D0     
      ASYDECR = 0.7D0
      ASYMIN = 0.01D0
      ASYMAX = 10.D0
      CALL UN(EEEN,N,1)
      CALL UN(EEEM,M,1)
      CALL ZERO(ZERON,N,1)
 
C     Mise a jour des asymptotes low etd upp
      IF (ITER.LT.3) THEN
          LOW = XVAL - ASYINIT*(XMAX-XMIN)
          UPP = XVAL + ASYINIT*(XMAX-XMIN)
      ELSE
          ZZZ = (XVAL-XOLD1)*(XOLD1-XOLD2)
          FACTOR = EEEN
          DO I=1,N
              IF (ZZZ(I,1).GT.0.D0) THEN
                  FACTOR(I,1) = ASYINCR
              ELSEIF (ZZZ(I,1).LT.0.D0) THEN
                  FACTOR(I,1) = ASYDECR
              ENDIF
          ENDDO
          LOW = XVAL - FACTOR*(XOLD1 - LOW)
          UPP = XVAL + FACTOR*(UPP - XOLD1)
          LOWMIN = XVAL - ASYMAX*(XMAX-XMIN)
          LOWMAX = XVAL - ASYMIN*(XMAX-XMIN)
          UPPMIN = XVAL + ASYMIN*(XMAX-XMIN)
          UPPMAX = XVAL + ASYMAX*(XMAX-XMIN)
          LOW = MAX(LOW,LOWMIN)
          LOW = MIN(LOW,LOWMAX)
          UPP = MIN(UPP,UPPMAX)
          UPP = MAX(UPP,UPPMIN)
      ENDIF
 
C     Calcul des bornes alfa et beta
      ZZZ1 = LOW + ALBEFA*(XVAL-LOW)
      ZZZ2 = XVAL - MOVE*(XMAX-XMIN)
      ZZZ  = MAX(ZZZ1,ZZZ2)
      ALFA = MAX(ZZZ,XMIN)
      ZZZ1 = UPP - ALBEFA*(UPP-XVAL)
      ZZZ2 = XVAL + MOVE*(XMAX-XMIN)
      ZZZ  = MIN(ZZZ1,ZZZ2)
      BETA = MIN(ZZZ,XMAX)
 
C     Calcul de p0, q0, P, Q et b
      XMAMI = XMAX - XMIN
      XMAMIEPS = 0.00001D0*EEEN
      XMAMI = MAX(XMAMI,XMAMIEPS)
      XMAMIINV = EEEN/XMAMI
      UX1 = UPP - XVAL
      UX2 = UX1*UX1
      XL1 = XVAL - LOW
      XL2 = XL1*XL1
      UXINV = EEEN/UX1
      XLINV = EEEN/XL1
      P0 = ZERON
      Q0 = ZERON
      P0 = MAX(DF0DX,0.D0)
      Q0 = MAX((-1.D0)*DF0DX,0.D0)
      PQ0 = 0.001D0*(P0 + Q0) + RAA0*XMAMIINV
      P0 = P0 + PQ0
      Q0 = Q0 + PQ0
      P0 = P0*UX2
      Q0 = Q0*XL2
      CALL ZERO(P,M,N)
      CALL ZERO(Q,M,N)
      P = MAX(DFDX,0.0D0)
      Q = MAX((-1.D0)*DFDX,0.D0)
      PQ = 0.001D0*(P + Q) + RAA0*MATMUL(EEEM,TRANSPOSE(XMAMIINV))
      P = P + PQ
      Q = Q + PQ
      CALL MULMV(TEMP,P,UX2,M,N)
      P = TEMP
      CALL MULMV(TEMP,Q,XL2,M,N)
      Q = TEMP
      B = MATMUL(P,UXINV) + MATMUL(Q,XLINV) - FVAL
 
C     Resolution du sous probleme par une methode de Newton primale-duale
      MN = MIN(M,N)
      CALL SUBSOLVE(M,N,MN,EPSIMIN,LOW,UPP,ALFA,BETA,P0,Q0,P,Q,
     &              A0,A,B,C,D,
     &              XMMA,YMMA,ZMMA,LAM,XSI,ETA,MU,ZET,S)
 
      END