Верификация DES-движка: независимая проверка по теории массового обслуживания
Верификация DES-движка: независимая проверка по теории массового обслуживания
Методика независимой математической верификации результатов дискретно-событийной симуляции Storm BPM. Позволяет убедиться, что движок выдаёт корректные результаты, сравнивая их с точными аналитическими формулами из теории массового обслуживания (queueing theory).
Для кого этот гайд
- Аналитики и инженеры, которые хотят убедиться в корректности результатов симуляции
- Руководители проектов, которым нужно обосновать точность движка перед заказчиком
- Исследователи, которые хотят самостоятельно воспроизвести и проверить результаты
Суть подхода
Мы сравниваем результат симуляции с точным аналитическим решением для пяти классических моделей массового обслуживания. Для каждой модели существуют формулы, дающие точные значения метрик (среднее время в системе, длина очереди, утилизация и т.д.). Если симуляция корректна, её результаты должны совпадать с этими формулами в пределах статистической погрешности.
Пять сценариев верификации
| # | Модель | Описание | Что проверяет |
|---|---|---|---|
| 1 | M/M/1 | 1 сервер, пуассоновские приходы, экспоненциальное обслуживание | Базовая корректность движка |
| 2 | M/M/1/K | Конечный буфер + отказы | Механика буферов и потерь |
| 3 | M/M/c | Несколько серверов (формула Эрланга C) | Многоресурсная обработка |
| 4 | Тандем | Две очереди последовательно (сеть Джексона) | Многошаговые процессы |
| 5 | M/D/1 | Детерминированное обслуживание (P-K формула) | Влияние вариативности |
Четыре уровня верификации
| Уровень | Что проверяем | Критерий прохождения |
|---|---|---|
| 1. Точечная проверка | Совпадение при одном значении нагрузки | Отклонение < 10% |
| 2. Чувствительность | Совпадение по всему диапазону нагрузки (ρ = 0.3..0.9) | Отклонение < 10% (< 20% при ρ ≥ 0.85) |
| 3. Доверительные интервалы | 95% CI из 10 прогонов содержит аналитическое значение | Статистический тест |
| 4. Проверка распределений | Генерируемые времена соответствуют заявленному распределению | KS-тест (α = 0.05) |
Почему именно такие допуски?
Важно понимать
Дискретно-событийная симуляция — это стохастический метод. Каждый прогон использует генератор случайных чисел, поэтому результат НИКОГДА не совпадёт с формулой с точностью до десятого знака. Это не баг, а фундаментальное свойство метода Монте-Карло.
Источники расхождения
Конечное число заявок. При ~5000-10000 заявок за прогон выборочное среднее отклоняется от теоретического на 3-8%. Это нормальная статистическая дисперсия.
Разогрев (warm-up). В начале симуляции система не находится в стационарном режиме. Первые 100-500 заявок проходят через пустую систему, что занижает средние. Чем дольше симуляция, тем меньше этот эффект.
Хвост прогона. При остановке симуляции часть заявок находится «в полёте» (в очереди или обработке). Они не учитываются в статистике, что создаёт небольшой сдвиг.
Высокая нагрузка (ρ ≥ 0.85). При утилизации > 85% очереди становятся очень длинными, а сходимость к стационарному режиму — очень медленной. Формально для M/M/1 при ρ → 1 время сходимости → ∞. Поэтому для ρ = 0.9 мы допускаем 20% вместо 10%.
Почему 10% — это нормально для единичного прогона
По центральной предельной теореме, стандартная ошибка среднего ≈ σ/√n. Для экспоненциального распределения σ = μ (среднее = стандартное отклонение), поэтому:
- При n = 5000 заявок: SE ≈ μ/√5000 ≈ 1.4% от среднего
- 95% доверительный интервал: ±2.8%
- С учётом warm-up и хвоста: реально 5-10%
Именно поэтому уровень 3 (доверительные интервалы) — самый строгий. Единичный прогон — это грубая проверка «не сломалось ли что-то», а 10 прогонов с CI — настоящая статистическая валидация.
Аналитические формулы
M/M/1 — одиночный сервер
Самая базовая модель. Приходы по Пуассону (экспоненциальные интервалы), обслуживание экспоненциальное, 1 сервер.
Обозначения:
- λ — интенсивность потока (заявок/сек)
- μ — интенсивность обслуживания (заявок/сек)
- ρ = λ/μ — утилизация (должна быть < 1)
Формулы:
| Метрика | Формула | Смысл |
|---|---|---|
| L = ρ/(1−ρ) | Среднее число заявок в системе | |
| Lq = ρ²/(1−ρ) | Среднее число заявок в очереди | |
| W = 1/(μ−λ) | Среднее время в системе | |
| Wq = ρ/(μ(1−ρ)) | Среднее время ожидания в очереди |
Закон Литтла: L = λ × W (связывает все метрики)
M/M/1/K — конечный буфер
Как M/M/1, но система вмещает максимум K заявок. Лишние отклоняются (потери).
Вероятности состояний:
p_n = (1−ρ) × ρⁿ / (1−ρ^(K+1)) для ρ ≠ 1
Вероятность отказа: P_reject = p_K
Эффективная интенсивность: λ_eff = λ × (1 − P_reject)
M/M/c — Эрланг C (несколько серверов)
c одинаковых серверов, одна общая очередь FIFO.
Формула Эрланга C (вероятность ожидания):
C(c, a) = [a^c / c! × 1/(1−ρ)] / [Σ(k=0..c-1) a^k/k! + a^c/c! × 1/(1−ρ)]
где a = λ/μ (предложенная нагрузка в Эрлангах), ρ = a/c.
Среднее время ожидания: Wq = C(c,a) / (cμ − λ)
Тандем (сеть Джексона)
Две очереди последовательно. По теореме Джексона (1957), в открытой сети M/M/c-очередей с пуассоновским входом каждая очередь ведёт себя как независимая M/M/c.
W_total = W₁ + W₂ = 1/(μ₁−λ) + 1/(μ₂−λ)
M/D/1 — детерминированное обслуживание
Приходы по Пуассону, но время обслуживания постоянное. Частный случай M/G/1.
По формуле Полачека-Хинчина:
- Wq(M/D/1) = ρ / (2μ(1−ρ))
- Ключевой инсайт: Wq(M/D/1) = Wq(M/M/1) / 2
Устранение вариативности обслуживания вдвое сокращает среднее время ожидания.
Python-скрипт верификации
Скрипт не использует внешних библиотек — только стандартную библиотеку Python (math, json, sys). Все формулы реализованы «с нуля» и подробно прокомментированы.
Установка и запуск
# Никаких зависимостей не требуется — только Python 3.6+
python verify_des.py --compute-only mm1 --rho 0.7
Режимы работы
| Команда | Описание |
|---|---|
--compute-only <scenario> | Вычислить аналитические значения без данных DES |
<scenario> <file.json> | Сравнить JSON-результат DES с формулами |
--sweep <scenario> <dir> | Проверить серию файлов (разные ρ) |
--ci <scenario> <dir> | Проверить доверительные интервалы из нескольких прогонов |
--ks <samples.json> | Тест Колмогорова-Смирнова на распределение |
Полный исходный код
Скачать verify_des.py — полный скрипт верификации (~500 строк, без внешних зависимостей).
Развернуть код verify_des.py
#!/usr/bin/env python3
"""
Storm BPM -- DES Engine Verification & Validation Suite
========================================================
Независимая математическая верификация результатов симуляции.
Все формулы выведены из первых принципов теории массового обслуживания.
Ссылки:
[1] Gross, Shortle et al. "Fundamentals of Queueing Theory" 5th ed.
[2] Kleinrock, "Queueing Systems, Vol. 1: Theory" (Wiley)
[3] Law, "Simulation Modeling and Analysis" 5th ed. (McGraw-Hill)
"""
import math
import json
import sys
import os
# ================================================================
# АНАЛИТИЧЕСКИЕ ФОРМУЛЫ
# ================================================================
def mm1_analytics(lambda_rate, mu_rate):
"""
M/M/1: один сервер, пуассоновские приходы, экспоненциальное обслуживание.
rho = lambda/mu (должно быть < 1)
L = rho/(1-rho) -- среднее число в системе
Lq = rho^2/(1-rho) -- среднее число в очереди
W = 1/(mu-lambda) -- среднее время в системе
Wq = rho/(mu*(1-rho)) -- среднее время ожидания
"""
rho = lambda_rate / mu_rate
if rho >= 1.0:
raise ValueError(f"M/M/1 requires rho < 1, got rho={rho:.4f}")
L = rho / (1 - rho)
Lq = rho ** 2 / (1 - rho)
W = 1.0 / (mu_rate - lambda_rate)
Wq = rho / (mu_rate * (1 - rho))
return {
"rho": rho, "L": L, "Lq": Lq, "W": W, "Wq": Wq,
"service_mean": 1.0 / mu_rate,
"lambda": lambda_rate, "mu": mu_rate,
"throughput": lambda_rate,
}
def mm1k_analytics(lambda_rate, mu_rate, K):
"""
M/M/1/K: конечная ёмкость K, лишние заявки отклоняются.
p_n = (1-rho)*rho^n / (1-rho^(K+1)) для rho != 1
P_reject = p_K
lambda_eff = lambda * (1 - P_reject)
"""
rho = lambda_rate / mu_rate
if abs(rho - 1.0) < 1e-10:
p = [1.0 / (K + 1)] * (K + 1)
else:
denom = 1.0 - rho ** (K + 1)
p = [(1.0 - rho) * rho ** n / denom for n in range(K + 1)]
P_reject = p[K]
L = sum(n * p[n] for n in range(K + 1))
lambda_eff = lambda_rate * (1.0 - P_reject)
W = L / lambda_eff if lambda_eff > 1e-15 else float('inf')
Lq = L - (1.0 - p[0])
Wq = Lq / lambda_eff if lambda_eff > 1e-15 else float('inf')
return {
"rho": rho, "K": K, "P_reject": P_reject,
"L": L, "Lq": Lq, "W": W, "Wq": Wq,
"lambda_eff": lambda_eff,
"lambda": lambda_rate, "mu": mu_rate,
}
def mmc_analytics(lambda_rate, mu_rate, c):
"""
M/M/c: c серверов, формула Эрланга C.
C(c,a) = [a^c/c! * 1/(1-rho)] / [sum(a^k/k!, k=0..c-1) + a^c/c! * 1/(1-rho)]
Wq = C(c,a) / (c*mu - lambda)
"""
a = lambda_rate / mu_rate
rho = a / c
if rho >= 1.0:
raise ValueError(f"M/M/c requires rho < 1, got rho={rho:.4f}")
sum_terms = sum(a ** k / math.factorial(k) for k in range(c))
last_term = (a ** c / math.factorial(c)) * (1.0 / (1.0 - rho))
C_erlang = last_term / (sum_terms + last_term)
Wq = C_erlang / (c * mu_rate - lambda_rate)
W = Wq + 1.0 / mu_rate
L = lambda_rate * W
Lq = lambda_rate * Wq
return {
"rho": rho, "c": c, "C_erlang": C_erlang,
"Wq": Wq, "W": W, "L": L, "Lq": Lq,
"lambda": lambda_rate, "mu": mu_rate,
}
def tandem_analytics(lambda_rate, mu1_rate, mu2_rate):
"""
Тандем (сеть Джексона): две очереди последовательно.
По теореме Джексона каждая очередь — независимая M/M/1.
W_total = 1/(mu1-lambda) + 1/(mu2-lambda)
"""
q1 = mm1_analytics(lambda_rate, mu1_rate)
q2 = mm1_analytics(lambda_rate, mu2_rate)
return {
"rho1": q1["rho"], "rho2": q2["rho"],
"W1": q1["W"], "W2": q2["W"],
"W_total": q1["W"] + q2["W"],
"L_total": q1["L"] + q2["L"],
"lambda": lambda_rate, "mu1": mu1_rate, "mu2": mu2_rate,
}
def md1_analytics(lambda_rate, mu_rate):
"""
M/D/1: детерминированное обслуживание.
Формула Полачека-Хинчина с Var[S]=0:
Wq = rho / (2*mu*(1-rho)) -- ровно в 2 раза меньше чем M/M/1
"""
rho = lambda_rate / mu_rate
if rho >= 1.0:
raise ValueError(f"M/D/1 requires rho < 1, got rho={rho:.4f}")
Wq = rho / (2.0 * mu_rate * (1.0 - rho))
W = Wq + 1.0 / mu_rate
Wq_mm1 = rho / (mu_rate * (1.0 - rho))
return {
"rho": rho, "W": W, "Wq": Wq,
"Wq_mm1": Wq_mm1, "Wq_ratio": 0.5,
"lambda": lambda_rate, "mu": mu_rate,
}
# ================================================================
# ТЕСТ КОЛМОГОРОВА-СМИРНОВА
# ================================================================
def ks_test_exponential(samples, expected_rate):
"""
KS-тест: H0: samples ~ Exp(rate)
D_n = max|F_n(x) - F(x)|, F(x) = 1 - exp(-rate*x)
Критическое значение (alpha=0.05): D_crit = 1.358/sqrt(n)
PASS если D_n < D_crit
"""
n = len(samples)
if n == 0:
return {"D": float('inf'), "D_crit": 0, "passed": False, "n": 0}
sorted_samples = sorted(samples)
D = 0.0
for i, x in enumerate(sorted_samples):
F_x = 1.0 - math.exp(-expected_rate * x)
D = max(D, abs((i + 1) / n - F_x), abs(i / n - F_x))
D_crit = 1.358 / math.sqrt(n)
return {"D": D, "D_crit": D_crit, "passed": D < D_crit, "n": n}
# ================================================================
# ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ
# ================================================================
# t-критические значения для 95% CI (df = n-1)
T_CRITICAL = {
1: 12.706, 2: 4.303, 3: 3.182, 4: 2.776, 5: 2.571,
6: 2.447, 7: 2.365, 8: 2.306, 9: 2.262, 10: 2.228,
14: 2.145, 19: 2.093, 29: 2.045, 49: 2.010, 99: 1.984,
}
def compute_ci_95(values):
"""95% CI = x_bar +/- t(alpha/2, n-1) * s / sqrt(n)"""
n = len(values)
if n < 2:
return {"mean": values[0] if values else 0, "ci_low": float('-inf'),
"ci_high": float('inf'), "n": n}
mean = sum(values) / n
std = math.sqrt(sum((x - mean) ** 2 for x in values) / (n - 1))
t_crit = T_CRITICAL.get(n - 1, 1.96)
hw = t_crit * std / math.sqrt(n)
return {"mean": mean, "std": std, "ci_low": mean - hw,
"ci_high": mean + hw, "half_width": hw, "n": n}
def check_ci_contains(values, analytical_value, metric_name="metric"):
"""Проверка: содержит ли 95% CI аналитическое значение?"""
ci = compute_ci_95(values)
passed = ci["ci_low"] <= analytical_value <= ci["ci_high"]
return passed, {
"metric": metric_name, "analytical": analytical_value,
"mean": ci["mean"], "ci_low": ci["ci_low"], "ci_high": ci["ci_high"],
"passed": passed,
}
# ================================================================
# СРАВНЕНИЕ МЕТРИК
# ================================================================
SERVICE_MEAN = 100.0 # 1/mu = 100 секунд
MU = 1.0 / SERVICE_MEAN
def default_params(scenario, rho=0.7, K=5, c=3, rho2=0.5):
"""Параметры по умолчанию для каждого сценария."""
if scenario == "mm1":
return {"lambda": rho * MU, "mu": MU}
elif scenario == "mm1k":
return {"lambda": rho * MU, "mu": MU, "K": K}
elif scenario == "mmc":
return {"lambda": rho * c * MU, "mu": MU, "c": c}
elif scenario == "tandem":
lam = 0.5 * MU
return {"lambda": lam, "mu1": MU, "mu2": lam / rho2}
elif scenario == "md1":
return {"lambda": rho * MU, "mu": MU}
else:
raise ValueError(f"Unknown scenario: {scenario}")
def extract_metrics(result_json):
"""Извлечение метрик из JSON-ответа SimulationResultV2."""
ov = result_json.get("overview", {})
val = result_json.get("validation", {})
res = result_json.get("resources", {})
stats = res.get("statistics", [])
rho_obs = stats[0].get("occupancy", {}).get("mean", 0) if stats else 0
total = ov.get("totalCreated", 0)
dropped = ov.get("totalDropped", 0)
return {
"W": ov.get("avgTimeInSystemSeconds", 0),
"Wq": ov.get("avgTotalWaitTimeSeconds", 0),
"ES": ov.get("avgProcessingTimeSeconds", 0),
"rho": rho_obs,
"P_reject": dropped / total if total > 0 else 0,
"total_created": total,
"total_processed": ov.get("totalProcessed", 0),
}
def compare_metric(name, observed, expected, tolerance=0.10):
"""Сравнение одной метрики с допуском."""
if expected == 0:
return abs(observed) < tolerance, {"metric": name, "obs": observed,
"exp": expected, "passed": abs(observed) < tolerance}
rel_err = abs(observed - expected) / abs(expected)
passed = rel_err <= tolerance
return passed, {
"metric": name, "observed": observed, "expected": expected,
"rel_error_pct": rel_err * 100, "tolerance_pct": tolerance * 100,
"passed": passed,
}
def verify_scenario(scenario, result_json, params, tolerance=0.10):
"""Верификация одного сценария."""
m = extract_metrics(result_json)
checks = []
if scenario == "mm1":
exp = mm1_analytics(params["lambda"], params["mu"])
checks.append(compare_metric("W", m["W"], exp["W"], tolerance))
checks.append(compare_metric("Wq", m["Wq"], exp["Wq"], tolerance))
checks.append(compare_metric("rho", m["rho"], exp["rho"], tolerance))
elif scenario == "mm1k":
exp = mm1k_analytics(params["lambda"], params["mu"], params["K"])
checks.append(compare_metric("W", m["W"], exp["W"], tolerance))
checks.append(compare_metric("P_reject", m["P_reject"], exp["P_reject"], tolerance))
elif scenario == "mmc":
exp = mmc_analytics(params["lambda"], params["mu"], params["c"])
checks.append(compare_metric("W", m["W"], exp["W"], tolerance))
checks.append(compare_metric("Wq", m["Wq"], exp["Wq"], tolerance))
elif scenario == "tandem":
exp = tandem_analytics(params["lambda"], params["mu1"], params["mu2"])
checks.append(compare_metric("W_total", m["W"], exp["W_total"], tolerance))
elif scenario == "md1":
exp = md1_analytics(params["lambda"], params["mu"])
checks.append(compare_metric("W", m["W"], exp["W"], tolerance))
checks.append(compare_metric("Wq", m["Wq"], exp["Wq"], tolerance))
all_passed = all(c[0] for c in checks)
return all_passed, [c[1] for c in checks]
# ================================================================
# ВЫВОД РЕЗУЛЬТАТОВ
# ================================================================
def print_results(scenario, params, checks):
"""Печать таблицы результатов."""
print(f"\n{'='*60}")
print(f" Сценарий: {scenario.upper()}")
print(f" Параметры: {params}")
print(f"{'='*60}")
print(f" {'Метрика':<20} {'Наблюд.':>10} {'Ожидаем.':>10} {'Ошибка%':>8} {'Статус':>8}")
print(f" {'-'*20} {'-'*10} {'-'*10} {'-'*8} {'-'*8}")
for c in checks:
status = "PASS" if c["passed"] else "FAIL"
print(f" {c['metric']:<20} {c.get('observed',0):>10.4f} "
f"{c.get('expected',0):>10.4f} "
f"{c.get('rel_error_pct',0):>7.1f}% {status:>8}")
all_ok = all(c["passed"] for c in checks)
print(f"\n Итог: {'ВСЕ ПРОВЕРКИ ПРОЙДЕНЫ' if all_ok else 'ЕСТЬ НЕПРОЙДЕННЫЕ'}")
# ================================================================
# ТОЧКА ВХОДА
# ================================================================
def main():
args = sys.argv[1:]
if not args or args[0] in ("-h", "--help"):
print(__doc__)
return
if args[0] == "--compute-only":
scenario = args[1] if len(args) > 1 else "mm1"
rho = 0.7
i = 2
while i < len(args):
if args[i] == "--rho":
rho = float(args[i + 1]); i += 2
else:
i += 1
params = default_params(scenario, rho=rho)
if scenario == "mm1":
a = mm1_analytics(params["lambda"], params["mu"])
elif scenario == "mmc":
a = mmc_analytics(params["lambda"], params["mu"], params["c"])
elif scenario == "md1":
a = md1_analytics(params["lambda"], params["mu"])
else:
a = params
print(f"Аналитические значения для {scenario.upper()}:")
for k, v in a.items():
if isinstance(v, float):
print(f" {k} = {v:.6f}")
else:
print(f" {k} = {v}")
else:
scenario = args[0]
json_path = args[1] if len(args) > 1 else f"scenario_{scenario}_output.json"
with open(json_path, 'r', encoding='utf-8') as f:
data = json.load(f)
params = data.get("params", default_params(scenario))
result = data.get("result", data)
passed, checks = verify_scenario(scenario, result, params)
print_results(scenario, params, checks)
sys.exit(0 if passed else 1)
if __name__ == "__main__":
main()
Пример использования
1. Вычисление аналитических значений
$ python verify_des.py --compute-only mm1 --rho 0.7
Аналитические значения для MM1:
rho = 0.700000
L = 2.333333 (среднее число в системе)
W = 333.333333 (среднее время в системе, сек)
Wq = 233.333333 (среднее время ожидания, сек)
2. Проверка результата DES
$ python verify_des.py mm1 scenario_mm1_output.json
============================================================
Сценарий: MM1
Параметры: {'lambda': 0.007, 'mu': 0.01, 'rho': 0.7}
============================================================
Метрика Наблюд. Ожидаем. Ошибка% Статус
-------------------- ---------- ---------- -------- --------
W 341.2345 333.3333 2.4% PASS
Wq 240.1234 233.3333 2.9% PASS
rho 0.6983 0.7000 0.2% PASS
Итог: ВСЕ ПРОВЕРКИ ПРОЙДЕНЫ
3. Сравнение M/M/1 vs M/D/1
Один из самых наглядных результатов — сравнение M/M/1 и M/D/1 при одинаковой нагрузке:
$ python verify_des.py --compute-only mm1 --rho 0.7
Wq = 233.33s
$ python verify_des.py --compute-only md1 --rho 0.7
Wq = 116.67s (ровно в 2 раза меньше!)
Это подтверждает формулу Полачека-Хинчина: устранение вариативности обслуживания вдвое сокращает ожидание. Если ваш DES-движок это не показывает — в нём ошибка.
Формат JSON-вывода
Скрипт ожидает JSON-файл с результатом симуляции Storm BPM. Ключевые поля:
{
"params": {
"lambda": 0.007,
"mu": 0.01,
"rho": 0.7
},
"result": {
"overview": {
"totalCreated": 3500,
"totalProcessed": 3495,
"totalDropped": 0,
"avgTimeInSystemSeconds": 341.23,
"avgTotalWaitTimeSeconds": 240.12,
"avgProcessingTimeSeconds": 101.11,
"virtualTimeSeconds": 500000
},
"resources": {
"statistics": [{
"occupancy": { "mean": 0.698 }
}]
},
"validation": {
"littlesLaw": {
"observedL": 2.35,
"calculatedL": 2.39,
"isValid": true
}
}
}
}
Хотите проверить сами?
Мы с радостью поможем вам запустить верификацию на ваших данных и разобрать результаты. Обратитесь к менеджеру, с которым вы работаете — мы предоставим JSON-выводы симуляции и поможем прогнать скрипт шаг за шагом.
Что мы можем сделать вместе
- Прогнать все 5 сценариев на вашем инстансе и сравнить с формулами
- Настроить параметры (нагрузку, число серверов, ёмкость буфера) под ваш конкретный кейс
- Объяснить каждую формулу и как она соотносится с вашим бизнес-процессом
- Провести воркшоп по интерпретации результатов симуляции
Литература
- Gross, Shortle et al. "Fundamentals of Queueing Theory" 5th ed. — базовый учебник, главы 2-6.
- Kleinrock, L. "Queueing Systems, Vol. 1: Theory" (Wiley) — классика, формулы M/M/1, M/G/1.
- Law, A.M. "Simulation Modeling and Analysis" 5th ed. (McGraw-Hill) — методология V&V симуляций.
- Jackson, J.R. "Networks of Waiting Lines" (1957) — теорема Джексона для сетей очередей.