1. %-------------------------------------------------------------------------%
  2. %%----------------Hito 2: Simulación de tráfico--------------------------%%
  3. %-------------------------------------------------------------------------%
  4.  
  5. % ---- Parametros constantes ---- %
  6. qc = 3000;       % ----- flujo critico          
  7. vc = 120;        % ----- vel. maxima/crític        
  8. kc = 25;         % ----- dens. critica (qc/vc)          
  9. kmax = 100;      % ----- dens. máxima
  10. kmax_via = 60;   % ----- dens. max en la via
  11. L = 5;           % ----- longitud por celda
  12. numCeldas = 5;   % ----- numero de celdas (No confundir con longitud)
  13. At = 1/120;      % ----- incremento de tiempo (horas)
  14. t = 7:At:10;     % ----- vector de tiempo
  15. f = length(t);   % ----- longitud del vector t
  16.  
  17. % ---- Iniciar matrices ---- %
  18. q0 = zeros(1,f);         % ----- flujo auxiliar
  19. Q0 = zeros(1,f);         % ----- flujo auxiliar
  20. k0 = zeros(1,f);         % ----- dens. auxiliar
  21. n = zeros(numCeldas,f);  % ----- num vehiculos    
  22. k = zeros(numCeldas,f);  % ----- densidad vehiculos
  23. q = zeros(numCeldas,f);  % ----- flujo de vehiculos saliente (real)
  24. Q = zeros(numCeldas,f);  % ----- flujo de vehiculos (teorico)
  25.  
  26. %% --- Simulacion 1: sin semaforos --- %
  27. for j = 1:f                   % Avance del tiempo
  28.     for i = 1:numCeldas       % Avance en secciones de via
  29.         k(i,j) = n(i, j)/L;   % Calculo densidad
  30.         Q(i,j) = funcTrafico(k(i,j), kc, kmax, qc); % calculo flujo
  31.         if i < 5   % flujo en cada celda
  32.             q(i,j) = min(Q(i,j), (kmax_via*L - n(i+1,j)) / At); % tras imponer restricción de flujo por otras celdas, obtenemos el flujo real
  33.         else
  34.             q(i,j) = Q(i,j); % cel 5 no esta limitada por la cap. de otra celda, por lo que se simplifica a esto
  35.         end
  36.        
  37.         % para la celda auxiliar (la celda "cero") usaremos variables independientes de las otras usadas, como sugiere la doc.
  38.         if i == 1 % la celda ficticia alimenta a la 1a celda real, empezamos por esa celda
  39.             qentrada(j) = 4000*exp(-(t(j) - 8) ^2);  % función de demanda conocida como e(t) (se comporta en función de un pulso gaussiano)
  40.             if k0(j) > 0
  41.                 Q0(j) = qentrada(j) + k0(j)*L / At; % si esta celda ficticia acumula vehiculos, se incrementa el flujo para añadirlos al arco lo antes posible
  42.             else                              
  43.                 Q0(j) = qentrada(j); % si es insuficiente, se almacena
  44.             end
  45.            
  46.             q0(j) = min(Q0(j), (kmax_via*L - n(1,j)) / At); % se impone restriccion y se obtiene un flujo "real"
  47.             k0(j+1) = k0(j) + (qentrada(j) - q0(j))*At / L; % actualizamos el valor de la densidad
  48.             n(i,j+1) = n(i,j) + (q0(j) - q(i,j))*At; % finalmente se ejecuta la ecuación de la ev. del tráfico.
  49.         else
  50.             % ecuacion que rige la evolución del trafico, ejecutada al final del bloque (pero para el resto de celdas)
  51.             n(i,j+1) = n(i,j) + (q(i-1,j) - q(i,j))*At;
  52.         end  
  53.     end  
  54. end
  55.  
  56. nsum_nosem = sum(n(:,1:f), 1); % sacar cuantos vehiculos hay en total en cada celda, en func del tiempo  
  57. qarco_nosem = q(numCeldas, 1:f); % sacar flujo de salida del arco en func del tiempo
  58.  
  59. % reinicio de variables para ejecutar de nuevo la sim con el semaforo
  60. q0 = zeros(1,f);    
  61. Q0 = zeros(1,f);    
  62. k0 = zeros(1,f);    
  63. n = zeros(numCeldas,f);      
  64. k = zeros(numCeldas,f);      
  65. q = zeros(numCeldas,f);      
  66. Q = zeros(numCeldas,f);    
  67.  
  68. %% ------ Simulación con semáforo ------ %
  69. for j = 1:f
  70.     for i = 1:numCeldas
  71.         k(i, j) = n(i, j)/L;
  72.         Q(i,j) = funcTrafico(k(i, j), kc, kmax, qc);
  73.  
  74.         % Se implementa funcSemaforo. El semaforo se coloca en la ultima region. Impone nuevas restricciones en la ultima celda
  75.         if i < 5
  76.             q(i,j) = min(Q(i,j), (kmax_via*L - n(i+1, j)) / At);
  77.         else
  78.             q(i,j) = min(Q(i,j), funcSemaforo(t(j))); % en vez del caso por defecto se implementa la funcSemaforo que limita el flujo de la celda 5
  79.         end
  80.  
  81.         % de nuevo, el caso especial de celdas ficticias
  82.         if i == 1
  83.             qentrada(j) = 4000*exp(-(t(j) - 8) ^2);
  84.             if k0(j) > 0
  85.                 Q0(j) = qentrada(j) + k0(j)*L / At;
  86.             else
  87.                 Q0(j) = qentrada(j);
  88.             end
  89.             q0(j) = min(Q0(j), (kmax_via*L - n(1,j)) / At);
  90.             k0(j+1) = k0(j) + (qentrada(j) - q0(j))*At / L;
  91.             n(i,j+1) = n(i,j) + (q0(j) - q(i,j))*At;
  92.         else
  93.             n(i,j+1) = n(i,j) + (q(i-1,j) - q(i,j))*At;
  94.         end
  95.     end
  96. end
  97.  
  98. nsum_sem = sum(n(:,1:f), 1);  
  99. qarco_sem = q(numCeldas, 1:f);
  100.  
  101. %% ----- Plot de los datos ----- %
  102. figure(1)
  103. subplot(2,1,1)
  104. plot(t, nsum_nosem, 'b');
  105. title('Num vehiculos tot. en el arco en func. del tiempo');
  106. ylabel('Vehículos/hora');
  107. hold on
  108. subplot(2,1,2)
  109. plot(t, qarco_nosem, 'b');
  110. title('Flujo de salida del arco en func. del tiempo');
  111. ylabel('Vehículos/hora');
  112. hold off
  113. saveas(figure(1), 'figura1sinSemaforo.png');
  114.  
  115. figure(2)
  116. subplot(2,1,1)
  117. plot(t, nsum_nosem, 'b');
  118. hold on;
  119. plot(t, nsum_sem, 'r');
  120. title('Num vehiculos tot. en el arco en func. del tiempo');
  121. ylabel('Vehículos/hora');
  122. % --- juntos
  123. subplot(2,1,2)
  124. plot(t, qarco_nosem, 'b');
  125. hold on;
  126. plot(t, qarco_sem, 'r');
  127. title('Flujo de salida del arco en func. del tiempo');
  128. ylabel('Vehículos/hora');
  129. legend('Sin semáforo', 'Con semáforo');
  130. saveas(figure(2), 'figura2comparaSemaforo.png');