Télécharger la plaquette Cast3M


fr en

Télécharger mmasub.eso

Retour à la liste

Numérotation des lignes : 

mmasub

  1. C MMASUB    SOURCE    FD218221  26/10/07    21:15:08     12662          
  2.       SUBROUTINE MMASUB(M,N,ITER,XVAL,XMIN,XMAX,XOLD1,XOLD2,
  3.      &                  DF0DX,FVAL,DFDX,MOVE,ALBEFA,RAA0,
  4.      &                  ASYINIT,ASYINCR,ASYDECR,ASYMIN,ASYMAX,
  5.      &                  LOW,UPP,ALFA,BETA,P0,Q0,P,Q,B)
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8 (A-H,O-Z)
  8.  
  9. C ----------------------------------------------------------------------
  10. C Ce code est une adaptation en FORTRAN/ESOPE de l'algorithme de la MMA
  11. C (Methode des Asymptotes Mobiles) propose initialement en langage
  12. C Matlab par K. Svarberg (sources : https://www.smoptit.se/)
  13. C ----------------------------------------------------------------------
  14. C    Cette subroutine definit un sous-probleme d'optimisation convexe
  15. C    selon la methode des asymptotes mobiles (MMA)
  16. C
  17. C    Le probleme d'optimisation initial est le suivant :
  18. C    Minimiser : f_0(x) + a_0*z + SUM[ c_i*y_i + 0.5*d_i*(y_i)^2 ]
  19. C     soumis a : f_i(x) - a_i*z - y_i <= 0,     i = 1,...,m
  20. C                xmin_j <= x_j <= xmax_j,       j = 1,...,n
  21. C                z >= 0,   y_i >= 0,            i = 1,...,m
  22. C
  23. C    Le probleme est transforme en une suite de sous-problemes convexes :
  24. C    Minimiser :   SUM[ p0_j/(upp_j-x_j) + q0_j/(x_j-low_j) ] + a_0*z +
  25. C                + SUM[ c_i*y_i + 0.5*d_i*(y_i)^2 ]
  26. 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
  27. C                  alfa_j <= x_j <= beta_j,  j = 1,...,n
  28. C                  z >= 0,   y_i >= 0,       i = 1,...,m
  29. C
  30. C    Cette subroutine calcule les matrices/vecteurs decrivant un sous-probleme
  31. C    Elle necessite d'etre iteree pour approcher la solution du probleme initial
  32. C
  33. C  Entrees :
  34. C  ---------
  35. C  M         = Nombre de contraintes
  36. C  N         = Nombre de variables x_j
  37. C  ITER      = Numero de l'iteration courante (=1 la premiere fois que mmasub est appellee)
  38. C  XVAL      = Vecteur (N) des variables x_j
  39. C  XMIN/XMAX = Vecteurs (N) des bornes inferieures/superieures pour les variables x_j
  40. C  XOLD1     = XVAL, a 1 iteration  precedente  (a condition que iter>1)
  41. C  XOLD2     = XVAL, a 2 iterations precedentes (a condition que iter>2)
  42. C  DF0DX     = Vecteur (N) des valeurs du gradient de la fonction objectif f_0,
  43. C              par rapport aux variables x_j, calcules en XVAL
  44. C  FVAL      = Vecteur (M) des valeurs des fonctions contraintes f_i, calculees en XVAL
  45. C  DFDX      = Matrice (M x N) des valeurs du gradient des fonctions contraintes f_i,
  46. C              par rapport aux variables x_j, calcules en XVAL
  47. C              DFDX(i,j) = gradient de f_i par rapport a x_j
  48. C  MOVE      = Taille maximale de l'increment des variables
  49. C  ALBEFA    = Pourcentage d'ecart par rapport aux asymptotes alfa et beta
  50. C  RAA0      = Parametre de conservatisme / courbure de l'objectif
  51. C  ASYINIT   = Facteur d'initialisation des asymptotes LOW et UPP
  52. C  ASYINCR   = Facteur d'eloignement des asymptotes LOW et UPP
  53. C  ASYDECR   = Facteur de rapprochement des asymptotes LOW et UPP
  54. C  ASYMIN    = Limite inferieure de la distance entre les asymptotes LOW et UPP
  55. C  ASYMAX    = Limite superieure de la distance entre les asymptotes LOW et UPP
  56. C  LOW/UPP   = Vecteurs (N) des asymptotes inferieures/superieures du sous-probleme
  57. C  ALFA/BETA = Vecteurs (N) des bornes inferieures/superieures du sous-probleme
  58. C
  59. C  Sorties :
  60. C  ---------
  61. C  LOW/UPP   = Vecteurs (N) des asymptotes inferieures/superieures mises a jour
  62. C  ALFA/BETA = Vecteurs (N) des bornes inferieures/superieures mises a jour
  63. C  P0/Q0     = Vecteurs (N) des coefficients des asymptotes superieure/inferieure
  64. C              dans la fonction objectif
  65. C  P/Q       = Matrices (M x N) des coefficients des asymptotes superieure/inferieure
  66. C              des contraintes
  67. C  B         = Vecteur (M) des seconds membres des contraintes
  68. C ----------------------------------------------------------------------
  69.  
  70. C     Arguments
  71.       INTEGER M,N,ITER
  72.       REAL*8  MOVE,ALBEFA,RAA0,ASYINIT,ASYINCR,ASYDECR,ASYMIN,ASYMAX
  73.       REAL*8  XVAL(N),XMIN(N),XMAX(N),XOLD1(N),XOLD2(N)
  74.       REAL*8  DF0DX(N),FVAL(M),DFDX(M,N)
  75.       REAL*8  LOW(N),UPP(N),ALFA(N),BETA(N)
  76.       REAL*8  P0(N),Q0(N),P(M,N),Q(M,N),B(M)
  77.  
  78. C     Variables locales
  79.       INTEGER I,J
  80.       REAL*8  ASYMP,FACT,DF0,DFJ,LOWMIN,LOWMAX,UPMIN,UPMAX,XMAMI
  81.  
  82. C     Mise a jour des asymptotes low et upp
  83.       DO 10 I = 1, N
  84.         XMAMI = XMAX(I) - XMIN(I)
  85.         IF (ITER .LE. 2) THEN
  86.           LOW(I) = XVAL(I) - ASYINIT * XMAMI
  87.           UPP(I) = XVAL(I) + ASYINIT * XMAMI
  88.         ELSE
  89.           FACT = (XVAL(I) - XOLD1(I)) * (XOLD1(I) - XOLD2(I))
  90.           IF (FACT .GT. 0.0D0) THEN
  91.             ASYMP = ASYINCR
  92.           ELSEIF (FACT .LT. 0.0D0) THEN
  93.             ASYMP = ASYDECR
  94.           ELSE
  95.             ASYMP = 1.0D0
  96.           ENDIF
  97.           LOW(I) = XVAL(I) - ASYMP * (XOLD1(I) - LOW(I))
  98.           UPP(I) = XVAL(I) + ASYMP * (UPP(I) - XOLD1(I))
  99.           LOWMIN = XVAL(I) - ASYMAX * XMAMI
  100.           LOWMAX = XVAL(I) - ASYMIN * XMAMI
  101.           UPMIN  = XVAL(I) + ASYMIN * XMAMI
  102.           UPMAX  = XVAL(I) + ASYMAX * XMAMI
  103.           LOW(I) = MAX(LOW(I),LOWMIN)
  104.           LOW(I) = MIN(LOW(I),LOWMAX)
  105.           UPP(I) = MAX(UPP(I),UPMIN)
  106.           UPP(I) = MIN(UPP(I),UPMAX)
  107.         ENDIF
  108.  10   CONTINUE
  109.  
  110. C     Mise a jour des bornes alfa et beta
  111.       DO 20 I = 1, N
  112.         ALFA(I) = MAX(LOW(I) + ALBEFA * (XVAL(I) - LOW(I)),
  113.      &                XVAL(I) - MOVE * (XMAX(I) - XMIN(I)),
  114.      &                XMIN(I))
  115.         BETA(I) = MIN(UPP(I) - ALBEFA * (UPP(I) - XVAL(I)),
  116.      &                XVAL(I) + MOVE * (XMAX(I) - XMIN(I)),
  117.      &                XMAX(I))
  118.  20   CONTINUE
  119.  
  120. C     Calcul de p0 et q0
  121.       DO 30 I = 1, N
  122.         DF0 = DF0DX(I)
  123.         XMAMI = MAX(XMAX(I)-XMIN(I), 0.00001D0)
  124.         P0(I) = MAX( DF0,0.D0) + 0.001D0*ABS(DF0) + RAA0/XMAMI
  125.         Q0(I) = MAX(-DF0,0.D0) + 0.001D0*ABS(DF0) + RAA0/XMAMI
  126.         P0(I) = P0(I) * ((UPP(I) - XVAL(I))**2)
  127.         Q0(I) = Q0(I) * ((XVAL(I) - LOW(I))**2)
  128.  30   CONTINUE
  129.  
  130. C     Calcul de p, q et b
  131.       DO 40 J = 1, M
  132.         B(J) = -FVAL(J)
  133.         DO 50 I = 1, N
  134.           DFJ = DFDX(J,I)
  135.           XMAMI = MAX(XMAX(I)-XMIN(I), 0.00001D0)
  136.           P(J,I) = MAX( DFJ,0.D0) + 0.001D0*ABS(DFJ) + RAA0/XMAMI
  137.           Q(J,I) = MAX(-DFJ,0.D0) + 0.001D0*ABS(DFJ) + RAA0/XMAMI
  138.           P(J,I) = P(J,I) * ((UPP(I)-XVAL(I))**2)
  139.           Q(J,I) = Q(J,I) * ((XVAL(I)-LOW(I))**2)
  140.           B(J) = B(J) + P(J,I) / (UPP(I) - XVAL(I)) 
  141.      &                + Q(J,I) / (XVAL(I) - LOW(I))
  142.  50     CONTINUE
  143.  40   CONTINUE
  144.  
  145.       RETURN
  146.       END
  147.  
  148.  
  149.  
  150.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales