REAL*8 FUNCTION QUEUECOST(x,delta,nmove,nlane) CR Returns the total user cost of the queues IMPLICIT NONE REAL*8 x(28,17) INTEGER delta(28,17) INTEGER nmove,nlane INTEGER i,j REAL*8 tmp INTEGER maxln QUEUECOST = 0.D0 DO j = 1,nlane tmp = 0.D0 maxln = 0 DO i = 1,nmove tmp = tmp + x(i,j) IF (delta(i,j).GT.maxln) maxln = delta(i,j) ENDDO IF (maxln.GT.0) QUEUECOST = QUEUECOST + tmp**2/maxln ENDDO RETURN END REAL*8 FUNCTION DQUEUECOST(x,y,alpha,delta,nmove,nlane) CR Returns the derivative of the total user cost of the queues with CR respect to alpha at alpha IMPLICIT NONE REAL*8 x(28,17) REAL*8 y(28,17) REAL*8 alpha INTEGER delta(28,17) INTEGER nmove,nlane INTEGER i,j REAL*8 tmp,tmp2 INTEGER maxln DQUEUECOST = 0.D0 DO j = 1,nlane tmp = 0.D0 tmp2 = 0.D0 maxln = 0 DO i = 1,nmove IF (delta(i,j).GT.maxln) maxln = delta(i,j) tmp = tmp + (x(i,j)+alpha*(y(i,j)-x(i,j))) tmp2 = tmp2 + (y(i,j)-x(i,j)) ENDDO IF (maxln.GT.0) DQUEUECOST = DQUEUECOST + tmp*tmp2/maxln ENDDO DQUEUECOST = 2*DQUEUECOST RETURN END REAL*8 FUNCTION D2QUEUECOST(x,y,alpha,delta,nmove,nlane) CR Returns the derivative of the total user cost of the queues with CR respect to alpha at alpha IMPLICIT NONE REAL*8 x(28,17) REAL*8 y(28,17) REAL*8 alpha INTEGER delta(28,17) INTEGER nmove,nlane INTEGER i,j REAL*8 tmp INTEGER maxln D2QUEUECOST = 0.D0 DO j = 1,nlane tmp = 0.D0 maxln = 0 DO i = 1,nmove IF (delta(i,j).GT.maxln) maxln = delta(i,j) tmp = tmp + (y(i,j)-x(i,j)) ENDDO IF (maxln.GT.0) D2QUEUECOST = D2QUEUECOST + tmp**2/maxln ENDDO D2QUEUECOST = 2*D2QUEUECOST RETURN END SUBROUTINE FINDDIST(x,y,delta,nmove,nlane,d) CR Finds the best distribution of auxiliary flows given the current CR loading #include "dyna.inc" REAL*8 x(28,17) REAL*8 y(28,17) INTEGER delta(28,17) INTEGER nmove,nlane INTEGER d(28) REAL c(17) INTEGER i,j INTEGER minlane REAL*8 minval Cr - Clear y DO j = 1,nlane c(j) = 0 DO i = 1,nmove y(i,j) = 0 ENDDO ENDDO DO j = 1,nlane DO i = 1,nmove IF (delta(i,j).GT.0) THEN c(j) = c(j) + x(i,j)/delta(i,j) ENDIF ENDDO ENDDO DO i = 1,nmove minlane = 0 minval = 999999 DO j = 1,nlane IF (delta(i,j).NE.0.AND.c(j).LT.minval) THEN minval = c(j) minlane = j ENDIF ENDDO IF (minlane.EQ.0) THEN WRITE(ostr,'(A,I5,A)') + 'No minimum lane found for movement ',i,CHAR(0) CALL DYNA_ERROR(ostr + ,DYNA_FATAL_ERROR + ,DYNA_LOGIC_BUG + ,DYNA_NO_MIN_LANE) ENDIF y(i,minlane) = d(i) ENDDO RETURN END INTEGER FUNCTION LANEOPT(x,delta,nmove,nlane,d) CR Finds the optimal lane usage for a given demand IMPLICIT NONE REAL*8 x(28,17) REAL*8 y(28,17) INTEGER delta(28,17) INTEGER nmove,nlane INTEGER d(28) INTEGER i,j REAL*8 alpha REAL*8 z,zold REAL*8 dz,dz2 REAL*8 atol,ztol INTEGER niter,naiter INTEGER dsum REAL*8 QUEUECOST,DQUEUECOST,D2QUEUECOST atol = 1E-5 ztol = 1E-5 LANEOPT = 0 DO i = 1,nmove DO j = 1,nlane x(i,j) = 0 ENDDO ENDDO dsum = 0 DO i = 1,nmove dsum = dsum + d(i) ENDDO IF (dsum.LE.0) RETURN Cr - Initialize x CALL FINDDIST(x,y,delta,nmove,nlane,d) DO i = 1,nmove DO j = 1,nlane x(i,j) = y(i,j) ENDDO ENDDO Cr - Start optimization loop z = 999999.D0 zold = 1 niter = 0 DO WHILE(DABS(z-zold)/zold.GT.ztol) niter = niter + 1 cr DO i = 1,nmove cr WRITE(*,'(I5,17(:F7.1))') i,(x(i,j),j=1,nlane) cr ENDDO Cr - Direction finding (find y) CALL FINDDIST(x,y,delta,nmove,nlane,d) cr WRITE(*,*) 'Step direction' cr DO i = 1,nmove cr WRITE(*,'(I5,17(:F7.1))') i,(y(i,j),j=1,nlane) cr ENDDO Cr - Calculate step size by newton's method alpha = 0.347748D0 dz = atol+1 naiter = 0 DO WHILE(DABS(dz).GT.atol) dz = DQUEUECOST(x,y,alpha,delta,nmove,nlane) dz2 = D2QUEUECOST(x,y,alpha,delta,nmove,nlane) alpha = alpha - dz/dz2 naiter = naiter + 1 cr WRITE(*,*) ' alpha = ',alpha,' dz = ',dz ENDDO cr WRITE(*,*) '# alpha iter = ',naiter Cr - Update DO i = 1,nmove DO j = 1,nlane x(i,j) = (1-alpha)*x(i,j) + alpha*y(i,j) ENDDO ENDDO zold = z z = QUEUECOST(x,delta,nmove,nlane) cr WRITE(*,*) 'Current solution = ',z,' (alpha = ',alpha,')' cr DO i = 1,nmove cr WRITE(*,'(I5,17(:F7.1))') i,(x(i,j),j=1,nlane) cr ENDDO ENDDO LANEOPT = niter RETURN END C ---------------------------------------------------------------------- SUBROUTINE NODE_LANEOPT(sim,net,n) C ---------------------------------------------------------------------- C - Allocates current demand on approaches to node n to the available C - lanes using a UE formation C - INCLUDED FILES: #include "dyna.inc" #include "sim.inc" #include "network.inc" C - UNMODIFIED ARGUMENTS: RECORD /Sim_Data/ sim INTEGER n C - MODIFIED ARGUMENTS: RECORD /Road_Network/ net C - MODIFIED GLOBAL DATA: ! ostr (see dyna.inc) C - LOCAL VARIABLES: INTEGER i,j,k INTEGER l INTEGER np1,np2 INTEGER mv ! move pointer INTEGER mtype ! move type INTEGER mi ! UE move index INTEGER delta(4*MAX_MOVE_TYPE,MAX_LANE_TYPE) INTEGER x(4*MAX_MOVE_TYPE,MAX_LANE_TYPE) INTEGER d(4*MAX_MOVE_TYPE) C - FUNCTIONS CALLED: INTEGER MOVECNT REAL*8 QUEUECOST C - RETURN VALUE: ! Simulation status C ---------------------------------------------------------------------- np1 = net.fs.npoint(n) np2 = net.fs.npoint(n+1) - 1 C - Set up the delta, and d vectors for solution DO i = np1,np2 l = net.fs.ifwdarc(i,2) delta((i-np1)*MAX_MOVE_TYPE+1,1) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,2) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,3) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,4) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,5) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,6) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,7) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,8) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,9) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,10) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,11) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,12) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,13) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,14) = 1 delta((i-np1)*MAX_MOVE_TYPE+1,15) = 0 delta((i-np1)*MAX_MOVE_TYPE+1,16) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,1) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,2) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,3) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,4) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,5) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,6) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,7) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,8) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,9) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,10) = 0 delta((i-np1)*MAX_MOVE_TYPE+2,11) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,12) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,13) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,14) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,15) = 1 delta((i-np1)*MAX_MOVE_TYPE+2,16) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,1) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,2) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,3) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,4) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,5) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,6) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,7) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,8) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,9) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,10) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,11) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,12) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,13) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,14) = 0 delta((i-np1)*MAX_MOVE_TYPE+3,15) = 1 delta((i-np1)*MAX_MOVE_TYPE+3,16) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,1) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,2) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,3) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,4) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,5) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,6) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,7) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,8) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,9) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,10) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,11) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,12) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,13) = 0 delta((i-np1)*MAX_MOVE_TYPE+4,14) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,15) = 1 delta((i-np1)*MAX_MOVE_TYPE+4,16) = 1 delta((i-np1)*MAX_MOVE_TYPE+5,1) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,2) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,3) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,4) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,5) = 1 delta((i-np1)*MAX_MOVE_TYPE+5,6) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,7) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,8) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,9) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,10) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,11) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,12) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,13) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,14) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,15) = 0 delta((i-np1)*MAX_MOVE_TYPE+5,16) = 0 C - Zero all the d's DO j = 1,MAX_MOVE_TYPE mi = MAX_MOVE_TYPE*(i-np1)+j d(mi) = 0 ENDDO C - How many lanes of this type on on this approach DO k = 1,MAX_LANE_TYPE DO j = 1,MAX_MOVE_TYPE mi = MAX_MOVE_TYPE*(i-np1)+j delta(mi,k) = + delta(mi,k)*net.node(n).lanenum(i-np1+1,k) ENDDO ENDDO C - Set up the demand for this approach DO j = 1,net.link(l).numdslinks mv = net.link(l).dsmoveptr(j) mtype = net.movement(mv).type mi = MAX_MOVE_TYPE*(i-np1)+mtype d(mi) = MOVECNT(sim,net,l,mtype) ENDDO ENDDO CALL LANEOPT(x,delta,MAX_MOVE_TYPE,MAX_LANE_TYPE,d) WRITE(ostr,600) net.node(n).number,i,CHAR(0) 600 FORMAT('NODE: 'I5' # of iter = 'I5,A) CALL SIMMSG(STATUS,ostr) WRITE(ostr,601) QUEUECOST(x,delta,MAX_MOVE_TYPE,MAX_LANE_TYPE,d) + ,CHAR(0) 601 FORMAT('Solution cost = ',F10.2,A) CALL SIMMSG(STATUS,ostr) DO i = 1,MAX_MOVE_TYPE WRITE(ostr,'(I1,17(F5.0),A)') i,(x(i,j),j=1,MAX_LANE_TYPE) + ,CHAR(0) CALL SIMMSG(STATUS,ostr) ENDDO END c$$$ PROGRAM TEST c$$$ c$$$ IMPLICIT NONE c$$$ c$$$ REAL*8 x(28,17) c$$$ INTEGER delta(28,17) c$$$ INTEGER d(28) c$$$ INTEGER nmove,nlane c$$$ INTEGER i,j c$$$ c$$$ INTEGER LANEOPT c$$$ REAL*8 QUEUECOST c$$$ c$$$ READ(*,*) nmove,nlane c$$$ READ(*,*) ((delta(i,j),j=1,nlane),i=1,nmove) c$$$ READ(*,*) (d(i),i=1,nmove) c$$$ c$$$ i = LANEOPT(x,delta,nmove,nlane,d) c$$$ c$$$ WRITE(*,*) '# of iterations = ',i c$$$ WRITE(*,*) 'Solution cost = ',QUEUECOST(x,delta,nmove,nlane,d) c$$$ DO i = 1,nmove c$$$ WRITE(*,'(I5,17(:F7.1))') i,(x(i,j),j=1,nlane) c$$$ ENDDO c$$$ c$$$ END c$$$