Робастная оптимизация: компромисс оптимальности и валидности решения

В комментариях к предыдущим публикациям поступали запросы рассмотреть формулировки оптимизационных задач и подходы к моделированию, которые позволяют учесть волатильность данных. В этой публикации рассмотрю один из вариантов учета неопределенности в данных в концепции робастной оптимизации.

Рассуждение будет построено на основе классической задачи Диеты Стиглера, добавим немного неопределенности и рассмотрим, как с ней бороться. Обратим внимание на два противоборствующих фактора: затраты на диету и степень удовлетворенности ограничений при различных сценариях отдельно.

В чем проблема?

Неопределенность. Интуитивно человек/сотрудник приходит к своему механизму работы с неизвестностью. То, насколько он справляется с этой задачей, по большей части определяет его успех. Трейдинг — рынок акций/валюты/крипты — крайне волатилен и трудно предсказуем. Но есть люди/компании, которые учитывают эту неопределенность лучше остальных.

Что насчет алгоритмов? Рассмотрим некоторые варианты учета неопределенности в теории прикладного мат.моделирования и оптимизации.

Худший сценарий (worst-case). Еще Дейл Карнеги писал, чтобы избавиться от тревожности, нужно представить худший исход и принять его. Как ни странно, этот подход присутствует в литературе по исследованию операций (прикладная математика). Подход предполагает использовать для оценки наихудший сценарий, например, закладывать минимальную цену продажи товара или минимальный ожидаемый спрос.

С другой стороны, можно сказать, что стакан наполовину полон. Тогда в игру вступает наилучший сценарий. Он позволяет получить оценку сверху.

Имея на руках наихудший, наилучший и еще какие-нибудь промежуточные сценарии, появляется возможность сформировать средний или средне-взвешенный сценарий, который будет отражать некоторую ожидаемую картину. В то же время, ожидаемая картина может быть далека от каждого отдельно взятого сценария.

Наихудший, наилучший или ожидаемый сценарий освещают одну сторону проблемы. Результат решения задачи оптимизации, построенного на таких данных, будет генерировать аналогичные однобокие картины. Такие решения несут в себе некоторую пользу, но по большей части они неадекватно моделируют происходящий процесс.

В попытках приблизиться к реальности разработаны подходы, которые позволяют учитывать неопределенность непосредственно внутри оптимизационной модели. Речь пойдет о робастной оптимизации.

Отмечу, что моделировать неопределенность будем посредством набора дискретных сценариев с учетом вероятности их реализации. Анализ проведу на основе известной академической задачи — диета Стиглера.

Постановка задачи: Диета Стиглера

Несколько слов о задаче.

Старая — 1944г.;
Одна из первых задач, к решению которой применялся симплекс метод;
Находится в портфеле базовых Examples большинства LP-солверов.

Если считаете свою работу рутинной, после следующего факта, возможно, переосмыслите это:

В 1947 году Джек Ладерман использовал симплекс-метод (тогда еще недавнее изобретение!) для определения оптимального решения. Потребовалось 120 человеко-дней девяти клерков на настольных калькуляторах, чтобы получить ответ.

О чем. В 1943 году Национальным исследовательским советом США были предложены диетические нормы (RDA) для «активного» мужчины весом 70 кг. С другой стороны, Джордж Стиглер задался вопросом о минимальной стоимости продуктов, которые позволят удовлетворить предложенные нормы (проблема студента). Рассматривался список из 77 продуктов питания и 9 питательных веществ (с входными данными можно ознакомиться здесь). В данных приведено содержание питательных веществ в продуктах на 1$.

Математическая постановка задачи

$i \in I$ — разнообразие доступных продуктов ();

$j \in J$ — множество питательных веществ ();

$a_{ij}$ — содержание питательного вещества в единице продукта ;

— нормативная суточная доза питательного вещества ;

$x_{i}$ — вещественная переменная, затраты на покупку продукта ;

Ограничения:

Удовлетворение суточной потребности в каждом из питательных веществ:

$\sum_{i \in I} a_{ij} x_{i} \ge b_j, \quad \forall j \in J.$

Целевая функция:

Минимизация стоимости корзины продуктов:

$\min \sum_{i \in I}x_{i}$

Программную реализацию модели выполним на python в OR-Tools и решателем GLOP. Детали реализации опущу, т.к. более подробно разбирал построение и оптимизацию в других публикациях (задача о назначениях, задача планирования расписаний).

import pandas as pd import numpy as np from ortools.linear_solver import pywraplp  # Загрузка данных по нормативам и продуктам df_nutrients_norm = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/03_diet_nutrients_limits.csv", sep=";") df_products = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/03_diet_products.csv", sep=";", encoding="cp1251")  # Реализация модели # Инициализируем солвер solver = pywraplp.Solver.CreateSolver('GLOP')  # Бесконечность ) infinity = solver.infinity()  # Инициализация переменных в DataFrame df_products["var_id"] = np.ogrid[:df_products.shape[0]] df_products["var"] = df_products["var_id"].apply(lambda x: solver.NumVar(0.0, infinity, f"x_{x}"))  # Ограничение: Удовлетворение суточной потребности в каждом из питательных веществ for nutrient in df_nutrients_norm.itertuples():   var_sum = sum(df_products[nutrient.name] * df_products["var"])   solver.Add(var_sum >= nutrient.min)  # Целевая функция: Минимизация стоимости корзины продуктов solver.Minimize(sum(df_products["var"]))  # Решаем status = solver.Solve()  # Извлекаем результат if status == pywraplp.Solver.OPTIMAL:   daily_cost = solver.Objective().Value()   print(f'Затраты на продукты в день = {round(daily_cost, 3)}$')   print(f'Затраты на продукты в год = {round(daily_cost * 365, 3)}$')  df_products["sol"] = df_products["var"].apply(lambda x: x.solution_value())  # Обработка результата filt_diet = df_products["sol"] > 0 df_diet = df_products[filt_diet].copy()

В ценах 1939 года минимальная стоимость продуктов, которая удовлетворяет минимальным нормам питания (1940-ых гг.), составляет 0.109$ в день или 39.662$ в год. Эти результаты похожи на результаты, полученные другими исследователями (например, в этой работе). Что касается продуктов и совокупных затрат на их приобретение, получили следующее:

Commodity	Commodity (rus)	Daily, $	Annual, $
Wheat Flour (Enriched)	Пшеничная мука (обогащенная)	0.03	10.774
Liver (Beef)	Печень (Говядина)	0.002	0.691
Cabbage	Капуста	0.011	4.093
Spinach	Шпинат	0.005	1.828
Navy Beans, Dried	Фасоль темно-синяя, сушеная	0.061	22.275
	Total	0.109	39.662

Пункт о разнообразии питания в модель не закладывали, предполагается всю жизнь на сушеной фасоли сидеть. По продуктам есть одно расхождение. Наша модель выбрала «Печень», в то время как у других получилось «Сгущенное молоко». О предположениях возникновения такого расхождения можете написать в комментах.

Перейдем к целевой части статьи — неопределенность.

Добавим немного неопределенности

Концепцию робастной оптимизации и схему демонстрации на примере задачи Диета Стиглера позаимствую из этой работы.

Будем рассматривать два вектора робастности:

Решение, «близкое» к оптимальному решению каждого отдельно взятого сценария, будем называть робастным решением (solution robust).
Решение, которое «почти» для любого отдельно взятого сценария является допустимым, будем называть робастным на уровне модели (model robust).

Откуда берется неопределенность? В магазинах на полках присутствует несколько сортов яблок, например: «Гало», «Сезонные», «Гренни Смит», «Гольден», «Белый налив» и др. Каждый из этих сортов имеет свой состав, даже внутри одного сорта состав питательных веществ может отличаться (собран в начале лета или конце; рос на солнце или в тени). Какой состав питательных веществ яблока считать эталонным? С одной стороны, мы можем расширить список продуктов всеми сортами яблок, что не решит проблему с волатильностью веществ внутри сорта. Или создать несколько наиболее вероятных сценариев состава яблок и принимать решение на основе нескольких исходов в данных.

Расширяем условие задачи. В нашем примере с Диетой Стиглера мы видим доминирование сушеной фасоли в решении. Кроме того, сушеная фасоль вносит значительный вклад в восполнении кальция. В целях простой интерпретации результатов рассмотрим неопределенность содержания кальция в темно-синей фасоли.

Пусть имеем 13 равновероятных сценариев $\Omega = \{S_1, \dots, S_8, \dots, S_{13} \}=$ $\{9.5, 9.75, 10, 10.25, 10.5, 11, 11.25, 11.4, 11.5, 11.75, 12, 12.25, 12.5\}$ содержания кальция в сушеной фасоли за доллар. Отмечу, так как состав продуктов приведен на 1$, то сценарии, в том числе, отражают неопределенность цены (не только состав).

Если рассмотреть наихудший сценарий , тогда получим решение, удовлетворяющее всем ограничениям для любого из сценариев. Но будет ли оно самым экономным/оптимальным для каждого из остальных сценариев в отдельности? Что если мы допустим незначительные нарушения по некоторым из сценариев, но при этом сформируем более оптимальную диету, т.е. повышаем робастность решения за счет снижения робастности модели. Одна из реализаций такого компромисса следующая:

$s \in S$ — множество рассматриваемых сценариев;

$c \in J$ — индекс кальция в списке питательных веществ;

$N \in I$ — индекс сушеной фасоли в списке продуктов;

— вероятность реализации сценария ;

$a^s_{Nc}$ — содержание кальция в сушеной фасоли по сценарию ;

$e^c_{s}$ — вещественная переменная, размер дефицита содержания кальция в диете по сценарию ;

Ограничения:

1. Удовлетворение суточной потребности в каждом из питательных веществ, кроме кальция:

$\sum_{i \in I} a_{ij} x_{i} \ge b_j, \quad \forall j \in J: j \ne c;$

2. Оценка дефицита кальция по каждому из сценариев:

$\sum_{i \in I: i\ne N} a_{ic} x_i + a^s_{Nc} x_{N} + e^c_{s} \ge b_с, \quad \forall s \in S;$

3. Выполнение минимального содержания кальция, как минимум, для наилучшего сценария:

$\sum_{i \in I: i\ne N} a_{ic} x_i + a^{13}_{Nc} x_{N} \ge b_с.$

Целевая функция:

Минимизация стоимости корзины продуктов в компромиссе с ожидаемым дефицитом кальция в диете:

$\min \sum_{i \in I}x_{i} + w\sum_{s \in S} p_s e^c_{s}.$

Вес используется для балансирования размера ожидаемого дефицита и затрат. Если положить , то будем решать задачу наилучшего сценария. Если присвоить очень большое значение (сдвигаем акцент на выполнение ограничений), то будем решать задачу наихудшего сценария. Выбор данного параметра зависит, в первую очередь, от экспертной оценки и целей использования модели.

Ограничения типа (2) выполняют функцию оценки отклонения содержания кальция в диете от норматива по каждому из сценариев:

$e^c_{s} = \max(0, b_с - \sum_{i \in I: i\ne N} a_{ic} x_i + a^s_{Nc} x_{N}).$

В свою очередь, важным является ограничение (3). Если его не добавлять в модель и положить , то модель может безнаказанно исключить кальций из рассмотрения. Во-первых, дефицит не будет штрафоваться при достаточно малых , во-вторых, экономия на затратах за счет исключения питательного вещества из рассмотрения. Добавление ограничения (3) делает снова актуальным учет норм по кальцию.

Почему сценарий $s_{13}$ наилучший? Это сценарий, когда сушеная фасоль содержит больше всего кальция в сравнении с другими сценариями. Следовательно, меньшим кол-вом фасоли можно закрыть потребность в кальции, а это приводит к снижению стоимости. Наилучший сценарий берется для того, чтобы была возможность ухудшать решение по затратам и улучшать по валидности (удовлетворению нормам) в разрезе каждого из сценариев в отдельности.

import pandas as pd import numpy as np from ortools.linear_solver import pywraplp  # Загрузка данных по нормативам и продуктам df_nutrients_norm = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/03_diet_nutrients_limits.csv", sep=";") df_products = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/03_diet_products.csv", sep=";", encoding="cp1251")  # Сценарии содержания кальция в сушеной фасоли lst_scenarios = [9.5, 9.75, 10, 10.25, 10.5, 11, 11.25, 11.4, 11.5, 11.75, 12, 12.25, 12.5] dct_scenarios = {idx + 1: val for idx, val in  enumerate(lst_scenarios)}  # Уровень компромиса w = 0.0498 s_cnt = len(dct_scenarios)  # Инициализируем солвер solver = pywraplp.Solver.CreateSolver('GLOP')  # Бесконечность ) infinity = solver.infinity()  # Инициализация переменных в DataFrame df_products["var_id"] = np.ogrid[:df_products.shape[0]] df_products["var"] = df_products["var_id"].apply(lambda x: solver.NumVar(0.0, infinity, f"x_{x}"))  # Переменные оценки нарушений dct_slack_vars = {i: solver.NumVar(0.0, infinity, f"e_{i}") for i in dct_scenarios}  # Ограничение: Удовлетворение суточной потребности в каждом из питательных веществ for nutrient in df_nutrients_norm.itertuples():   if nutrient.name == "Calcium (g)":     # Фильтруем фасоль     filt_navie = df_products["Commodity (rus)"] == 'Фасоль темно-синяя, сушеная'     # Ограничения для каждого сценария     for s_id, volume in dct_scenarios.items():       # Подменяем значение из сценария       df_products["Calcium (g)"] = np.where(filt_navie, volume, df_products["Calcium (g)"])        var_sum = sum(df_products[nutrient.name] * df_products["var"])       solver.Add(var_sum + dct_slack_vars[s_id] >= nutrient.min)     solver.Add(var_sum >= nutrient.min)   else:     var_sum = sum(df_products[nutrient.name] * df_products["var"])     solver.Add(var_sum >= nutrient.min)  # Целевая функция: Минимизация стоимости корзины продуктов в компромиссе # с ожидаемым нарушением ограничения по содержанию кальция в диете solver.Minimize(sum(df_products["var"]) + w / s_cnt * sum(dct_slack_vars.values()))  # Решаем status = solver.Solve()  # Извлекаем результат if status == pywraplp.Solver.OPTIMAL:   df_products["sol"] = df_products["var"].apply(lambda x: x.solution_value())   daily_cost = df_products["sol"].sum()   print(f'Затраты на продукты в день = {round(daily_cost, 3)}$')   print(f'Затраты на продукты в год = {round(daily_cost * 365, 3)}$')    # Суммарные нарушения   dct_slack_val = {key: var.solution_value() for key, var in dct_slack_vars.items()}   total_slack = sum(dct_slack_val.values())   print(f'Ожидаемый дефицит = {total_slack / s_cnt}g')

Обратимся к результатам эксперимента с подставлением различных значений . На графиках ниже представлена ошибка или дефицит кальция по соответствующему сценарию для определенного уровня .

Видим, что при увеличении снижается кол-во случаев дефицита кальция — повышается робастность модели. При только у сценария $S_{13}$ нет нарушений и этот случай имеет наименьшие затраты — робастное решение.

Рассмотрим поведение общих затрат и размер дефицита кальция в зависимости от . Разброс по затратам составляет от 0.107$ до 0.113$ в день (~5%), ожидаемый дефицит кальция от 0г до 0.08г (~10% от нормы кальция 0.8г).

Как выбрать ? У измерительных приборов есть погрешность при замерах с которой приходится мириться. Аналогичной методикой можно воспользоваться и здесь, предварительно оценивая границы допустимого ожидаемого нарушения. Например, ограничить ожидаемую ошибку 0.02г (составляет 2,5% от нормы), это будет соответствовать . При таких допущениях, в бюджет на питание можно закладывать 0.1093$ в день (если вернуться в тот самый 1939г). Немного дороже, чем в исходной линейной постановке, зато надежно (робастно).

Заключение

Разобрали принцип работы робастной оптимизации на урезанном кейсе Диеты Стиглера для случая неопределенного содержания кальция в сушеной фасоли. Несложными манипуляциями можно масштабировать подход на все питательные вещества и на все продукты, а так же на другие задачи. Отмечу, что это один из подходов к учету и обработке неопределенности данных.

Ссылки

Справка по задаче Диета Стиглера;
Версия статьи в Jupyter notebook.

Предыдущие материалы:

Модель назначения машин такси на заявки;
Модель планирования расписаний сотрудников;
Задача мостов Кенигсберга.

ссылка на оригинал статьи https://habr.com/ru/articles/751226/