ЛР-03: Муниципальная готовность к ЧС

Worked example: civil 03

Это полностью разобранный example по duality и анализу чувствительности. Он показывает полный цикл от прямой модели до сценарных пересчётов.

1. Исходный кейс

Сценарий с ресурсами готовности к чрезвычайным ситуациям.

Ограничения ресурсов

Ресурс	Лимит
Бюджет	94
Трудозатраты	56
Операционная ёмкость	48

Программы

Программа	Эффект	Бюджет	Трудозатраты	Операционная ёмкость
Резерв мобильных генераторов	90	42	17	16
Запас питьевой воды	72	20	9	8
Обучение волонтёров	64	16	19	7
Полевые пункты размещения	98	52	26	23

import numpy as np
import pandas as pd
from scipy.optimize import linprog

def solve_primal(effects, A_ub, b_ub):
    result = linprog(-effects, A_ub=A_ub, b_ub=b_ub, bounds=[(0, 1)] * len(effects), method='highs')
    if not result.success:
        raise RuntimeError(result.message)
    shadow_prices = -result.ineqlin.marginals
    return result, shadow_prices

def solve_dual(effects, A_ub, b_ub):
    m, n = A_ub.shape
    c_dual = np.concatenate([b_ub, np.ones(n)])
    A_dual = -np.hstack([A_ub.T, np.eye(n)])
    b_dual = -effects
    result = linprog(c_dual, A_ub=A_dual, b_ub=b_dual, bounds=[(0, None)] * (m + n), method='highs')
    if not result.success:
        raise RuntimeError(result.message)
    return result

def rerun_with_resource_change(effects, A_ub, b_ub, resource_index, delta):
    new_b = b_ub.copy()
    new_b[resource_index] += delta
    result, _ = solve_primal(effects, A_ub, new_b)
    return new_b, result

program_names = ['Резерв мобильных генераторов', 'Запас питьевой воды', 'Обучение волонтёров', 'Полевые пункты размещения']
resource_names = ['Бюджет', 'Трудозатраты', 'Операционная ёмкость']
effects = np.array([90, 72, 64, 98], dtype=float)
A_ub = np.array([[42, 20, 16, 52], [17, 9, 19, 26], [16, 8, 7, 23]], dtype=float)
b_ub = np.array([94, 56, 48], dtype=float)

primal_result, shadow_prices = solve_primal(effects, A_ub, b_ub)
dual_result = solve_dual(effects, A_ub, b_ub)

plan_df = pd.DataFrame({
    'программа': program_names,
    'x*': np.round(primal_result.x, 4),
    'эффект на единицу': effects,
})
resources_df = pd.DataFrame({
    'ресурс': resource_names,
    'лимит': b_ub,
    'slack': np.round(primal_result.slack, 4),
    'shadow_price': np.round(shadow_prices, 4),
    'binding': np.isclose(primal_result.slack, 0.0),
})

print('Оптимальный эффект (primal):', round(-primal_result.fun, 4))
print('Оптимальное значение dual:', round(dual_result.fun, 4))
print()
print('Оптимальный план:')
display(plan_df)
print('Ресурсный разбор:')
display(resources_df)

Оптимальный эффект (primal): 256.1538
Оптимальное значение dual: 256.1538

Оптимальный план:

	программа	x*	эффект на единицу
0	Резерв мобильных генераторов	1.0000	90.0
1	Запас питьевой воды	1.0000	72.0
2	Обучение волонтёров	1.0000	64.0
3	Полевые пункты размещения	0.3077	98.0

Ресурсный разбор:

	ресурс	лимит	slack	shadow_price	binding
0	Бюджет	94.0	0.0000	1.8846	True
1	Трудозатраты	56.0	3.0000	0.0000	False
2	Операционная ёмкость	48.0	9.9231	0.0000	False

scenario_rows = []
for label, resource_index, delta in [
    ('Бюджет +3', 0, 3),
    ('Трудозатраты +2', 1, 2),
]:
    new_b, new_result = rerun_with_resource_change(effects, A_ub, b_ub, resource_index, delta)
    predicted = shadow_prices[resource_index] * delta
    actual = (-new_result.fun) - (-primal_result.fun)
    scenario_rows.append({
        'сценарий': label,
        'ресурс': resource_names[resource_index],
        'delta': delta,
        'прогноз по shadow price': round(predicted, 4),
        'факт после пересчёта': round(actual, 4),
        'разница': round(actual - predicted, 4),
    })

scenario_df = pd.DataFrame(scenario_rows)
display(scenario_df)

	сценарий	ресурс	delta	прогноз по shadow price	факт после пересчёта	разница
0	Бюджет +3	Бюджет	3	5.6538	5.6538	0.0
1	Трудозатраты +2	Трудозатраты	2	0.0000	0.0000	0.0

objective_label = 'Волонтёры +8 к эффекту'
objective_index = 2
objective_delta = 8

new_effects = effects.copy()
new_effects[objective_index] += objective_delta
objective_result, _ = solve_primal(new_effects, A_ub, b_ub)
objective_df = pd.DataFrame({
    'сценарий': [objective_label],
    'новый оптимальный эффект': [round(-objective_result.fun, 4)],
    'изменение эффекта': [round((-objective_result.fun) - (-primal_result.fun), 4)],
    'новое значение программы': [round(objective_result.x[objective_index], 4)],
})
display(objective_df)

	сценарий	новый оптимальный эффект	изменение эффекта	новое значение программы
0	Волонтёры +8 к эффекту	264.1538	8.0	1.0

2. Что важно проговорить в выводе

• какие ограничения оказались binding и почему именно они держат оптимум;

• какой ресурс имеет наибольшую теневую цену и что это означает содержательно;

• насколько хорошо совпал прогноз по shadow price с фактическим пересчётом;

• меняется ли структура плана при небольших изменениях b и c.