Почему скрипт «замирает» после отчета HiGHS в SciPy MILP и как это исправить

Пауза после отчета HiGHS в SciPy 1.15.2 (MILP) вызвана O(N^2)-постобработкой результата. Обновитесь до 1.15.3 или 1.16.0 — задержка исчезнет для задач MILP.

Когда модель MILP в SciPy выводит отчет решателя HiGHS, а затем скрипт Python «замирает» на минуты прежде чем вернуться, может показаться, что ограничение по времени проигнорировано или решатель завис. На деле это нередко вызвано конкретным замедлением в SciPy 1.15.2, которое проявляется уже после завершения решателя: этап формирования результата внутри SciPy способен занять львиную долю общего времени и создать длинную паузу между напечатанным отчетом и вашим следующим print.

Воспроизводимый пример, демонстрирующий задержку

Ниже — фрагмент, повторяющий типичную постановку MILP с использованием SciPy’s milp, но с другими именами переменных для наглядности. Он генерирует случайные данные, собирает переменные и линейные ограничения, решает задачу с лимитом по времени и печатает общую длительность, измеренную в Python. В SciPy 1.15.2 она может оказаться значительно больше, чем «Timing (total)», который сообщает HiGHS.

import pandas as pd
import numpy as np
from scipy.optimize import Bounds, milp, LinearConstraint
import time
import scipy
print(scipy.__version__)
cap_hi = 2
cap_lo = 1
fr_hi = 1
fr_lo = 1
N_total = 300
np.random.seed(42)
data_tbl = pd.DataFrame({'Col1': np.random.rand(N_total), 'Col2': np.random.rand(N_total)})
def run_and_collect(batch_df, max_wall):
    n_rows = len(batch_df)
    if n_rows < 10:
        return [], batch_df.copy(), False
    group_cnt = n_rows // 10
    cap_arr = batch_df['Col1'].astype(float).to_numpy()
    fr_arr  = batch_df['Col2'].astype(float).to_numpy()
    obj_vec = np.concatenate([np.zeros(n_rows*group_cnt), -np.ones(group_cnt), np.zeros(n_rows)])
    integrality_mask = np.ones(n_rows*group_cnt + group_cnt + n_rows, dtype=int)
    dom_limits = Bounds(0, 1)
    L1 = np.zeros((n_rows, n_rows*group_cnt + group_cnt + n_rows)); r1_up = np.zeros(n_rows)
    L2 = np.zeros((n_rows, n_rows*group_cnt + group_cnt + n_rows)); r2_up = np.zeros(n_rows)
    for i in range(n_rows):
        L1[i, i + np.arange(group_cnt)*n_rows] = 1
        L1[i, n_rows*group_cnt + group_cnt + i] = -1
        L2[i, i + np.arange(group_cnt)*n_rows] = -1
        L2[i, n_rows*group_cnt + group_cnt + i] =  1
    C2a = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs2a = np.zeros(group_cnt)
    C2b = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs2b = np.zeros(group_cnt)
    for j in range(group_cnt):
        C2a[j, j*n_rows:(j+1)*n_rows] = 1
        C2a[j, n_rows*group_cnt + j]  = -10
        C2b[j, j*n_rows:(j+1)*n_rows] = -1
        C2b[j, n_rows*group_cnt + j]  = 10
    C3 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs3 = np.zeros(group_cnt)
    C4 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs4 = np.zeros(group_cnt)
    C5 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs5 = np.zeros(group_cnt)
    C6 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs6 = np.zeros(group_cnt)
    for j in range(group_cnt):
        C3[j, j*n_rows:(j+1)*n_rows]   = cap_arr
        C3[j, n_rows*group_cnt + j]    = -cap_hi*10
        C4[j, j*n_rows:(j+1)*n_rows]   = -cap_arr
        C4[j, n_rows*group_cnt + j]    =  cap_lo*10
        C5[j, j*n_rows:(j+1)*n_rows]   = fr_arr
        C5[j, n_rows*group_cnt + j]    = -fr_hi*10
        C6[j, j*n_rows:(j+1)*n_rows]   = -fr_arr
        C6[j, n_rows*group_cnt + j]    =  fr_lo*10
    Csym = np.zeros((group_cnt-1, n_rows*group_cnt + group_cnt + n_rows)); rhs_sym = np.zeros(group_cnt-1)
    for j in range(group_cnt-1):
        Csym[j, n_rows*group_cnt + j]   = 1
        Csym[j, n_rows*group_cnt + j+1] = -1
    C7 = np.zeros((1, n_rows*group_cnt + group_cnt + n_rows))
    C7[0, n_rows*group_cnt + group_cnt:] = 1
    C7[0, n_rows*group_cnt:n_rows*group_cnt + group_cnt] = -10
    cons = [
        LinearConstraint(C2a, -np.inf, rhs2a),
        LinearConstraint(C2b, -np.inf, rhs2b),
        LinearConstraint(C3,  -np.inf, rhs3),
        LinearConstraint(C4,  -np.inf, rhs4),
        LinearConstraint(C5,  -np.inf, rhs5),
        LinearConstraint(C6,  -np.inf, rhs6),
        LinearConstraint(C7,  -np.inf, 0),
        LinearConstraint(L1,  -np.inf, r1_up),
        LinearConstraint(L2,  -np.inf, r2_up),
        LinearConstraint(Csym, rhs_sym, np.inf)
    ]
    t0 = time.time()
    ans = milp(
        c=obj_vec,
        integrality=integrality_mask,
        bounds=dom_limits,
        constraints=cons,
        options={'disp': True, 'time_limit': max_wall}
    )
    t1 = time.time()
    print(f'Duration: {t1 - t0:.2f}')
run_and_collect(data_tbl, 10)

Что на самом деле вызывает задержку

Такое поведение совпадает с известной проблемой SciPy 1.15.2, затрагивающей оптимизацию на базе HiGHS. Внутри SciPy в _highs_wrapper() чтение col_status реализовано циклом с трудоемкостью O(N^2). Эта работа выполняется после завершения решателя, когда SciPy формирует объект результата. Поскольку это постобработка, на нее не влияет параметр time_limit — отсюда быстро появляющийся отчет решателя и значительно более позднее возвращение управления в Python. На практике даже задачи малого и среднего размера могут пострадать, если число переменных достаточно велико, чтобы шаг O(N^2) стал доминирующим.

Решение: обновите SciPy

Практичный выход — обновиться до версии SciPy, где это замедление устранено. Переход с 1.15.2 на 1.15.3 или 1.16.0 убирает O(N^2)-обработку col_status и выравнивает время возврата с таймингом самого решателя. Менять вашу модель или код не требуется.

Если хотите проверить установленную версию во время выполнения, выведите ее так:

import scipy
print(scipy.__version__)

На той же бенчмарк-конфигурации наблюдаемая длительность резко сокращается после обновления. Например, в приведенном выше примере «Timing (total)» HiGHS составляет около 1.6 секунды, тогда как сквозная длительность в Python падает примерно с 48 секунд на 1.15.2 до около 1.7 секунды на 1.15.3.

Почему это важно

В рабочих процессах MILP часто задают ограничение времени через time_limit и оборачивают вызовы решателя таймерами. Когда узким местом становится упаковка результатов после решения, эти договоренности незаметно нарушаются: задания выглядят как выходящие за бюджет, хотя решатель завершился вовремя. Несоответствие приводит к запутанной телеметрии, хрупким пайплайнам и лишним затратам. Устранение накладных расходов возвращает ожидаемую связь между отчетом решателя и реальным временем приложения.

Выводы

Если в SciPy 1.15.2 вы наблюдаете большие и непредсказуемые паузы после отчета HiGHS, замедление, скорее всего, возникает на этапе подготовки результата в обертке, а не в самой MILP-задаче. Обновление до SciPy 1.15.3 или 1.16.0 решает проблему и не требует правок в моделировании. Отслеживайте печатаемую версию SciPy и сравнивайте собственные измерения «настенных» часов с «Timing (total)» решателя, чтобы подтвердить улучшение.