Теория очередей и модель M/M/1

Очередь в банке, заявки в колл-центре, пакеты на сетевом интерфейсе — всё это системы массового обслуживания. У них есть красивая математика, и простейшая её модель называется M/M/1. Сегодня мы и посчитаем её формулами, и проверим симуляцией.

M/M/1 — это система массового обслуживания с одним сервером, в которой клиенты приходят пуассоновским потоком (интервалы между приходами распределены экспоненциально со средним 1/lambda), а время обслуживания тоже экспоненциально (со средним 1/mu). Три символа кодируют: M — приходы без памяти, M — обслуживание без памяти, 1 — один сервер.

Зачем это нужно? Теория очередей отвечает на практические вопросы: сколько касс открыть, чтобы клиенты не ждали дольше пяти минут? Хватит ли одного сервера, или при росте нагрузки очередь «взорвётся»? M/M/1 — самая простая модель, на которой эти эффекты видно как на ладони.

Загрузка системы rho

Ключевая величина — загрузка (utilization) rho = lambda / mu. Это отношение скорости прихода к скорости обслуживания, то есть доля времени, которую сервер занят.

Если rho < 1 — сервер успевает обслуживать поток, система стабильна, очередь колеблется вокруг конечного среднего.
Если rho >= 1 — приходят быстрее, чем уходят, очередь растёт неограниченно. Система нестабильна.
Если rho стремится к 1 снизу (например 0.95, 0.99) — система формально стабильна, но очередь становится огромной.

Теоретические формулы

Для стабильной M/M/1 (rho < 1) теория даёт замечательно простые ответы:

Среднее число клиентов в системе (в очереди плюс на обслуживании): L = rho / (1 - rho).
Среднее время, проведённое клиентом в системе (ожидание плюс обслуживание): W = 1 / (mu - lambda).

Эти две формулы связаны законом Литтла L = lambda * W — фундаментальным результатом теории очередей. Посмотрите, как L ведёт себя у границы: при rho = 0.5 получаем L = 1; при rho = 0.9 уже L = 9; при rho = 0.99 — целых L = 99. Рост лавинообразный.

Симуляция M/M/1

Проверим формулы дискретно-событийной симуляцией. У нас два типа будущих событий: следующий приход (next_arr) и следующий уход (next_dep). На каждой итерации мы берём ближайшее из них, переводим часы и обрабатываем. Попутно копим площадь под графиком числа клиентов (area), чтобы усреднением получить L по времени.

import random, math
from collections import deque
random.seed(11)
lam, mu = 0.7, 1.0
clock = 0.0
next_arr = random.expovariate(lam)
next_dep = math.inf
in_system = 0
served = 0
total_wait = 0.0
arrivals = deque()
END = 20000.0
area = 0.0
last = 0.0
while clock < END:
    t = min(next_arr, next_dep)
    area += in_system * (t - last)
    last = t
    clock = t
    if next_arr <= next_dep:
        in_system += 1
        arrivals.append(clock)
        if in_system == 1:
            next_dep = clock + random.expovariate(mu)
        next_arr = clock + random.expovariate(lam)
    else:
        in_system -= 1
        a = arrivals.popleft()
        total_wait += clock - a
        served += 1
        next_dep = clock + random.expovariate(mu) if in_system > 0 else math.inf
rho = lam / mu
print(f"rho (загрузка) = {rho}")
print(f"Обслужено: {served}")
print(f"L: симуляция={area/END:.3f}, теория={rho/(1-rho):.3f}")
print(f"W: симуляция={total_wait/served:.3f}, теория={1/(mu-lam):.3f}")

Вывод:

rho (загрузка) = 0.7
Обслужено: 14085
L: симуляция=2.377, теория=2.333
W: симуляция=3.375, теория=3.333

Симуляция и теория совпали с точностью до сотых: L получилось 2.377 против 2.333, W — 3.375 против 3.333. Маленькое расхождение — это статистический шум конечного прогона; увеличив END, мы сделали бы его меньше. Главное, что симуляция, не зная формул, воспроизвела их результат — значит, и формулы, и симуляция описывают одно и то же.

Как работает под капотом

Разберём ключевые приёмы.

expovariate — экспоненциальные интервалы. Вызов random.expovariate(lam) возвращает случайный интервал до следующего прихода со средним 1/lam. Экспоненциальное распределение «без памяти»: сколько бы мы ни ждали, ожидаемое время до следующего события не меняется. Именно это свойство делает поток пуассоновским и даёт ту самую букву M.

Два указателя в будущее. Вместо полноценной кучи здесь хватает двух переменных: next_arr и next_dep. Сервер один, поэтому уход в любой момент времени запланирован максимум один. Когда сервер пуст, next_dep = math.inf — «ухода не предвидится», и min всегда выберет приход.

Площадь под графиком. Число клиентов in_system — это ступенчатая функция времени. Чтобы усреднить её по времени, мы на каждом шаге прибавляем in_system * (t - last) — площадь прямоугольника шириной от прошлого события до текущего. В конце делим накопленную площадь на длину прогона и получаем среднее по времени число клиентов — это и есть L.

число клиентов
 3 |          ____
 2 |     ____|    |__
 1 |  __|           |____
 0 |_|                   |___
   +--------------------------- время
   площадь под этой ступенькой / END = L

Очередь приходов. deque хранит времена приходов клиентов, ещё не покинувших систему. Когда клиент уходит, мы достаём из её левого конца время самого старого прихода (popleft) и считаем, сколько он провёл в системе — это вклад в total_wait. В конце total_wait / served даёт среднее W.

Что будет при перегрузке

Поставьте мысленно lam = 0.95, оставив mu = 1.0. Тогда rho = 0.95, и теория обещает L = 0.95 / 0.05 = 19 клиентов в системе и W = 1 / 0.05 = 20 единиц времени. Сравните с нашими 2.33 и 3.33 при rho = 0.7! Стоит подобраться к rho = 0.99 — и L взлетает до 99. Вот почему перегруженные системы (касса в час пик, сервер под нагрузкой) «встают»: при загрузке у самой единицы любая случайная пачка приходов раскачивает очередь до огромных значений, и рассосаться она не успевает. Запас по производительности — не роскошь, а условие устойчивости.

Частые ошибки

Думать, что при rho < 1 очереди не будет. Будет: при rho = 0.9 среднее число в системе уже 9. Стабильность означает лишь, что среднее конечно, а не мало.
Считать L и W линейными по нагрузке. Они растут гиперболически: знаменатель (1 - rho) или (mu - lambda) стремится к нулю, и значения взрываются у границы.
Усреднять число клиентов по событиям, а не по времени. L — это среднее по времени. Нужно взвешивать каждое значение длительностью интервала (площадь под графиком), иначе результат смещён.
Забыть про math.inf для пустого сервера. Если не «отключать» уход, когда система пуста, min может выбрать несуществующее завершение обслуживания и сломать логику.
Доверять одному короткому прогону. Симуляция шумит. Чем дольше END, тем ближе оценки к теории; по одному короткому прогону судить о системе нельзя.

Итоги

M/M/1 — очередь с одним сервером, пуассоновскими приходами (rate lambda) и экспоненциальным обслуживанием (rate mu).
Загрузка rho = lambda/mu — доля занятости сервера; при rho < 1 система стабильна, при rho >= 1 очередь растёт без предела.
Теория даёт L = rho/(1-rho) и W = 1/(mu-lambda), связанные законом Литтла L = lambda*W.
DES-симуляция с двумя событиями (приход/уход) и площадью под графиком воспроизводит эти формулы с точностью до статистического шума.
У границы rho стремится к 1 величины L и W взрываются — поэтому перегруженные системы лавинообразно «встают», и нужен запас по производительности.