# Другие рекомендации (Прогнозируй и рекомендуй). Часть 2
## Интро
Уже достаточно продолжительное время различные команды занимаются разработкой и внедрением алгоритмов, "под капотом" которых находятся модели машинного обучения и оптимизационные алгоритмы, в промышленности: черная и цветная металлургия, нефтяная отрасль и др.

Не буду сравнивать сложность внедрение решений на базе ML в промышленности и, например, в ритейле или банкинге - об этом можно послушать и почитать в других материалах. Но, очевидно, что любое промышленное производство - это очень сложная система со всех точек зрения: технология, IT, оборудование, логистика и пр.

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

**Глоссарий**
- плавка
- конвертер
- отдача
- сталевар
- доводка

**Описание процесса (поверхностно)**
В конвертерном цехе на сталелитейном производстве выполняется часть обработки плавки следующих агрегатов: конвертер, установка печь ковш (УПК) и установка ковшевого вакуумирования (УКВ). В плавку, состоящую из чугуна, лома, извести и флюсов добавляют ферросплавы для получения стали с требуемым содержанием химических элементов: $Si, Mn, C, N, Al$ и др. Разные требования по химическому составу соответствуют разным маркам стали, т.е. в некоторых видах изделий требуется высокое содержание $Si$, а в некоторых - нет.
!добавить ссылку на описание процесса

**Пролематика**
Стоимость ферросплавов высокая, в процессе производства стали в конвертерном цехе требуется минимизировать расходы на ферросплавы при достижении требуемого химического состава.
Избыточное использование ферросплавов может быть вызвано:
- сложностью процесса и, как следствие, проблематичностью учета всевозможных факторов;
- опытом сталеваров;
- требованиями по производительности.

**Задача**
Сформировать рекомендации об объемах отдачи различных ферросплавов на разных этапах/агрегатах конвертерного цеха.

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

## Моделирование
**Исходные данные**
Чуть подробнее рассмотрим ключевые данные:
- добавляемые материалы: $ferro_i$ -- **управляемые параметры**, ранее обозначенные как $u_i, i=\overline{0, m-1}$
Например: $ferro_0$ - $FeMn\; 78\%$
- технологические параметры: $feature_i$ -- **фиксированные параметры**, ранее обозначенные как $x_i, i=\overline{0, n-1}$.
Например: $feature_2$ -- температура плавки на агрегате, например, на УКП
- целевые параметры: $target_i$
- $target_1$ - изменение между текущим и предыдущим химическим анализ (далее ХА) по $Mn$
- $target_2$ - изменение между текущим и предыдущим ХА по $Si$
- дополнительные параметры
- $target\_min\_bound_i, i\in\{1, 2\}$ - нижняя граница по требованиями ХА для $Mn$ и $Si$
- $target\_max\_bound_i, i\in\{1, 2\}$ - верхняя граница по требованиями ХА для $Mn$ и $Si$
- $initial\_val\_target_1$ - предыдущий ХА плавки по $Mn$
- $initial\_val\_target_2$ - предыдущий ХА плавки по $Si$

все вышеприведенные параметры имеют размерность $\mathrm{R}^n$, где $n$ -- количество наблюдений в данных.

**Прогнозные модели**
Строим следующие прогнозные модели
\begin{equation}
pred\_target_i^k = f_i(feature_0^k, \dots, feature_{n-1}^k, ferro_0^k, \dots, ferro_{m-1}^k), \; i \in \{1, 2\},
\end{equation} где $k$ -- номер наблюдения, $k=\overline{1, n},\; pred\_target_i^k\in \mathrm{R}$ -- прогноз целевой переменной для $k$-го наблюдения при заданных фиксированных и управляемых параметрах.

**Оптимизационная задача**
Перед нами стоит задача минимизировать затраты за счет более экономного использования ферросплавов, но при этом ни в коем случае это не должно влиять на качество изготавливаемой продукции. Следовательно, из этого вытекает следующая критериальная функция: \begin{equation}
\sum_{i=\overline{0, m-1}}c_i\cdot ferro_i^k \to \min\limits_{ferro_i^k \in U}
\end{equation} где $c_i$ - цена за 1 тонну $i$-го ферросплава. При минимизации стоимости плавки требуется соблюдать требования по химии, т.е. могут быть наложены следующие ограничения на $U$
\begin{gather*}
target\_min\_bound_i \leq pred\_target_i^k \leq target\_max\_bound_i, \; i \in \{1, 2\}, \\ g_j(pred\_target_1^k,\; pred\_target_2^k) \leq 0, \; j=\overline{1, l}
\end{gather*} через ограничения $g_j$ могут устанавливаться различные дополнительные технологические ограничения, например, ограничения на соотношение между ключевыми элементами или более узкие диапазоны для $pred\_target_i^k$, для того, чтобы попадать в определенную цель и стабилизировать процесс доводки стали.

Альтернативным примером критериальной функции может быть
\begin{equation}
|pred\_target_i^k - goal\_target_i^k| \to \min\limits_{ferro_i \in U},
\end{equation}где $goal\_target_i^k$ является некоторой целью по содержанию химического элемента, в которую необходимо попасть. Чаще всего люди с производства имеют регламент с заданными целями или владеют этим ценным знанием в своей голове.

**Заметка**
- Критериальные функции можно комбинировать, но нужно быть аккуратным, чтобы не складывать тонны с рублями;
- Выше приведена одна из возможных постановок задачи. Может меняться практически каждый компонент в зависимости от ваших предпочтений и требований людей, которые в дальнейшем будут эксплуатировать систему;
- Далее по ходу статьи постановка задачи будет незначительно модифицироваться в зависимости от подхода к решению.

## Подходы

При решении задачи оптимизации мы используем внутреннюю библиотеку, которая поддерживает различные методы и постановки задачи. Для демонстрации части подходов будут использоваться методы из библиотеки.

### Перебор по сетке
В данном подходе используется подход декомпозиции исходной задачи, который в данном случае заключается в следующем:
1. Строится модель для прогноза $Mn$ на части ферросплавов, которые используются для доводки $Mn$
2. Формируется рекомендация ферросплавов для доводки $Mn$
3. Строится модель для прогноза $Si$ на ферросплавах, которые используются для доводки $Si$ и $Mn$
4. Формируется рекомендация ферросплавов для доводки $Si$, но с учетом того, что в п.2 при рекомендации ферросплавы также содержат $Si$ необходимо учесть, сколько из уже добавленных ферросплавов будет прирост по $Si$.

Для решения задачи необходимо определить следующие пункты:
1. Сетка для поиска

search_space = {
    'ferro_0': np.linspace(ferro_0_min, ferro_0_max, 100),
    'ferro_1': np.linspace(ferro_1_min, ferro_1_max, 100)
}

2. Критериальная функция - абсолютное отклонение от цели
(здесь и далее префикс $gs$ расшифровывается как $GridSearch$)

def gs_opt_obj(x, models, target_point):
    # Target variable prediction by model
    target_pred = models['target_1']['model'].predict(x[models['target_1']['features']].values)
    
    # Diff between true and predicted target variable value
    target_loss = abs(target_point['target_1_point'] - target_pred)
        
    return target_loss

3. Прогнозная модель - градиентный бустинг

gs_optim = gs_opt.GSOpt(
    models= {
        'target_1': {
            'model': gb_models['target_1']['model'],
            'features': gb_models['target_1']['features']
        }
    }
)

4. Запуск поиска по сетке

# set goal_target
gs_obj_params = [{'target_point': data.loc[idx, ['target_1_point']].to_dict()}]

opt_res, opt_hist=gs_optim.run_gs_opt(
    x=data.loc[idx], # set point for optimization
    objective_params=gs_obj_params, # set params 
    objectives=[gs_opt_obj], # set objective
    search_space=search_space, # set search space
    opt_direction=[0, 'minimize'] # set direction
)

После выполнения поиска результат выглядит так:

ferro_0                       0.378000
ferro_1                       0.000000
ferro_2                       0.000000
feature_0                     0.180000
feature_1                    -0.007000
feature_2                  1560.000000
feature_3                    70.100000
feature_4                   295.763000
feature_5                     4.219158
initial_val_target_1          1.360000
target_1_min_bound            1.450000
target_1_max_bound            1.600000
target_1                      0.080000
target_1_point                0.100000
ferro_0_opt                   0.104848
ferro_1_opt                   0.468808
target_1_opt                  0.099997
dtype: float64

В данном случае эксперимент проводился на исторических данных, поэтому у нас есть возможность сравнить рекомендуемые результаты с фактическими:
- **Факт**
- Добавляемые материалы:
- Ферросплав №0: 378 кг.
- Ферросплав №1 и №2: 0 кг.
- Прирост по $Mn$: 0.08%
- **Рекомендации**
- Цель (target_1_point) прирост по $Mn$: 0.1%
- Оценка прироста моделью (target_1_opt): 0.099997%
- Добавляемые материалы:
- Ферросплав №0: 104 кг.
- Ферросплав №1: 468 кг.
- Ферросплав №2: 0 кг.
- **Резюме по рекомендациям**
- по сетке удается найти точку, которая будет близка к установленной цели;
- в используемой критериальной функции не учитывается стоимость, поэтому полученное перераспределение не говорит об экономии.

Анализ вида критериальной функции возможно выполнить по графикам зависимости целевой функции от агрументов Рисунок 2 (добавить подпись картинке) и от номера итерации Рисунок 3 (добавить подпись картинке), по которому заметно, что задача имеет не единственное решение, и в качестве прогнозной моделеи используется градиентный бустинг.

### Линейная модель на ферросплавах
Для разработки линейной прогнозной модели используем робастные к выбросам алгоритмы, такие как $HuberRegressor$ из библиотеки $sklearn$. Важно, что линейная модель строится только признаках, содержащих в себе информацию о добавленных ферроспалавах, т.е. в данном случае на основании исторических данных мы пробуем подобрать оптимальные коэффициенты усвояемости материалов. Возможные подходы к подбору оптимальных коэффициентов - вопрос отдельной статьи, здесь предложен наиболее простой вариант.

Задачу оптимизации определим следующим образом:
\begin{equation}
1880*ferro_0^k+2000*ferro_1^k+3600*ferro_2^k \to \min\limits_{ferro_0^k, ferro_1^k, ferro_2^k} \end{equation} ограничения на прогнозы моделей
\begin{equation}
target\_min\_bound_i \leq pred\_target_i^k \leq target\_max\_bound_i, \; i \in \{1, 2\},
\end{equation}
допустимые границы ферросплавов
\begin{equation}
min\_ferro_j^k \leq ferro_j^k \leq min\_ferro_j^k, \; j \in \{0, 1, 2\},
\end{equation}
дополнительные ограничения на ферросплавы
\begin{gather*}
ferro_0^k + ferro_1^k \leq 0.4 \\
ferro_0^k + ferro_2^k \leq 0.5 \\
ferro_1^k + ferro_2^k \leq 0.4 \\
\end{gather*}

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

ferro_cost = {
    'ferro_0': 1880,
    'ferro_1': 2000,
    'ferro_2': 3600
}

2. Ограничения на ферросплавы

upper_bound = data.loc[idx,['tar_delta_1_upper_bound', 'tar_delta_2_upper_bound']]
lower_bound = data.loc[idx, ['tar_delta_1_lower_bound', 'tar_delta_2_lower_bound']]

A = np.array([[1, 1, 0],
              [1, 0, 1],
              [0, 1, 1]])

b=np.array([0.4, 0.5, 0.4])

3. Установить коэффициенты прогнозной модели и критериальной функции

lin_optim = lin_opt.LinOpt(
    objective_coefs=ferro_cost,
    model_coefs = model_coefs
)

5. Запустить оптимизатор

opt_res, opt_hist = lin_optim.run_linprog_opt(
    data=data.loc[idx],
    lower_bound=lower_bound,
    upper_bound=upper_bound,
    feat_bounds=ferro_bounds,
    A_ub=A,
    b_ub=b
)

**Заметка**
Очевидно, что линейной моделью проблематично описать влияние различных факторов технологического процесса. Но, с другой стороны, линейная модель очень удобна в оптимизации и хорошо интерпретируема. По этой причине мы используем комбинированный подход, состоящий из следующих шагов
1. Обучение линейной модели на добавках $\to lin\_pred\_target_i^k$
2. Расчет ошибки линейной модели
3. Обучение градиентного бустинга на ошибку линейной модели $\to error\_pred\_target_i^k$, причем из обучения исключаются управляемые параметры, т.е. итоговый прогноз модели состоит из двух компонент
\begin{equation}
pred\_target_i^k = lin\_pred\_target_i^k + error\_pred\_target_i^k
\end{equation}

Данный подход позволяет:
- сохранить возможность использовать методы для задачи линейного программирования;
- повышает точность прогноза;
- интерпритируемый прогноз.

## Байесовская оптимизация
### Смешивание критериев/ограничений
Как ранее говорилось, при использовании байесовской оптимизации нельзя в явном виде учитывать ограничения на прогноз модели и одновременно ставить задачу о снижении стоимости рассходуемых материалов. Решением данной задачи может быть учет нескольких критериев и ограничений через критериальную функцию, но это сразу же усложняет настройку алгоритма, а найденная точка необязательно будет удовлетворять всем ограничениям.
\begin{gather*}
L_1 = \frac{1}{C_{cost}}(1880*ferro_0^k+2000*ferro_1^k+3600*ferro_2^k), \\
L_2 = \sum_{i \in \{1, 2\}}\frac{1}{C_{target}^i}|pred\_target_i^k - goal\_target_i^k|, \\
L = L_1 + k \cdot L_2 \to \min\limits_{ferro_0^k, ferro_1^k, ferro_2^k}
\end{gather*} допустимые границы ферросплавов
\begin{equation}
min\_ferro_j^k \leq ferro_j^k \leq min\_ferro_j^k, \; j \in \{0, 1, 2\},
\end{equation}где
- $C_{target}^i$ -- нормировочный коэффициент для $i$-го таргета;
- $C_{cost}$ -- нормировочный коэффициент для стоимостной добавки;
- $k$ -- коэффициент взвешивания компонент критериальной функции.

Для применения данного подхода необходимо выполнить следующие шаги:
1. Ограничения на ферросплавы

space = {
    'ferro_0': hp.quniform('ferro_0', ferro_0_min, ferro_0_max, 0.01),
    'ferro_1': hp.quniform('ferro_1', ferro_1_min, ferro_1_max, 0.01),
    'ferro_2': hp.quniform('ferro_2', ferro_2_min, ferro_2_max, 0.01)
}

2. Расчет нормировочных коэффициентов

# Normalization constant for target variables
targets_norm = {
    'target_1': data['target_1'].abs().mean(),
    'target_2': data['target_2'].abs().mean()
}

# Normalization constant for summary cost of additives
median_cost = (
    data['ferro_0'].mean() * ferro_cost['ferro_0'] +
    data['ferro_1'].mean() * ferro_cost['ferro_1'] +
    data['ferro_2'].mean() * ferro_cost['ferro_2']
)

# Coefficient for part of summary cost of additives in total loss
cost_coeff = 0.25

3. Определить критериальную функцию

def hyper_opt_obj(x, models, target_point, targets_norm=targets_norm, ferro_cost=ferro_cost, median_cost=median_cost):
    total_loss = 0
    cost = 0
    for target_name in models.keys():
        # Target variable prediction by model
        target_pred = models[target_name]['model'].predict(x[models[target_name]['features']].values.reshape(1, -1))[0]
        
        # Diff between true and predicted target variable value
        target_loss = abs(target_point[f'{target_name}_point'] - target_pred) / targets_norm[target_name]
        
        # Total loss (sum of loss for each target variable)
        total_loss += target_loss
        
        # Loss by cost of additives
    for ferro in ferro_cost.keys():
        cost += x[ferro] * ferro_cost[ferro] 
    
    cost_loss = cost / median_cost
    total_loss += cost_loss * cost_coeff
    
    return total_loss

4. Запустить оптимизатор

hyper_optim = hyp_opt.HypOpt(
    models = gb_models
)

hyp_obj_params = {'target_point': data.loc[idx, ['target_1_point', 'target_2_point']].to_dict()}

opt_res, opt_hist = hyper_optim.run_hyper_opt(
    x=data.loc[idx],
    objective_params=hyp_obj_params,
    objective=hyper_opt_obj,
    search_space=space,
    max_evals=500,
)

### Optuna
Интресной фичей библиотеки Optuna является возможность решить задачу многокритериальной оптимизации, т.е. нет необходимости изобретать комбинированную целевую функцию, как это было в предыдущем разделе. Без дополнительных преобразований задача может быть записана следующим образом:
\begin{gather*}
L_1 \to \min\limits_{ferro_0^k, ferro_1^k, ferro_2^k}, \;L_2 \to \min\limits_{ferro_0^k, ferro_1^k, ferro_2^k}
\end{gather*} допустимые границы ферросплавов
\begin{equation}
min\_ferro_j^k \leq ferro_j^k \leq min\_ferro_j^k, \; j \in \{0, 1, 2\}.
\end{equation}

## Общие советы
По результатам прочтения статьи у вас должно было сложиться впечатление (но не точно), что при создании рекомендательного сервиса есть очень много нюансов практичных и математических, которые не сразу же бросаются в глаза. Несколько советов, о которых стоит помнить:
- выбранные параметры должны быть физически реализуемы
- например, на производстве редко есть возможность регулировать подачу ферросплавов с точностью до 1 кг.
- управляемый параметр, который вычисляется через какие-то статистики достаточно использовать в реальной жизни;
- выбранные исходные данные должны быть доступны в момент формирования рекомендации, этому может помешать:
- ручной ввод данных
- задержки во время передачи данных
- и много других уникальных ситуаций.
- данные в CSV и данные с прода редко совпадают;
- технические специалисты, которые в дальнейшем будут эксплуатировать систему, должны понимать и доверять ей, т.е. важна точность и интерпритация моделей.
- прогнозную модель на базе градиентного бустинга из коробки сложно сделать физичной, т.е. в случаях, где должна явно наблюдаться монотонная линейная зависимость, а бустинг ее не обеспечивает. В данном случае может помочь переход к комбинированной линейной модели и модели градиентного бустинга;
- в некоторых случаях можно предложить совсем иной подход к рекомендательному сервису.
Например, не строить прогнозную модель и не решать задчу оптимизации, а искать похожие плавки в истории и рекомендовать аналогичный исторический режим обработки. Причем это применимо во многих задачах и не только ферросплавов.
- и это только верхушка айсберга

# Итоги

**Ссылки на статьи**