Продолжаю рассказывать широкой общественности об интересном ML проекте, результаты которого внедряются в реальный технологический процесс. В Первой части разобрались, что такое глушение и почему важно уметь рассчитывать объемы жидкостей глушения. В этой части будет непосредственно все то, как мы решали эту задачу с помощью МЛ:

• Построили двухконтурную систему: офлайн-обучение на XGBoost и CatBoost — и лeгкий онлайн-инференс через Flask.

• Вместо одного .fit() с дефолтным лоссом внедрили K‑method — асимметричную функцию потерь. Теперь модель «боится» недолить жидкость сильнее, чем перелить, потому что в реальности эти ошибки стоят по-разному.

• CatBoost лучше по удобству работы с «сырыми» категориями, XGBoost потребовал кастомного кодирования, но дал сравнимые метрики.

• На малых данных (~350 строк) случайное разбиение творит хаос: метрики скачут от сида к сиду. Поэтому отбираем топ‑20 лучших random_state, а гиперпараметры усредняем частотным методом.

• Весь пайплайн — от подбора параметров до прогона K‑сетки — завернули в Airflow ради повторяемости, а все эксперименты и логи складываем в MLflow.

Дисклеймер: Из-за NDA в статье нет реальных названий месторождений, конкретных значений давлений и других признаков, а также точных формул блокирующих составов. Все результаты (, MAE и др.) приведены приблизительно. Весь показанный в статье код — это общие инженерные паттерны (кастомные loss-функции для градиентного бустинга, отбор random_state, частотное агрегирование гиперпараметров). Они не содержат доменных и проектных секретов и могут быть переиспользованы в любых регрессионных задачах с асимметричной стоимостью ошибки.

1. Почему просто «прогноз объeмов» не всегда работает и как “улучшить” таргет

Напомним кратко из первой части: стандартные гидродинамические модели плохо работают на трещиноватых коллекторах с аномально низким давлением. Технологи используют «интуицию + эмпирику», но это дорого (блокирующие составы стоят недешево). В процессе работы мы поняли: предсказывать абсолютный объeм V сложно, а вот предсказывать изменение относительно предыдущего глушения (delta) — гораздо стабильнее. Особенно для скважин, которые глушат повторно. Мы это объясняем тем, что по историческим данным в большинстве сложных глушений расчет необходимого объема основывался как раз на том, сколько использовалось блокирующих составов в предыдущее глушение, с некоторой корректировкой, учитывающей нынешнее состояние пласта/скважины. Также такое преобразование позволяет получить “более” нормальное распределение таргета. Для модели теперь таргет “-30” может быть одинаково интерпретируемым и для глушений с абсолютными объемами в 100 и 600 .

def add_delta_targets(df, liquids=("жидкость 1", "жидкость 2", "жидкость 3")):
    for liquid in liquids:
        df[f"разность_объемов_{liquid}"] = (
            df[f"объем_{liquid}"] - df[f"объем_{liquid} предыдущего глушения"]
        )
    return df

Поэтому в проекте существуют два независимых режима обучения:
V — предсказываем абсолютный объeм (используется как фолбэк для первичных глушений, где исторического prev_volume просто нет);
delta — предсказываем «сколько добавить/убавить» к предыдущему глушению (основной режим для повторных операций).

На инференсе для delta мы делаем простое действие: final_volume = prev_volume + predicted_delta.

2. Архитектура обучения: почему этапов много?

Обычно в ML-туториалах всe ограничивается train_test_split → fit → predict. В промышленности этого катастрофически мало. У нас данные:

• Маленькие (всего ~250-400 глушений на каждую модель после очистки).

• С шумом (ошибки в данных, человеческий фактор, см. Первую часть).

• С сильной зависимостью от разбиения (два случайных разбиения могут давать или ).

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

[Данные] ➔ [Этап 0: Подбор гиперпараметров (50 split)] ➔ [Частотное агрегирование]
               ➔ [Этап 1: Обучение на 50 random_state] ➔ [Этап 2: Отбор топ-20 моделей]
               ➔ [Этап 3: Применение K-loss для асимметрии]

Этап 0. Поиск гиперпараметров (частотный метод)

Проблема классического подхода

В стандартном ML мы делаем так:

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
search = RandomizedSearchCV(model, param_grid, cv=5)
search.fit(X_train, y_train)
best_params = search.best_params_

Что здесь плохо для наших данных:

1. Зависимость от одного разбиения. Так как выборка мала (~350), если random_state=42 неудачно разделил данные (например, все сложные глушения ушли в тест, а простые — в трейн), мы получим гиперпараметры, которые переобучены под «лeгкую» часть данных.

2. Нестабильность выбора параметров. На соседних random_state (41 и 43) best_params могут скакать: сегодня depth=7, завтра depth=4. Какой брать?

Наше решение: частотное агрегирование

Мы запускаем поиск не один раз, а для 50 разных разбиений. Каждое разбиение — свой random_state (от 1 до 50). Для каждого получаем свой словарь best_params. 50 — эмпирически найденная точка насыщения: при меньшем числе результаты нестабильны, при большем — время растeт без существенного прироста качества.

def _most_common_params(param_list: list[dict]) -> dict:
    """Возвращает моду (наиболее частое значение) для каждого параметра."""
    if not param_list:
        return {}
    
    param_frequencies: dict[str, dict[str, int]] = {}
    for params in param_list:
        if not params:
            continue
        for param, value in params.items():
            param_frequencies.setdefault(param, {})
            value_str = str(value)
            param_frequencies[param][value_str] = \
                param_frequencies[param].get(value_str, 0) + 1
    
    most_common_params = {}
    for param, frequencies in param_frequencies.items():
        # Берeм самое частое значение
        value_str, _ = max(frequencies.items(), key=lambda x: x[1])
        # Пробуем преобразовать обратно в число
        try:
            value = int(value_str)
        except ValueError:
            try:
                value = float(value_str)
            except ValueError:
                value = value_str
        most_common_params[param] = value
    return most_common_params

Почему это робастно:

• Если параметр действительно важен (например, max_depth влияет на переобучение), он будет часто выбираться во многих разбиениях.

• Если параметр нейтрален (например, min_child_weight в широком диапазоне даeт одинаковый loss), разные разбиения будут выбирать случайные значения — мода всe равно укажет на «среднюю» точку, не переобученную под шум.

• Этот метод автоматически отбраковывает выбросы: если в 48 разбиениях выбрано learning_rate=0.05, а в двух — 0.5 (переобучились под шум), мы возьмeм 0.05.

Этап 1. Перебор random_state

Фиксируем найденные гиперпараметры. Обучаем 50 моделей, меняя только random_state разбиения. Сохраняем каждую в отдельный .pkl. Выбор 50 обусловлен тем же принципом, что и на этапе поиска гиперпараметров: это минимальное количество разбиений, при котором топ-20 лучших random_state перестают меняться при перезапуске процедуры отбора. Мы проверяли это через Jaccard similarity между составами топ-20 для разного числа экспериментов — на 50 разбиениях Jaccard достигает 0.9.

Зачем? Чтобы потом отобрать устойчивые разбиения, которые дают стабильно высокий composite‑score (R² + MAE).

Этап 2. Отбор топ-20 разбиений

Вводим композитную метрику:

def calculate_composite_score(r2, mae, r2_weight=0.4, mae_weight=0.6) -> float:
    normalized_mae_score = 1 / (1 + mae) if mae > 0 else 1.0
    normalized_r2 = max(0, r2)
    return normalized_r2 * r2_weight + normalized_mae_score * mae_weight

• R² показывает, насколько модель объясняет дисперсию.

• MAE — средняя абсолютная ошибка в кубометрах.

• Веса (0.4 / 0.6) подобраны эмпирически: MAE важнее, потому что перелив/недолив в 3 м³ — это реальные деньги.

Почему именно композитный скор, а не просто R² или MAE?

Важный нюанс: отбор лучших моделей по валидационным метрикам — это не «утечка данных» (data leakage). Каждая модель обучается строго независимо на своем подмножестве. Мы просто отбираем те разбиения, которые оказались наиболее репрезентативными для генеральной совокупности.

Чтобы убедиться, что топ-20 моделей не переобучились под конкретные легкие примеры, мы проверяем равномерность покрытия всей выборки валидационными фолдами:

# Псевдокод проверки покрытия
coverage_matrix = pd.DataFrame(index=all_well_ids, columns=top_random_states)

for rs in top_random_states:
    train_idx, test_idx = train_test_split(df, random_state=rs)
    for well_id in test_idx:
        coverage_matrix.loc[well_id, rs] = True

# Проверяем: каждое ли глушение хотя бы раз попало в тест?
wells_not_in_test = coverage_matrix.any(axis=1).value_counts()

Вывод: топ-20 разбиений дают равномерное покрытие сложных примеров, поэтому модель видит в трейне достаточно «трудных» скважин и учится на них, а не просто запоминает легкие.

Подводим промежуточные итоги: для имеющихся данных мы строим ансамбль из 20 лучших моделей, гиперпараметры которых мы также отлично настроили. 20 лучших моделей соответствуют 20 различным разбиениям глушений на трейн и тест. Итоговым решением является усредненный объем по этим моделям.

Этап 3. K‑метод: учитываем асимметричную стоимость ошибки

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

Мы решили наказывать модель за недолив сильнее, чем за перелив, модифицировав стандартную MSE.

Математика K-метода

Стандартная MSE:

Градиент:

Гессиан (для XGBoost): (константа)

Наша модификация (асимметричная MSE):

Градиент:

Гессиан остаeтся (константа), что хорошо для XGBoost, которому для построения деревьев критически важно иметь стабильный гессиан.

Реализация для CatBoost

CatBoost позволяет передать произвольную функцию потерь через obj и метрику через eval_metric:

class RmseObjective:
    """Кастомная loss-функция для CatBoost"""

    def __init__(self, K: int):
        self.K = int(K)

    def calc_ders_range(self, approxes, targets, weights):
        assert len(approxes) == len(targets)
        if weights is not None:
            assert len(weights) == len(approxes)

        result = []
        for index in range(len(targets)):
            if approxes[index] < targets[index]:
                der1 = targets[index] - approxes[index] + self.K
                der2 = -1
            else:
                der1 = targets[index] - approxes[index]
                der2 = -1

            if weights is not None:
                der1 *= weights[index]
                der2 *= weights[index]

            result.append((der1, der2))
        return result


class RmseMetric:
    """Кастомная метрика для CatBoost"""

    def __init__(self, K: int):
        self.K = int(K)

    def get_final_error(self, error, weight):
        return np.sqrt(error / (weight + 1e-38))

    def is_max_optimal(self):
        return False

    def evaluate(self, approxes, target, weight):
        assert len(approxes) == 1
        assert len(target) == len(approxes[0])

        approx = approxes[0]

        error_sum = 0.0
        weight_sum = 0.0

        for i in range(len(approx)):
            w = 1.0 if weight is None else weight[i]
            if approx[i] < target[i]:
                weight_sum += w
                error_sum += w * ((approx[i] - target[i]) ** 2) + self.K * (target[i] - approx[i])
            else:
                weight_sum += w
                error_sum += w * ((approx[i] - target[i]) ** 2)

        return error_sum, weight_sum


# Обучаем модель с кастомным лоссом
model = CatBoostRegressor(
    **base_params, 
    loss_function=RmseObjective(K), 
    eval_metric=RmseMetric(K)
)

Обратим внимание, что calc_ders_range возвращает антиградиент der1 и отрицательный гессиан der2. Мы лишь следуем документации CatBoost, хотя этот момент может сбить с толку.

Реализация для XGBoost

А вот XGBoost требует функцию, возвращающую именно (grad, hess):

class XGBoostCustomLoss:
    """Кастомная loss-функция для XGBoost"""

    def __init__(self, K: int):
        self.K = int(K)

    def custom_loss(self, y_pred, dtrain):
        y_true = dtrain.get_label()
        grad = np.zeros_like(y_true)
        hess = np.ones_like(y_true)

        if self.K == 0:
            grad = y_pred - y_true
        else:
            for i in range(len(y_true)):
                if y_pred[i] < y_true[i]:
                    grad[i] = y_pred[i] - y_true[i] - self.K
                else:
                    grad[i] = y_pred[i] - y_true[i]

        return grad, hess


class XGBoostCustomMetric:
    """Кастомная метрика для XGBoost"""

    def __init__(self, K: int):
        self.K = int(K)

    def custom_metric(self, y_pred, dtrain):
        y_true = dtrain.get_label()

        error_sum = 0.0
        weight_sum = 0.0

        for i in range(len(y_pred)):
            w = 1.0
            if y_pred[i] < y_true[i]:
                weight_sum += w
                error_sum += w * ((y_pred[i] - y_true[i]) ** 2) + self.K * (y_true[i] - y_pred[i])
            else:
                weight_sum += w
                error_sum += w * ((y_pred[i] - y_true[i]) ** 2)

        final_error = np.sqrt(error_sum / (weight_sum + 1e-38))
        return "custom_rmse", final_error


# Пример использования
custom_loss = XGBoostCustomLoss(K)
custom_metric = XGBoostCustomMetric(K)
model = xgb.train(
        params,
        dtrain,
        num_boost_round=1000,
        evals=[(dtrain, "train"), (dtest, "test")],
        early_stopping_rounds=100,
        obj=custom_loss.custom_loss,
        custom_metric=custom_metric.custom_metric,
        verbose_eval=False,
    )

Сравнение K-метода с квантильной регрессией

Квантильная регрессия — это классический способ получить интервальные оценки. Но у нас не просто «нужен верхний квантиль», а штраф за ошибку разный для недолива и перелива.

Сравнение на реальных данных ( vs квантиль 0.8, указанные результаты приблизительны):

Вывод: K-метод при дал ту же защиту от недолива (18% < 20%), но с меньшим переливом (3.0 м³ против 6.0 м³). Потому что модель не просто «сдвинула все предсказания вверх», а научилась смещать только те случаи, где это действительно нужно (по признакам).

Подводим промежуточные итоги: теперь для всех 20 лучших моделей мы реализуем К-метод (К может варьироваться по-разному, но чтобы не плодить много похожих моделей, меняется с шагом 20 от 0 до, например, 600). Для каждого имеем усредненное значение по 20 моделям. Осталось определить оптимальный , чтобы минимизировать количество неуспешных глушений (случаи, когда модель прогнозирует заведомо низкий объем) и минимизировать количество глушений, для которых мы спрогнозировали заведомо завышенный объем.

3. Подбор оптимального K (после обучения)

Для наглядности построим два графика: график уменьшения доли неуспешных глушений (недоливов) с ростом и график увеличения доли переливов с ростом .

Как строятся эти графики: рассмотрим все тесты из 20 лучших разбиений. Мы знаем, что все глушения попали в эти тесты (в среднем по 7 раз). Для для каждого глушения находим средний прогнозный объем, определяем долю недолива (количество глушений, у которых средний прогнозный объем меньше фактического на 1 ) и долю перелива (количество глушений, у которых средний прогнозный объем больше фактического на 2 ). 1 , 2 - это пороги, они задаются по-разному, в зависимости от месторождения и обычно обсуждаются непосредственно с теми, кто занимается глушениями.

Теперь для каждого мы можем посчитать долю недолива и долю перелива и наконец-то построить красивые графики.

Внимательно посмотрим на кривые: классическая задача trade-off. Можно реализовать какой-то скрипт для задачи оптимизации, можно построить график в терминах перелив-недолив, который позволит найти , при котором недолив падает, а перелив практически не растет. Все это мы сделали, но математическая оптимизация здесь часто выдавала пограничные значения . И тогда мы вспомнили, что наша основная задача - это обеспечить долю неуспешных глушений . Значение (допустимая доля неуспешных глушений) определяется совместно с технологическим отделом. Например, если исторически 25% глушений требовали повторных операций, мы можем поставить цель снизить этот показатель до 15%. Тогда = 15% — и мы находим , который это обеспечивает. Если заказчик говорит: “недолив недопустим”, мы берeм = 0% (или 1-2% на статистический шум). Такой подход прозрачен и легко согласуется с бизнесом. В этом случае все упрощается: просто найдем , при котором доля недоливов равна . Почему это лучшее решение? Потому что это интерпретируемо, объяснить заказчику такой подход становится очень легко. “Мы используем такой метод и такой , потому что это гарантирует долю неуспешных глушений в , при этом возможен минимальный перерасход жидкостей. Если глушение раньше было сложным, то можем слегка поднять для гарантии успеха.” На опытно-промышленных испытаниях такой подход показал себя отлично.

4. CatBoost vs XGBoost: честная битва на наших данных

Внимательный читатель мог заметить, что выше мы приводили примеры кода для обеих библиотек. Так вышло не сразу. До последнего момента мы использовали CatBoost как единственную модель машинного обучения, так как изначально проверили метрики на одном random_state, и CatBoost показал наилучшие и MAE. И конечно, некоторое количество категориальных признаков сыграло в пользу выбора этой модели. Однако в финальной части проекта мы решили удостовериться, что наш выбор действительно оправдан. Для этого провели исследование: проанализировали метрики (и composite score тоже) топ-20 разбиений для Random Forest, LightGBM, XGBoost. Random Forest и LightGBM показали заметно хуже метрики на таких малых данных. А вот XGBoost показал MAE на примерно 10% меньше, чем CatBoost. Да, требуется дополнительная обработка категориальных признаков, например такая:

bin_enc = ce.BinaryEncoder(cols=["номер скважины", "куст"])
ohe_enc = ce.OneHotEncoder(cols=["месторождение", "пласт"])
X_encoded = ohe_enc.fit_transform(X)
X_encoded = bin_enc.fit_transform(X_encoded)

Не забываем сохранять модель и кодирование категориальных признаков, чтобы на этапе инференса все работало правильно. После этого реализовали полноценный К-метод для XGBoost и на одном графике сравнили результаты CatBoost и XGBoost. Видим, что с ростом К CatBoost стремится минимизировать риск любой ценой, что на малых данных приводит к стратегии «перестраховки». А XGBoost демонстрирует более консервативное и робастное поведение.

Что это дает на практике? У XGBoost MAE примерно на 10% ниже, чем CatBoost, что означает меньший перерасход жидкостей в среднем. Однако CatBoost лучше контролирует хвосты распределения — количество катастрофических ошибок (сильный недолив) у него ниже. Поэтому на практике мы используем оба метода: CatBoost — для сложных скважин (где цена ошибки высока), XGBoost — для типовых глушений (где важна экономия).

5. Инференс и генерация отчeта (Flask + визуализация)

Пайплайн инференса в веб-интерфейсе выглядит так:

1. Пользователь (технолог) заполняет параметры скважины в UI.

2. Система извлекает исторический контекст (параметры и объемы предыдущего глушения этой скважины из БД).

3. Строка признаков прогоняется через ансамбль из 20 сохраненных моделей для сетки .

4. На экран выводится интерактивный график зависимости объема от коэффициента , а также Feature Importance для конкретного предсказания (локальная интерпретируемость SHAP/усредненные веса).

[Интерфейс Flask] ➔ [Запрос истории глушений] ➔ [Расчет Ансамбля 20 моделей] ➔ [Вывод графика объема от K]

Такой подход дает инженеру гибкость: он видит базовый прогноз () и может оценить, сколько кубометров состава добавит модель, если он захочет подстраховаться (выбрав K > 0). Для удобства мы завернули решение в легковесное Flask-приложение с возможностью дообучения моделей при загрузке свежих данных.

6. Airflow + MLflow

Чтобы проект не превратился в «набор Jupyter-ноутбуков на коленке», все этапы обучения и валидации оркеструются в Apache Airflow:

• DAG автоматически запускает переобучение при накоплении критической массы новых данных по глушениям.

• Задачи подбора параметров для V и delta (find_hyperparams) выполняются параллельно.

• Далее последовательно выполняются шаги: train_50_models select_top_20 generate_k_curves.

• Все метрики (, MAE, параметры лучших моделей, графики trade-off) автоматически логируются в реестр моделей MLflow.

Заключение

Подходы, реализованные в проекте, позволили нам:

• Уйти от предсказания абсолютов к дельтам, стабилизировав таргет.

• Внедрением кастомного -loss решить проблему асимметричной стоимости ошибок.

• Построить устойчивый ансамбль на топ-20 разбиениях, решив проблему дефицита и шумности данных.

• Объединить сильные стороны XGBoost и CatBoost в единую систему принятия решений.

Главный итог: метод признан успешным по итогам опытно-промышленных испытаний. Внедрение системы позволило снизить затраты на глушение за счёт более точного прогнозирования объёмов и контролируемого баланса между риском недолива и перерасходом жидкостей. Результаты проекта тиражируются на другие месторождения.

Вопросы в комментариях приветствуются — постараюсь ответить, не нарушая NDA.