function d = run_lab_2(ctrl) dbstop if error d.p.ctrl = ctrl; d = build_setup(d); for t = 1:d.p.t_final d = compute_speed(d,t); d = evolve_states(d,t); end d = compute_headways(d); % calculate mean bus speed for the whole simulation d.p.mean_v = mean(d.s.v(:)); end function d = build_setup(d) rng(1) % choose controller %d.p.ctrl = 1; % no control %d.p.ctrl = 2; % I-control %d.p.ctrl = 3; % PI-control %d.p.ctrl = 4; % holding control % controller parameters d.p.K_I = 1; d.p.K_P = 10; d.p.H_th = -50; % simulation length (in number of time steps) d.p.t_final = 12*60*60; % number of buses d.p.K_b = 4; % number of stops d.p.K_s = 8; % timestep (s) d.p.T = 1; % dwell time bounds (in number of time steps) d.p.dt_min = 10; d.p.dt_max = 50; d.s.dt_mean = 30; d.s.dt_stdev = 20; pd_dt = makedist('Normal',... 'mu',d.s.dt_mean,... 'sigma',d.s.dt_stdev); d.p.tpd_dt = truncate(pd_dt,d.p.dt_min,d.p.dt_max); % distance between two consecutive stops (m) d.p.D_stop = 2000; % vector of stop positions (m) d.p.L_j = (0:d.p.D_stop:(d.p.K_s-1)*d.p.D_stop)'; % bus loop length (m) d.p.L_max = d.p.L_j(end) + d.p.D_stop; % circular distances of stops to stop 1 [m] d.p.L_j_to_1 = [d.p.L_max;d.p.L_j(2:end)]; % initial distance between two consecutive buses (m) d.p.D_bus = d.p.L_max/d.p.K_b; % vector of initial bus positions (m) d.p.x_0 = (0:d.p.D_bus:(d.p.K_b-1)*d.p.D_bus)'; % preallocate memory d.s.e = NaN(d.p.K_b,d.p.t_final); d.s.s_r = NaN(d.p.K_b,d.p.t_final); d.s.s_f = NaN(d.p.K_b,d.p.t_final); d.s.x = NaN(d.p.K_b,d.p.t_final+1); d.s.dt = zeros(d.p.K_b,d.p.t_final); d.s.dt_record = []; d.s.f = NaN(d.p.K_b,d.p.t_final+1); d.s.v = NaN(d.p.K_b,d.p.t_final); d.s.h = zeros(d.p.K_b,d.p.t_final); % set initial state d.s.x(:,1) = d.p.x_0+100; d.s.f(:,1) = ones(d.p.K_b,1); % lower bounds on bus speed due to traffic (m/s) d.p.v_min_l = 4; d.p.v_min_u = 8; d.s.v_min_traffic_mean = 6; d.s.v_min_traffic_stdev = 1; pd_v_min = makedist('Normal',... 'mu',d.s.v_min_traffic_mean,... 'sigma',d.s.v_min_traffic_stdev); tpd_v_min = truncate(pd_v_min,d.p.v_min_l,d.p.v_min_u); d.s.v_min_traffic = random(tpd_v_min,d.p.K_b,d.p.t_final); % upper bounds on bus speed due to traffic (m/s) d.p.v_max_l = 10; d.p.v_max_u = 18; d.s.v_max_traffic_mean = 14; d.s.v_max_traffic_stdev = 2; pd_v_max = makedist('Normal',... 'mu',d.s.v_max_traffic_mean,... 'sigma',d.s.v_max_traffic_stdev); tpd_v_max = truncate(pd_v_max,d.p.v_max_l,d.p.v_max_u); d.s.v_max_traffic = random(tpd_v_max,d.p.K_b,d.p.t_final); end function d = compute_speed(d,t) d = compute_spacing_error(d,t); switch d.p.ctrl case 1 % no control d.s.v(:,t) = d.s.v_max_traffic(:,t); case 2 % I control d = compute_I_control(d,t); case 3 % PI control d = compute_PI_control(d,t); case 4 % holding control d.s.v(:,t) = d.s.v_max_traffic(:,t); d = compute_H_control(d,t); end end function d = compute_spacing_error(d,t) s_f = NaN(d.p.K_b,1); s_r = NaN(d.p.K_b,1); for i = 1:d.p.K_b if i == d.p.K_b s_f(i) = d.s.x(1,t) - d.s.x(i,t); if s_f(i) < 0 s_f(i) = d.s.x(1,t) + d.p.L_max - d.s.x(i,t); end else s_f(i) = d.s.x(i+1,t) - d.s.x(i,t); if s_f(i) < 0 s_f(i) = d.s.x(i+1,t) + d.p.L_max - d.s.x(i,t); end end end for i = 1:d.p.K_b if i == 1 s_r(i) = s_f(d.p.K_b); else s_r(i) = s_f(i-1); end end d.s.s_f(:,t) = s_f; d.s.s_r(:,t) = s_r; d.s.e(:,t) = s_f - s_r; end function d = evolve_states(d,t) for i = 1:d.p.K_b if d.s.f(i,t) == 1 % bus i is cruising if t == 1 % evolve position of bus i d.s.x(i,t+1) = d.s.x(i,t) + d.p.T*d.s.v(i,t); % keep flag at cruising d.s.f(i,t+1) = 1; else if d.s.f(i,t-1) == 0 % bus i was stopping during previous time step % evolve position of bus i d.s.x(i,t+1) = d.s.x(i,t) + d.p.T*d.s.v(i,t); % keep flag at cruising d.s.f(i,t+1) = 1; else % bus i was cruising during previous time step if d.s.x(i,t) > d.p.L_j(end) % bus i is cruising to stop 1 if d.s.x(i,t) + d.p.T*d.s.v(i,t) > d.p.L_max % bus i will reach stop 1 at time t % correct bus speed d.s.v(i,t) = (d.p.L_max-d.s.x(i,t))/d.p.T; % snap position to stop 1 d.s.x(i,t+1) = 0; % change flag to stopping d.s.f(i,t+1) = 0; else % bus i will not reach stop 1 at time t % evolve position of bus i d.s.x(i,t+1) = d.s.x(i,t) + d.p.T*d.s.v(i,t); % keep flag at cruising d.s.f(i,t+1) = 1; end else % bus i is cruising to a stop different than stop 1 % find the index of next stop j_next = find((d.s.x(i,t) > d.p.L_j) == 0, 1, 'first'); if d.s.x(i,t) + d.p.T*d.s.v(i,t) > d.p.L_j(j_next) % bus i will reach stop j_next at time t % correct bus speed d.s.v(i,t) = (d.p.L_j(j_next)-d.s.x(i,t))/d.p.T; % snap position to stop j_next d.s.x(i,t+1) = d.p.L_j(j_next); % change flag to stopping d.s.f(i,t+1) = 0; else % bus i will not reach stop j_next at time t % evolve position of bus i d.s.x(i,t+1) = d.s.x(i,t) + d.p.T*d.s.v(i,t); % keep flag at cruising d.s.f(i,t+1) = 1; end end end end else % bus i is stopping % keep bus position d.s.x(i,t+1) = d.s.x(i,t); d.s.v(i,t) = 0; if d.s.f(i,t-1) == 1 % bus i was cruising during previous time step % create dwell time for bus i d.s.dt(i,t) = ceil(random(d.p.tpd_dt,1,1)); d.s.dt_record(end+1) = d.s.dt(i,t); % keep flag at stopping d.s.f(i,t+1) = 0; else % bus i was stopping during previous time step if d.s.dt(i,t) == 0 && d.s.h(i,t) == 0 % dwell time counter reached 0 % AND % holding command is zero (no holding command given) % change flag to cruising d.s.f(i,t+1) = 1; else % decrease the dwell time counter by 1 d.s.dt(i,t) = d.s.dt(i,t-1) - 1; % keep flag at stopping d.s.f(i,t+1) = 0; end end end end end function d = compute_headways(d) for h = 1:d.p.K_s d.r.hw_times{h} = []; for t = 2:d.p.t_final for i = 1:d.p.K_b if d.s.x(i,t-1) < d.p.L_j_to_1(h) && d.s.x(i,t) >= d.p.L_j_to_1(h) d.r.hw_times{h}(end+1) = t; end end end end for h = 1:d.p.K_s d.r.hws{h} = []; for idx = 2:length(d.r.hw_times{h}) d.r.hws{h}(end+1) = (d.p.T/60)*... (d.r.hw_times{h}(idx) - d.r.hw_times{h}(idx-1)); end end d.r.hws_nw = []; for h = 1:d.p.K_s d.r.hws_nw = [d.r.hws_nw d.r.hws{h}]; end d.p.mH = mean(d.r.hws_nw); d.p.sH = std(d.r.hws_nw); end