optBathtub/continuous/f77/massn.f

      PROGRAM MASS
C This program computes optimum mass entries in the bathtub model
C 
C The main program and subroutine are in this file. Routines from
C reading dat input and displaying output are in put in file
C "auxiliary.f".
C 
C The main program call a first subroutines to read data input and then
C calls the subroutine that computes the optimum entry rates.
C
C A SUMMARY OF THE SUBROUTINES 
C
C UTIL      : computes the objective function
C D_UTIL    : computes the gradient of the objective function
C CONS      : computes the value of the constraint. This version 
C             has one constraint (the sum of the entry rates is equal to
C             one
C D_CONS    : computes the gradient of the constraint
C SPLITLINE : for input, it takes the current approximation line with
C             respect variable m and breaks and evenly subdivides the largest
C             intervals into two intervals. Bu doing so we intend to get
C             a smoother approximation of the solution
C
C For the optimization step we use the fortran version IPOPT (the large-scale 
C optimization package): 
C
C https://www.coin-or.org/Ipopt/ipopt-fortran.html
C 
C TO RUN THIS PROGRAM

C - install IPOPT (it requires LAPACK and BLAS)
C - compile through a command like 
C       gfortran -o massn massn.f auxiliary.f -llapack -lblas -lipopt  
C   by providing the correct path to the libraries
C - set data in file "params.dat"
C - run the program
C       ./massn
C
C The process is little tedious, in particular for IPOPT.
C
C
C KNOWN ISSUES AND POSSIBLE IMPROVEMENTS
C
C 1. This version does not explicitly set constraint on :
C    - last arrival time should be smaller than TMAX
C    - traffic density should be smaller than road jam capacity
C    and relies on the optimization to carry on implicitly these
C    constraints. This does not occur in all cases and fine tuning of
C    parameter M0, MBAR and ITMAX is frequently required for convergence.
C
C 2. M0 and MBAR are fixed. It is suggested to provide a large window in
C    the start. The solution will entail a first phase where without
C    entries. Similarly it will entail a last phase where entries are
C    equal to zero. A new of the computation can then be considered by
C    zooming a new window with higher M0 and smaller MBAR. This process
C    can be repeated until satisfactory values are obtained.
C
C 3. The Hessian is not provided to the optimizer. It uses finite
C    difference to approximate it (IQUASI = 6, in the options passed to
C    IPOPT).
C
C 4. The program does not produce any graphical output. For instance,
C    output table can be exported to a text file (redirect the output
C    using the unix standard "./massn > outfile") and then use gnuplot
C    or another program to make the plots (after deleting useless
C    lines).
C
C These issues will be addressed in the updates
C
C
C NOTATION THROUGH AN EXAMPLE   
C  - number of entry (and exit) points: N = 10
C  - number of entry points before within a cycle:  NL = 2
C  - total number of breakpoints: N + NL = 12
C
C (0)    t1     t2    t3     t4    t5     t6     t7    t8     t9   t10    t11      
C (1)  P----S-------P----S-------P----S-------P-----S-------P----S-----P--------S--- 
C
C (2)  0    Mb-4L   L    Mb-3L   2L   Mb-2L   3L    Mb-L    4L   Mb    (4+1)L   Mb+L   
C
C (3)  1    2       3    4       5    6       7     8       9    10    11       12 
C
C (4)  e1   e2      e3   e4      e5   e6      e7    e8      e9   e10   
C
C (5)               e1   e2      e3   e4      e5    e6      e7   e8    e9       e10   
C 
C (6)  T1   T2      T3   T4      T5   T6      T7    T8      T9   T10   T11      T12
C
C
C  (0) t_i clock-time from m_i to m_(i+1)
C  (1) Primary (P) and secondary (S) breakpoints
C  (3) indices for variable m
C  (4) entries at each breakpoint m_i (at each breakpoint)
C  (5) exits at each breakpoint m_i (at each breakpoint)
C  (6) T_i clock-time at mile m_i 

      IMPLICIT NONE
      
      DOUBLE PRECISION UT
      DOUBLE PRECISION R0, R1, S0, S1, VF, K, L, P, TMAX 
      DOUBLE PRECISION M0, MBAR
      INTEGER          ITMAX
      COMMON /MODELPARAM/  R0, R1, S0, S1, VF, K, L, P, TMAX 
      COMMON /OTHERPARAM/  M0, MBAR, ITMAX

C Read parameters from file "params.dat"
      CALL READPARAM("params.dat")
      PRINT '(11F9.2)', R0,R1,S0,S1,VF,K,L,P,TMAX,M0,MBAR
C Call subroutine UTILOPT to compute the optimum entry rates
      CALL UTILOPT ( M0, MBAR , UT, ITMAX)

C Print the value of the objective function
      PRINT*, UT

      END


      SUBROUTINE UTILOPT (M0,MBAR,F,ITMAX)
C Computes entry masses at points m_i that maximize aggregate
C utility given that M0 and MBAR are fixed
      IMPLICIT NONE

      DOUBLE PRECISION M0, MBAR
      INTEGER          I, NCYC, NL, NMAX, IT, ITMAX
      PARAMETER        (NMAX=10000)
      DOUBLE PRECISION U

C Parameters of the problem
      DOUBLE PRECISION R0, R1, S0, S1, VF, K, L, P, TMAX 

C Variable definition for IPOPT
      INTEGER          LRW,  LIW
      PARAMETER        ( LRW = 10000, LIW = 10000)
      INTEGER          IW(LIW)
      DOUBLE PRECISION RW(LRW)
      DOUBLE PRECISION IPOPT
      INTEGER          N, M
      DOUBLE PRECISION LAM(NMAX)
      DOUBLE PRECISION C(NMAX)
      DOUBLE PRECISION G(NMAX)
      DOUBLE PRECISION E(NMAX)
      INTEGER          NLB, NUB
      INTEGER          ILB(NMAX), IUB(NMAX)
      DOUBLE PRECISION BNDSL(NMAX), BNDSU(NMAX)
      DOUBLE PRECISION V_L(NMAX), V_U(NMAX)
      DOUBLE PRECISION DAT  (NMAX)
      INTEGER          IDAT (NMAX)
      INTEGER          NARGS
      DOUBLE PRECISION ARGS (50)
      CHARACTER*20     CARGS(50)
      INTEGER          IERR, ITER
      DOUBLE PRECISION F

      EXTERNAL       UTIL, CONS, D_UTIL, D_CONS
      EXTERNAL       EV_H_DUMMY,EV_HLV_DUMMY, EV_HOV_DUMMY, EV_HCV_DUMMY
C End of definition of IPOPT variables

      COMMON /MODELPARAM/  R0, R1, S0, S1, VF, K, L, P, TMAX
 
      NCYC= FLOOR ( (MBAR-M0) / L )  
      N = 2 * (  NCYC+ 1 )
      IF(MBAR.GT.L) THEN
        NL = 2
      ELSE
        NL =1
      ENDIF
      IDAT(1) = NCYC
      IDAT(2) = NL

C Compute the breakpoints. This is the first approximation that will be
C used.
      DO 10, I = 1,  NCYC+ NL
C Primary breakpoints
        DAT(2*I-1) = M0   + DBLE( I - 1 ) * L
C Secondary breakpoints. Notice that : 
C       MBAR + L -  NCYCL + (I-1) * L = MBAR - (NCYC-I) * L 
        DAT(2*I)   = MBAR - DBLE (  NCYC- I + 1 ) * L
 10   CONTINUE

C Compute a simple initial solution
      DO 20, I = 1, N
        E(I) = ( DAT(I+1) - DAT(I) ) / ( (  NCYC + 1 ) * L ) 
 20   CONTINUE

C Set algorithmic parameters (for IPOPT)
       NARGS   = 1
       ARGS(1) = 1.D-8
       CARGS   = 'dtol' 


C MAIN LOOP
C It starts by optimizing over the initial breakpoints, then 
C uses the obtained solution to increase the number of approximation
C points. Variable ITMAX indicates the number of iterations (how many
C time the approximation line) is subdivided.

      IT = 1
 100  CONTINUE

C Number of variable with lower bounds
      NLB = N 
C Indices of variables with lower bounds
      DO 1, I = 1, N
      ILB(I) = I
 1    CONTINUE
C Values of lower bounds
      DO 2, I = 1, N
      BNDSL(I) = 1.0D-6
 2    CONTINUE
C Number of variables with upper bounds
      NUB = 0
C i.e. no uppers bounds on the entry rates

C There is only one constraint on the problem (the sum of all entry
C rates is one)
      M = 1
      PRINT* , N, NL, IT

C  Call IPOPT (the optimizer)
      F = - IPOPT(N, E, M, NLB, ILB, BNDSL, NUB, IUB, BNDSU, V_L, V_U,
     1   LAM, C, LRW, RW, LIW, IW, ITER, IERR, UTIL, CONS, D_UTIL,
     1   D_CONS, EV_H_DUMMY, EV_HLV_DUMMY, EV_HOV_DUMMY, EV_HCV_DUMMY,
     1   DAT, IDAT, NARGS, ARGS, CARGS)

      PRINT '(A20,F9.3)' , 'Value of objective function : ', F 

       PRINT '(F9.5)', SUM(E) 

      IT = IT + 1 

      IF(IT.LE.ITMAX) THEN
        CALL SPLITLINE(N, DAT,NL,E,0.80D0)
        IDAT(2) = NL
        GOTO 100
      ENDIF

C Print the last solution obtained
      CALL STATDES(N, E, DAT, IDAT)

      RETURN
      END

C ======================================
C Computation of the objective function
C (aggregate utility)
C ======================================
      SUBROUTINE UTIL(N, E, U, DAT, IDAT)
      IMPLICIT NONE
      INTEGER          I, N, NL, NMAX
      PARAMETER        (NMAX=10000)
      DOUBLE PRECISION U, E(N)
      DOUBLE PRECISION DAT(NMAX)
      INTEGER          IDAT(NMAX)
C Variables used for intermediate computations
      DOUBLE PRECISION V(N+IDAT(2)-1), T(N+IDAT(2)-1), TT(N+IDAT(2))
      DOUBLE PRECISION UE, UX
      DOUBLE PRECISION SUME
C Parameters of the problem
      DOUBLE PRECISION R0, R1, S0, S1, VF, K, L, P, TMAX 
      COMMON /MODELPARAM/  R0, R1, S0, S1, VF, K, L, P, TMAX

      NL = IDAT(2) 

C Compute travel speed, travel time and clock-time on each interval
      I=1
      TT(I) = DAT(I) / VF
      SUME=E(I)
      V(I) = VF * ( 1.0D0 - P * SUME / K )
      T(I) = ( DAT(I+1) - DAT(I) ) / V(I)
C     PRINT*, V(I), VF, P, K, SUME,T(I),TT(I),DAT(2)- DAT(1)
      TT(I+1) = TT(I) + T(I)
      DO 1, I = 2, N + NL - 1
        IF(I.LE.NL) THEN
          SUME = SUME + E(I)
        ELSEIF (I.LE.N) THEN
          SUME = SUME + E(I) - E(I-NL)
        ELSE
          SUME = SUME        - E(I-NL)
        ENDIF
        V(I) = MAX(0.D0,VF * ( 1.0D0 - P * SUME / K ))
        T(I) = ( DAT(I+1) - DAT(I) ) / V(I)
        TT(I+1) = TT(I) + T(I)
C       PRINT '(I5,5F9.3)', I,DAT(I),V(I), T(I), TT(I+1), SUME
 1    CONTINUE 

C Compute entry and exit sub-utilities
      UE = 0.0D0 
      UX = 0.0D0 
      DO 20, I = 1, N
         UE = UE + R0 * DLOG( R1 *          TT(I     ) ) * E(I) * P
         UX = UX + S0 * DLOG( S1 * ( TMAX - TT(I + NL))) * E(I) * P
C     PRINT '(2I5,4F9.3)', I, NL,TMAX,TT(I+NL),UE, UX
 20   CONTINUE


C AGREGATE UTILITY 
C Notice that we consider the opposite of the utility since we set the
C problem as a minimization program
      U = - ( UE + UX )  


      RETURN
      END


C =====================================================
C Computation of the gradient of the objective function
C =====================================================
      SUBROUTINE D_UTIL( N, E, G, DAT, IDAT)
      IMPLICIT NONE

      INTEGER          N, NL, NMAX
      PARAMETER        (NMAX=10000)
      DOUBLE PRECISION G(N), E(N)
      DOUBLE PRECISION DAT(NMAX)
      INTEGER          IDAT(NMAX)
      DOUBLE PRECISION D_T(N+IDAT(2)-1), D_TT(N+IDAT(2),N)
      DOUBLE PRECISION D_UE(N,N), D_UX(N,N)
      DOUBLE PRECISION V(N+IDAT(2)-1), T(N+IDAT(2)-1), TT(N+IDAT(2))
      DOUBLE PRECISION SUME
      INTEGER          I, J, II
C Parameters of the problem
      DOUBLE PRECISION R0, R1, S0, S1, VF, K, L, P, TMAX 
      COMMON /MODELPARAM/  R0, R1, S0, S1, VF, K, L, P, TMAX
      
      NL = IDAT(2)
C Compute travel speed on each interval
C This step is replicated from UTIL subroutine.
      I = 1
      TT(I)  = DAT(I) / VF
      SUME   = E(I)
      V(I)   = MAX(0.0D0,VF * ( 1.D0 - P * SUME / K ))
      T(I)   = ( DAT(I+1) - DAT(I) ) / V(I)
      TT(I+1) = TT(I) + T(I)
C     PRINT*, V(I), VF, P, K, SUME,T(I),TT(I),DAT(2)- DAT(1)
      DO 1, I = 2, N + NL - 1
        IF(I.LE.NL) THEN
          SUME = SUME + E(I)
        ELSEIF (I.LE.N) THEN
          SUME = SUME + E(I) - E(I-NL)
        ELSE
          SUME = SUME        - E(I-NL)
        ENDIF
        V(I) = VF * ( 1.D0 - P * SUME / K )
        T(I) = ( DAT(I+1) - DAT(I) ) / V(I)
        TT(I+1) = TT(I) + T(I)
C     PRINT*, V(I), VF, P, K, SUME,T(I),TT(I),DAT(I+1), DAT(I)
 1    CONTINUE 

C Compute the derivative
C
C   d t_i    m_(i+1) - m_i     vf * P
C   ----- =  -------------  *  ------
C   d e_i        v_i^2           K
C
      DO 10, I = 1, N + NL - 1
      D_T(I) = (DAT(I+1) - DAT(I))/V(I)**2 * VF*P/K
 10   CONTINUE 
      
C Compute the derivative
C
C                  /
C                 |   0                            if  j <= i
C                 |      
C                 |
C      d T_j      |  d t_i           d t_(j-1)
C      -----  =  /   -----  + ... +  ---------     if  i < j <=i+NL 
C      d e_i     \   d e_i           d e_(j-1)
C                 |      
C                 |
C                 |  d T_(i+NL)    
C                 |  ----------                    if  j > i + NL
C                 |   d e_i         
C                  \
C 
      
      D_TT = 0.0D0
      DO 20, I = 1, N
      DO 20, J = I+1, N+NL
        IF (J.LE.I) THEN
          D_TT(J,I) = 0.0D0
        ELSEIF (J.LE.I+NL) THEN
          DO 22, II = I, J-1
          D_TT(J,I) = D_TT(J,I) + D_T(II)
22        CONTINUE
        ELSE
          D_TT(J,I) = D_TT(I+NL,I)
        ENDIF 
20    CONTINUE 

C IMPACTS ON ENTRY SUBUTILITIES
C An increase of the entry rate at T_i increases the size of the group
C itself and delays the entry time for all groups that enter at T_j with
C j > i.
C Users who enter at T_N do not delay the entry of any other group. The
C structure of the loop automatically takes into account this case since.
C Indeed, for I=N, J varies varies from N+1 to N and, then, the loop is
C not executed.
      D_UE = 0.0D0
      DO 30, I =   1, N 
      D_UE(I, I) = R0 * DLOG( R1 * TT(I) ) * P
      DO 30, J = I+1, N
      D_UE(J, I) = R0 * D_TT(J,I) / TT(J) * E(J) * P 
30    CONTINUE 

C IMPACTS ON EXIT SUBUTILITIES
C An increase of the size of the group of users entering at T_i
C has two impacts: (1) it delays the arrival of all groups entering at 
C T_{max(1,i-(NL-1))}, including the group itself, and (2) it increases
C the size of the group of those entering at T_i. 
      D_UX = 0.0D0 
      DO 40, I  = 1, N 
      DO 40, J  = MAX(I-(NL-1), 1), N
      D_UX(J, I) = - S0 * D_TT(J+NL,I) / (TMAX - TT(J+NL)) * E(J) * P 
40    CONTINUE
      DO 45, I  = 1, N
      D_UX(I, I) = D_UX(I, I) + S0 * DLOG( S1 * (TMAX - TT(I+NL)) ) * P
45    CONTINUE 

C TOTAL IMPACTS (the gradient)
C Notice that we compute the gradient of - UTIL since we set
C the problem as a minimization program.
      G = 0.0D0
      DO 50, I = 1, N
      DO 50, J = 1, N 
      G(I) = G(I) - D_UE(J,I) - D_UX(J,I)
C     PRINT '(7F9.3)', G(I), E(1:6) 
 50   CONTINUE 

      RETURN
      END 


C =========================================
C Computation of the equality constraints
C In this problem, there is only one constraint
C stating that the sum of all entry rates 
C is equal to one :
C
C      e_1 + e_2 + .. + e_n = 1 
C
C ========================================= 
      SUBROUTINE CONS(N, E, M, C, DAT, IDAT)
      IMPLICIT NONE

      INTEGER          N, M, NMAX
      PARAMETER (NMAX=10000)
      DOUBLE PRECISION C(M), E(N)
      DOUBLE PRECISION DAT(NMAX)
      INTEGER          IDAT(NMAX)


      DOUBLE PRECISION SUME
      INTEGER          I

      SUME = 0.D0
      DO 10, I = 1, N
        SUME = SUME + E(I)
 10   CONTINUE

      C(1) = SUME - 1.D0

      RETURN
      END
 

C =============================================
C Computation of the Jacobian of the constraint
C =============================================
      SUBROUTINE D_CONS(TASK, N, E, NZ, A, ACON, AVAR, DAT, IDAT)
      IMPLICIT NONE
      INTEGER          TASK, N, NZ
      INTEGER          NMAX, I
      PARAMETER (NMAX=10000)
      DOUBLE PRECISION E(N), A(NZ)
      INTEGER          ACON(NZ), AVAR(NZ)
      DOUBLE PRECISION DAT(NMAX)
      INTEGER          IDAT(NMAX)
      
      IF(TASK.EQ.0) THEN
        NZ = N
      ELSE
        DO 10, I = 1, N
        AVAR(I) = I
        ACON(I) = 1
10      CONTINUE 
        DO 20, I = 1, N
        A(I) = 1.D0
20      CONTINUE
      ENDIF

      RETURN
      END


C =============================================
C Split interval
C =============================================
      SUBROUTINE SPLITLINE(N, M, NL, E, THR)
      IMPLICIT NONE
C Increases the number of approximation points by breaking the largest
C intervals.
C N   the number of entry rates, or mile-points where entry occurs
C M   table containing the approximation points
C NL  index the last element between0 and L: we have M(NL+1) = L  
C THR intervals with length larger than THR * X are swplit into two
C     intervals, where X is the the length of the largest interval

      INTEGER          NL, N, I, J, NMAX
      PARAMETER        (NMAX=10000) 
      DOUBLE PRECISION M(NMAX), E(NMAX)
      DOUBLE PRECISION THR
      DOUBLE PRECISION M2(NMAX), E2(NMAX)
      INTEGER          N2,NPTS,NL2,NTOT
      DOUBLE PRECISION Y(NMAX), YMAX

      IF(N.LT.2) THEN
        PRINT*, "A table of legth 2 at least is expected"
        STOP
      ENDIF
      IF(THR.LE.0.0D0.OR.THR.GE.1.0D0) THEN
        PRINT*, "THR should be between 1/2 and 1."
        STOP
      ENDIF 


C Compute the length of the largest interval
      YMAX = M(2) - M(1)
      DO 10, I = 1, NL
        Y(I) = M(I+1) - M(I)
        IF(Y(I).GT.YMAX) YMAX = Y(I)
 10   CONTINUE 

C Store some values
      N2 = N
      NL2 = NL
      NTOT = N + NL
      E2 = 0.0D0

C Split intervals with length higher than
C THR * YMAX and update indexes
      J = 0
      DO 20, I = 1, N+NL-1 
        IF((M(I+1)-M(I)).GT.(THR*YMAX)) THEN
          M2(I+J  ) = M(I)
          M2(I+J+1) = M(I) + (M(I+1) - M(I)) / 2.0D0
          J  = J + 1
          IF(I.LE.NL) THEN
            E2(I+J  ) = E(I)
            E2(I+J-1) = 0.0D0
            NL2  = NL2 + 1
            N2   = N2  + 1
            NTOT = NTOT + 1
          ELSEIF(I.LT.N) THEN
            E2(I+J  ) = E(I)
            E2(I+J-1) = 0.0D0
            N2   = N2  + 1
            NTOT = NTOT + 1
          ELSE
            NTOT = NTOT + 1
          ENDIF
        ELSE
          M2(I+J) = M(I)
          E2(I+J) = E(I)
        ENDIF 
 20   CONTINUE 
      E2(N2  ) = E(N)
      M2(NTOT) = M(N+NL)

C Set valriables to new values
      M  = M2
      E  = E2
      N  = N2
      NL = NL2


      RETURN
      END