aoi_fcfs_md1.m

function [meanAoI, varAoI, peakAoI] = aoi_fcfs_md1(lambda, d)
%AOI_FCFS_MD1 Mean, variance, and peak AoI for M/D/1 FCFS queue
%
% [meanAoI, varAoI, peakAoI] = aoi_fcfs_md1(lambda, d)
%
% Computes the Age of Information metrics for an M/D/1 queue with
% First-Come First-Served (FCFS) discipline.
%
% M/D/1: Poisson arrivals with rate lambda, deterministic service time d.
%
% Parameters:
%   lambda (double): Arrival rate (Poisson arrivals)
%   d (double): Deterministic service time
%
% Returns:
%   meanAoI (double): Mean (average) Age of Information
%   varAoI (double): Variance of Age of Information
%   peakAoI (double): Mean Peak Age of Information
%
% Formulas (from Inoue et al., IEEE Trans. IT, 2019):
%   Uses M/GI/1 FCFS results with deterministic service.
%   For M/D/1:
%     E[H] = d, E[H^2] = d^2 (deterministic)
%     rho = lambda * d
%     E[W] = lambda * d^2 / (2*(1-rho))  (Pollaczek-Khinchine)
%     E[T] = E[W] + d
%     E[A] = d*(1/2 + 1/(2*(1-rho)) + ((1-rho)/rho)*exp(rho))  (exact)
%
% Reference:
%   Y. Inoue, H. Masuyama, T. Takine, T. Tanaka, "A General Formula for
%   the Stationary Distribution of the Age of Information and Its
%   Application to Single-Server Queues," IEEE Trans. Information Theory,
%   vol. 65, no. 12, pp. 8305-8324, 2019.
%
% See also: aoi_fcfs_mm1, aoi_fcfs_mgi1, aoi_lcfspr_md1
 
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
% Validate inputs
if lambda <= 0
    line_error(mfilename, 'Arrival rate lambda must be positive');
end
if d <= 0
    line_error(mfilename, 'Service time d must be positive');
end
 
% Compute utilization
rho = lambda * d;
 
% Check stability
if rho >= 1
    line_error(mfilename, 'System unstable: rho = lambda*d = %.4f >= 1', rho);
end
 
% Service time moments (deterministic)
E_H = d;
E_H2 = d^2;
 
% Mean waiting time (Pollaczek-Khinchine for M/G/1)
% E[W] = lambda * E[H^2] / (2 * (1 - rho))
E_W = lambda * E_H2 / (2 * (1 - rho));
 
% Mean system time (sojourn time)
E_T = E_W + E_H;
 
% Mean interarrival time
E_Y = 1 / lambda;
 
% Mean AoI for M/D/1 FCFS (exact, Inoue et al. 2019)
% E[A] = d * (1/2 + 1/(2*(1-rho)) + ((1-rho)/rho) * exp(rho))
% accounts for the negative correlation between interarrival and waiting
% times (validated against simulation)
meanAoI = d * (0.5 + 1 / (2 * (1 - rho)) + ((1 - rho) / rho) * exp(rho));
 
% Mean Peak AoI
% E[Apeak] = E[T] + E[Y]
peakAoI = E_T + E_Y;
 
% Variance of AoI
% For M/D/1, the variance depends on the second moment of AoI.
% Using the general M/G/1 result with deterministic service:
% Var[A] can be computed from LST derivatives, but for deterministic
% service, we use:
% Var[Y] = 1/lambda^2 (exponential interarrival)
% Var[H] = 0 (deterministic service)
%
% A simpler approximation for variance:
% For now, compute using the relationship with second moments
% This is an approximation; exact formula requires LST analysis
E_Y2 = 2 / lambda^2;  % Second moment of exponential
E_T2 = E_H2 + 2*E_H*E_W + E_W^2 + lambda * E_H2 / (1 - rho); % Approx
 
% Use relationship: E[A^2] approx from paper results
% For deterministic service, variance is lower than exponential case
varAoI = E_Y2 - E_Y^2 + 2*E_W*E_H / (1-rho) + E_H2 * rho / (1-rho)^2;
 
% Ensure non-negative variance (numerical safety)
if varAoI < 0
    varAoI = 0;
end
 
end