Capítulo 5


Introdução à Teoria das Filas



“The quiet statisticians have changed our world – not by discovering new facts or technical developments but by changing the ways we reason, experiment and form our opinions about it.”



– Ian Hacking, 1984



5.1. Introdução


As filas de espera são fenômenos inerentes a sistemas nos quais, em determinados momentos, a demanda por serviços excede a capacidade de atendimento disponível, resultando na formação de filas físicas ou conceituais. Tais sistemas são compostos essencialmente por um processo de chegada de clientes, um ou mais servidores responsáveis pela prestação do serviço e uma disciplina que define a ordem de atendimento. A aleatoriedade associada tanto aos tempos de chegada quanto aos tempos de serviço confere a esses sistemas uma natureza estocástica, cuja dinâmica é diretamente influenciada pela relação entre a taxa de chegada e a taxa de atendimento. Quando a primeira supera a segunda, ocorre o congestionamento, obrigando os clientes a aguardarem até que haja capacidade disponível para o atendimento. Esse fenômeno é observado em diversos contextos, como serviços de saúde, instituições financeiras, centrais de atendimento, sistemas logísticos, transporte urbano e, de forma geral, em qualquer ambiente no qual recursos de serviço são finitos frente a uma demanda variável e aleatória.


5.2. Elementos de Uma Fila


A modelagem de sistemas de filas fundamenta-se na caracterização de quatro elementos básicos: processo de chegada, processo de serviço, disciplina de atendimento e configuração dos servidores (ou canais de serviço).


  • (I) Processo de Chegada: O processo de chegada descreve o mecanismo estocástico pelo qual os clientes acessam o sistema. É formalmente caracterizado pela taxa média de chegada \((\lambda)\), expressa em número médio de chegadas por unidade de tempo, e pela distribuição de probabilidade dos tempos entre chegadas. Na maioria dos modelos clássicos, assume-se que as chegadas seguem um processo de Poisson, o que implica tempos entre chegadas distribuídos segundo uma exponencial com média \(1/\lambda\). Além disso, verifica-se se o processo é estacionário (isto é, se \(\lambda\) se mantém constante ao longo do tempo) e se as chegadas são independentes entre si.


  • (II) Processo de Serviço: O processo de serviço representa a dinâmica estocástica associada ao atendimento dos clientes. É definido pela taxa média de serviço \((\mu)\), correspondente ao número médio de atendimentos realizados por unidade de tempo, por servidor. Pressupõe-se que os tempos de serviço são aleatórios, geralmente modelados por uma distribuição exponencial com média \(1/\mu\), embora outros modelos sejam possíveis dependendo do contexto. Este processo deve ser analisado quanto à sua estacionariedade e independência em relação às chegadas e aos demais atendimentos. A capacidade de atendimento está diretamente relacionada ao número de servidores \(s\), operando em paralelo, e à taxa \(\mu\). A eficiência do serviço é avaliada por meio de métricas como o tempo médio de espera na fila, o tempo médio no sistema e a taxa de utilização.


  • (III) Disciplina de Atendimento: A disciplina de atendimento estabelece a regra de ordenamento dos clientes na fila, definindo como são selecionados para o serviço. As mais recorrentes na literatura são:

    • FIFO (First-In, First-Out): ordem cronológica de chegada.
    • LIFO (Last-In, First-Out): o último a chegar é o primeiro a ser atendido.
    • SIRO (Service In Random Order): a ordem de atendimento é aleatória.
    • Atendimento com Prioridade: clientes são atendidos conforme uma hierarquia predefinida.
    • Jockeying: permite que clientes transitem entre filas paralelas, buscando reduzir seus tempos de espera.
    • Fila Única: todos os clientes aguardam em uma única fila, sendo alocados ao primeiro servidor disponível.


  • (IV) Configuração dos Servidores (Canais de Serviço): Este componente refere-se ao número de servidores, denotado por \(s\), e à sua disposição estrutural no sistema, que pode assumir diferentes configurações conforme a natureza do atendimento. Os servidores podem estar organizados em paralelo, possibilitando o atendimento simultâneo de múltiplos clientes; em série, quando o atendimento ocorre em etapas sequenciais obrigatórias; ou em estruturas mais complexas, como redes de filas, que interligam múltiplas estações de serviço. A capacidade agregada do sistema relaciona-se ao número de servidores e às respectivas taxas individuais de atendimento, representada pelo produto \(s \mu\) no caso de servidores paralelos. A taxa de utilização do sistema é expressa por \(\rho = \lambda/(s \mu)\), indicando a fração do tempo em que, em média, os servidores permanecem ocupados. A estabilidade do sistema exige que \(\rho < 1\), condição necessária para evitar o crescimento indefinido da fila de espera. A definição da configuração estrutural e do número de servidores demanda um equilíbrio entre os custos operacionais e os níveis de desempenho desejados, os quais são avaliados por meio de métricas como o tempo médio de espera, o tamanho da fila e os níveis de serviço, considerando as particularidades do sistema analisado.


5.3. Notação de Um Sistema de Fila


A notação mais adotada na representação de sistemas de filas foi proposta por Kendall (1953), estabelecendo-se como padrão consolidado na literatura de teoria das filas. Este formalismo adota a forma \(A/B/m/k/M\), na qual cada elemento descreve uma característica estrutural ou operacional fundamental do sistema. Especificamente, o termo \(A\) caracteriza o processo estocástico que governa os tempos entre chegadas dos clientes, definido por sua distribuição de probabilidade e pelas propriedades dos seus incrementos (como independência, estacionaridade e memória). Analogamente, \(B\) especifica o processo estocástico associado aos tempos de serviço, incluindo sua distribuição e sua dinâmica probabilística. Os valores atribuídos aos termos \(A\) e \(B\) indicam as classes dos processos estocásticos envolvidos, sendo os mais comuns:

  • M (Markovian): interchegadas ou tempos de serviço seguem uma distribuição exponencial, associada a um processo de Poisson no caso das chegadas, ou a um processo Markoviano no serviço (memória nula e incrementos independentes e estacionários).
  • D (Deterministic): tempos constantes e fixos, tanto nas chegadas quanto nos serviços.
  • G (General): qualquer distribuição de probabilidade genérica, com média e variância finitas, sem suposição de forma específica.
  • \(E_k\) (Erlang de ordem k): distribuição Erlang, que resulta da soma de \(k\) variáveis exponenciais independentes, representando processos com menor variabilidade que o exponencial.
  • \(GI\) (General Independent): processo com tempos entre chegadas ou de serviço provenientes de distribuições quaisquer, independentes entre si, mas não necessariamente idênticas.


Já o parâmetro \(m\) indica o número de servidores ou canais de atendimento em operação simultânea, enquanto \(k\) representa a capacidade máxima do sistema, contabilizando tanto os clientes em atendimento quanto aqueles em espera na fila; na ausência de restrição explícita, considera-se \(k \to \infty\). Por fim, o termo \(M\) define a disciplina de atendimento, responsável por determinar a ordem em que os clientes são selecionados da fila para o serviço, sendo a disciplina FIFO (First-In, First-Out) a mais frequentemente adotada. Alguns exemplos práticos de modelos que ilustram esta notação incluem:


  • \(M/M/1\): descreve um sistema com chegadas segundo um processo de Poisson, tempos de serviço exponenciais, um único servidor, capacidade ilimitada (\(k \to \infty\)) e disciplina FIFO. Este modelo caracteriza situações de atendimento isolado, como caixas bancários, postos de serviços e guichês.

  • \(M/M/s\): generaliza o anterior para \(s\) servidores em operação paralela, mantendo as mesmas premissas para chegadas e tempos de serviço. Este modelo aplica-se a ambientes como centrais de atendimento telefônico, check-in em aeroportos, unidades hospitalares e agências bancárias.

  • \(M/M/1/k\): incorpora uma capacidade finita \(k\), de modo que qualquer cliente que chega quando o sistema está cheio é bloqueado e não recebe atendimento. Este modelo representa ambientes com recursos limitados, como estacionamentos, buffers de comunicação ou sistemas de transmissão de dados com capacidade restrita.

  • \(M/D/1\): adota chegadas segundo um processo de Poisson e tempos de serviço determinísticos, com duração fixa para cada atendimento, um único servidor e disciplina FIFO. Este modelo representa sistemas industriais automatizados e processos produtivos com tempo de operação constante por unidade processada.

  • \(M/G/1\): assume chegadas Poisson e tempos de serviço com distribuição arbitrária (geral), mantendo um único servidor. Este modelo é apropriado quando os tempos de serviço apresentam variabilidade não compatível com o modelo exponencial, como operações de manutenção, atendimento médico especializado ou processos administrativos complexos.

  • \(GI/M/1\): admite uma distribuição geral e independente (GI) para os tempos entre chegadas, mantendo tempos de serviço exponenciais, com um único servidor. Este modelo é adequado para situações em que o padrão de chegada não segue um processo Poisson, como fluxos com dependências temporais, agendamentos ou demanda sazonal.

  • \(E_k/M/1\): adota uma distribuição Erlang de ordem \(k\) para os tempos entre chegadas, tempos de serviço exponenciais e um único servidor. Este modelo descreve cenários nos quais os tempos entre chegadas apresentam variabilidade reduzida, típica de sistemas com controle rígido sobre o fluxo de entrada, como linhas de produção balanceadas, transporte regular ou serviços com agendamento fixo.

  • \(G/G/1\): constitui a formulação mais geral, permitindo qualquer distribuição para os tempos entre chegadas e para os tempos de serviço, ambos independentes. Este modelo representa sistemas operacionais altamente complexos, sujeitos a grande variabilidade tanto na demanda quanto na prestação do serviço, como operações logísticas, redes hospitalares ou processos produtivos sob demanda flutuante.


5.4. Conceitos Fundamentais em Filas


Considere um sistema de fila com servidor único, no qual as chegadas de clientes ocorrem segundo um processo estocástico com taxa média \(\lambda\) (clientes por unidade de tempo) e cada atendimento possui tempo médio de serviço \(\bar{x}\). Para que o sistema seja estável — ou seja, para que a fila não cresça indefinidamente — é necessário que a taxa média de chegada seja inferior à taxa máxima de atendimento. Neste caso, a condição de estabilidade é definida por:

\[\begin{align}\\ \lambda < \frac{1}{\bar{x}}\\\\ \end{align}\]

Quando essa condição não é satisfeita, o sistema é dito instável, o comprimento da fila diverge para infinito e qualquer análise em regime de equilíbrio perde a validade. Este requisito de estabilidade conduz naturalmente ao conceito de utilização do servidor.


Definição 5.1 (Nível de Utilização - Bhat 2008; Haviv 2013) O nível de utilização ou intensidade de tráfego do servidor é definido como \(\rho = \lambda \bar{x}\), onde \(\rho\) representa a fração do tempo em que o servidor está ocupado. De forma complementar, \(1 - \rho\) corresponde à fração do tempo em que o servidor permanece ocioso. Note que, à medida que \(\rho\) se aproxima de 1, o sistema opera próximo da sua saturação, o que resulta em aumento significativo dos tempos médios de espera e do comprimento médio da fila.


No cenário em que não há variabilidade nos processos de chegada e de serviço, não ocorre formação de fila. De fato, considere que, ao surgir a primeira chegada, o servidor realiza o atendimento por um intervalo de duração \(\bar{x}\). Em seguida, o servidor permanece ocioso durante um período igual a \(1/\lambda - \bar{x}\) até a chegada do próximo cliente, o qual será atendido por um tempo \(\bar{x}\), e assim sucessivamente. Dessa forma, a probabilidade de o servidor estar ocupado é função do tempo e, em particular, não converge para um valor limite estacionário. Por outro lado, quando há variabilidade nos processos de chegada e/ou de serviço, a existência de uma probabilidade limite torna-se possível. Tal hipótese não apenas reflete as características observadas na maioria dos sistemas reais de filas, mas também é fundamental para viabilizar a análise teórica dos modelos. Sob essa premissa, denotamos por \(\rho\) a probabilidade limite — isto é, no horizonte temporal assintoticamente longo — de que o servidor esteja ocupado, correspondendo, portanto, à fração do tempo em que o sistema encontra-se em operação ativa.


5.4.1. Lei de Little



Ao analisar sistemas em regime de equilíbrio, evidencia-se uma relação universal entre três quantidades fundamentais: a taxa média de chegada, o número médio de unidades presentes no sistema e o tempo médio de permanência dessas unidades. Nesse contexto, considere, por exemplo, um sistema ao qual chegam clientes de forma contínua e indefinida, que permanecem por um intervalo finito antes de partirem, assumindo-se que nenhum cliente permanece no sistema de forma permanente. Denotemos por \(\lambda\) a taxa média limite de chegadas — que coincide com a taxa média limite de saídas —, por \(W\) o tempo médio limite de permanência de cada cliente no sistema e por \(L\) o número médio limite de clientes presentes em um dado instante, pressupondo a existência desses limites. Essa relação conduz-nos diretamente à formulação da lei de Little.


Teorema 5.1 (Lei de Little - Cobham 1954; Little 1961). Considere um sistema em regime estacionário, para o qual são definidas as seguintes quantidades:


  • \(L\): número médio de unidades no sistema (incluindo fila e serviço);
  • \(\lambda\): taxa média de chegadas (e saídas);
  • \(W\): tempo médio que uma unidade passa no sistema (tempo de espera mais tempo de serviço).


Nessas condições, a relação fundamental, conhecida como lei de Little, que caracteriza o sistema é dada por \(L = \lambda W\).


Demonstração. Considere um sistema em regime estacionário, no qual as unidades (clientes, tarefas, etc.) chegam continuamente, permanecem por um tempo finito no sistema e depois partem, sem jamais ficarem retidas permanentemente. Seja \(N(t)\) o número de unidades no sistema no instante \(t\), e seja \(T_i\) o tempo total que a \(i\)-ésima unidade permanece no sistema (tempo de espera na fila somado ao tempo de atendimento). Defina:


  • \(A(t)\) como o número total de unidades que chegaram ao sistema até o instante \(t\);
  • \(S(t)\) como o número total de unidades que saíram do sistema até o instante \(t\).


Como o sistema está em regime estacionário e estável, as taxas médias de chegada e saída coincidem, ou seja, existe o limite:

\[\begin{align}\\ \lambda = \lim_{t \to \infty} \frac{A(t)}{t} = \lim_{t \to \infty} \frac{S(t)}{t}\\\\ \end{align}\]

Além disso, o número médio de unidades no sistema é definido como o valor esperado da quantidade média de unidades presentes, isto é,

\[\begin{align}\\ L = \lim_{T \to \infty} \frac{1}{T} \int_0^T N(t) \, dt\\\\ \end{align}\]

e o tempo médio de permanência de uma unidade no sistema é:

\[\begin{align}\\ W = \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^n T_i\\\\ \end{align}\]

Note que a integral \(\int_0^T N(t) dt\) representa a “área” total acumulada do número de unidades no sistema ao longo do intervalo \([0,T]\). Em outras palavras, é a soma dos tempos de permanência de todas as unidades que estiveram no sistema durante esse intervalo, pois cada unidade contribui com seu tempo de permanência para essa área. Suponha que, durante o intervalo \([0,T]\), tenham chegado \(A(T)\) unidades, com tempos de permanência \(T_1, T_2, \ldots, T_{A(T)}\). Então, a soma dos tempos de permanência dessas unidades é:

\[\begin{align}\\ \sum_{i=1}^{A(T)} T_i\\\\ \end{align}\]

Como cada unidade permanece no sistema por um tempo \(T_i\), o total acumulado de unidades no sistema integrado sobre o tempo deve ser igual a essa soma:

\[\begin{align}\\ \int_0^T N(t) dt = \sum_{i=1}^{A(T)} T_i\\\\ \end{align}\]

Assim, dividindo ambos os lados por \(T\), tem-se que:

\[\begin{align}\\ \frac{1}{T} \int_0^T N(t) dt = \frac{1}{T} \sum_{i=1}^{A(T)} T_i = \frac{A(T)}{T} \cdot \frac{1}{A(T)} \sum_{i=1}^{A(T)} T_i\\\\ \end{align}\]

Logo, tomando o limite quando \(T \to \infty\) e utilizando a definição de \(L\), \(\lambda\) e \(W\), obtém-se que:

\[\begin{align}\\ L = \lim_{T \to \infty} \frac{1}{T} \int_0^T N(t) dt = \lim_{T \to \infty} \frac{A(T)}{T} \cdot \frac{1}{A(T)} \sum_{i=1}^{A(T)} T_i \\\\ \end{align}\]

Como:


  • \(\lim_{T \to \infty} \frac{A(T)}{T} = \lambda\) (taxa média de chegadas);
  • \(\lim_{T \to \infty} \frac{1}{A(T)} \sum_{i=1}^{A(T)} T_i = W\) (tempo médio de permanência).


Concluímos que \(L = \lambda W\), como queríamos demonstrar.

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

A lei de Little foi originalmente justificada por meio de um argumento heurístico relativamente simples. Contudo, sua relevância advém da ampla flexibilidade que proporciona na definição do conceito de “sistema”. Por exemplo, considere-se um fluxo de chegadas a uma fila constituído pela união de dois subfluxos distintos, como o de indivíduos do sexo masculino e do sexo feminino. Independentemente da disciplina de atendimento adotada — que pode tanto favorecer quanto desfavorecer um desses grupos —, o número médio de unidades do subgrupo correspondente será igual ao produto entre sua taxa média de chegada e o tempo médio de permanência desse subgrupo no sistema. Assim, neste contexto, o “sistema” pode ser compreendido exclusivamente como o conjunto de unidades do subgrupo analisado.


Analogamente, uma linha de produção complexa, composta por múltiplos servidores e filas, pode ser concebida como um sistema único. Por exemplo, ao se delimitar pontos lógicos de início e término em uma linha de montagem, a quantidade de trabalho em processo — mensurada pelo número de unidades que podem ser produzidas a partir de determinado estágio, caso não ocorra nova produção — corresponde ao produto entre a taxa média de produtos finalizados por unidade de tempo e o tempo médio transcorrido desde o início da produção até a saída do produto do sistema.


Sua interpretação também pode ser intuída por meio de um argumento econômico: se cada cliente paga uma taxa de 1 unidade monetária por unidade de tempo em que permanece no sistema, o faturamento por unidade de tempo será \(\lambda W\). Alternativamente, se o operador cobra 1 unidade monetária por unidade de tempo para cada cliente presente no sistema, o faturamento será \(L\). Como ambos os métodos devem gerar, em equilíbrio, a mesma receita, conclui-se que \(L = \lambda W\).


5.4.2. Tempo Residual de Serviço



Suponha que os tempos de serviço consecutivos sejam independentes e identicamente distribuídos. Denote-se por \(\overline{x^2}\) o valor esperado do quadrado do tempo de serviço. Considere, neste contexto, exclusivamente o processo de serviço. Em particular, imagine um relógio que permanece pausado sempre que o servidor está ocioso e avança apenas durante os períodos em que o servidor se encontra ocupado. Sob essa contagem de tempo, o processo de serviço configura-se como um processo de renovação, no qual cada ciclo corresponde a um atendimento completo. Em particular, se um cliente é selecionado aleatoriamente enquanto está em atendimento, então — da definição de processos de renovação — a sua idade média de serviço \((A)\), o seu tempo residual médio de serviço (\(R\)) e o seu tempo total de serviço (\(S\)) são, respectivamente,

\[\begin{align}\\ A =\dfrac{\overline{x^2}}{2\bar{x}}, \quad R =\dfrac{\overline{x^2}}{2\bar{x}}, \quad S= \dfrac{\overline{x^2}}{\bar{x}}\\\\ \end{align}\]

Dentre esses parâmetros, o tempo residual médio de serviço assume particular relevância, pois quantifica o tempo médio restante até a conclusão do atendimento em curso, refletindo diretamente quanto tempo ainda será necessário até que o servidor se torne disponível para atender o próximo cliente.


Definição 5.2 (Tempo Residual Médio de Serviço - Haviv 2013) Seja \(X\) a variável aleatória que representa o tempo de serviço, com média \(\bar{x}\) e segundo momento \(\overline{x^2} = \mathrm{E}[X^2]\). O tempo residual médio de serviço, denotado por \(R\), condicionado ao fato de que há um cliente em atendimento no instante considerado, é dado por: \[\begin{align}\\ R = \dfrac{\overline{x^2}}{2\bar{x}}\\\\ \end{align}\] representando o tempo médio restante até a conclusão do serviço em andamento. Este valor corresponde à esperança do tempo residual de um ciclo observado em um instante aleatório no qual o servidor se encontra ocupado.


Considerando que a probabilidade de o servidor estar ocupado é \(\rho = \lambda \bar{x}\) — e que, naturalmente, a quantidade média de trabalho no servidor, condicionada ao fato de que este se encontra ocioso, é igual a zero —, a quantidade média incondicional de trabalho acumulado no servidor, denotada por \(W_s\), é dada por:

\[\begin{align}\\ W_s = (1-\rho)\cdot 0 + \rho \cdot R = \rho \cdot \dfrac{\overline{x^2}}{2\bar{x}} = \lambda \cdot \dfrac{\overline{x^2}}{2}\\\\ \end{align}\]

onde \(R = \overline{x^2}/2\bar{x}\) representa o tempo residual médio de serviço. Observe que, no caso específico em que os tempos de serviço seguem uma distribuição exponencial — cujo segundo momento é \(\overline{x^2} = 2\bar{x}^2\) —, obtém-se:

\[\begin{align}\\ W_s = \lambda \cdot \dfrac{2\bar{x}^2}{2} = \lambda \bar{x}^2 = \rho\bar{x}\\\\ \end{align}\]

Importa destacar que essa expressão não é insensível, pois depende explicitamente do segundo momento dos tempos de serviço, \(\overline{x^2}\), e não apenas da média \(\bar{x}\), como seria requerido para que a métrica fosse insensível à forma da distribuição de serviço. Por fim, ressalta-se que este resultado é válido sob a única hipótese de que a disciplina de atendimento seja work-conserving, isto é, o servidor nunca permanece ocioso enquanto houver clientes aguardando serviço.


5.4.3. Tempo Virtual de Espera



Em sistemas de filas operando sob a disciplina First-In, First-Out (FIFO), compreender o estado instantâneo do sistema é importante para avaliar o desempenho. Uma das métricas mais informativas para esse propósito é o tempo virtual de espera, que indica o tempo que uma unidade hipotética necessitaria aguardar caso chegasse naquele exato instante. Diferentemente do tempo médio de espera — que é uma medida agregada no tempo —, o tempo virtual de espera oferece uma visão pontual da carga acumulada no sistema, refletindo tanto os clientes já enfileirados quanto o atendimento em andamento. Além disso, esse tempo está diretamente relacionado à quantidade total de trabalho presente na fila, isto é, à soma dos tempos de serviço que ainda precisam ser executados antes que o cliente hipotético possa ser atendido. Esse trabalho inclui duas parcelas: (i) o tempo residual de serviço do cliente atualmente em atendimento (se houver) e (ii) os tempos completos de serviço dos clientes na fila.


Definição 5.3 (Tempo Virtual de Espera - Haviv 2013). Seja \(V_q(t)\) o tempo virtual de espera no instante \(t\), definido como a quantidade total de trabalho acumulada na fila naquele momento. Em um sistema FIFO e work-conserving, \(V_q(t)\) corresponde à soma do tempo residual de serviço do cliente em atendimento (se houver) e dos tempos completos de serviço dos clientes aguardando na fila. Assim, assumindo regime estacionário, a relação média entre o tempo virtual de espera \(V_q\), o número médio de clientes na fila \(L_q\), o tempo médio de serviço \(\bar{x}\), a utilização do servidor \(\rho\), e o tempo residual médio de serviço \(R\) é dada por:

\[\begin{align}\\ V_q = L_q \bar{x} + \rho R\\\\ \end{align}\]

em que \(L_q \bar{x}\) representa o trabalho total dos clientes que aguardam na fila (todos com tempos de serviço médios completos \(\bar{x}\)); e \(\rho R\) corresponde ao trabalho residual médio do cliente em atendimento, ponderado pela probabilidade de o servidor estar ocupado. Aplicando a lei de Little na fila, obtém-se:

\[\begin{align}\\ V_q = \lambda W_q \bar{x} + \rho R = \lambda W_q \bar{x} + \rho \frac{\overline{x^2}}{2 \bar{x}}\\\\ \end{align}\]

que pode ser expressa diretamente em função de \(W_q\) como: \(V_q = \rho \left(W_q + R\right)\).


A formulação descrita na Definição 5.3 evidencia que o tempo virtual de espera incorpora dois componentes fundamentais: (i) o trabalho acumulado devido aos clientes aguardando na fila e (ii) o tempo residual médio do serviço em andamento, ponderado pela utilização do servidor. Assim, \(V_q\) representa a carga instantânea que uma unidade hipotética enfrentaria ao ingressar no sistema naquele momento, sob as premissas da disciplina FIFO e da propriedade work-conserving.


5.4.4. Distribuições nos Instantes de Chegada e Saída



A análise do número de clientes no sistema, particularmente em modelos de filas do tipo \(G/G/1\), pode ser conduzida sob distintas perspectivas temporais, cada uma capturando um aspecto específico da dinâmica estocástica do sistema. Em termos gerais, a caracterização da distribuição do número de clientes pode ser realizada (i) em um instante arbitrário de tempo, (ii) no instante imediatamente anterior a uma chegada ou (iii) no instante imediatamente posterior a uma partida. Cada uma dessas abordagens oferece uma visão particular sobre o estado do sistema e, consequentemente, conduz, em geral, a distribuições não equivalentes.


Em sistemas que operam, por exemplo, sob condição de estabilidade — isto é, quando a taxa média de chegada \(\lambda\) é estritamente inferior à taxa média de serviço \(\mu\), garantindo \(\rho = \lambda \bar{x} < 1\) —, torna-se importante compreender as relações entre essas distribuições. Tal compreensão não apenas aprimora a interpretação dos indicadores de desempenho, mas também permite construir modelos mais aderentes à realidade observacional, dependendo da perspectiva adotada (monitoramento contínuo, análise nas épocas de chegada ou nas épocas de partida). Para formalizar essa distinção, definem-se as seguintes variáveis aleatórias associadas ao número de clientes no sistema:


  • \(Q(t)\): número de clientes no sistema em um instante arbitrário de tempo \(t\);
  • \(Q_a(n)\): número de clientes imediatamente antes da \(n\)-ésima chegada;
  • \(Q_d(n)\): número de clientes imediatamente após a \(n\)-ésima partida.


As sequências \(\{Q_a(n)\}_{n \ \geqslant \ 1}\) e \(\{Q_d(n)\}_{n \ \geqslant \ 1}\) constituem cadeias de Markov em tempo discreto, associadas, respectivamente, aos processos de chegada e de partida. Por outro lado, o processo \(\{Q(t)\}_{t \ \geqslant \ 0}\) é um processo estocástico em tempo contínuo, que descreve a evolução dinâmica do sistema entre os eventos de chegada e de partida.


Definição 5.4 (Distribuições nas Épocas de Chegada e de Partida - Bhat 2008; Haviv 2013). Sejam \(Q_a(n)\) o número de clientes no sistema imediatamente antes da \(n\)-ésima chegada e \(Q_d(n)\) o número de clientes imediatamente após a \(n\)-ésima partida. Diz-se que o sistema se encontra em regime estacionário quando as distribuições de \(Q_a(n)\) e \(Q_d(n)\) convergem, respectivamente, para distribuições limite, denotadas por \(Q_a\) e \(Q_d\), tais que

\[\begin{align}\\ \lim_{n \to \infty} \mathrm{P}(Q_a(n) = k) &= \mathrm{P}(Q_a = k)\\\\ \lim_{n \to \infty} \mathrm{P}(Q_d(n) = k) &= \mathrm{P}(Q_d = k), \quad \forall k \in \mathbb{N}_0\\\\ \end{align}\]

Sob a hipótese de estabilidade e assumindo que a disciplina de atendimento seja work-conserving — isto é, o servidor jamais permanece ocioso quando há clientes na fila —, verifica-se que as distribuições nas épocas de chegada e de partida coincidem no regime estacionário, ou seja,

\[\begin{align}\\ \mathrm{P}(Q_a = k) = \mathrm{P}(Q_d = k)\\\\ \end{align}\]

para todo \(k \in \mathbb{N}_0\). Essa equivalência reflete diretamente o princípio de conservação de fluxo no sistema: no equilíbrio, a estatística do número de clientes observada imediatamente antes de uma chegada é idêntica à estatística observada imediatamente após uma partida. Isso decorre do fato de que, em regime estacionário, a taxa média de chegadas deve, necessariamente, igualar-se à taxa média de saídas.


É importante destacar que a distribuição limite — comum às épocas de chegada e de partida — formalizada na Definição 5.4 não coincide, em geral, com a distribuição do número de clientes observada em um instante arbitrário de tempo, representada por \(Q(t)\). Esta diferença decorre de um fenômeno conhecido como viés de tamanho (size-biasing). Intuitivamente, estados nos quais o sistema se encontra mais carregado tendem a ter durações médias mais longas. Consequentemente, um observador que realiza amostragens em tempo contínuo possui maior probabilidade de encontrar o sistema em estados mais congestionados, quando comparado a um observador que efetua amostragens restritas aos eventos discretos de chegada ou de partida.


Essa distinção entre amostragem contínua no tempo e amostragem nos pontos discretos dos eventos é, portanto, de natureza estrutural e não acessória. De fato, salvo em casos muito específicos — notadamente no modelo \(M/M/1\), no qual a propriedade de Poisson das chegadas, combinada à falta de memória da distribuição exponencial dos tempos de serviço, elimina completamente o efeito de viés —, verifica-se que:

\[\begin{align}\\ \mathrm{P}(Q(t) = k) \neq \mathrm{P}(Q_a = k) = \mathrm{P}(Q_d = k), \quad \forall k\in \mathbb{N}_0\\\\ \end{align}\]

Essa distinção não constitui uma sutileza meramente teórica, mas possui consequências práticas diretas na avaliação de desempenho de sistemas reais. Por exemplo, métricas calculadas sob a perspectiva de um cliente que observa o sistema ao chegar — como o tempo médio de espera ou a probabilidade de enfrentar uma fila vazia — podem ser significativamente distintas das métricas obtidas por um operador que realiza monitoramento contínuo do sistema ao longo do tempo, particularmente no que diz respeito à ocupação média, à probabilidade de congestionamento e ao tempo médio em sistema.


5.4.5. Propriedade ASTA e Fórmula de Khintchine–Pollaczek



Vimos, na Definição 5.3, que, em geral, o tempo virtual de espera na fila, denotado por \(V_q\), e o tempo médio de espera na fila observado por um cliente na sua chegada, denotado por \(W_q\), não são, em princípio, idênticos, uma vez que essa diferença decorre diretamente da natureza das médias envolvidas: \(W_q\) representa uma média condicional, calculada nos instantes aleatórios de chegada dos clientes, enquanto \(V_q\) corresponde a uma média temporal, avaliada sobre o estado do sistema em instantes arbitrários, escolhidos independentemente do processo de chegada.


Contudo, apesar dessa distinção conceitual fundamental, existem classes de modelos de filas para as quais esses dois valores coincidem exatamente, sendo essa coincidência formalizada por meio de uma propriedade estrutural denominada Arrival Sees Time Averages (ASTA), que é uma generalização do princípio clássico PASTA (Poisson Arrivals See Time Averages) para contextos nos quais o processo de chegada não necessariamente segue uma distribuição Poisson, mas em que certas condições de independência e estacionaridade são satisfeitas. Em essência, a propriedade ASTA assegura que a distribuição do sistema observada por um cliente no instante de sua chegada coincide exatamente com a distribuição do estado do sistema em um instante arbitrário no tempo, implicando, consequentemente, que o tempo médio de espera na fila na perspectiva do cliente que chega (\(W_q\)) coincide com o tempo virtual de espera (\(V_q\)), que reflete o trabalho médio presente na fila no instante da observação.


Definição 5.5 (Propriedade ASTA - Haviv 2013). Diz-se que um sistema de filas satisfaz a propriedade Arrival Sees Time Averages (ASTA) se, no regime estacionário, o tempo médio de espera na fila observado por um cliente na sua chegada coincide exatamente com o tempo virtual de espera na fila, isto é, se \(W_q = V_q\).


Sistemas que satisfazem a propriedade ASTA admitem formulações analíticas mais simples e elegantes para diversas métricas de desempenho. Em particular, no contexto de sistemas de servidor único, a validade de ASTA permite expressar o tempo médio de espera na fila por meio da fórmula de Khintchine–Pollaczek (K-P), a qual estabelece uma relação direta entre o tempo médio de espera, a intensidade de tráfego no sistema e a variabilidade dos tempos de serviço.


Teorema 5.2 (Fórmula de Khintchine–Pollaczek - Pollaczek 1930). Considere um sistema de filas de servidor único, estável e work-conserving, que satisfaz a propriedade ASTA. Seja \(X\) a variável aleatória que representa o tempo de serviço, com média \(\bar{x} = \mathrm{E}[X]\) e segundo momento \(\overline{x^2} = \mathrm{E}[X^2]\). Seja \(\lambda\) a taxa média de chegadas e \(\rho = \lambda \bar{x}\) a intensidade de tráfego, com \(\rho < 1\). Nessas condições, o tempo médio de espera na fila, denotado por \(W_q\), e, por consequência, o tempo virtual de espera, denotado por \(V_q\) (em virtude da propriedade ASTA), são dados por:

\[\begin{align}\\ W_q = V_q = \frac{\lambda \overline{x^2}}{2(1 - \rho)}\\\\ \end{align}\]

Dessa relação, tem-se que tempo médio de espera na fila cresce de forma hiperbólica à medida que \(\rho \to 1\), ilustrando não apenas o impacto da carga do sistema, mas também a variabilidade dos tempos de serviço — capturada pelo segundo momento \(\overline{x^2}\). Isto é, quanto maior a variabilidade, maior tende a ser o tempo médio de espera, mesmo sob uma mesma carga média \(\rho\).


Demonstração. Considere um sistema de filas de servidor único, work-conserving, operando sob condições de estabilidade (\(\rho = \lambda \bar{x} < 1\)) e que satisfaz a propriedade ASTA. Por definição, o tempo virtual de espera \(V_q\) é composto por dois termos: o trabalho acumulado na fila, proveniente dos clientes que aguardam atendimento; e o tempo residual de serviço do cliente em atendimento (se houver), ponderado pela probabilidade de que o servidor esteja ocupado, que é \(\rho\). Isto é, \(V_q\) é dado por:

\[\begin{align}\\ V_q = L_q \bar{x} + \rho R\\\\ \end{align}\]

em que:


  • \(L_q\) é o número médio de clientes na fila (excluindo o atendimento em andamento);
  • \(\bar{x}\) é o tempo médio de serviço;
  • \(R\) é o tempo residual médio de serviço, dado por:


Aplicando a lei de Little na fila \(L_q = \lambda W_q\), obtém-se que:

\[\begin{align}\\ V_q = \lambda W_q \bar{x} + \frac{\rho\overline{x^2}}{2 \bar{x}}\\\\ \end{align}\]

conforme descrito na Definição 5.3. Agora, dada a hipótese da propriedade ASTA, tem-se que \(V_q = W_q\). Então,

\[\begin{align}\\ W_q = \lambda W_q \bar{x} + \frac{\rho\overline{x^2}}{2 \bar{x}}\\\\ \end{align}\]

Isolando \(W_q\) nessa equação, obtém-se que:

\[\begin{align}\\ W_q = \dfrac{\lambda \overline{x^2}}{2(1 - \rho)}\\\\ \end{align}\]

como queríamos demonstrar.

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

5.5. Modelos Markovianos de Filas


Sistemas de filas representam, na essência, processos estocásticos nos quais entidades — tais como clientes, pacientes ou requisições — chegam aleatoriamente, aguardam serviço e, após serem atendidas, deixam o sistema. Sob certas hipóteses sobre os mecanismos de chegada e de atendimento, é possível modelar uma ampla classe desses sistemas como processos markovianos de nascimento e morte em tempo contínuo, cuja dinâmica é descrita por transições entre estados discretos (\(\mathcal{E_1, \mathcal{E}_2},\ldots\)) ao longo do tempo. Nessa formulação, o estado do sistema é representado pelo número de clientes presentes no sistema em um dado instante, incluindo aqueles em atendimento e, eventualmente, na fila. O processo evolui por meio de dois tipos de eventos estocásticos:


  • Chegadas, que correspondem a nascimentos, incrementando o número de clientes no sistema em uma unidade;
  • Serviços, que correspondem a mortes, reduzindo o número de clientes em uma unidade.


Seja \(n\) o número de clientes no sistema no instante \(t\). As dinâmicas de chegada e de atendimento são descritas pelas funções de taxa:


  • \(\lambda_n\) como a taxa de chegada, quando o sistema se encontra no estado \(n\);
  • \(\mu_n\) como a taxa de serviço, também dependente do estado atual \(n\).


Sob essas premissas, as regras que governam a dinâmica do processo podem ser formalizadas da seguinte maneira:


  • Mecanismo de Chegadas: No o intervalo infinitesimal \((t, t + \Delta t]\), a probabilidade de que nenhuma chegada ocorra é descrita por: \[\begin{align}\\ 1 - \lambda_n \Delta t + o(\Delta t)\\\\ \end{align}\] Já a probabilidade de que ocorra exatamente uma chegada é de \(\lambda_n \Delta t + o(\Delta t)\), enquanto a probabilidade de ocorrência de duas ou mais chegada é descrita por \(o(\Delta t)\). Isto é, \[\begin{align}\\ P(\text{uma chegada}) &= \lambda_n \Delta t + o(\Delta t) \\\\ P(\text{nenhuma chegada}) &= 1 - \lambda_n \Delta t + o(\Delta t) \\\\ P(\text{mais de uma chegada}) &= o(\Delta t)\\\\ \end{align}\] em que \(o(\Delta t)\) é tal que \(o(\Delta t)/\Delta t \to 0\) quando \(\Delta t \to 0\). Adicionalmente, assume-se que os tempos entre chegadas são exponencialmente distribuídos e independentes, de modo que o instante da próxima chegada é independente do tempo já decorrido desde a última chegada — isto é, o mecanismo de chegada satisfaz a propriedade de falta de memória.


  • Mecanismo de Serviço: De modo análogo ao mecanismo de chegada, tem-se que, no intervalo infinitesimal \((t, t + \Delta t]\), a probabilidade de que nenhuma conclusão de serviço ocorra é descrita por: \[\begin{align}\\ 1 - \mu_n \Delta t + o(\Delta t)\\\\ \end{align}\] Por outro lado, a probabilidade de ocorrer exatamente uma conclusão de serviço é \(\mu_n \Delta t + o(\Delta t)\), enquanto a probabilidade de que haja duas ou mais conclusões de serviço ocorra é \(o(\Delta t)\). Isto é, \[\begin{align}\\ P(\text{uma conclusão de serviço}) &= \mu_n \Delta t + o(\Delta t) \\\\ P(\text{nenhuma conclusão de serviço}) &= 1 - \mu_n \Delta t + o(\Delta t) \\\\ P(\text{mais de uma conclusão de serviço}) &= o(\Delta t)\\\\ \end{align}\] em que \(o(\Delta t)\) é tal que \(o(\Delta t)/\Delta t \to 0\) quando \(\Delta t \to 0\). Analogamente, os tempos de serviço também são exponencialmente distribuídos e independentes, garantindo que o instante de conclusão de um serviço não dependa do tempo decorrido desde seu início — isto é, o mecanismo de serviço também satisfaz a propriedade de falta de memória.


Seja \(Q(t)\) o número de clientes no sistema no instante \(t\). Defina:

\[\begin{align}\\\ P_{i,n}(t) = P[Q(t) = n \mid Q(0) = i]\\\\ \end{align}\]

como a probabilidade de que o sistema esteja no estado \(n\) no tempo \(t\), dado que começou no estado \(i\) no instante inicial. Logo, incorporando as probabilidades de transição durante o intervalo infinitesimal \((t, t + \Delta t]\), conforme as leis de evolução descritas anteriormente, obtém-se:

\[\begin{align}\\ P_{n,n+1}(\Delta t) &= \lambda_n \Delta t + o(\Delta t), & n = 0, 1, 2, \ldots \\\\ P_{n,n-1}(\Delta t) &= \mu_n \Delta t + o(\Delta t), & n = 1, 2, 3, \ldots \\\\ P_{n,n}(\Delta t) &= 1 - \lambda_n \Delta t - \mu_n \Delta t + o(\Delta t), & n = 1, 2, 3, \ldots \\\\ P_{n,j}(\Delta t) &= o(\Delta t), & j \neq n-1, n, n+1\\\\ \end{align}\]

A partir dessas equações de transição, constrói-se, então, o gerador infinitesimal — denotado por \(A\), e também denominado matriz de taxas infinitesimais — associado ao processo markoviano de nascimento e morte que modela o sistema de filas. Essa matriz \(A\), em particular, possui a seguinte estrutura tridiagonal:

\[\begin{align}\\ A = \begin{bmatrix} -\lambda_0 & \lambda_0 & 0 & 0 & \cdots \\ \mu_1 & -(\lambda_1 + \mu_1) & \lambda_1 & 0 & \cdots \\ 0 & \mu_2 & -(\lambda_2 + \mu_2) & \lambda_2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}\\\\ \end{align}\]

A partir de \(A\), pode-se, então, formular o sistema de equações diferenciais de Kolmogorov prospectivas, que descrevem a evolução temporal das probabilidades de ocupação dos estados. Estas equações, também denominadas equações de equilíbrio transitório, assumem, neste caso, a seguinte forma:

\[\begin{align}\\ P_0'(t) &= -\lambda_0 P_0(t) + \mu_1 P_1(t)\\\\ P_n'(t) &= -(\lambda_n + \mu_n) P_n(t) + \lambda_{n-1} P_{n-1}(t) + \mu_{n+1} P_{n+1}(t), \quad n = 1, 2, \ldots\\\\ \end{align}\]

onde \(P_n(t) = P_{i,n}(t)\) representa a probabilidade de o processo estar no estado \(n\) no tempo \(t\), dado que começou no estado inicial \(i\). Assim, para determinar \(P_n(t)\) (equivalente a \(P_{i,n}(t)\)), as equações acima devem ser resolvidas juntamente com a condição inicial:

\[\begin{align}\\ P_n(0) = \begin{cases} 1 & \text{se} \quad n = 0\\ 0 & \text{se} \quad n \neq 0 \end{cases}\\\\ \end{align}\]

Mesmo em configurações simples, como no caso em que \(\lambda_n = \lambda\) e \(\mu_n = \mu\) para todo \(n = 0, 1, 2, \dots\) — isto é, quando as chegadas seguem um processo de Poisson e os tempos de serviço são exponenciais (modelo \(M/M/1\)) —, a obtenção de uma expressão analítica fechada para \(P_n(t)\) é um processo complexo. Na maioria das aplicações, no entanto, o comportamento transitório — isto é, a evolução temporal da distribuição de probabilidade — não constitui o principal foco de interesse. A análise de maior relevância, tanto sob a perspectiva teórica quanto prática, concentra-se no comportamento assintótico do processo, isto é, na caracterização da distribuição estacionária, a qual é obtida considerando o limite das equações diferenciais de Kolmogorov quando \(t \to \infty\). No regime estacionário, o comportamento do processo torna-se independente tanto do tempo quanto do estado inicial. Isto é,

\[\begin{align}\\ \lim_{t \to \infty} P_{i,n}(t) = \pi_n, \quad n = 0, 1, 2, \ldots\\\\ \end{align}\]

o que implica que:

\[\begin{align}\\ P_n'(t) \to 0 \quad \text{quando} \quad t \to \infty\\\\ \end{align}\]

Substituindo essa condição limite nas equações diferenciais de Kolmogorov, obtém-se o seguinte sistema de equações de equilíbrio:

\[\begin{align}\\ 0 &= -\lambda_0 \pi_0 + \mu_1 \pi_1, \\\\ 0 &= -(\lambda_n + \mu_n) \pi_n + \lambda_{n-1} \pi_{n-1} + \mu_{n+1} \pi_{n+1}, \quad n=1,2,\ldots\\\\ \end{align}\]

Esse Sistema, em particular, pode ser resolvido de forma recursiva. Isto é, a partir da primeira equação, tem-se:

\[\begin{align}\\ \pi_1 = \frac{\lambda_0}{\mu_1} \pi_0\\\\ \end{align}\]

Para \(n = 1\), a segunda equação resulta em:

\[\begin{align}\\ (\lambda_1 + \mu_1) \pi_1 = \lambda_0 \pi_0 + \mu_2 \pi_2\\\\ \end{align}\]

Usando o resultado para \(\pi_1\), tem-se que:

\[\begin{align}\\ \pi_2 = \frac{\lambda_1 \lambda_0}{\mu_2 \mu_1} \pi_0\\\\ \end{align}\]

Prosseguindo por indução para \(n = 2, 3, \dots\), tem-se a relação de recorrência:

\[\begin{align}\\ \mu_n \pi_n = \lambda_{n-1} \pi_{n-1}\\\ \end{align}\]

que conduz, de forma geral, à solução:

\[\begin{align}\\ \pi_n = \frac{\lambda_0 \lambda_1 \cdots \lambda_{n-1}}{\mu_1 \mu_2 \cdots \mu_n} \pi_0\\\\ \end{align}\]

Aplicando a condição de normalização, \(\sum_{n \ \in \ \mathcal{E}} p_n = 1\), obtém-se que:

\[\begin{align}\\ \pi_0 = \left[ 1 + \sum_{n=1}^\infty \frac{\lambda_0 \lambda_1 \cdots \lambda_{n-1}}{\mu_1 \mu_2 \cdots \mu_n} \right]^{-1}\\\\ \end{align}\]

Portanto, a distribuição estacionária do número de clientes no sistema é dada por \(\{\pi_n, n = 0, 1, 2, \ldots \}\). É importante destacar que essa distribuição é propriamente definida — ou seja, os valores de \(\pi_n\) são finitos e somam 1 — se, e somente se, a seguinte condição de convergência é satisfeita:

\[\begin{align}\\ 1 + \sum_{n=1}^\infty \frac{\lambda_0 \lambda_1 \cdots \lambda_{n-1}}{\mu_1 \mu_2 \cdots \mu_n} < \infty\\\\ \end{align}\]

Essa condição estabelece, de forma rigorosa, o critério de estabilidade do sistema, garantindo que o processo de nascimento e morte associado ao modelo de filas não exploda ao longo do tempo e, consequentemente, que uma distribuição de equilíbrio exista.


5.5.1. Modelo M/M/1



O modelo \(M/M/1\) é o mais simples dentro da teoria de filas e serve como base para a compreensão de sistemas mais complexos. Este modelo representa um sistema com uma única fila de espera e um único servidor, no qual as chegadas ocorrem de acordo com um processo de Poisson e os tempos de serviço seguem uma distribuição exponencial. As principais suposições do modelo \(M/M/1\) são:


  • Processo de Chegada: As chegadas ocorrem segundo um processo de Poisson com taxa \(\lambda\), o que implica que o número de chegadas no intervalo \((0, t]\) é uma variável aleatória \(N(t)\) com distribuição de Poisson de parâmetro \(\lambda t\): \[\begin{align}\\ \mathrm{P}[N(t) = j] = e^{-\lambda t} \frac{(\lambda t)^j}{j!}, \quad j = 0, 1, 2, \ldots\\\\ \end{align}\] Além disso, os tempos entre chegadas consecutivas são independentes e identicamente distribuídos segundo uma distribuição exponencial com parâmetro \(\lambda\), cuja função densidade de probabilidade é dada por: \[\begin{align}\\ a(x) = \lambda e^{-\lambda x}, \quad x > 0\\\\ \end{align}\]

  • Processo de Serviço: Os tempos de serviço são independentes e seguem uma distribuição exponencial com taxa \(\mu\), com função densidade de probabilidade definida por: \[\begin{align}\\ b(x) = \mu e^{-\mu x}, \quad x > 0\\\\ \end{align}\]

  • Disciplina de Atendimento: FIFO (First-In, First-Out).

  • Capacidade da Fila: Ilimitada, ou seja, não há perdas.

  • População de Clientes: Infinita, de modo que a taxa de chegada não depende do número de clientes presentes.


Com essas suposições, tem-se que o tempo médio entre chegadas e o tempo médio de serviço são descritos, respectivamente, como:

\[\begin{align}\\ E[\text{tempo entre chegadas}] &= \frac{1}{\lambda} = \frac{1}{\text{taxa de chegada}}\\\\ E[\text{tempo de serviço}] &= \frac{1}{\mu} = \frac{1}{\text{taxa de serviço}}\\\\ \end{align}\]

O parâmetro que governa o comportamento do sistema é a intensidade de tráfego (\(\rho\)) (ou fator de utilização), definido por:

\[\begin{align}\\ \rho = \frac{\text{taxa de chegada}}{\text{taxa de serviço}} = \frac{\lambda}{\mu}\\\\ \end{align}\]

Para que o sistema atinja o equilíbrio — isto é, não ocorra crescimento indefinido do número de clientes na fila —, é necessário que \(\rho < 1\), ou seja, a taxa de serviço deve exceder a taxa de chegada.


(I) Equações de Kolmogorov


O processo estocástico \(\{Q(t)\}_{t \geqslant 0}\), que representa o número de clientes no sistema em um instante \(t\), que descreve o modelo \(M/M/1\) é um caso especial do processo de nascimento e morte, com taxas constantes, \(\lambda_n = \lambda\), \(\mu_n = \mu \,\, (n\geqslant 1)\), \(\mu_0 = 0\), e espaço de estados \(\mathcal{E} = \{0, 1, 2, \ldots \}\). Para este processo, as equações de Kolmogorov prospectivas, que descrevem a evolução das probabilidades de transição \(P_n(t) = \mathrm{P}[Q(t) = n]\), são:

\[\begin{align}\\ P_0'(t) &= -\lambda P_0(t) + \mu P_1(t)\\\\ P_n'(t) &= -(\lambda + \mu) P_n(t) + \lambda P_{n-1}(t) + \mu P_{n+1}(t), \quad n = 1, 2, \ldots\\\\ \end{align}\]

com \(P_n(0) = 1\) quando \(n = i\) e \(P_n(0)= 0\) caso contrário. Para uma solução dessas equações diferenciais, é necessário o uso de transformadas de Laplace, ou de métodos computacionais, seja pela integração numérica direta das equações diferenciais, seja por meio de simulação estocástica do processo \(\{Q(t)\}_{t \ \geqslant \ 0}\).


(II) Distribuição Estacionária


Assumindo \(\rho < 1\), tem-se que as probabilidades limite são denotadas por \(\pi_n = \lim_{t \to \infty} P_n(t)\) e satisfazem as equações de equilíbrio:

\[\begin{align}\\ \lambda \pi_0 &= \mu \pi_1, \qquad\qquad\qquad \,\,\,\, n = 0\\\\ (\lambda + \mu) \pi_n &= \lambda \pi_{n-1} + \mu \pi_{n+1}, \qquad n = 1, 2, 3, \ldots\\\\ \end{align}\]

A solução geral desse sistema, considerando a condição de normalização \(\sum_{n=0}^\infty \pi_n = 1\), é:

\[\begin{align}\\ \pi_n = (1 - \rho) \rho^n, \quad n = 0, 1, 2, \ldots\\\\ \end{align}\]

em que \(\rho = \lambda / \mu < 1\). Portanto, a distribuição estacionária do número de clientes no sistema para o modelo \(M/M/1\) é uma distribuição geométrica. Uma consequência direta dessa distribuição é a probabilidade de ocupação do servidor, isto é, a probabilidade de que o sistema não se encontre no estado vazio, que é expressa por:

\[\begin{align}\\ \text{P(servidor ocupado)} &= 1 - \pi_0 \\\\ &= 1 - (1 - \rho) \\\\ &= \rho\\\\ \end{align}\]

que é igual a intensidade de tráfego, que, no modelo \(M/M/1\), é interpretada como a fração do tempo, no equilíbrio, em que o servidor permanece ocupado, além de representar a razão entre a taxa média de chegada \(\lambda\) e a taxa média de serviço \(\mu\).


(III) Número Médio de Clientes


Recorde, agora, que definimos \(Q(t)\) como o número de clientes no sistema no instante \(t\). Então, no regime estacionário, defina \(Q = Q(\infty)\) e seja \(Q_q\) o número de clientes na fila, excluindo aquele em serviço. Duas quantidades nos interessam:


  • \(L = E[Q]\), que representa o número médio de clientes no sistema.
  • \(L_q = E[Q_q]\), que representa o número médio de clientes na fila.


A partir da distribuição estacionária \(\pi_n = (1 - \rho)\rho^n\), \(n = 0, 1, 2, \dots\), obtém-se que o número médio de clientes no sistema é descrito por:

\[\begin{align}\\ L = E[Q] = \sum_{n=1}^\infty n (1-\rho) \rho^n = \frac{\rho}{1-\rho}\\\\ \end{align}\]

o que pode ser alternativamente expresso em função das taxas de chegada e de serviço como:

\[\begin{align}\\ L = E[Q] = \frac{\lambda}{\mu - \lambda}\\\\ \end{align}\]

De forma análoga, o número médio de clientes na fila, \(L_q\), é obtido como:

\[\begin{align}\\ L_q = E[Q_q] &= \sum_{n=1}^\infty (n-1) \pi_n \\\\ &= \sum_{n=1}^\infty n \pi_n - \sum_{n=1}^\infty \pi_n \\\\ &= L - \rho = \frac{\rho^2}{1-\rho} \\\\ &= \frac{\lambda^2}{\mu(\mu - \lambda)}\\\\ \end{align}\]

Observe que a intensidade de tráfego \(\rho\) representa, no regime estacionário, a probabilidade de que o servidor esteja ocupado, ou equivalentemente, o número esperado de clientes em serviço. Dessa forma, a decomposição do número médio total de clientes no sistema se dá pela relação:

\[\begin{align}\\ L = L_q +\rho\\\\ \end{align}\]

que expressa, de forma direta, que o número médio total de clientes é a soma do número médio de clientes na fila com o número médio em serviço. Além dos valores médios, uma medida adicional relevante na análise de desempenho do sistema é a variância do número de clientes no sistema, denotada por \(V(Q) = \mathrm{Var}(Q)\). Da distribuição estacionária, obtemos que a variância é dada por:

\[\begin{align}\\ V(Q) = \frac{\rho}{(1-\rho)^2} = \frac{\lambda \mu}{(\mu - \lambda)^2}\\\\ \end{align}\]

(IV) Tempo de Espera do Cliente


Do ponto de vista do cliente, duas métricas de interesse são importantes na análise de desempenho: o tempo total no sistema e o tempo de espera na fila. No regime estacionário, denote por \(T\) o tempo total que um cliente passa no sistema (incluindo espera na fila e serviço) e por \(T_q\) o tempo de espera na fila até o início do serviço. Sob a disciplina FIFO, o tempo de espera na fila \(T_q\) corresponde ao tempo necessário para completar o atendimento de todos os clientes que estão no sistema no momento da chegada. Já o tempo total no sistema \(T\) é simplesmente a soma de \(T_q\) com o tempo de serviço. Como os tempos de serviço são assumidos independentes e distribuídos exponencialmente com taxa \(\mu\), o tempo total de serviço de \(n\) clientes possui distribuição Erlang de ordem \(n\) com função densidade de probabilidade definida por:

\[\begin{align}\\ f_n(x) = e^{-\mu x} \frac{\mu^n x^{n-1}}{(n-1)!}\\\\ \end{align}\]

Seja \(F_q(t) = \mathrm{P}(T_q \leqslant t)\) a função distribuição acumulada do tempo de espera na fila. Note, inicialmente, que:

\[\begin{align}\\ F_q(0) &= \mathrm{P}(T_q = 0) \\\\ &= \mathrm{P}(Q = 0) \\\\ &= 1 - \rho\\\\ \end{align}\]

ou seja, o cliente não espera na fila se, no instante de sua chegada, o sistema está vazio. Este salto na função distribuição é consequência direta da natureza discreta do processo de chegadas e da ausência de espera quando não há clientes no sistema. Agora, devido à propriedade de falta de memória da distribuição exponencial, o tempo residual do serviço em andamento (quando há) é também exponencial com parâmetro \(\mu\). Para \(t > 0\), a probabilidade de que o tempo de espera \(T_q\) esteja no intervalo \((t, t + dt]\) é dada por:

\[\begin{align}\\ dF_q(t) = \sum_{n=1}^\infty p_n e^{-\mu t} \frac{\mu^n t^{n-1}}{(n-1)!} dt \\\\ \end{align}\]

onde o termo Erlang corresponde ao tempo necessário para atender \(n\) clientes (o cliente em serviço, se houver, mais os \(n-1\) na fila). Substituindo a distribuição estacionária \(\pi_n = (1 - \rho) \rho^n\), obtemos que:

\[\begin{align}\\ dF_q(t) = (1-\rho) \sum_{n=1}^\infty \rho^n e^{-\mu t} \frac{\mu^n t^{n-1}}{(n-1)!} dt\\\\ \end{align}\]

Usando a expansão em série de Taylor da função exponencial, a expressão acima se simplifica para:

\[\begin{align}\\ dF_q(t) = \lambda (1-\rho) e^{-\mu (1-\rho) t} dt\\\\ \end{align}\]

Devido à descontinuidade em 0 na distribuição de \(T_q\), temos, via integral, a função distribuição acumulada do tempo de espera na fila é:

\[\begin{align}\\ F_q(t) = \mathrm{P}(T_q = 0) + \int_0^t dF_q(s) = 1 - \rho e^{-\mu (1-\rho) t}\\\\ \end{align}\]

A partir dessa função distribuição, calcula-se facilmente a esperança e a variância de \(T_q\). De fato, o valor esperado do tempo de espera na fila é dado por:

\[\begin{align}\\ W_q = E(T_q) = \frac{\rho}{\mu (1-\rho)} = \frac{\lambda}{\mu (\mu - \lambda)}\\\\ \end{align}\]

enquanto a variância é dada por:

\[\begin{align}\\ V(T_q) = \frac{\rho (2 - \rho)}{\mu^2 (1-\rho)^2}\\\\ \end{align}\]

Por outro lado, dado que o tempo total no sistema \(T\) é a soma do tempo de espera na fila \(T_q\) e do tempo de serviço, que é exponencial com parâmetro \(\mu\), obtemos que que o valor esperado de \(T\) é descrito por:

\[\begin{align}\\ W &= E[T]\\\\ &= E[T_q] + E[\text{serviço}] \\\\ &= \frac{\lambda}{\mu (\mu - \lambda)} + \frac{1}{\mu} \\\\ &= \dfrac{1}{\mu(1 - \rho)}\\\\ \end{align}\]

enquanto a variância de \(T\) é dada por:

\[\begin{align}\\ V(T) &= V(T_q) + V(\text{serviço}) \\\\ &= \frac{\rho(2 - \rho)}{\mu^2 (1 - \rho)^2} + \frac{1}{\mu^2} \\\\ &= \dfrac{1}{\mu^2(1 - \rho)^2}\\\\ \end{align}\]

(V) Período de Ocupação


O período de ocupação é definido como o intervalo de tempo durante o qual o servidor permanece continuamente ocupado, isto é, o tempo decorrido desde a transição do sistema do estado vazio (zero clientes) para qualquer estado ocupado, até o momento em que o sistema retorna novamente ao estado vazio. A cada término de um período de ocupação, inicia-se um período de ociosidade, correspondente ao tempo em que o sistema permanece vazio até a próxima chegada. No modelo \(M/M/1\), como as chegadas ocorrem segundo um processo de Poisson com taxa \(\lambda\), o período de ociosidade é uma variável aleatória exponencial com parâmetro \(\lambda\).


Diversos métodos podem ser empregados para a derivação da distribuição do período de ocupação no modelo \(M/M/1\), embora nenhum deles seja trivial. Neste texto, será adotada uma abordagem baseada nas equações de Kolmogorov prospectivas associadas ao processo de Markov subjacente. Nesta abordagem, o período de ocupação corresponde ao tempo até a primeira visita ao estado \(0\), condicionado a que o processo tenha sido iniciado no estado \(1\). Trata-se, portanto, de um problema clássico de tempo até absorção, em que o estado \(0\) é transformado em um estado absorvente, enquanto os demais estados formam uma classe transitória irredutível.


Para formalizar este problema, constrói-se, então, uma versão modificada da matriz geradora infinitesimal do processo, na qual as transições que partem do estado \(0\) são eliminadas (tornando-o absorvente), e as dinâmicas entre os estados \(n \geqslant 1\) são preservadas. A matriz geradora modificada assume, portanto, a seguinte forma:

\[\begin{align}\\ A = \begin{bmatrix} 0 & 0 & 0 & 0 & \cdots \\ \mu & -(\lambda + \mu) & \lambda & 0 & \cdots \\ 0 & \mu & -(\lambda + \mu) & \lambda & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}\\\\ \end{align}\]

As equações de Kolmogorov prospectivas associadas às probabilidades \(P_n(t) = \mathrm{P}[Q(t) = n]\) são:

\[\begin{align}\\ P_0'(t) &= \mu P_1(t)\\\\ P_1'(t) &= -(\lambda + \mu) P_1(t) + \mu P_2(t)\\\\ P_n'(t) &= -(\lambda + \mu) P_n(t) + \lambda P_{n-1}(t) + \mu P_{n+1}(t), \quad n = 2, 3, \ldots\\\\ \end{align}\]

com condição inicial \(P_1(0) = 1\) e \(P_n(0) = 0\) para \(n \neq 1\). Para determinar a solução deste sistema diferencial, é necessário o uso transformadas de Laplace. Neste caso, seja \(\gamma(\theta)\) a transformada de Laplace da função densidade do período de ocupação, definida como:

\[\begin{align}\\ \gamma(\theta) = \int_0^\infty e^{-\theta t} P_0'(t) \, dt, \quad \mathrm{Re}(\theta) > 0\\\\ \end{align}\]

que, após operações adequadas, se reduz a expressão:

\[\begin{align}\\ \gamma(\theta) = \frac{1}{2\lambda} \left[ \theta + \lambda + \mu - \sqrt{(\theta + \lambda + \mu)^2 - 4 \lambda \mu} \right]\\\\ \end{align}\]

Note que a inversão de \(\gamma(\theta)\) fornece, em particular, a função densidade de probabilidade do período de ocupação, isto é,

\[\begin{align}\\ P_0'(t) = e^{-(\lambda + \mu) t} \sqrt{\frac{\mu}{\lambda}} \, t \, I_1\left( 2 \sqrt{\lambda \mu} t \right)\\\\ \end{align}\]

em que \(I_1(\cdot)\) é a função de Bessel modificada de primeira ordem (Abramowitz e Stegun 1965), definida por:

\[\begin{align}\\ I_j(x) = \sum_{n=0}^\infty \frac{(x/2)^{2n + j}}{n! (n + j)!}\\\\ \end{align}\]

Uma formulação alternativa para \(P_0'(t)\), obtida por argumentos combinatórios, é descrita por:

\[\begin{align}\\ P_0'(t) = e^{-(\lambda + \mu) t} \sum_{n=1}^\infty \frac{\lambda^{n-1} \mu^n t^{2n-2}}{n! (n-1)!}\\\\ \end{align}\]

Agora, defina \(B\) como a variável aleatória que representa a duração do período de ocupação no modelo \(M/M/1\), isto é, o tempo decorrido desde a transição do estado vazio (0 clientes) para o estado ocupado, até o retorno ao estado vazio. Note que a função densidade de probabilidade de \(B\), neste caso, é precisamente a expressão \(P_0'(t)\), obtida anteriormente. Assim, a partir da transformada de Laplace de \(B\), é possível determinar seus momentos utilizando as propriedades clássicas da transformada de Laplace, dadas por:

\[\begin{align}\\ E(B) &= -\gamma'(0)\\\\ E(B^2) &= \gamma''(0)\\\\ \end{align}\]

A aplicação direta dessas propriedades sobre \(\pi_0(\theta)\) conduz às seguintes expressões para a média e a variância do período de ocupação:

\[\begin{align}\\ E[B] &= \frac{1}{\mu - \lambda}\\\\ V[B] &= \frac{1 + \rho}{\mu^2 (1 - \rho)^3}\\\\ \end{align}\]

Dessa análise, conclui-se que, quando o período de ocupação se inicia com \(i\) clientes no sistema, sua duração pode ser formalmente interpretada como a soma de \(i\) períodos de ocupação independentes, cada um correspondente à trajetória do processo de \(k \to k-1\), para \(k = i, i-1, \dots, 1\). Essa decomposição decorre diretamente da propriedade markoviana do processo subjacente, uma vez que, após cada transição de \(k\) para \(k-1\), o processo reinicia com a mesma dinâmica, mas com um cliente a menos. Logo, como a duração de \(B_i\) corresponde à soma de \(i\) cópias independentes e identicamente distribuídas de \(B\), segue diretamente da linearidade da esperança e da propriedade aditiva da variância em somas de variáveis independentes que:

\[\begin{align}\\ E[B_i] &= i E[B] = \frac{i}{\mu - \lambda}\\\\ V[B_i] &= i V[B] = \frac{i (1 + \rho)}{\mu^2 (1 - \rho)^3}\\\\ \end{align}\]

Por fim, a análise do comportamento assintótico da transformada de Laplace \(\pi_0(\theta)\) quando \(\theta \to 0\) permite uma interpretação probabilística para a caracterização do período de ocupação. De fato, tem-se que:

\[\begin{align}\\ \lim_{\theta \to 0} \pi_0(\theta) &= \lim_{\theta \to 0} \frac{(\theta + \lambda + \mu) - \sqrt{(\theta + \lambda + \mu)^2 - 4 \lambda \mu}}{2 \lambda} \\\\ &= \frac{1}{2\lambda} \left[\lambda + \mu - |\mu - \lambda|\right] \\\\ &= \begin{cases} 1, & \mu \geqslant \lambda, \\\\ \dfrac{\lambda}{\mu}, & \mu < \lambda, \end{cases}\\\\ \end{align}\]

o que equivale a

\[\begin{align}\\ \lim_{\theta \to 0} \pi_0(\theta) = \begin{cases} 1, & \rho < 1, \\\\ \dfrac{\lambda}{\mu}, & \rho > 1. \end{cases}\\\\ \end{align}\]

Mas \(\lim_{\theta \to 0} \pi_0(\theta) = \int_0^\infty P_0'(t) dt\), onde \(P_0'(t)\) denota a densidade da distribuição do período de ocupação. Portanto, conclui-se que a distribuição do período de ocupação é própria se, e somente se, \(\rho \leqslant 1\). Para \(\rho > 1\), a distribuição é imprópria, o que significa que há uma probabilidade positiva — dada por \(1 - \rho^{-1}\) — do período de ocupação não terminar. Este resultado, em particular, expressa a condição de estabilidade do modelo \(M/M/1\): o sistema é estável e recorrente positivo quando \(\rho \leqslant 1\), retornando ao estado zero em tempo finito com probabilidade igual a um. Por outro lado, para \(\rho > 1\), o sistema torna-se instável, caracterizado pelo crescimento ilimitado do número de clientes e pela ausência de retorno ao estado vazio.


(VI) Processo de Saída


O processo de saída de uma fila \(M/M/1\) corresponde ao resultado combinado dos processos de chegada e serviço. Sempre que o servidor permanece continuamente ocupado, o processo de saída coincide com o próprio processo de serviço. No entanto, em momentos de ociosidade — isto é, quando o sistema se encontra vazio — ocorrem interrupções nas saídas. Apesar disso, quando o sistema está em regime estacionário (isto é, para \(\rho < 1\)), é possível caracterizar rigorosamente as propriedades do processo de saída sem a necessidade de fazer referência direta às dinâmicas individuais dos processos de chegada e de serviço. Para aprofundar esta ideia, considere \(t_1, t_2, \dots\) os instantes de saída dos clientes do sistema, e defina \(T_n = t_{n+1} - t_n\) como o intervalo entre saídas consecutivas. No regime estacionário, denote essa variável por \(T\). Seja ainda \(Q(x)\) o número de clientes no sistema \(x\) unidades de tempo após uma saída, e defina:

\[\begin{align}\\ F_n(x) = \mathrm{P}[Q(x) = n, T > x]\\\\ \end{align}\]

Na fila \(M/M/1\), a distribuição estacionária do número de clientes no sistema, denotada por \(\pi_n = (1 - \rho) \rho^n\), é válida não apenas para instantes arbitrários de tempo, mas também condicionada a pontos específicos da trajetória, como tempos de chegada ou de saída. Portanto, para qualquer \(x \geqslant 0\), tem-se que:

\[\begin{align}\\ \mathrm{P}[Q(x) = n] = (1 - \rho) \rho^n, \quad n = 0, 1, 2, \ldots, \quad (x \geq 0)\\\\ \end{align}\]

Por outro lado, a função

\[\begin{align}\\ F(x) = \mathrm{P}(T > x) = \sum_0^\infty F_n(x)\\\\ \end{align}\]

representa a probabilidade de que o intervalo entre duas saídas consecutivas seja maior que \(x\). Observe que, devido à propriedade markoviana do processo subjacente, o valor de \(T\) depende apenas do número de clientes no sistema imediatamente após uma saída, e não dos intervalos anteriores. Para determinar \(F(x)\), analisamos a evolução do sistema no intervalo infinitesimal \((x, x + \Delta x]\). A probabilidade \(F_0(x + \Delta x)\) corresponde à situação em que não há clientes no sistema após um tempo \(x + \Delta x\), e nenhuma saída ocorre até esse instante. Como não há clientes, a única dinâmica possível é a chegada de novos clientes, de modo que:

\[\begin{align}\\ F_0(x + \Delta x) = F_0(x)[1 - \lambda \Delta x] + o(\Delta x)\\\\ \end{align}\]

Por outro lado, para \(n \geqslant 1\), a dinâmica contempla tanto chegadas quanto saídas, e a evolução é descrita por:

\[\begin{align}\\ F_n(x + \Delta x) = F_n(x)[1 - \lambda \Delta x - \mu \Delta x] + F_{n-1}(x) \lambda \Delta x + o(\Delta x), \quad n = 1, 2, \ldots\\\\ \end{align}\]

Subtraindo \(F_n(x)\) de ambos os lados, dividindo por \(\Delta x\) e tomando o limite quando \(\Delta x \to 0\), obtemos o sistema de equações diferenciais:

\[\begin{align}\\ F_0'(x) &= -\lambda F_0(x)\\\\ F_n'(x) &= -(\lambda + \mu) F_n(x) + \lambda F_{n-1}(x), \quad n = 1, 2, \ldots\\\\ \end{align}\]

com as condições iniciais:

\[\begin{align}\\ F_n(0) = \mathrm{P}[Q(0) = n] = \pi_n\\\\ \end{align}\]

A equação diferencial para \(F_0(x)\), em particular, possui solução defina por:

\[\begin{align}\\ F_0(x) = \pi_0 e^{-\lambda x} = (1 - \rho)e^{-\lambda x}\\\\ \end{align}\]

A solução para \(F_n(x)\), para \(n \geqslant 1\), é obtida por indução. Supondo a validade da solução para \(F_{n-1}(x)\), a equação diferencial resulta em:

\[\begin{align}\\ F_n'(x) + (\lambda + \mu) F_n(x) = \lambda \pi_{n-1} e^{-\lambda x}\\\\ \end{align}\]

Multiplicando ambos os lados por \(e^{(\lambda + \mu) x}\), integrando, e utilizando a condição inicial, obtemos:

\[\begin{align}\\ F_n(x) = \pi_n e^{-\lambda x}, \quad n = 1, 2, 3, \ldots\\\\ \end{align}\]

Assim, a função \(F(x)\), que representa a sobrevivência do tempo entre saídas, é

\[\begin{align}\\ F(x) = \sum_{n=0}^\infty \pi_n e^{-\lambda x} = e^{-\lambda x}\\\\ \end{align}\]

uma vez que \(\sum_{n=0}^\infty \pi_n = 1\). Portanto, a distribuição do tempo entre saídas é exponencial com parâmetro \(\lambda\), idêntica à distribuição dos tempos entre chegadas. Este resultado confirma que o processo de saída, em equilíbrio, é um processo de Poisson com taxa \(\lambda\), exatamente igual ao processo de chegada. Esse é um resultado fundamental na teoria de filas, conhecido como o Teorema de Burke (Burke 1956), que estabelece que, para a fila \(M/M/1\) em equilíbrio, o processo de saída é um processo de Poisson com a mesma taxa de chegada, e que o número de clientes no sistema após uma saída tem distribuição estacionária \(\pi_n\), independentemente da história passada. Consequentemente, o número esperado de clientes atendidos em um intervalo de tempo \(t\), quando o sistema se encontra em equilíbrio, é simplesmente \(\lambda t\).

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 5.1 (Modelo \(M/M/1\) Aplicado a uma Pista de Pouso). Considere um aeroporto que dispõe de uma única pista para pouso de aeronaves. Estudos operacionais indicaram que os aviões chegam segundo um processo de Poisson com taxa média de \(\lambda = 15\) aviões por hora. O tempo necessário para a realização de cada pouso segue uma distribuição exponencial, com tempo médio de 3 minutos por pouso, o que corresponde a uma taxa de serviço de \(\mu = 20\) pousos por hora. Sob essas condições, o sistema pode ser adequadamente modelado como uma fila \(M/M/1\). A seguir, são calculadas as principais medidas de desempenho operacionais desse sistema.


(a) Utilização da Pista


A utilização da pista, ou fator de utilização, é a fração do tempo em que a pista permanece ocupada no regime estacionário. É dada por:

\[\begin{align}\\ \rho = \frac{\lambda}{\mu} = \frac{15}{20} = \frac{3}{4} = 0.75\\\\ \end{align}\]

Isso significa que a pista está, em média, ocupada 75% do tempo, permanecendo ociosa nos 25% restantes.


(b) Número Médio de Aeronaves Aguardando na Fila (Não Incluindo a que Está Pousando)


A quantidade média de aviões em espera, isto é, na fila, é calculada por:

\[\begin{align}\\ L_q = \frac{\rho^2}{1 - \rho} = \frac{(0.75)^2}{0.25} = 2.25\\\\ \end{align}\]

Em média, aproximadamente 2 aeronaves estão aguardando na fila para pousar.


(c) Tempo Médio de Espera na Fila (Antes do Início do Pouso)


O tempo médio que uma aeronave permanece aguardando na fila, antes do início do pouso, é dado, pela lei de Little, por:

\[\begin{align}\\ W_q = \frac{\lambda}{\mu (\mu - \lambda)} = \frac{15}{20 (20 - 15)} = \frac{3}{20} \text{ hora} = 9 \text{ minutos}\\\\ \end{align}\]

Portanto, uma aeronave espera, em média, 9 minutos na fila antes de iniciar o procedimento de pouso.


(d) Probabilidade de Determinadas Situações de Espera


A distribuição do tempo de espera \(T_q\) possui uma massa de probabilidade em zero (quando o avião não espera) e uma componente contínua para tempos positivos. As probabilidades de interesse são:


  • Probabilidade de não haver espera: \[\begin{align}\\ P(\text{sem espera}) = P(T_q = 0) = 1 - \rho = 0.25\\\\ \end{align}\] Ou seja, 25% das aeronaves iniciam o pouso imediatamente ao chegar.


  • A probabilidade de que o tempo de espera seja superior a \(t\) minutos é dada por:: \[\begin{align}\\ P(T_q > t) = \rho e^{-\mu (1-\rho) t}\\\\ \end{align}\] onde o tempo \(t\) é expresso em minutos. Então, para \(t = 5\) minutos, tem-se que: \[\begin{align}\\ P(T_q > 5 \text{ minutos}) = \frac{3}{4} e^{-20 (1 - \frac{3}{4}) \frac{5}{60}} = \frac{3}{4} e^{-\frac{25}{60}} = 0.4944\\\\ \end{align}\] e, para \(t = 10\) minutos, tem-se que: \[\begin{align}\\ P(T_q > 10 \text{ minutos}) = \frac{3}{4} e^{-\frac{50}{60}} = 0.3259\\\\ \end{align}\]


Portanto:


  • A probabilidade de uma aeronave esperar mais de 5 minutos é aproximadamente 49.44%.
  • A probabilidade de esperar mais de 10 minutos é aproximadamente 32.59%.


(e) Número Esperado de Pousos em um Intervalo de 20 Minutos


O processo de saída, por ser um processo de Poisson com taxa \(\lambda = 15\) pousos por hora, gera uma quantidade esperada de pousos proporcional à duração do intervalo. Para um intervalo de 20 minutos (\(1/3\) hora):

\[\begin{align}\\ E[\text{número de pousos}] = 15 \times \frac{20}{60} = 5\\\\ \end{align}\]

Portanto, espera-se que 5 aeronaves realizem pouso nesse intervalo.


O Código 1 apresenta a implementação, em linguagem R, dos cálculos analíticos e da visualização gráfica dos principais indicadores de desempenho do sistema de filas \(M/M/1\) aplicado a um aeroporto com uma única pista de pouso considerado neste exemplo. As chegadas de aeronaves seguem um processo de Poisson com taxa \(\lambda = 15\) por hora, e os tempos de pouso seguem uma distribuição exponencial com tempo médio de 3 minutos, correspondente a uma taxa de serviço \(\mu = 20\) por hora.


A Figura 1 ilustra a função de sobrevivência do tempo de espera na fila, isto é, a probabilidade de que o tempo de espera antes do pouso exceda um determinado valor; e, também, apresenta a distribuição da quantidade de pousos esperada em um intervalo de 20 minutos, que segue uma distribuição de Poisson com média igual a 5.


Código 1. Cálculo e visualização das medidas de desempenho do sistema \(M/M/1\) aplicado a uma pista de pouso de aeroporto, com \(\lambda = 15\) e \(\mu = 20\) (tempo médio de serviço igual a 3 minutos).

# ---------------------------------------------------------
# Análise do Modelo de Filas M/M/1 Aplicado a um Aeroporto
# ---------------------------------------------------------

# --- 0. Pacotes Necessários ---

library(dplyr)
library(ggplot2)
library(tidyr)
library(patchwork)

# --- 1. Configurações Iniciais ---

lambda            <- 15               # Taxa de Chegada (Aviões por Hora)
mu                <- 60 / 3           # Taxa de Serviço (Aviões por Hora), 3 Minutos por Pouso
t_int             <- 20/60            # Intervalo de 20 Minutos (Em Horas)
rho               <- lambda / mu      # Fator de Utilização

# --- 2. Cálculo das Métricas de Desempenho ---


L                 <- rho / (1 - rho)                      # Número Médio no Sistema
Lq                <- rho^2 / (1 - rho)                    # Número Médio na Fila
W                 <- 1 / (mu - lambda)                    # Tempo Médio no Sistema (horas)
Wq                <- rho / (mu * (1 - rho))               # Tempo Médio na Fila (horas)

P_zero            <- 1 - rho                              # Probabilidade de Não Esperar
P_wait_5          <- rho * exp(-mu * (1 - rho) * (5/60))  # Espera > 5 min
P_wait_10         <- rho * exp(-mu * (1 - rho) * (10/60)) # Espera > 10 min

exp_land_20min    <- lambda * t_int                       # Número Esperado de Pousos em 20 Minutos

# --- 3. Organização dos Resultados ---

resultados        <- tibble(
                        Métrica = c(
                          "Fator de Utilização (ρ)",
                          "Nº Médio no Sistema (L)",
                          "Nº Médio na Fila (Lq)",
                          "Tempo Médio no Sistema (W) [min]",
                          "Tempo Médio na Fila (Wq) [min]",
                          "P(Sem Espera)",
                          "P(Espera > 5 min)",
                          "P(Espera > 10 min)",
                          "Nº Esperado de Pousos em 20 min"
                        ),
                        Valor = c(
                          round(rho, 4),
                          round(L, 4),
                          round(Lq, 4),
                          round(W * 60, 4),   # Convertido para Minutos
                          round(Wq * 60, 4),  # Convertido para Minutos
                          round(P_zero, 4),
                          round(P_wait_5, 4),
                          round(P_wait_10, 4),
                          round(exp_land_20min, 4)
                        )
                      )

print(resultados)
# A tibble: 9 × 2
  Métrica                           Valor
  <chr>                             <dbl>
1 Fator de Utilização (ρ)           0.75 
2 Nº Médio no Sistema (L)           3    
3 Nº Médio na Fila (Lq)             2.25 
4 Tempo Médio no Sistema (W) [min] 12    
5 Tempo Médio na Fila (Wq) [min]    9    
6 P(Sem Espera)                     0.25 
7 P(Espera > 5 min)                 0.494
8 P(Espera > 10 min)                0.326
9 Nº Esperado de Pousos em 20 min   5    
# --- 4. Função de Sobrevivência do Tempo de Espera ---

tempo_min         <- seq(0, 30, by = 0.5)

df_surv           <- tibble(Tempo = tempo_min,
                            Sobrevivencia = rho * exp(-mu * (1 - rho) * Tempo / 60))

grafico1          <- ggplot(df_surv, aes(x = Tempo, y = Sobrevivencia)) +
                        geom_line(color = "darkblue", linewidth = 1) +
                        geom_point(color = "darkblue", size = 2) +
                        labs(
                          title = "",
                          x = "Tempo (minutos)",
                          y = expression(P(T[q]>t))
                        ) +
                        theme_minimal() +
                        theme(
                          axis.title.x = element_text(margin = margin(t = 0.5, unit = "lines")),
                          axis.title.y = element_text(margin = margin(r = 0.5, unit = "lines"))
                        )

# --- 5. Distribuição do Número de Pousos em 20 Minutos ---

x                 <- 0:15

df_poisson        <- tibble(x = x, Probabilidade = dpois(x, lambda * t_int))

grafico2          <- ggplot(df_poisson, aes(x = x, y = Probabilidade)) +
                        geom_col(fill = "darkred") +
                        geom_text(aes(label = round(Probabilidade, 2)), vjust = -0.5, size = 3) +
                        labs(
                          title = "",
                          x = "Nº de Pousos",
                          y = "Probabilidade"
                        ) +
                        theme_minimal() +
                        theme(
                          axis.title.x = element_text(margin = margin(t = 0.5, unit = "lines")),
                          axis.title.y = element_text(margin = margin(r = 0.5, unit = "lines"))
                        )

# --- 6. Visualização dos Gráficos ---

grafico1 + grafico2


Figura 5.1.. Ilustração da função de sobrevivência do tempo de espera na fila e da distribuição da quantidade de pousos esperada em um intervalo de 20 minutos.

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

5.5.2. Modelo M/M/s



O modelo \(M/M/s\) é uma generalização direta do modelo \(M/M/1\), adequado para a representação de sistemas com uma única fila de espera e \(s\) servidores paralelos (\(s \in \mathbb{N}\)), onde as chegadas ocorrem segundo um processo de Poisson e os tempos de serviço seguem uma distribuição exponencial. As principais suposições do modelo \(M/M/s\) são:


  • Processo de Chegada: As chegadas ocorrem segundo um processo de Poisson com taxa \(\lambda\), o que implica que o número de chegadas no intervalo \((0, t]\) é uma variável aleatória \(N(t)\) com distribuição de Poisson de parâmetro \(\lambda t\), e os tempos entre chegadas consecutivas são independentes e identicamente distribuídos segundo uma distribuição exponencial com parâmetro \(\lambda\).

  • Processo de Serviço: Os tempos de serviço são independentes e distribuídos exponencialmente com taxa \(\mu\). O sistema possui \(s\) servidores operando em paralelo. Quando há \(n \leqslant s\) clientes no sistema, todos são atendidos simultaneamente; quando \(n > s\), \(s\) estão em serviço e \((n-s)\) aguardam na fila.

  • Disciplina de Atendimento: FIFO (First-In, First-Out).

  • Capacidade da Fila: Ilimitada, sem perdas.

  • População de Clientes: Infinita.


O parâmetro que governa o comportamento do sistema é a intensidade de tráfego por servidor (\(\rho\)) (ou fator de utilização), definido por:

\[\begin{align}\\ \rho = \frac{\lambda}{s\mu}\\\\ \end{align}\]

Para que o sistema atinja o equilíbrio — isto é, não ocorra crescimento indefinido do número de clientes —, é necessário que \(\rho < 1\), ou equivalentemente, que \(\lambda < s\mu\).


(I) Equações de Kolmogorov


O processo estocástico \(\{Q(t)\}_{t \geqslant 0}\), que representa o número de clientes no sistema em um instante \(t\), que descreve o modelo M/M/s é também um caso especial do processo de nascimento e morte, com taxas definidas por:

\[\begin{align}\\ \lambda_n &= \lambda\\\\ \mu_n &= \begin{cases} n\mu \,\, & n < s\\ s\mu \,\, & n \geqslant s \end{cases}\\\\ \end{align}\]

e espaço de estados \(\mathcal{E} = \{0, 1, 2, \ldots \}\). Para este processo, as equações de Kolmogorov prospectivas, que descrevem a evolução das probabilidades de transição \(P_n(t) = \mathrm{P}[Q(t) = n]\), são:

\[\begin{align}\\ P_0'(t) &= -\lambda P_0(t) + \mu P_1(t)\\\\ P_n'(t) &= -(\lambda + n\mu) P_n(t) + \lambda P_{n-1}(t) + (n+1)\mu P_{n+1}(t), \quad 1\leqslant n < s\\\\ P_n'(t) &= -(\lambda + s\mu) P_n(t) + \lambda P_{n-1}(t) + s\mu P_{n+1}(t), \quad n \geqslant s\\\\ \end{align}\]

com \(P_n(0) = 1\) quando \(n = i\) e \(P_n = 0\) caso contrário. Para determinar uma solução dessas equações diferenciais, também é necessário o uso de transformadas de Laplace, ou de métodos computacionais, seja pela integração numérica direta das equações diferenciais, seja por meio de simulação estocástica do processo \(\{Q(t)\}_{t \ \geqslant \ 0}\).


(II) Distribuição Estacionária


A distribuição estacionária \(\pi_n = \lim_{t \to \infty} P_n(t)\) é obtida a partir das equações de equilíbrio, que, por balanço de fluxos, são:

\[\begin{align}\\ \lambda \pi_0 &= \mu \pi_1, \qquad\qquad\qquad\qquad\quad\,\,\,\,\, n = 0\\\\ (\lambda + n\mu) \pi_n &= \lambda \pi_{n-1} + (n+1)\mu \pi_{n+1}, \qquad 1\leqslant n < s\\\\ (\lambda + s\mu) \pi_n &= \lambda \pi_{n-1} + s\mu \pi_{n+1}, \qquad\qquad\,\,\,\,\, n\geqslant s\\\\ \end{align}\]

Utilizando um procedimento recursivo, como no caso M/M/1, obtém-se que:

\[\begin{align}\\ \pi_n &= \dfrac{1}{n!} \left(\dfrac{\lambda}{\mu}\right)^n \pi_0, \qquad\qquad 0 \leqslant n \leqslant s\\\\ \pi_{n+r} &= \left(\dfrac{\lambda}{s\mu}\right)^r \pi_s, \qquad\qquad\quad r = 0, 1, \ldots\\\\ \pi_{n} &= \left(\dfrac{\lambda}{s\mu}\right)^{n-s} \pi_s, \qquad\qquad n = s, s+1, \ldots\\\\ \end{align}\]

Sendo \(\rho = \lambda/s\mu\), tem-se que:

\[\begin{align}\\ \pi_n &= \dfrac{1}{n!} (s\rho)^n \pi_0, \qquad\quad\,\,\,\, 0 \leqslant n \leqslant s\\\\ \pi_n &= \dfrac{1}{s!} (s\rho)^s \rho^{n-s} \pi_0, \qquad s \leqslant n < + \infty\\\\ \end{align}\]

Para determinar a constante \(\pi_0\), utiliza-se a condição de normalização \(\sum_{n=0}^{\infty} \pi_n = 1\). Neste caso, tem-se que:

\[\begin{align}\\ \pi_0 &= \left[\sum_{r=0}^{s-1} \dfrac{(s\rho)^r}{r!}+\dfrac{(s\rho)^s}{s!(1-\rho)}\right]^{-1}\\\\ \end{align}\]

Assim,

\[\begin{align}\\ \pi_n &= \rho^{n-s} \pi_s \qquad n \geqslant s\\\\ \end{align}\]

isto é, quando o número de clientes no sistema é maior ou igual a \(s\), o sistema se comporta como um modelo \(M/M/1\) com taxa de serviço \(s\mu\). Por conveniência, pode-se também escrever \(\alpha = \lambda/\mu\) tal que \(\alpha/s = \rho\). Neste caso,

\[\begin{align}\\ \pi_0 &= \left[\sum_{r=0}^{s-1} \dfrac{\alpha^r}{r!}+\dfrac{\alpha^s}{s!}\left(1 - \dfrac{\alpha}{s}\right)^{-1}\right]^{-1}\\\\ \pi_n &= \dfrac{\alpha^n}{n!} \pi_0, \qquad\qquad\quad\,\,\, 0 \leqslant n \leqslant s\\\\ \pi_n &= \dfrac{\alpha^s}{s!} \left(\dfrac{\alpha}{s}\right)^{n-s} \pi_0, \qquad s \leqslant n < + \infty\\\\ \end{align}\]

Note que os clientes precisarão aguardar pelo serviço apenas se o número de clientes no sistema for maior ou igual a \(s\). A probabilidade desse evento no modelo \(M/M/s\) é chamada de probabilidade de atraso, e é definida por:

\[\begin{align}\\ \text{P(atraso)} = C(s,\alpha) &= \sum_{n = s}^\infty \pi_n = \dfrac{\alpha^s}{s!}\left[1 - \dfrac{\alpha}{s}\right]^{-1} \pi_0 \\\\ \end{align}\]

A fórmula para \(C(s, \alpha)\) é conhecida na literatura como fórmula de atraso de Erlang (Erlang 1917) ou segunda fórmula de Erlang, sendo também denotada por \(E_{2,s}(\alpha)\). Antes do surgimento dos computadores, a indústria telefônica utilizava tabelas e gráficos de \(C(s, \alpha)\), construídos para diferentes combinações de \(s\) e \(\alpha\).


(III) Número Médio de Clientes


Seja \(Q = Q(\infty)\) o número de clientes no sistema no regime estacionário, e \(Q_q\) o número de clientes na fila, excluindo aqueles em serviço. Define-se:


  • \(L = E[Q]\), que representa o número médio de clientes no sistema.
  • \(L_q = E[Q_q]\), que representa o número médio de clientes na fila.


A partir das equações da distribuição estacionária \(\pi_n\) e escrevendo, por conveniência, \(\alpha = s\rho\), tem-se que:

\[\begin{align}\\ \sum_{n=1}^\infty n\pi_n &= \pi_0\left[\sum_{n=1}^s n \dfrac{\alpha^n}{n!} + \sum_{n = s + 1}^{\infty} n\rho^{n-s} \dfrac{\alpha^s}{s!}\right]\\\\ &= \pi_0\left[\alpha \sum_{r=0}^{s-1} \dfrac{\alpha^r}{r!} + \dfrac{\alpha^s}{s!}\left\{\dfrac{\rho}{(1-\rho)^2} + \dfrac{s\rho}{1-\rho}\right\}\right]\\\\ &= \dfrac{\rho\alpha^s\pi_0}{s!(1-\rho)^2} + \alpha \pi_0 \left[\sum_{r=0}^{s-1} \dfrac{\alpha^r}{r!} + \dfrac{\alpha^s}{s!\rho(1-\rho)}\right]\\\\ &= \dfrac{\rho\alpha^s\pi_0}{s!(1-\rho)^2} + \alpha \pi_0 \pi_{0}^{-1}\\\\ &= \alpha + \dfrac{\rho\alpha^s\pi_0}{s!(1-\rho)^2}\\\\ \end{align}\]

Isso implica que,

\[\begin{align}\\ L = \alpha + \dfrac{\rho\alpha^s\pi_0}{s!(1-\rho)^2} = \alpha +\dfrac{\rho\pi_s}{(1-\rho)^2}\\\\ \end{align}\]

De forma análoga, o número médio de clientes na fila, \(L_q\), é obtido como:

\[\begin{align}\\ L_q &= \sum_{n=s + 1}^\infty (n-s) \dfrac{\alpha^s}{s!}\rho^{n-s}\pi_0\\\\ &= \dfrac{\alpha^s}{s!}\pi_0\sum_{n-s = 1}^\infty (n-s)\rho^{n-s}\\\\ &= \dfrac{\alpha^s}{s!}\pi_0\sum_{r = 1}^\infty r\rho^{r}\\\\ &= \frac{\rho\alpha^s\pi_0}{s!(1-\rho)^2}\\\\ &= \dfrac{\rho\pi_s}{(1-\rho)^2}\\\\ \end{align}\]

Comparando as expressões para \(L\) e \(L_q\), pode-se deduzir que \(s\rho\) (isto é, \(\alpha\)) representa o número esperado de servidores ocupados. Este valor também pode ser interpretado como a contribuição do fator de utilização em um sistema com \(s\) servidores. Por exemplo, a utilização individual de um servidor pode ser expressa como:

\[\begin{align}\\ \sum_{n=1}^{s-1} \dfrac{n}{s} \pi_n + \sum_{n=s}^\infty \pi_n\\\\ \end{align}\]

Ao substituir as expressões para \(\pi_n\) na equação acima e realizar as simplificações, obtém-se que o fator de utilização individual de cada servidor é igual a \(\rho\). Ou seja, no longo prazo, a probabilidade (ou a fração de tempo) de que um servidor esteja ocupado é \(\rho\).


(IV) Tempo de Espera do Cliente


Para a análise dos tempos de espera dos clientes, note que, quando o número de clientes no sistema é maior ou igual a \(s\), os tempos entre saídas consecutivas são variáveis exponenciais com parâmetro \(s\mu\). Então, seja \(T_q\) o tempo de espera na fila de um cliente no regime estacionário (isto é, quando \(t \to \infty\)), e denote sua função distribuição acumulada por \(F_q(t) = \mathrm{P}[T_q \leqslant t]\). É evidente que:

\[\begin{align}\\ F_q(0) &= \mathrm{P}[T_q = 0] = \mathrm{P}(Q < s) = \pi_0 \sum_{n = 0}^{s-1} \frac{\alpha^n}{n!}\\\\ \end{align}\]

onde \(\alpha = \lambda / \mu\). Da primeira equação de equilíbrio, tem-se que:

\[\begin{align}\\ \sum_{n = 0}^{s-1} \frac{\alpha^n}{n!} = \frac{1}{\pi_0} - \frac{\alpha^s}{s!(1 - \rho)}\\\\ \end{align}\]

o que nos permite escrever:

\[\begin{align}\\ F_q(0) = 1 - \frac{\alpha^s \pi_0}{s!(1 - \rho)}\\\\ \end{align}\]

Além disso, seguindo um raciocínio análogo ao adotado para o modelo \(M/M/1\), obtém-se, no caso multisservidor para \(n \geqslant s\), que:

\[\begin{align}\\ dF_q(t) &= \sum_{n = s}^{\infty} \pi_n e^{-s\mu t} \frac{(s\mu t)^{n-s}}{(n-s)!} s\mu \, dt \\\\ &= \pi_s e^{-s\mu t} \sum_{n = s}^{\infty} \frac{\{\rho(s\mu t)\}^{n-s}}{(n-s)!} s\mu \, dt\\\\ \end{align}\]

Usando a expansão em série de Taylor da função exponencial, a expressão acima se simplifica para:

\[\begin{align}\\ dF_q(t) = s\mu \pi_s e^{-s\mu(1 - \rho)t} \, dt\\\\ \end{align}\]

Alternativamente, expressando \(\pi_s\) em função de \(\pi_0\), pode-se reescrever:

\[\begin{align}\\ dF_q(t) = \frac{s\mu \alpha^s}{s!} \pi_0 e^{-s\mu(1 - \rho)t} \, dt\\\\ \end{align}\]

Observe que o termo \(F_q(0)\) não contribui para o valor esperado de \(T_q\). Então, o tempo médio de espera na fila é calculado como:

\[\begin{align}\\ W_q &= \int_0^{\infty} t \, dF_q(t) \\\\ &= \int_0^{\infty} s\mu \pi_s t e^{-s\mu(1 - \rho)t} \, dt\\\\ &= \frac{\pi_s}{s\mu(1 - \rho)^2}\\\\ \end{align}\]

Expressando, novamente, \(\pi_s\) em função de \(\pi_0\), tem-se que:

\[\begin{align}\\ W_q = \frac{\alpha^s \pi_0}{s! \, s\mu (1 - \rho)^2}\\\\ \end{align}\]

A função distribuição acumulada \(F_q(t)\) do tempo de espera na fila, considerando os resultados anteriores, é dada por:

\[\begin{align}\\ F_q(t) &= F_q(0) + \int_0^t \frac{s\mu \alpha^s}{s!} \pi_0 e^{-s\mu(1 - \rho)x} \, dx \\\\ &= 1 - \dfrac{\alpha^s\pi_0}{s!(1-\rho)} + \dfrac{\alpha^s\pi_0}{s!(1-\rho)}\int_{0}^{t}s\mu(1-\rho)e^{-s\mu(1-\rho)x} \, dx\\\\ &= 1 - \frac{\alpha^s \pi_0}{s!(1 - \rho)} e^{-s\mu(1 - \rho)t}\\\\ \end{align}\]

(V) Período de Ocupação


O período de ocupação no modelo \(M/M/s\) é formalmente definido como o intervalo de tempo durante o qual todos os \(s\) servidores permanecem simultaneamente ocupados, iniciando no instante em que o sistema transita do estado \(n = s-1\) para o estado \(n = s\), e encerrando quando o número de clientes no sistema retorna para \(n \leqslant s-1\). Assim, trata-se do tempo até a primeira ocorrência de uma saída que reduza o número de clientes no sistema de \(s\) para \(s-1\), condicionado ao fato de que o processo foi iniciado no estado \(n = s\). A formulação completa da distribuição deste período de ocupação no modelo \(M/M/s\) pode ser desenvolvida a partir das equações de Kolmogorov associadas a um processo de Markov em tempo contínuo, no qual o estado \(n = s-1\) é tornado absorvente, enquanto os estados \(n \geqslant s\) formam uma classe transitória. Essa abordagem é conceitualmente análoga àquela empregada para o modelo \(M/M/1\), porém apresenta maior complexidade matemática, uma vez que, no modelo \(M/M/s\), a taxa de serviço total é variável — igual a \(s\mu\) enquanto \(n \geqslant s\), e reduz-se proporcionalmente à medida que o número de clientes no sistema diminui para valores inferiores a \(s\).


Por outro lado, se definirmos o período de ocupação como o intervalo de tempo durante o qual ao menos um dos servidores está ocupado, isto é, o sistema não está vazio, torna-se necessária uma formulação substancialmente distinta. A análise desse conceito exige o desenvolvimento de modelos que considerem a evolução do sistema a partir da transição do estado \(n = 0\) para qualquer estado \(n \geqslant 1\), encerrando-se no retorno ao estado vazio. A obtenção da distribuição e dos momentos associados a esse tipo de período de ocupação requer resultados teóricos adicionais, envolvendo uma formulação mais elaborada das equações diferenciais associadas, e, por essa razão, ultrapassa o escopo deste texto.


(VI) Processo de Saída


De acordo com o teorema de Burke (Burke 1956), o processo de saída de uma fila \(M/M/s\), quando o sistema opera em regime estacionário (\(\rho < 1\)), é um processo de Poisson com taxa \(\lambda\), rigorosamente igual ao processo de chegada. Isso significa que, no equilíbrio, os intervalos de tempo entre partidas consecutivas são independentes e seguem uma distribuição exponencial com parâmetro \(\lambda\), independentemente do número de clientes no sistema ou da quantidade de servidores ocupados no instante considerado. Consequentemente, o número esperado de clientes que concluem o atendimento em um intervalo de duração \(t\) é simplesmente \(\lambda t\). Este resultado decorre diretamente da estrutura Markoviana do modelo, da falta de memória dos tempos exponenciais de chegada e de serviço e da propriedade de reversibilidade associada a processos de nascimento e morte balanceados em regime estacionário.


Essa propriedade não é exclusiva do modelo \(M/M/1\). Conforme discutido na análise do processo de saída do modelo \(M/M/1\), o procedimento ali desenvolvido se estende de forma natural ao modelo \(M/M/s\). De fato, as equações diferenciais que descrevem o processo de saídas podem ser generalizadas para acomodar taxas de serviço que variam de acordo com o número de clientes no sistema. De forma mais formal, seja \(F_n(x)\) a função de distribuição do tempo até a próxima partida, condicionado ao fato de que o sistema se encontra no estado \(n\) no instante inicial. Dado que o sistema está em equilíbrio, a probabilidade de o sistema estar no estado \(n\) é \(\pi_n\), e a probabilidade de nenhuma chegada ocorrer até o tempo \(x\) é \(e^{-\lambda x}\), pois as chegadas seguem um processo de Poisson com taxa \(\lambda\). Assim, a função de distribuição condicional é dada por:

\[\begin{align}\\ F_n(x) = P(T\leqslant x \mid N = n) = \pi_n e^{-\lambda x}, \quad n = 0, 1, 2, 3, \ldots\\\\ \end{align}\]

A função de distribuição acumulada marginal do tempo até a próxima saída, , denotada por \(F(x)\), é obtida como a soma das distribuições condicionais ponderadas pelas respectivas probabilidades estacionárias, isto é,

\[\begin{align}\\ F(x) = \sum_{n=0}^\infty F_n(x) = e^{-\lambda x} \sum_{n=0}^\infty \pi_n = e^{-\lambda x}\\\\ \end{align}\]

uma vez que \(\sum_{n=0}^\infty \pi_n = 1\). Portanto, o tempo até a próxima saída no equilíbrio é uma variável aleatória exponencial com taxa \(\lambda\), o que confirma que o processo de saída é um processo de Poisson de mesma taxa do processo de chegada. Este é um resultado mostra que, apesar da dinâmica interna complexa associada à presença de múltiplos servidores e da variabilidade no número de clientes no sistema, o fluxo de saída preserva a simplicidade de um processo de Poisson com taxa \(\lambda\).

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 5.2 (Modelo \(M/M/s\) Aplicado a Duas Pistas de Pouso). Considere novamente o aeroporto do Exemplo 5.1, mas agora equipado com duas pistas de pouso operando em paralelo, mantendo-se as mesmas taxas de chegada e de serviço: as aeronaves chegam segundo um processo de Poisson com taxa média de \(\lambda = 15\) aviões por hora e o tempo de pouso de cada aeronave segue uma distribuição exponencial com média de 3 minutos, isto é, taxa de serviço \(\mu = 20\) por hora por pista. Sob essas condições, o sistema é modelado por uma fila \(M/M/s\) com \(s = 2\). A seguir, são calculadas as principais medidas de desempenho operacionais desse sistema.


(a) Utilização de Cada Pista


O fator de utilização de cada pista, isto é, a fração do tempo em que uma pista está ocupada, é dado por:

\[\begin{align}\\ \rho = \frac{\lambda}{s\mu} = \frac{15}{2 \times 20} = \frac{3}{8} = 0.375\\\\ \end{align}\]

Isso significa que Cada pista está ocupada, em média, 37.5% do tempo, permanecendo ociosa nos 62.5% restantes.


(b) Número Médio de Aeronaves Aguardando na Fila (Não Incluindo as que Estão Pousando)


O cálculo do número médio de aeronaves aguardando na fila para o modelo \(M/M/s\) é dado por:

\[\begin{align}\\ L_q = \frac{\rho \alpha^s \pi_0}{s! (1 - \rho)^2}\\\\ \end{align}\]

onde \(\alpha = s\rho = 2 \times 0.375 = 0.75\). Para determiná-lo, deve-se, primeiramente, calcular \(\pi_0\) (probabilidade de que não haja nenhuma aeronave no sistema), isto é,

\[\begin{align}\\ \pi_0 &= \left[\sum_{r=0}^{1} \dfrac{\alpha^r}{r!}+\dfrac{\alpha^s}{s!(1 - \rho)}\right]^{-1}\\\\ &= \left[1 + \frac{3}{4} + \frac{(3/4)^2}{2} \cdot \left(1 - \frac{3}{8} \right)^{-1}\right]^{-1}\\\\ &= \left(1 + 0.75 + 0.45\right)^{-1} \\\\ &= 0.4545 \end{align}\]

Agora, calculando \(L_q\), tem-se que:

\[\begin{align}\\ L_q = \frac{0.375 \times (0.75)^2 \times 0.4545}{2 \times (5/8)^2} = 0.1227\\\\ \end{align}\]

Em média, aproximadamente 0.12 aeronaves estão aguardando na fila para pousar, ou seja, o tempo de espera é praticamente inexistente na maioria das situações.


(c) Tempo Médio de Espera na Fila (Antes do Início do Pouso)


O tempo médio que uma aeronave permanece aguardando na fila, antes do início do pouso, é dado, pela lei de Little, por:

\[\begin{align}\\ W_q = \frac{L_q}{\lambda} = \frac{0.1227}{15} = 0.00818 \text{ hora} = 0.49 \text{ minutos}\\\\ \end{align}\]

Portanto, uma aeronave espera, em média, 0.49 minutos (ou, aproximadamente, menos de 30 segundos) na fila antes de iniciar o procedimento de pouso.


(d) Probabilidades Associadas ao Tempo de Espera


A distribuição do tempo de espera \(T_q\) possui uma massa de probabilidade em zero (quando o avião não espera) e uma componente contínua para tempos positivos. As probabilidades de interesse são:


  • Probabilidade de não haver espera: \[\begin{align}\\ P(\text{sem espera}) = F_q(0) = 1 - \dfrac{\alpha^s \pi_0}{s! (1-\rho)} = 1 - \dfrac{(0.75)^2 \times 0.4545}{2 \times (5/8)} = 0.7955\\\\ \end{align}\] Ou seja, em 79.5% dos casos, a aeronave inicia o pouso imediatamente ao chegar.


  • A probabilidade de que o tempo de espera seja superior a \(t\) minutos é dada por:

\[\begin{align}\\ P(T_q > t) = \dfrac{\alpha^s \pi_0}{s! (1-\rho)} e^{-s\mu(1-\rho)t}\\\\ \end{align}\] onde o tempo \(t\) é expresso em minutos. Então, para \(t = 5\) minutos, tem-se que: \[\begin{align}\\ P(T_q > 5 \text{ minutos}) = 0.2045 e^{-2.0833} = 0.0255\\\\ \end{align}\] e, para \(t = 10\) minutos, tem-se que: \[\begin{align}\\ P(T_q > 10 \text{ minutos}) = 0.2045 e^{-4.16667} = 0.0032\\\\ \end{align}\]


Portanto:


  • A probabilidade de esperar mais de 5 minutos é de aproximadamente 2.55%.
  • A probabilidade de esperar mais de 10 minutos é de aproximadamente 0.32%.


(e) Número Esperado de Pousos em um Intervalo de 20 Minutos


Como o processo de saída é um processo de Poisson com taxa \(\lambda = 15\) pousos por hora, o número esperado de pousos em um intervalo de 20 minutos é dada por:

\[\begin{align}\\ E[\text{número de pousos}] = 15 \times \frac{20}{60} = 5\\\\ \end{align}\]

Portanto, espera-se que 5 aeronaves realizem pouso nesse intervalo, o que é idêntico ao caso de uma única pista, pois o processo de saída depende apenas da taxa de chegada \(\lambda\) em equilíbrio, conforme o Teorema de Burke.


O Código 2 apresenta a implementação, em linguagem R, dos cálculos analíticos e da visualização gráfica dos principais indicadores de desempenho do sistema de filas \(M/M/s\) aplicado ao aeroporto com duas pistas de pouso.


A Figura 2 ilustra a função de sobrevivência do tempo de espera na fila, bem como a distribuição da quantidade de pousos esperada em um intervalo de 20 minutos, que segue uma distribuição de Poisson com média igual a 5.


Código 2. Cálculo e visualização das medidas de desempenho do sistema \(M/M/s\), com \(s = 2\), aplicado a duas pistas de pouso, com \(\lambda = 15\) e \(\mu = 20\) (tempo médio de serviço igual a 3 minutos).

# ---------------------------------------------------------------------
# Análise do Modelo de Filas M/M/s (com s = 2) Aplicado a um Aeroporto
# ---------------------------------------------------------------------

# --- 0. Pacotes Necessários ---

library(dplyr)
library(ggplot2)
library(tidyr)
library(patchwork)

# --- 1. Configurações Iniciais ---

lambda            <- 15                 # Taxa de Chegada (Aviões por Hora)
mu                <- 20                 # Taxa de Serviço por Pista (Aviões por Hora)
s                 <- 2                  # Número de Pistas (Servidores)
rho               <- lambda / (s * mu)  # Fator de Utilização por Servidor
alpha             <- s * rho            # Parâmetro Auxiliar

# --- 2. Cálculo de π_0 (Probabilidade do Sistema Vazio) ---

sum1              <- sum(alpha^(0:(s-1)) / factorial(0:(s-1)))
sum2              <- (alpha^s) / (factorial(s) * (1 - rho))
pi_0              <- 1 / (sum1 + sum2)

# --- 3. Cálculo das Métricas de Desempenho ---

# Número Médio na Fila

Lq                <- (rho * alpha^s * pi_0) / (factorial(s) * (1 - rho)^2)

# Tempo Médio de Espera na Fila (Em Horas)

Wq                <- Lq / lambda                                            

# Probabilidade de Não Esperar

Fq0               <- 1 - (alpha^s * pi_0) / (factorial(s) * (1 - rho))      

# Probabilidades de Esperar Mais Que t Minutos

t1                <- 5   # 5 Minutos
t2                <- 10  # 10 Minutos
expo_rate         <- s * mu * (1 - rho)  # Taxa do Expoente

P_wait_5          <- (alpha^s * pi_0) / (factorial(s) * (1 - rho)) * exp(-expo_rate * t1 / 60)
P_wait_10         <- (alpha^s * pi_0) / (factorial(s) * (1 - rho)) * exp(-expo_rate * t2 / 60)

# Número Esperado de Pousos em 20 Minutos

t_int             <- 20/60
exp_land          <- lambda * t_int

# --- 4. Organização dos Resultados ---

resultados <- tibble(
                  Métrica = c(
                    "Fator de Utilização (ρ)",
                    "Nº Médio na Fila (Lq)",
                    "Tempo Médio na Fila (Wq) [min]",
                    "P(Sem Espera)",
                    "P(Espera > 5 min)",
                    "P(Espera > 10 min)",
                    "Nº Esperado de Pousos em 20 min"
                  ),
                  Valor = c(
                    round(rho, 4),
                    round(Lq, 4),
                    round(Wq * 60, 4),
                    round(Fq0, 4),
                    round(P_wait_5, 4),
                    round(P_wait_10, 4),
                    round(exp_land, 4)
                  )
                )

print(resultados)
# A tibble: 7 × 2
  Métrica                          Valor
  <chr>                            <dbl>
1 Fator de Utilização (ρ)         0.375 
2 Nº Médio na Fila (Lq)           0.123 
3 Tempo Médio na Fila (Wq) [min]  0.491 
4 P(Sem Espera)                   0.796 
5 P(Espera > 5 min)               0.0255
6 P(Espera > 10 min)              0.0032
7 Nº Esperado de Pousos em 20 min 5     
# --- 5. Função de Sobrevivência do Tempo de Espera ---

tempo_min         <- seq(0, 30, by = 0.5)
surv_rate         <- (alpha^s * pi_0) / (factorial(s) * (1 - rho))

df_surv           <- tibble(Tempo = tempo_min,
                            Sobrevivencia = surv_rate * exp(-expo_rate * Tempo / 60))


grafico1          <- ggplot(df_surv, aes(x = Tempo, y = Sobrevivencia)) +
                        geom_line(color = "darkblue", linewidth = 1) +
                        geom_point(color = "darkblue", size = 2) +
                        labs(
                          title = "",
                          x = "Tempo (minutos)",
                          y = expression(P(T[q]>t))
                        ) +
                        theme_minimal() +
                        theme(
                          axis.title.x = element_text(margin = margin(t = 0.5, unit = "lines")),
                          axis.title.y = element_text(margin = margin(r = 0.5, unit = "lines"))
                        )

# --- 6. Distribuição do Número de Pousos em 20 Minutos ---

x                 <- 0:15

df_poisson        <- tibble(x = x, Probabilidade = dpois(x, exp_land))

grafico2          <- ggplot(df_poisson, aes(x = x, y = Probabilidade)) +
                        geom_col(fill = "darkred") +
                        geom_text(aes(label = round(Probabilidade, 2)), vjust = -0.5, size = 3) +
                        labs(
                          title = "",
                          x = "Nº de Pousos",
                          y = "Probabilidade"
                        ) +
                        theme_minimal() +
                        theme(
                          axis.title.x = element_text(margin = margin(t = 0.5, unit = "lines")),
                          axis.title.y = element_text(margin = margin(r = 0.5, unit = "lines"))
                        )

# --- 7. Visualização dos Gráficos ---

grafico1 + grafico2


Figura 5.2.. Ilustração da função de sobrevivência do tempo de espera na fila e da distribuição da quantidade de pousos esperada em um intervalo de 20 minutos para o modelo M/M/s, com s = 2.

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

5.5.3. Modelo M/G/1



Considere um sistema de filas no qual os clientes chegam segundo um processo de Poisson com taxa \(\lambda\) e são atendidos por um único servidor. Suponha que o tempo de serviço desses clientes seja uma sequência de variáveis aleatórias i.i.d., \(\{S_n, n = 1, 2, 3, \dots\}\), com função de distribuição acumulada \(F_S(x) = P(S_n \leqslant x) = B(x)\) para \(x \geqslant 0\), valor esperado \(\mathrm{E}(S_n) = b\) e variância \(\mathrm{Var}(S_n) = \sigma_S^2\). Assume-se aqui que \(S_n\) é o tempo de serviço do \(n\)-ésimo cliente.


Defina \(Q(t)\) como o número de clientes no sistema no instante \(t\) com \(t_0 = 0, t_1, t_2, \dots\) sendo os instantes de partida dos clientes. Por construção, nesses instantes, os tempos restantes de serviço dos clientes são zero. Seja \(Q_n = Q(t_n + 0)\) o número de clientes no sistema logo após a \(n\)-ésima partida. Este processo, em particular, é uma cadeia de Markov. De fato, seja \(X_n\) o número de clientes que chegam durante o tempo de serviço \(S_n\). Com a suposição de Poisson para o processo de chegadas, tem-se que:

\[\begin{align}\\ k_j &= \mathrm{P}(X_n = j) \\\\ &=\int_{0}^{\infty}P(X_n = j\mid S_n) P(t < S_n \leqslant t + dt)\\\\ &= \int_0^\infty e^{-\lambda t}\frac{(\lambda t)^j}{j!} \, dB(t), \quad j = 0, 1, 2, \dots\\\\ \end{align}\]

Considere, agora, o relacionamento entre \(Q_n\) e \(Q_{n+1}\). Com base em \(X_n\), tem-se que:

\[\begin{align}\\ Q_{n+1} = \begin{cases} Q_n + X_{n+1} - 1, & \text{se } Q_n > 0, \\ X_{n+1}, & \text{se } Q_n = 0. \end{cases}\\\\ \end{align}\]

Como pode ser visto na equação acima, \(Q_{n+1}\) pode ser expresso em termos de \(Q_n\) e de uma variável aleatória \(X_{n+1}\), que não depende de nenhum evento anterior a \(t_n\). Como \(X_{n+1}\) é i.i.d., ele também não depende de \(Q_n\). Portanto, a propriedade de dependência de um passo das cadeias de Markov é satisfeita. Assim, \(\{Q_n, n = 0, 1, 2, \dots\}\) é uma cadeia de Markov. Seu espaço de parâmetros é formado pelos pontos de partida, e o espaço de estados \(S\) é dado pelo número de clientes no sistema, isto é, \(\mathcal{S} = \{0, 1, 2, \dots\}\). Devido à natureza do espaço de parâmetros, \(\{Q_n, n = 0, 1, 2, \dots\}\) é conhecida como uma cadeia de Markov embutida nos instantes de saída.


Seja, agora, \(P_{ij}^{(n)} = \mathrm{P}(Q_{n} = j \mid Q_0 = i)\) para \(i,j\in S\) tal que \(P_{ij}^{(1)} \equiv P_{ij}\). A partir das relações anteriores e da definição de \(k_j\), tem-se que:

\[\begin{align}\\ P_{ij} = P(Q_{n+1} = j \mid Q_n = i) = \begin{cases} k_{j-i+1}, & \text{se } i > 0, \\ k_j, & \text{se } i = 0. \end{cases}\\\\ \end{align}\]

A matriz de transição \(\mathbf{P}\), neste caso, é dada por:

\[\begin{align}\\ \mathbf{P} = \begin{bmatrix} 0 & 1 & 2 & \ldots\\ \hline k_0 & k_1 & k_2 & \dots \\ k_0 & k_1 & k_2 & \dots \\ & k_0 & k_1 & \dots \\ & & k_0 & \dots \\ & & & \ddots \end{bmatrix}\begin{matrix} \\ 0 \\ 1 \\ 2 \\ 3 \\ \vdots \end{matrix}\\\\\\ \end{align}\]

Observação 5.1. Para que a cadeia de Markov seja irredutível (ou seja, o espaço de estados tenha uma única classe de comunicação), as seguintes duas condições devem ser satisfeitas: \(k_0 > 0\) e \(k_0 + k_1 < 1\). É fácil perceber que, se \(k_0 = 0\), com uma ou mais chegadas para cada partida, não há possibilidade de o sistema atingir estabilidade, e o número de clientes no sistema aumentará indefinidamente com o tempo. Se \(k_0 + k_1 = 1\), apenas os dois estados, \(\{0, 1\}\), são possíveis no sistema — Se o sistema começa com \(i > 1\) clientes, uma vez que atinja 0 ou 1, permanecerá em \(\{0, 1\}\). A classificação adicional dos estados depende de \(\mathrm{E}(X_n)\), o número esperado de clientes que chegam durante um tempo de serviço.


Agora, defina a transformada de Laplace–Stieltjes (Harchol-Balter 2010) da distribuição do tempo de serviço como:

\[\begin{align}\\ \psi(\theta) = \int_0^\infty e^{-\theta t} \, dB(t), \quad \mathrm{Re}(\theta) > 0\\\\ \end{align}\]

e a função geradora de probabilidades (PGF) do número de chegadas durante um tempo de serviço como:

\[\begin{align}\\ K(z) = \sum_{j=0}^\infty k_j z^j, \quad |z| \leqslant 1. \end{align}\]

Utilizando as propriedades clássicas das transformadas de Laplace–Stieltjes e das PGFs, obtém-se, a partir da equação para \(k_j\), que:

\[\begin{align}\\ K(z) = \int_0^\infty e^{-\lambda t} \sum_{j=0}^\infty \frac{(\lambda t z)^j}{j!} \, dB(t) = \int_0^\infty e^{-(\lambda - \lambda z)t} \, dB(t) = \psi(\lambda - \lambda z)\\\\ \end{align}\]

Portanto,

\[\begin{align}\\ K'(z) &= -\lambda \psi'(\lambda - \lambda z)\\\\ K'(1) &= -\lambda \psi'(0) = \lambda b \\\\ \end{align}\]

Note que \(\lambda b = (\text{taxa de chegada}) \times (\text{tempo médio de serviço})\). Esta quantidade, em particular, é chamada de intensidade de tráfego do modelo \(M/G/1\), denotada por \(\rho\). O valor de \(\rho\), assim como nos modelos anteriores, determina se o sistema está em equilíbrio (atinge estado estacionário) quando o parâmetro de tempo \(n\) (relacionado aos instantes \(t_n\)) tende ao infinito.


(I) Distribuição Estacionária


Seja \(\boldsymbol{\pi} = (\pi_0, \pi_1, \pi_2, \dots)\) a distribuição limite da cadeia embutida \(\{Q_n, n = 0, 1, 2, \dots\}\) do modelo \(M/G/1\). Fazendo o uso matriz de probabilidades de transição na equação \(\boldsymbol{\pi}P = \boldsymbol{\pi}\), obtemos:

\[\begin{align}\\ k_0 \pi_0 + k_0 \pi_1 &= \pi_0\\\\ k_1 \pi_0 + k_1 \pi_1 + k_0 \pi_2 &= \pi_1\\\\ k_2 \pi_0 + k_2 \pi_1 + k_1 \pi_2 + k_0 \pi_3 &= \pi_2\\\\ \vdots\\\\ \end{align}\]

A solução dessas equações, naturalmente, deve ser obtida por meio de procedimentos computacionais. Para isso, define-se \(\nu_0 \equiv 1\) e \(\nu_i = \dfrac{\pi_i}{\pi_0}\), reescrevendo-se as equações em termos de \(\nu_i\) (\(i = 1, 2, \dots\)) da seguinte maneira:

\[\begin{align}\\ \nu_1 &= \frac{1 - k_0}{k_0}\\\\ \nu_2 &= \frac{1 - k_1}{k_0}\nu_1 - \frac{k_1}{k_0}\\\\ &\vdots\\\\ \nu_j &= \frac{1 - k_1}{k_0}\nu_{j-1} - \frac{k_2}{k_0} \nu_{j-2} - \cdots - \frac{k_{j-1}}{k_0}\nu_{1} - \frac{k_{j-1}}{k_0}\\\\ &\vdots\\\\ \end{align}\]

Essas equações podem ser resolvidas de forma recursiva para determinar os valores de \(\nu_i\) (\(i = 1, 2, \dots\)). Ademais, as probabilidades limite \((\pi_0, \pi_1, \pi_2, \dots)\) são funções monótonas e côncavas em relação ao estado, de modo que, para valores suficientemente elevados de \(i\), essas probabilidades se tornam extremamente pequenas. Consequentemente, os valores de \(\nu_i = \pi_i / \pi_0\) também apresentam esse mesmo comportamento. Essa propriedade permite, do ponto de vista computacional, a definição de um valor de corte apropriado para o espaço de estados, restringindo-o a uma dimensão finita sem perda significativa de precisão. Assim, para recuperar os valores de \(\pi_i\) a partir de \(\nu_i\), observa-se que:

\[\begin{align}\\ \sum_{i=0}^\infty \nu_i = 1 + \sum_{i=1}^\infty \frac{\pi_i}{\pi_0} = \frac{1}{\pi_0} \sum_{i=0}^\infty \pi_i = \frac{1}{\pi_0}\\\\ \end{align}\]

uma vez que a condição de normalização é dada por \(\sum_{i=0}^\infty \pi_i = 1\). Logo, obtém-se que:

\[\begin{align}\\ \pi_0 = \left(\sum_{i=0}^\infty \nu_i \right)^{-1} = \left(1 + \sum_{i=1}^\infty \nu_i \right)^{-1} \\\\ \end{align}\]

e,

\[\begin{align}\\ \pi_i = \frac{\nu_i}{1 + \sum_{i=1}^\infty \nu_i} = \nu_i \cdot \pi_0\\\\ \end{align}\]

Do ponto de vista analítico, a distribuição limite \((\pi_0, \pi_1, \pi_2, \dots)\) pode ser determinada por meio da resolução das equações anteriores utilizando funções geradoras. Entretanto, a obtenção de expressões explícitas para as probabilidades individuais requer, em geral, a inversão da função geradora de probabilidades resultante, o que nem sempre é trivial. No caso do modelo \(M/G/1\) pode-se definir:

\[\begin{align}\\ \Pi(z) = \sum_{j=0}^\infty \pi_j z^j, \quad |z| \leqslant 1\\\\ \end{align}\]

e,

\[\begin{align}\\ K(z) = \sum_{j=0}^\infty k_j z^j, \quad |z| \leqslant 1\\\\ \end{align}\]

Multiplicando as equações definidas no sistema de equações \(\boldsymbol{\pi}P = \boldsymbol{\pi}\) por potências apropriadas de \(z\) e somando-as, pode-se reescrever \(\Pi(z)\) em função de \(K(z)\) da seguinte forma:

\[\begin{align}\\ \Pi(z) &= \pi_0 K(z) + \pi_1 K(z) + \pi_2 z K(z) + \cdots\\\\ &= \pi_0 K(z) + \dfrac{K(z)}{z} \left[\pi_1 z + \pi_2 z^2 + \pi_3 z^3 + \cdots \right]\\\\ &= \pi_0 K(z) + \dfrac{K(z)}{z} [\Pi(z) - \pi_0]\\\\ \end{align}\]

isto é,

\[\begin{align}\\ \Pi(z) \left[1 - \frac{K(z)}{z}\right] &= \pi_0 K(z) \left[1 - \frac{1}{z}\right] \quad \implies \quad \Pi(z) = \frac{\pi_0 K(z)(z - 1)}{z - K(z)}\\\\ \end{align}\]

Nesta equação, a quantidade desconhecida \(\pi_0\) pode ser determinada utilizando a condição de normalização \(\sum_{j=0}^\infty \pi_j = 1\). Para tal, considere:

\[\begin{align}\\ \Pi(1) = \sum_{j=0}^\infty \pi_j = 1\\\\ \end{align}\]

Então, fazendo \(z \to 1\) em \(\Pi(z)\) e aplicando a regra de L’Hôpital, obtém-se que:

\[\begin{align}\\ 1 = \frac{\lim_{z \to 1}\pi_0 [K(z) - (z - 1) K'(z)]}{\lim_{z \to 1}[1 - K'(z)]}\\\\ \end{align}\]

Como \(K(1) = 1\) e \(K'(1) = \rho\), tem-se que que a constante \(\pi_0\) é definida por:

\[\begin{align}\\ \pi_0 = 1 - \rho \end{align}\]

e, portanto,

\[\begin{align}\\ \Pi(z) = \frac{(1 - \rho)(z - 1) K(z)}{z - K(z)}\\\\ \end{align}\]

As expressões explícitas para as probabilidades \({\pi_j, j = 0, 1, 2, \dots}\) podem ser obtidas, em casos especiais, por meio da expansão em série de potências da função geradora \(\Pi(z)\). De fato, uma formulação alternativa de \(\Pi(z)\) mais conveniente para a realização dessa expansão é descrita como:

\[\begin{align}\\ \Pi(z) = \frac{(1 - \rho)K(z)}{[z - K(z)]/[z - 1]} = \frac{(1 - \rho)K(z)}{1 - [1 - K(z)]/[1 - z]}\\\\ \end{align}\]

Note que \(\sum_{j=0}^\infty z^j(k_{j+1} + k_{j+2} + \ldots)\) pode ser escrito como:

\[\begin{align}\\ \sum_{j=0}^\infty z^j(k_{j+1} + k_{j+2} + \ldots) = \frac{1 - K(z)}{1 - z} = C(z)\\\\ \end{align}\]

Para \(|z| \leqslant 1\), tem-se que:

\[\begin{align}\\ |C(z)| = \left| \frac{1 - K(z)}{1 - z} \right| < 1 \quad \text{ se } \quad \rho < 1\\\\ \end{align}\]

Agora, utilizando a expansão em série geométrica, pode-se escrever \(\Pi(z)\) como:

\[\begin{align}\\ \Pi(z) = (1 - \rho)K(z) \sum_{j=0}^\infty [C(z)]^j\\\\ \end{align}\]

Neste caso, a expressão explícita para \(\pi_j\) é obtida expandindo o lado direito da equação acima como uma série em potências de \(z\) e selecionando o coeficiente de \(z^j\). Além disso, a função geradora \(\Pi(z)\) também fornece a distribuição limite \(\lim_{t \to \infty} Q(t)\). De fato, a expressão explícita para a distribuição limite de \(Q(t)\) (e também da cadeia embutida \(Q_n\)) é dada por:

\[\begin{align}\\ \pi_0 &= 1 - \rho\\\\ \pi_j &= (1 - \rho) \int_{0}^\infty e^{-\lambda t} \sum_{n=0}^\infty \left[ \frac{ (\lambda t)^{n + j - 1}}{(n + j - 1)!} - \frac{ (\lambda t)^{n + j}}{(n + j)!} \right] dB_n(t)\\\\ \end{align}\]

para \(\rho < 1\), onde \(B_n(t)\) é a \(n\)-convolução de \(B(t)\) com ela mesma.


(II) Número Médio de Clientes


Seja \(Q = Q(\infty)\) o número de clientes no sistema no regime estacionário, e \(Q_q\) o número de clientes na fila, excluindo aquele que está eventualmente em serviço. Define-se:


  • \(L = \mathrm{E}[Q]\), que representa o número médio de clientes no sistema;
  • \(L_q = \mathrm{E}[Q_q]\), que representa o número médio de clientes na fila.


A partir da função geradora de probabilidades associada à distribuição estacionária de \(Q\), denotada por:

\[\begin{align}\\ \Pi(z) = \sum_{j=0}^\infty \pi_j z^j, \quad |z| \leqslant 1\\\\ \end{align}\]

é possível determinar tanto o valor esperado como a variância do número de clientes no sistema, utilizando técnicas diferenciais aplicadas à PGF. De fato, lembrando que a função geradora \(\Pi(z)\) pode ser escrita como:

\[\begin{align}\\ \Pi(z) = (1 - \rho) \cdot \dfrac{(z - 1) K(z)}{z - K(z)}\\\\ \end{align}\]

onde \(K(z)\) é a função geradora da distribuição de número de chegadas durante um tempo de serviço genérico, definida como:

\[\begin{align}\\ K(z) = \sum_{j=0}^\infty k_j z^j\\\\ \end{align}\]

com \(k_j = \mathrm{P}(\text{chegam } j \text{ clientes durante um tempo de serviço})\), a média e a variância do número de clientes no sistema são obtidas por meio das derivadas primeira e segunda da PGF avaliadas em \(z = 1\). Formalmente,

\[\begin{align}\\ L &= \mathrm{E}[Q] = \Pi'(1)\\\\ \mathrm{Var}(Q) &= \Pi''(1) + \Pi'(1) - \left[\Pi'(1)\right]^2\\\\ \end{align}\]

Então, calculando a derivada de \(\Pi(z)\), obtém-se:

\[\begin{align}\\ \Pi'(z) = \dfrac{(1 - \rho)}{[z - K(z)]^2} \Big\{ [z - K(z)]\cdot\big[K(z) + (z - 1)K'(z)\big] - (z - 1)\cdot[1 - K'(z)]\cdot K(z) \Big\}\\\\ \end{align}\]

Aplicando a regra de L’Hôpital duas vezes quando \(z \to 1\), tem-se:

\[\begin{align}\\ \Pi'(1) = \dfrac{2K'(1)\cdot[1 - K'(1)] + K''(1)}{2[1 - K'(1)]}\\\\ \end{align}\]

Agora, dado que \(K'(1) = \rho\) e \(K''(1) = \lambda^2 \mathrm{E}[S^2]\), onde \(\mathrm{E}[S^2]\) é o segundo momento do tempo de serviço, a média do número de clientes no sistema resulta em:

\[\begin{align}\\ L = \mathrm{E}[Q] = \dfrac{\rho + \dfrac{\lambda^2 \mathrm{E}[S^2]}{2}}{1 - \rho} = \rho + \dfrac{\lambda^2 \mathrm{E}[S^2]}{2(1 - \rho)}\\\\ \end{align}\]

enquanto a variância do número de clientes no sistema é dada por:

\[\begin{align}\\ \mathrm{Var}(Q) = \rho(1 - \rho) + \dfrac{\lambda^2 \mathrm{E}[S^2]}{2(1 - \rho)}\cdot\left[3 - 2\rho + \dfrac{\lambda^2 \mathrm{E}[S^2]}{2(1 - \rho)}\right] + \dfrac{\lambda^3 \mathrm{E}[S^3]}{3(1 - \rho)}\\\\ \end{align}\]

A partir da relação \(\mathrm{Var}(S) = \mathrm{E}[S^2] - [\mathrm{E}[S]]^2\), é possível expressar \(L\) em termos da média e da variância do tempo de serviço. De fato, considerando que \(\rho = \lambda \mathrm{E}[S]\), a fórmula de \(L\) pode ser reescrita como:

\[\begin{align}\\ L = \rho + \dfrac{\rho^2}{2(1 - \rho)} + \dfrac{\lambda^2 \mathrm{Var}(S)}{2(1 - \rho)}\\\\ \end{align}\]

Essa formulação, em particular, evidencia que o tamanho médio da fila cresce não apenas com a carga média \(\rho\), mas também com a variabilidade dos tempos de serviço, refletida no termo proporcional à variância \(\mathrm{Var}(S)\). De fato, quando \(\mathrm{Var}(S) = 0\), isto é, quando o tempo de serviço é constante, tem-se que:

\[\begin{align}\\ L = \rho + \dfrac{\rho^2}{2(1 - \rho)} = \dfrac{\rho}{1 - \rho} \left[1 - \dfrac{\rho}{2}\right]\\\\ \end{align}\]

Por outro lado, quando a distribuição do tempo de serviço é uma distribuição Erlang com média \(1/\lambda\) e parâmetro de escala \(k\) tal que \(\mathrm{Var}(S) = 1/(k\lambda^2)\), tem-se que:

\[\begin{align}\\ L = \rho + \dfrac{\rho^2}{2(1 - \rho)} + \dfrac{\rho^2}{2k(1 - \rho)} = \rho + \dfrac{\rho^2(1 + k)}{2k(1-\rho)}\\\\ \end{align}\]

A partir dessas relações, tem-se que a média do número de clientes na fila, denotada por \(L_q\), é dada pela diferença \(L_q = L - \rho\).


(III) Tempo de Espera do Cliente


Seja \(T\) o tempo total que um cliente permanece no sistema, composto pelo tempo de espera na fila (\(T_q\)) e o tempo de serviço (\(S\)). Então, considerando o regime estacionário, defina \(W = \mathrm{E}[T]\) como o tempo médio no sistema e \(W_q = \mathrm{E}[T_q]\) como o tempo médio de espera na fila. Defina, também, \(F(\cdot)\) como a função de distribuição de \(T\), cuja transformada de Laplace–Stieltjes é dada por:

\[\begin{align}\\ \Phi(\theta) = \int_0^\infty e^{-\theta t} \, dF(t), \quad \text{Re}(\theta) > 0\\\\ \end{align}\]

Suponha que um cliente que acaba de deixar o sistema, após ter permanecido nele por um tempo total \(T\). Se, no instante de sua partida, esse cliente deixa \(n\) outros clientes no sistema, isso significa que esses \(n\) clientes chegaram durante o intervalo de duração \(T\). Portanto, a probabilidade de que haja \(n\) clientes no sistema, no regime estacionário, é dada por:

\[\begin{align}\\ \mathrm{P}(Q = n) = \int_0^\infty e^{-\lambda t} \dfrac{(\lambda t)^n}{n!} \, dF(t), \quad n \geqslant 0\\\\ \end{align}\]

Neste caso, a partir da definição de função geradora de probabilidade, tem-se que:

\[\begin{align}\\ \Pi(z) &= \sum_{n=0}^\infty \mathrm{P}(Q = n) z^n \\\\ &= \sum_{n=0}^\infty z^n \int_{0}^\infty e^{-\lambda t} \dfrac{(\lambda t)^n}{n!} \, dF(t)\\\\ &= \int_{0}^\infty e^{-\lambda t} \sum_{n=0}^\infty \dfrac{(\lambda t z)^n}{n!} \, dF(t)\\\\ &= \Phi(\lambda - \lambda z)\\\\ \end{align}\]

Comparando com a função gerada \(\Pi(z)\) obtida anteriormente, obtém-se que:

\[\begin{align}\\ \dfrac{(1 - \rho)(z - 1) K(z)}{z - K(z)} = \Phi(\lambda - \lambda z)\\\\ \end{align}\]

Agora, dado que \(K(z) = \psi(\lambda - \lambda z)\), onde \(\psi(\cdot)\) é a transformada de Laplace–Stieltjes da distribuição dos tempos de serviço, tem-se que:

\[\begin{align}\\ \Phi(\lambda - \lambda z) = \dfrac{(1 - \rho)(z - 1) \psi(\lambda - \lambda z)}{z - \psi(\lambda - \lambda z)}\\\\ \end{align}\]

Defina \(\theta = \lambda - \lambda z\), ou equivalentemente \(z = 1 - \dfrac{\theta}{\lambda}\). Então, a transformada de Laplace–Stieltjes do tempo no sistema pode ser reescrita como:

\[\begin{align}\\ \Phi(\theta) = \dfrac{(1 - \rho) \, \theta \, \psi(\theta)}{\theta - \lambda [1 - \psi(\theta)]}\\\\ \end{align}\]

Como o tempo no sistema \(T\) é a soma do tempo de espera na fila \(T_q\) com o tempo de serviço \(S\), a transformada de Laplace–Stieltjes do tempo de espera na fila, denotada por \(\Phi_q(\theta)\), é obtida a partir da relação funcional:

\[\begin{align}\\ \Phi(\theta) = \psi(\theta) \cdot \Phi_q(\theta)\\\\ \end{align}\]

a qual, por comparação com a expressão previamente obtida para \(\Phi(\theta)\), fornece:

\[\begin{align}\\ \Phi_q(\theta) = \dfrac{(1 - \rho) \, \theta}{\theta - \lambda [1 - \psi(\theta)]} = (1 - \rho) \sum_{n=0}^\infty \left[ \dfrac{\lambda}{\theta} \{1 - \psi(\theta)\} \right]^n\\\\ \end{align}\]

válida sob a condição de convergência \(\left| \lambda/\theta \{1 - \psi(\theta)\} \right| < 1\), a qual se satisfaz sempre que \(\rho < 1\). Assim, ao calcular a inversa dessa transformada, obtém-se a função de distribuição do tempo de espera na fila, \(T_q\), dada por a qual se satisfaz sempre que \(\rho < 1\). Assim, ao calcular a inversa dessa transformada, obtém-se a função de distribuição do tempo de espera na fila, \(T_q\), dada por:

\[\begin{align}\\ F_q(t) = (1 - \rho) \sum_{n=0}^\infty \rho^n R^{(n)}(t)\\\\ \end{align}\]

em que \(R^{(n)}(t)\) representa a \(n\)-ésima convolução da função de distribuição do tempo residual de serviço (ou forward recurrence time) consigo mesma. A partir dessa distribuição e da aplicação da lei de Little, que estabelece \(L = \lambda W\) e \(L_q = \lambda W_q\), obtém-se que os tempos médios, no sistema e de espera na fila, são dados, respectivamente, por:

\[\begin{align}\\ W = \mathrm{E}[S] + \dfrac{\lambda \, \mathrm{E}[S^2]}{2(1 - \rho)}\\\\ \end{align}\]

e

\[\begin{align}\\ W_q = \dfrac{\lambda \, \mathrm{E}[S^2]}{2(1 - \rho)}\\\\ \end{align}\]

(IV) Período de Ocupação


No contexto de uma cadeia de Markov embutida, a duração do período de ocupação é mensurada pelo número de transições realizadas pela cadeia até que o estado \(0\) seja visitado pela primeira vez, a partir de um estado inicial distinto de zero. Mais precisamente, seja \(B_i\) o número de transições realizadas pela cadeia antes de atingir o estado \(0\) pela primeira vez, dado que ela iniciou no estado \(i\). Defina:

\[\begin{align}\\ g_i^{(n)} = \mathrm{P}(B_i = n), \quad n = 1, 2, \dots\\\\ \end{align}\]

em que \(g_i^{(n)}\) representa a distribuição de probabilidade do número de transições no período de ocupação iniciado no estado \(i\). Uma propriedade fundamental de \(B_i\) é que ele pode ser interpretado como a soma de \(i\) variáveis aleatórias independentes, cada uma com a mesma distribuição de \(B_1\). Isso equivale a afirmar que a transição do estado \(i\) para o estado \(0\) ocorre em \(i\) etapas sucessivas:

\[\begin{align}\\ i \to i-1, i-1 \to i-2, \ldots, 1 \to 0\\\\ \end{align}\]

Isso se justifica pelo fato de que as transições descendentes podem ocorrer apenas um passo por vez. Como todas essas transições são estruturalmente idênticas, conclui-se que \(B_i\) é a soma de \(i\) variáveis aleatórias independentes, cada uma distribuída como \(B_1\). Consequentemente, \(g_i^{(n)}\) corresponde à \(i\)-ésima convolução de \(g_1^{(n)}\) consigo mesma. Assim, a função geradora de probabilidades associada a \(g_i^{(n)}\) é dada por:

\[\begin{align}\\ G_i(z) = \sum_{n=i}^\infty g_i^{(n)} z^n = [G_1(z)]^i\\\\ \end{align}\]

Seja, agora, \(k_r\) a probabilidade de que ocorram \(r\) chegadas durante um período de serviço. Como o período de ocupação não pode se encerrar antes da \(i\)-ésima transição da cadeia, estabelecem-se as seguintes relações de recorrência para a distribuição de \(B_i\):

\[\begin{align}\\ g_i^{(i)} &= k^{(i)}_0\\\\ g_i^{(n)} &= \sum_{r=1}^{n-1} k^{(i)}_r \, g_r^{(n-i)}\\\\ \end{align}\]

em que \(k^{(i)}_r\) representa a \(i\)-ésima convolução de \(k_r\). No caso particular de \(i = 1\), as equações se reduzem a:

\[\begin{align}\\ g_1^{(1)} &= k_0\\\\ g_1^{(n)} &= \sum_{r=1}^{n-1} k_r \, g_r^{(n-1)} \end{align}\]

Multiplicando ambos os lados dessas equações pelas potências adequadas de \(z\) e somando-as, obtém-se a seguinte equação para a função geradora de probabilidades \(G_1(z)\):

\[\begin{align}\\ G_1(z) = zk_0 + z \sum_{r=1}^\infty k_r \sum_{n=r+1}^\infty z^{n-1}g_r^{(n-1)} = zK[G_1(z)]\\\\ \end{align}\]

em que \(K(z)\) é a função geradora de probabilidades do número de chegadas durante um único tempo de serviço, que, conforme previamente definido, é dada por:

\[\begin{align}\\ K(z) = \psi(\lambda - \lambda z)\\\\ \end{align}\]

sendo \(\psi(\theta)\) a transformada de Laplace–Stieltjes da distribuição dos tempos de serviço. Assim, a função geradora de probabilidades \(G_1(z) \equiv G(z)\) satisfaz a seguinte relação:

\[\begin{align}\\ \omega = z \cdot \psi\big(\lambda - \lambda\omega\big)\\\\ \end{align}\]

Dessa equação, tem-se que a média do número de transições no período de ocupação pode ser obtida por diferenciação implícita. Para isso, seja:

\[\begin{align}\\ G(z) = z \cdot \psi\big(\lambda - \lambda G(z)\big)\\\\ \end{align}\]

Derivando ambos os lados em relação a \(z\), tem-se:

\[\begin{align}\\ G'(z) = \psi\big[\lambda - \lambda G(z)\big] + z \cdot \psi'\big[\lambda - \lambda G(z)\big][-\lambda G'(z)]\\\\ \end{align}\]

Fazendo \(z \to 1\), obtém-se:

\[\begin{align}\\ G'(1) &= \psi(0) - \lambda \psi'(0) G'(1) \quad \implies \quad G'(1)= \dfrac{1}{1 + \lambda \psi'(0)}\\\\ \end{align}\]

Dado que \(\lambda \psi'(0) = -\rho\), segue que a média do número de transições durante o período de ocupação iniciado por um único cliente é dada por:

\[\begin{align}\\ E[B_1] = \dfrac{1}{1 - \rho}\\\\ \end{align}\]

De forma análoga, considerando que o período de ocupação iniciado por \(i\) clientes é a soma de \(i\) períodos independentes distribuídos como \(B_1\), segue que:

\[\begin{align}\\ E[B_i] = \dfrac{i}{1 - \rho}\\\\ \end{align}\]

Até este ponto, a contagem foi realizada em número de transições; assim, o comprimento médio do período de ocupação, em unidades de tempo, é obtido multiplicando-se esse valor pelo tempo médio de serviço \(E[S]\), ou seja,

\[\begin{align}\\ \text{Comprimento médio do período de ocupação} = \dfrac{E[S]}{1 - \rho}\\\\ \end{align}\]

Por fim, considerando que um ciclo completo de operação é composto por um período de ocupação seguido de um período de ociosidade — cujo tempo médio, no modelo \(M/G/1\) com taxa de chegada \(\lambda\), é \(1/\lambda\) —, o comprimento médio de um ciclo completo é dado por:

\[\begin{align}\\ \text{Comprimento médio do ciclo} = \dfrac{E[S]}{1 - \rho} + \dfrac{1}{\lambda} = \dfrac{1}{\lambda(1 - \rho)}\\\\ \end{align}\]

Este resultado evidencia que, embora o tempo médio de ocupação dependa da distribuição dos tempos de serviço (via \(E[S]\)), o comprimento médio do ciclo é determinado apenas pelos parâmetros \(\lambda\) e \(\rho\).

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 5.3 (Oficina Mecânica com Chegadas Poissonianas e Tempos de Serviço Genéricos - Modelo \(M/G/1\)). Em uma oficina mecânica operada por um único mecânico, os registros mantidos pelo proprietário indicam a seguinte distribuição para o número de veículos que chegam durante o tempo de serviço de um veículo:


  • \(\mathrm{P}(\text{nenhuma chegada}) = 0{.}5\)
  • \(\mathrm{P}(\text{uma chegada}) = 0{.}3\)
  • \(\mathrm{P}(\text{duas chegadas}) = 0{.}2\)


Assumindo que as chegadas de veículos seguem um processo de Poisson, este sistema pode ser modelado como uma fila do tipo \(M/G/1\), mesmo na ausência de informações sobre a distribuição dos tempos de serviço. Dessa forma, com base nos dados fornecidos, os coeficientes da função geradora do número de chegadas durante um tempo de serviço são:

\[\begin{align}\\ k_0 = 0.5,\quad k_1 = 0.3,\quad k_2 = 0.2\\\\ \end{align}\]

A intensidade de tráfego \(\rho\) corresponde à média do número de chegadas durante um tempo de serviço, sendo calculada por:

\[\begin{align}\\ \rho = E(\text{número de chegadas}) = 0 \cdot k_0 + 1 \cdot k_1 + 2 \cdot k_2 = 0.7\\\\ \end{align}\]

Como a distribuição dos tempos de serviço não é especificada, o método computacional baseado na determinação dos coeficientes \(\nu_j\) torna-se adequado para obter a distribuição estacionária do número de clientes no sistema. Esses coeficientes são calculados recursivamente segundo as expressões:

\[\begin{align}\\ \nu_1 &= \frac{1 - k_0}{k_0}\\\\ \nu_2 &= \frac{1 - k_1}{k_0}\nu_1 - \frac{k_1}{k_0}\\\\ &\vdots\\\\ \nu_j &= \frac{1 - k_1}{k_0}\nu_{j-1} - \frac{k_2}{k_0} \nu_{j-2} - \cdots - \frac{k_{j-1}}{k_0}\nu_{1} - \frac{k_{j-1}}{k_0}\\\\ &\vdots\\\\ \end{align}\]

Neste exemplo, considerando que \(k_n = 0\) para \(n \geqslant 3\), as fórmulas simplificam-se para:

\[\begin{align}\\ \nu_1 &= \frac{1 - k_0}{k_0}\\\\ \nu_2 &= \frac{1 - k_1}{k_0}\nu_1 - \frac{k_1}{k_0}\\\\ \nu_3 &= \frac{1 - k_1}{k_0}\nu_{2} - \frac{k_2}{k_0} \nu_{1} - \frac{k_{2}}{k_0}\\\\ \nu_j &= \frac{1 - k_1}{k_0}\nu_{j-1} - \frac{k_2}{k_0} \nu_{j-2}, \qquad j \geqslant 4\\\\ \end{align}\]

Baseando-se nessas expressões recursivas, considerou-se aqui o cálculo dos coeficientes \(\nu_j\) até \(j = 10\), pois essa escolha proporciona um equilíbrio adequado entre precisão numérica e complexidade computacional, uma vez que os termos posteriores apresentam valores negligenciáveis, garantindo assim a convergência e uma aproximação confiável da distribuição estacionária. Assim, os valores dos coeficientes \(\nu_j\), para \(j = 0, \ldots, 10\), são dados por:

\[\begin{align}\\ \begin{matrix} \nu_0 = 1 &&& \nu_1 = 1 &&& \nu_2 = 0.8 \\ \nu_3 = 0.32 &&& \nu_4 = 0.128 &&& \nu_5 = 0.051 \\ \nu_6 = 0.021 &&& \nu_7 = 0.008 &&& \nu_8 = 0.003 \\ \nu_9 = 0.001 &&& \nu_{10} = 0.001&&&\\ \end{matrix}\\\\ \end{align}\]

O somatório total dos coeficientes é dado por:

\[\begin{align}\\ \sum_{i=0}^{10} \nu_i = 3.333\\\\ \end{align}\]

Assim, a probabilidade de não haver clientes no sistema é descrita como:

\[\begin{align}\\ \pi_0 = \frac{1}{\sum_{i=0}^{10} \nu_i} = 0.300\\\\ \end{align}\]

e as demais probabilidades são obtidas de acordo com a relação:

\[\begin{align}\\ \pi_i = \nu_i \cdot \pi_0\\\\ \end{align}\]

que resulta em:

\[\begin{align}\\ \begin{matrix} \pi_0 = 0.300 &&& \pi_1 = 0.300 &&& \pi_2 &= 0.240 \\ \pi_3 = 0.096 &&& \pi_4 = 0.038 &&& \pi_5 &= 0.015 \\ \pi_6 = 0.006 &&& \pi_7 = 0.002 &&& \pi_8 &= 0.001 \\ \pi_9 = 0.000 &&& \pi_{10} = 0.000 &&&\\ \end{matrix}\\\\ \end{align}\]

Neste caso, o número médio de veículos no sistema é dado por:

\[\begin{align}\\ L = E(Q) = \sum_{i=0}^{10} i \cdot \pi_i = 1.366\\\\ \end{align}\]

E, por fim, aplicando a lei de Little, o tempo médio no sistema é:

\[\begin{align}\\ W = \frac{L}{\lambda} = \frac{1.366}{0.7} = 1.951\\\\ \end{align}\]

Desses resultados, a intensidade de tráfego \(\rho = 0.7\) indica que o sistema está operando com 70% da capacidade total do mecânico, o que caracteriza um sistema estável, porém suscetível à formação de filas. A probabilidade \(\pi_0 = 0.3\) mostra que o sistema estará vazio — isto é, sem veículos em atendimento ou esperando — em 30% do tempo, evidenciando uma alternância entre períodos de ociosidade e atendimento ativo. As probabilidades \(\pi_i\) apresentam uma queda rápida conforme \(i\) aumenta, indicando que a ocorrência de muitos veículos no sistema é pouco frequente, o que justifica a negligência dos estados mais elevados. Ademais, o valor médio \(L = 1.366\) revela que, em média, há aproximadamente um veículo no sistema, somando os que aguardam e os que estão em atendimento. Por fim, o tempo médio no sistema \(W = 1.951\) unidades de tempo indica que, em média, um veículo permanece quase o dobro do tempo médio de serviço, considerando o tempo de espera e o atendimento.


O Código 3 apresenta a implementação, em linguagem R, dos cálculos analíticos dos principais indicadores de desempenho do sistema de filas \(M/G/1\) aplicado à oficina mecânica, cujas chegadas seguem um processo de Poisson e os tempos de serviço genéricos, considerando a distribuição empírica do número de chegadas durante um tempo de serviço.


Código 3. Implementação, em linguagem R, dos cálculos analíticos dos principais indicadores de desempenho do sistema de filas \(M/G/1\) aplicado à oficina mecânica com chegadas Poissonianas e tempos de serviço genéricos, considerando a distribuição empírica do número de chegadas durante um tempo de serviço.

# -----------------------------------------------------------------
# Análise do Modelo de Filas M/G/1 Aplicado a uma Oficina Mecânica
# -----------------------------------------------------------------

# --- 0. Pacotes Necessários ---

# --- Não são necessários pacotes ---

# --- 1. Parâmetros Iniciais ---

k0          <- 0.5
k1          <- 0.3
k2          <- 0.2

lambda      <- 0* k0 + 1 * k1 + 2 * k2  # Intensidade de Tráfego
n_max       <- 10                       # Valor Máximo de j

# --- 2. Cálculo dos Coeficientes ν_j ---

nu          <- numeric(n_max + 1)
nu[1]       <- 1                                                      # ν_0
nu[2]       <- (1 - k0) / k0                                          # ν_1
nu[3]       <- (1 - k1) / k0 * nu[2] - (k1 / k0)                      # ν_2
nu[4]       <- (1 - k1) / k0 * nu[3] - (k2 / k0) * nu[2] - (k2 / k0)  # ν_3

# Para ν_j Com j >= 4

for (j in 4:n_max)
{
  nu[j + 1] <- (1 - k1) / k0 * nu[j] - (k2 / k0) * nu[j - 1]
}

# --- 3. Cálculo das Probabilidades Estacionárias π_i ---

sum_nu      <- sum(nu)
pi_0        <- 1 / sum_nu
pi          <- nu * pi_0

# --- 4. Resultados de ν_j e π_i ---

df_nu       <- data.frame(j = 0:n_max,
                          nu = round(nu, 3),
                          pi = round(pi, 3))

print(df_nu)
    j    nu    pi
1   0 1.000 0.300
2   1 1.000 0.300
3   2 0.800 0.240
4   3 0.320 0.096
5   4 0.128 0.038
6   5 0.051 0.015
7   6 0.020 0.006
8   7 0.008 0.002
9   8 0.003 0.001
10  9 0.001 0.000
11 10 0.001 0.000
# --- 4. Métricas do Sistema ---

rho         <- lambda
L           <- sum((0:n_max) * pi)   # Número Médio no Sistema
W           <- L / rho               # Tempo Médio no Sistema (Lei de Little)

resultados  <- data.frame(
                      Indicador = c(
                            "Intensidade de Tráfego (ρ)",
                            "Probabilidade do Sistema Vazio (π0)",
                            "Número Médio no Sistema (L)",
                            "Tempo Médio no Sistema (W)"
                          ),
                          Valor = c(
                            round(rho, 3),
                            round(pi_0, 3),
                            round(L, 3),
                            round(W, 3)
                          )
)

print(resultados)
                            Indicador Valor
1          Intensidade de Tráfego (ρ) 0.700
2 Probabilidade do Sistema Vazio (π0) 0.300
3         Número Médio no Sistema (L) 1.366
4          Tempo Médio no Sistema (W) 1.951

\[\small \begin{align} \tag*{$\blacksquare$}\\\\ \end{align}\]