6.6 Cenni alla teoria delle code
I processi a stati discreti che abbiamo introdotto sopra hanno applicazioni in tantissimi ambiti. In questa sezione mostriamo come semplici modelli possano essere utilizzati per studiare la teoria delle code, ossia delle linee d’attesa che si possono formare in situazioni realistiche, ad esempio quando più persone vogliono accedere ad un servizio (entrare in un negozio, o pagare alla cassa), oppure dei veicoli si presentano ad un casello autostradale, o ancora delle istanze di calcolo che devono essere eseguite da una o più processori in un computer. Lo studio delle code permette di individuare strategie per migliorare l’esperienza di chi è in attesa (ridurre i tempi) rendendone più efficiente il servizio (e quindi eventualmente ridurre i costi). La teoria delle code è un campo molto esteso e noi ne presentiamo solamente i modelli più semplici come esempi interessanti di processi di Markov a salti.
Usiamo il termine clienti (in inglese si usa anche il termine jobs) per indicare genericamente le persone, le auto, i processi ecc. che nello specifico esempio di coda devono essere serviti da uno o più serventi (in inglese servers).
Gli aspetti fondamentali che si vogliono modellizzare di una coda sono l’ingresso di uno o più clienti, l’attesa (eventualmente nulla) che un servente prenda in carico il compito richiesto, e infine l’uscita dalla coda quando il compito è svolto. Una volta introdotto un modello di coda, è di interesse calcolare quantità come il tempo medio di attesa, il numero medio di clienti in coda, ma anche ovviamente stimare i parametri di un modello sulla base di quantità osservate in una coda reale.
Per classificare i vari modelli di code studiati in letteratura, Kendall propose una notazione abbreviata16: si usano due lettere e un numero (\(A/S/c\)), in cui la prima lettera (\(A\)) indica un “processo” di arrivo dei clienti, la seconda (\(S\)) la legge del tempo di servizio per ciascun cliente, e il numero \(c\) il numero dei serventi.
In questa sezione consideriamo solamente i modelli \(M/M/c\), in cui gli arrivi e i servizi sono Markoviani a tempi continui, più precisamente con tempi esponenziali di due parametri (\(\lambda\) per il tasso di arrivo dei clienti e \(\pi\) per il servizio), e sono quindi formalmente definiti come processi di Markov a salti nell’insieme degli stati \(E = \mathbb{N}\). Lo stato \(n \in E\) rappresenta infatti la situazione in cui vi siano \(n\) clienti in servizio oppure in attesa di essere serviti. Una volta che un cliente è servito, esso “scompare” dalla coda, che quindi passa dallo stato \(n\) allo stato \(n-1\). L’arrivo di un cliente è invece rappresentato con una transizione dallo stato \(n\) allo stato \(n+1\) (non supporremo mai che due o più clienti arrivino oppure lascino la coda allo stesso istante). A seconda del numero di serventi \(c\) definiamo una matrice delle intensità di salto diversa.
6.6.1 Processo di Poisson
Il modello più semplice rappresenta la situazione in cui non vi siano serventi (o meglio si è interessati solo al processo di arrivo dei clienti): potrebbe essere classificato come \(M/M/0\), anche se più comunemente è detto processo di Poisson di intensità \(\lambda>0\) (il tasso di ingresso dei clienti in coda). Le uniche transizioni avvengono da uno stato \(n\) a uno stato \(n+1\), e si pone, per ogni \(n \in \mathbb{N}\), \[ L_{n \to n+1} = \lambda, \quad L_{n \to n} = - \lambda\] e \(L_{n \to k} = 0\) se \(k \neq n\), \(k \neq n+1\).
Figura 6.7: grafo associato al processo di Poisson
Ogni stato è quindi transitorio, e non esiste una distribuzione invariante. Infatti, se \(\pi\) fosse invariante, allora \[ 0 = (\pi L)_0 = \pi_0 L_{0 \to 0} = - \pi_0 \lambda \quad \text{e quindi} \quad \pi_0 = 0,\] mentre \[ 0 = (\pi L)_1 = \pi_0 L_{0 \to 1} + \pi_1 L_{1\to 1} = - \pi_1 \lambda \quad \text{e quindi} \quad \pi_1 = 0,\] e similmente si ottiene, per ogni \(n \ge 1\), \(\pi_n = 0\). Non è possibile quindi che \(\pi\) sia una densita discreta di probabilità (non può essere \(\sum_{n} \pi_n = 1\)).
C’è un legame preciso tra il processo di Poisson e la densità discreta di Poisson definita nell’Esempio 3.5. Infatti, se definiamo, per ogni \(t \ge 0\), \[ \pi^t_n \propto \frac{ (t \lambda)^n}{n!} =\frac{ (t \lambda)^n}{n!} \exp\bra{-t\lambda}.\] la densità Poisson di parametro \(t \lambda\), allora \(\pi^t\) è la densità della marginale \(X_t\) di un processo di Poisson tale che \(X_0 = 0\) (la densità \(\pi^0\) vale infatti \(1\) nel valore \(0\)). Basta verificare che valga la master equation, per ogni \(n \in \mathbb{N}\), \(t \ge 0\), \[ \frac{d}{dt} \pi^t_n = (\pi^t L)_n = \begin{cases} - \lambda \pi^t_0 &\text{se $n=0$,} \\ \lambda( \pi^t_{n-1} - \pi^t_n ) & \text{se $n \ge 1$.}\end{cases} \] Calcoliamo quindi \[ \frac {d}{dt} \exp\bra{-t\lambda}\frac{ (t \lambda)^n}{n!} = \begin{cases} - \lambda \exp\bra{- t \lambda} & \text{se $n=0$,}\\ \frac{ n t^{n-1} \lambda^n}{n!} \exp\bra{- t \lambda} - \lambda \frac{(t\lambda)^n}{n!} \exp\bra{- t \lambda} & \text{ se $n \ge 1$.} \end{cases} \] Per concludere nel caso \(n \ge 1\) basta notare che \[\frac{ n t^{n-1} \lambda^n}{n!} \exp\bra{- t \lambda} = \lambda \frac{ (t\lambda)^{n-1}}{(n-1)!} \exp\bra{- t \lambda} = \lambda \pi^t_{n-1}.\] Consideriamo infine il problema di stimare il parametro \(\lambda\) a partire dall’osservazione di un cammino \(\gamma = (x_0 \to x_1 \to \ldots \to x_{n})\) con tempi di permanenza \(t_1\) (nello stato \(x_0\)), \(t_2\) (in \(x_1\)), …., \(t_{n+1}\). Osserviamo che, poiché i salti avvengono solo tra uno stato \(n\) e il successivo \(n+1\), deve essere \(x_1 = x_0+1\), \(x_2=x_0+2\), ecc., quindi la formula per la verosimiglianza in questo caso diventa \[ L( \Lambda = \lambda; X = \gamma) = \prod_{k=1}^n \exp\bra{- \lambda t_k} \lambda = \lambda^n \exp\bra{- \lambda T}. \] dove abbiamo supposto per semplicità che fosse noto a priori che \(X_0 = x_0\) e abbiamo indicato con \(T = \sum_{k=1}^n t_k\). La stima di massima verosimiglianza \(\lambda_{MLE}\) si trova passando al logaritmo e imponendo che la derivata si annulli. Si trova \[ \frac{n}{\lambda_{MLE}} - T = 0 \quad \text{quindi} \quad \lambda_{MLE} = \frac{n}{T}.\]
Esempio 6.20 In un intervallo di tempo \(T = 5\) minuti si osservano entrare \(n = 10\) persone in un supermercato. Si può introdurre quindi un processo di Poisson con intensità \(\lambda = 2\) (persone/minuto) per modellizzare gli ingressi.
Per un approccio bayesiano, in cui i calcoli siano particolarmente semplici si può introdurre una densità a priori per la variabile \(\Lambda\) del tipo Gamma, ossia \[ p(\Lambda = \lambda) \propto \lambda^{\alpha-1 } \exp\bra{ - \beta \lambda },\] dove \(\alpha\), \(\beta>0\) sono parametri (si scrive anche \(\Gamma(\alpha, \beta)\)). Il valor medio di \(\Lambda\) (a priori) si può calcolare ed è dato da \(\alpha/\beta\), mentre la moda è \((\alpha-1)/\beta\) (per \(\alpha\ge 1\)). Nel caso \(\alpha=1\) essa coincide con una densità esponenziale di parametro \(\beta\). Questa densità rappresenta in modo preciso una possibile informazione nota sul parametro \(\lambda\), ad esempio informalmente che \(\lambda \approx \alpha/\beta\).
La densità a posteriori diventa quindi \[ p(\Lambda = \lambda| X = \gamma) \propto \lambda^{n+\alpha-1 } \exp\bra{ -(\beta+T) \lambda },\] ossia una densità \(\Gamma(n+\alpha, \beta+T)\). La moda della densità a posteriori è quindi (se \(n +\alpha \ge 1\)) \[ \lambda_{MAP} = \frac{ n+\alpha}{T +\beta}.\]
6.6.2 Code \(M/M/1\)
Consideriamo ora la situazione in cui vi sia un solo servente (\(M/M/1\)), e che il tempo di servizio per ciascun cliente sia una variabile esponenziale di parametro \(\mu\) (ogni cliente sia indipendente dagli altri). Per modellizzare la coda con un processo di Markov a salti, conviene considerare prima il caso in cui non vi siano arrivi. In tal caso si osserveranno solamente salti da uno stato \(n\) verso \(n-1\) (se \(n \ge 1\)) con dei tempi di permanenza esponenziali di parametro \(\mu\). Pertanto, si avrà (se \(n\ge 1\)) \[L_{n \to n-1} = \mu.\] Nel caso in cui vi siano arrivi con un processo di Poisson di intensità \(\lambda\), poniamo, per \(n \ge 0\), \[ L_{n \to n+1} = \lambda,\] e di conseguenza \[ L_{ n \to n} = \begin{cases} -\lambda & \text{se $n = 0$}\\ -(\lambda + \mu) & \text{se $n \ge 1$,}\end{cases}\] avendo posto \(L_{ n \to k} = 0\) se \(k \notin \cur{n-1, n, n+1}\).
Figura 6.8: grafo associato ad una coda \(M/M/1\)
Ogni stato è ricorrente, ma essendo infiniti stati non è ovvio che esista una distribuzione invariante. Vi è infatti una competizione tra il tasso di arrivo \(\lambda\) e di uscita \(\mu\). Analogamente a quanto fatto nel caso di Poisson, si può risolvere l’equazione \(\pi L = 0\) e ottenere per \(n \ge 0\), \[\pi_n \propto \bra{ \frac{\lambda}{\mu}}^n.\] È un semplice esercizio verificare che \(\pi\) soddisfa l’equazione, ossia \[ 0 = (\pi L)_n = \pi_{n-1} L_{n-1 \to n} + \pi_n L_{n\to n} +\pi_{n+1}L_{n+1\to n} = \lambda \pi_{n-1} -(\lambda+\mu) \pi_n + \mu \pi_{n+1}\] (per \(n \ge 1\), mentre per \(n=0\) bisogna porre \(L_{-1 \to 0} = 0\)).
Dovendo garantire che tale \(\pi\) sia una densità di probabilità, bisogna che \[\sum_{n} \bra{ \frac\lambda \mu}^n < \infty,\] ma tale serie (geometrica) converge se e solo se \(\lambda<\mu\). In altri termini, esiste un equilibrio per la coda se e solo se il tasso di arrivo è strettamente minore di quelo di uscita, altrimenti il numero di persone in coda cresce (più lentamente del processo di Poisson, ma comunque in modo inarrestabile).
La distribuzione invariante se \(\lambda<\mu\) è quindi una densità geometrica di parametro \(1-\lambda/\mu\). In particolare, il valor medio del numero di clienti nella coda (in attesa o in servizio) \(N\) in regime stazionario (ossia se la marginale del processo a salti è \(\pi\)) vale \[ \E{N} = \sum_{n} n \pi_n = \frac{\lambda}{\mu -\lambda},\] una quantità che diverge al tendere di \(\lambda\) verso \(\mu\).
deltal <- 0.001
l <- seq(0.5, 0.99, by = deltal)
mu <- 1
plot(l, l/(1 - l), type = "l", col = miei_colori[2],
lwd = 3, xlab = "tasso di ingresso",
ylab = "numero medio di clienti")
abline(h = 50, col = miei_colori[1], lwd = 3)
Figura 6.9: grafico di \(\E{N}\) per \(\mu=1\) in funzione di \(\lambda\) (in rosso) e una soglia massima di possibili persone in coda (in azzurro)
Tale divergenza è problematica se \(\lambda\) si trova molto vicino a \(\mu\) e per qualche motivo il tasso di ingresso aumenta, anche di poco, portando a superare il limite massimo possibile di clienti in coda (che nella realtà esiste sempre).
Esempio 6.21 Prima della pandemia gli ospedali operavano in modo da usare tutti o quasi i posti letto disponibili – usando quindi efficientemente tutte le risorse di personale e strutture. Rappresentando un ospedale come una coda, erano quindi molto vicini al limite massimo di possibili “clienti” (i pazienti) in coda. L’arrivo del nuovo coronavirus ha avuto l’effetto di aumentare il tasso di ingresso, causando un aumento notevole di \(N\) con conseguenze potenzialmente catastrofiche.
Veniamo ora alla stima dei parametri \((\lambda, \mu)\) sulla base dell’osservazione di un cammino \(\gamma = (n_0 \to n_1 \to \ldots \to n_\ell)\) con i soliti tempi di permanenza \(t_1\), \(t_2\), …, \(t_\ell\) e poniamo pure \(T = \sum_{k=1}^\ell t_i\). Supponiamo inoltre che il cammino osservato non passi mai per lo stato \(0\) (quindi c’è sempre almeno un cliente in coda). Usando l’espressione (6.7), si ottiene la verosimiglianza \[ L( \lambda, \mu; X=\gamma) = \exp\bra{-(\lambda+\mu)T} \lambda^{\gamma_+} \mu^{\gamma_-},\] dove \(\gamma_+\) indica il numero di arrivi osservati in \(\gamma\) (ossia transizioni da uno stato \(n\) a \(n+1\)), mentre \(\gamma_-\) il numero di uscite. La stima di massima verosimiglianza è quindi \[ \lambda_{MLE} = \frac{\gamma_+}{T}, \quad \mu_{MLE} = \frac{\gamma_-}{T}.\] Se invece il cammino trascorre un tempo \(T_0\) nello stato \(0\), l’espressione per la verosimiglianza cambia (al posto di \(-(\lambda+\mu)T\) si trova \(-\lambda T -\mu (T-T_0)\) e di conseguenza \(\lambda_{MLE}\) non cambia, ma \[ \mu_{MLE} = \frac{\gamma_-}{T-T_0}.\] L’interpetazione è che il tempo trascorso con la coda vuota non può essere utile alla stima del tasso di uscita dei clienti dalla coda, e quindi va sottratto.
Esempio 6.22 In un intervallo di \(10\) minuti si osservano \(5\) persone arrivare alla cassa di un supermercato e \(3\) persone uscirne. Supponendo che la cassa non sia mai senza lavoro si stimano i parametri \(\lambda = 1/2\) persone al minuto, \(\mu = 3/10\) persone al minuto. Se invece la cassa è rimasta priva di persone in coda per \(4\) minuti, si stima \(\mu = 3/6 =1/2\) persone al minuto.
Tralasciamo l’approccio bayesiano, che è simile al caso del Poisson (supponendo ad esempio \(\lambda\), \(\mu\) indipendenti a priori).
6.6.3 Code \(M/M/\infty\)
Consideriamo infine la situazione opposta, in cui vi sono un numero arbitrariamente grande, idealmente infinito, di serventi (\(M/M/\infty\)). Supponiamo ancora che il tempo di servizio per ciascun cliente sia una variabile esponenziale di parametro \(\mu\) (e ogni cliente sia indipendente dagli altri). Per capire quali intensità di salto definire, conviene considerare ancora il caso in cui non vi siano arrivi. In tal caso si osserveranno solamente salti da uno stato \(n\) verso \(n-1\) (se \(n \ge 1\)) con dei tempi di permanenza dati dal minimo di \(n\) variabili aleatorie \(T_1\), \(T_2\), , \(T_n\) esponenziali indipendenti tra loro (infatti, la transizione avviene appena il cliente che impega meno tempo tra gli \(n\) in servizo lascia la coda). Possiamo allora affermare che tale tempo è una variabile esponenziale \(T\), di parametro \(n\mu\): infatti si calcola la funzione di sopravvivenza (per \(t\ge 0\)) \[\begin{equation*}\begin{split} \SUR_T(t) & = P( \min \cur{T_1, T_2, \ldots, T_k} > t )\\ & = P( T_1>t, T_2>t, \ldots, T_n >t) \\ & = P(T_1>t) P(T_2>t) \cdot \ldots\cdot P(T_n>t) \\ & = e^{-\mu t} \cdot e^{-\mu t} \cdot \ldots \cdot e^{-\mu t}\\ & = e^{-n\mu t}, \end{split} \end{equation*}\] e derivando si ottiene la densità esponenziale.
Pertanto, si avrà (se \(n\ge 1\)) \[L_{n \to n-1} = n \mu\] Nel caso in cui vi siano arrivi con un processo di Poisson di intensità \(\lambda\), poniamo, per \(n \ge 0\), \[ L_{n \to\ n+1} = \lambda,\] e di conseguenza \[ L_{ n \to n} = \begin{cases} -\lambda & \text{se $n = 0$}\\ -(\lambda + n\mu) & \text{se $n \ge 1$,}\end{cases}\] avendo posto \(L_{ n \to k} = 0\) se \(k \notin \cur{n-1, n, n+1}\).
Figura 6.10: grafo associato ad una coda \(M/M/\infty\)
Come nel caso \(M/M/1\), ogni stato è ricorrente, ma essendo infiniti stati non è ovvio che esista una distribuzione invariante. Vi è ancora una competizione tra il tasso di arrivo \(\lambda\) e di uscita \(\mu\), ma decisamente “smorzata” dal fatto che per \(n\) abbastanza grande si avrà comunque \(\lambda <n\mu\). Questo suggerisce che una distribuzione invariante esista sempre. Infatti, si può risolvere l’equazione \(\pi L = 0\) e ottenere per \(n \ge 0\), \[\pi_n \propto \frac{1}{n!} \bra{ \frac{\lambda}{\mu}}^n,\] ossia una densità Poisson di parametro \(\lambda/\mu\). Lasciamo per esercizio di verificare che \(\pi\) soddisfa l’equazione \(\pi L = 0\). Il valor medio del numero di clienti nella coda è quindi \(\E{N} = \lambda/\mu\) che cresce linearmente al crescere di \(\lambda\) (non presenta asintoti).
Infine, la stima dei parametri \((\lambda, \mu)\) sulla base dell’osservazione di un cammino \(\gamma = (n_0 \to n_1 \to \ldots \to n_{\ell-1})\) con i tempi di permanenza \(t_1\), \(t_2\), …, \(t_\ell\) si può effettuare tramite il metodo di massima verosimiglianza, usando l’espressione (6.7): \[ L( \lambda, \mu; X=\gamma) \propto \exp\bra{-\lambda T - \mu T_\gamma )} \lambda^{\gamma_+} \mu^{\gamma_-},\] dove \(T = \sum_{k=1}^\ell t_i\), \(\gamma_+\) e \(\gamma_-\) sono come nel caso \(M/M/1\) e infine \[ T_\gamma = \sum_{k=1}^\ell t_i n_i,\] (il tempo totale trascorso da tutti i clienti osservati nella coda). La stima di massima verosimiglianza è quindi \[ \lambda_{MLE} = \frac {\gamma_-}{T}, \quad \mu_{MLE} = \frac{\gamma_-}{T_\gamma}.\]
Esempio 6.23 Un ipermercato dispone di un numero molto grande di casse, e al mattino è poco frequentato cosicché ogni cliente trova sempre una cassa libera. Si osserva che per \(4\) minuti consecutivi tutte le casse erano libere, per \(3\) minuti \(1\) sola cassa era occupata, poi si è liberata per \(1\) minuto tutte le casse erano di nuovo libere, e infine per \(2\) minuti in un intervallo di \(10\) minuti minuti tutte le casse erano libere, per \(3\) minuti una sola cassa era occupata, e per i rimanenti \(3\) minuti \(5\) casse erano occupate.
Remark. Il caso generale \(M/M/c\) con \(2 \le c < \infty\) è intermedio tra i due estremi che abbiamo considerato. In particolare una distribuzione invariante esiste se e solo se \(\lambda < c\mu\).