% Deterministic queueing simulation of an airport facility
% Considers demand and capacity as continuous (time-dependent)  functions
% 
% Integrates the demand and supply flow rates to get
% cumulative supply and demand curves

% Given inputs to the problem:
%
% demand = demand function (lambda) as a function of time (t)
% capacity = Supply = Capacity function (mu) as a function of time (t)
% time = Time vector (units of time up to you and defined by the problem context)

% Outputs of the problem:

% p(1) = L(t) = queue length as a function of time
% p(2) = A(t) = area under the queue length function (~ total delay).
% p(2) is the Integral of L(t)

% Units of L(t) are entities (aircraft, passengers, etc.)
% Units of A(t) are entities-time (aircraft-hours, passenger-hours, etc.)
%
% Basic state equations =>

% p(1) = L(t) = L(t-1) + [demand(t) - supply(t)] (dt)           % queue length 
% p(2) = A(t) = A(t-1) + [L(t)] (dt)                                      % area under L(t) curve
%
%  Demand / Supply
%		| 
%		|________			Demand 	
%		|   	 \			__________
%		|_________\_______/___________\______ Supply	
%     |			-----			   \
%		|		   
%		|		    
%		|___________________________________
%			    			Time
%
%
% Programmer: Toni Trani
% Date: November 10, 2023
% 
% Calls: fqueue_deterministic_mod.m
%
%  fqueue_deterministic_mod.m computes the integrals of the demand and supply rates
%  of change

clear
%
% Define global variables

global demand capacity time

% Enter demand function as an array of values over time

% time = time vector containing discrete points in time
% demand = vector with values of demand rate  (entities/time) observed in the system
% capacity = vector with values of capacity rate (entities/time) for the system

	% time = [0 1 1.00001 2 2.00001 3 4 ];                           % values of time (time vector)
	% demand = [1500 1500 500 500 500 500 500];             % values of demand over time
	% capacity = [1000 1000 1000 1000 1500 1500 1500 ];  % values of capacity over time

   % Alternative demand function (to illustrate the time-depent nature of the solution)
   %
   % Example in course notes for two runways used in segregated mode
   time = 0.5:1:23.5;                                                                        % values of time (time vector in vector form)
   demand = [7 10 4 3 10 21 22 33 40 38 32 20 32 23 37 40 25 23 20 37 36 29 20 13];             % arrival values of demand over time
    capacity = [32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32]; 

% Compute min and maximum values for proper scaling in plots
mintime 	= min(time);
maxtime 	= max(time+3);
maxd	 	= max(demand);
maxc 		= max(capacity);
mind	 	= min(demand);
minc 		= min(capacity);
scale		= round(.2 *(maxc+maxd)/2)
minplot = min(minc,mind) - scale;
maxplot	= max(maxc,maxd) + scale;

po = [0 0 0 0 0];			% intial number of vehicles in the queue (assumed to be zero)
to = mintime;
tf = maxtime;
tspan = [to tf];

% where:
% to is the initial time to solve this equation
% tf is the final time
% tspan is the time span to solve the simulation

% Call the ODE solver in Matlab and call function fqueue_deterministic that
% contains the rates of change of demand and capacity

[t,p] = ode23('fqueue_deterministic_mod',tspan,po);

% Compute statistics

Ltmax 		= max(p(:,1));
tdelay 		= max(p(:,2));

  clc
  disp([blanks(5),'Deterministic Queueing Model '])
  disp(' ')
  disp(' ')
  disp([blanks(5),'  Simulation Period (hours) = ', num2str(maxtime)])

  disp(' ')
  
  disp(' ')
  disp([blanks(5),'  Total delay (aircraft-hour) = ', num2str(tdelay)])
  disp([blanks(5),'  Max queue length (aircraft) = ', num2str(Ltmax)])
  disp(' ')


% Plot the demand and supply functions

figure
set(gca,'fontsize', 20) 

subplot(2,1,1)
plot(time,demand,'b',time,capacity,'--k')
xlabel('Time (hours)','fontsize',20)
ylabel('Demand or Capacity Rates (Entities/hour)','fontsize',20)
axis([mintime maxtime minplot maxplot])
legend('Demand (\lambda)','Capacity (\mu)')
grid

subplot(2,1,2)
plot(t,p(:,1),'b')
xlabel('Time (Hours)','fontsize',20)
ylabel('Queue Length (L(t))','fontsize',20)
grid

% pause

% Plot the results of the numerical integration procedure

figure
set(gca,'fontsize', 20) 

subplot(2,1,1)
plot(t,p(:,1),'b')
xlabel('Time (hours)','fontsize',20)
ylabel('Queue Length (L(t)) - Entities','fontsize',20)
grid

subplot(2,1,2)
plot(t,p(:,2),'k')
xlabel('Time (hours)','fontsize',20)
ylabel('Total Delay (Entities-hour)','fontsize',20)
grid


% New figure to plot the integration results of demand  and supply (capacity) in the same
% plot
figure

plot(t,p(:,4),'o-k',t,p(:,5),'^--r')
set(gca,'fontsize', 20) 
xlabel('Time (hours)','fontsize',20)
ylabel('Cumulative Number of Entities','fontsize',20)
legend('Cumulative Demand (\int_{0}^{t} \lambda(t) dt)','Cumulative Capacity (\int_{0}^{t} \mu(t) dt)','fontsize',20)
grid


