Оптимизация расхода рабочего тела в промконтуре маслоохладителей АЭС по частоте вращения насоса#

 На момент публикации в практическом руководстве представлен всего лишь один метод оптимизации, который является глобальным (Differential Evolution). Данный метод является избыточным для многих практических задач на платформе REPEAT VISION (количество оптимизационных параметров слишком мало, пространство параметров является простым по форме, целевая функция имеет единственный экстремум и т.д.). Поэтому стоит задача исследовать другие оптимизационные методы в библиотеке Scipy. Для тестирования методов была предоставлена модель промконтура маслоохладителей АЭС, в которой надо оптимизировать расход рабочего тела по частоте вращения насоса.

Благодарность

Выражаем благодарность автору данного руководства, Герасимову Степану Андреевичу, за вклад в развитие платформы REPEAT VISION.

1.  Краткий обзор модели на платформе REPEAT VISION#

 На рис. 1 предствлена полная модель промконтура маслоохладителей АЭС:

Модель

Рис. 1 Полная модель промконтура маслоохладителей АЭС#

 Для задачи оптимизации представляет интерес блок (группа 21AP001), внутри которого находится центробежный насос и источник угловой скорости. На рис. 2 видна переменная w, которую и надо оптимизировать для достижения целевого расхода.

Рис. 2 Переменная w

Рис. 2 Переменная w#

 Расход измеряется датчиком массового расхода на самом выходе из система (на рис. 1 - справа). Для данного расхода установлен блок Выходной порт с именем переменной m. Она будет использоваться в дальнейшей оптимизации как целевая переменная.

2.  Постановка оптимизационной задачи#

 Пусть \(m(w)\) — массовый расход на выходе промконтура при угловой скорости насоса \(w\). Требуется найти такое значение \(w^*\), при котором расход равен заданной величине (1) :

(1)#\[{m}(w^*) = M_{\text{target}}, \quad M_{\text{target}} = 326.6\ \text{кг/с}\]

 Решается уравнение (2) :

(2)#\[m(w) - 326.6 = 0\]

 Для применения численных оптимизаторов переформулируем это в задачу минимизации квадрата невязки (3) :

(3)#\[\min_{w} F(w), \quad F(w) = \bigl( {m}(w) - 326.6 \bigr)^2\]

 При монотонности функции \(m(w)\) (что характерно для центробежного насоса) \(F(w)\) имеет единственный глобальный минимум (равный 0) в точке \(w^*\).

 Физические ограничения на скорость вращения насоса (рад/с) (4) :

(4)#\[w_{\min} = 140, \quad w_{\max} = 200\]

 Они найдены эмпирически, т.к. функция \(m(w)\) монотонная.

3.  Выбор методов оптимизации#

3.1.  Метод Брента для поиска корня#

Метод root_scalar с параметром brentq реализует классический метод Брента поиска корня.

Принцип поиска корня#

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

(5)#\[m(w) - 326.6 = 0\]

Примечание

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

 Алгоритм не выбирает один метод, а динамически комбинирует сразу три в зависимости от ситуации:

  1. Обратная квадратичная интерполяция: Это главный «ускоритель» метода. Когда у нас есть три точки с их значениями, алгоритм строит через них параболу, но не от x к y, а наоборот, от y к x. Затем он находит, где y=0, и получает очень хорошее приближение корня, что позволяет методу сходиться очень быстро на гладких функциях. На рис. 3 показана формула для данной обратной интерполяции:

Обратная квадратичная интерполяция

Рис. 3 Формула для обратной интерполяции#

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

  2. Метод деления отрезка пополам: Это «страховочный метод» алгоритма. Если быстрые, но рискованные шаги (интерполяция или секущие) привели к точке, которая не гарантирует надежного схватывания корня, алгоритм делает один гарантированный шаг, просто деля текущий интервал с корнем пополам. Это гарантирует, что интервал с корнем будет уменьшаться, и алгоритм никогда не «потеряет» корень.

3.2.  Оптимизационный метод Брента#

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

 Разберёмся, почему это так, и рассмотрим сильные стороны каждого метода, чтобы вы могли сделать осознанный выбор.

 Этот метод является де-факто стандартом в большинстве вычислительных библиотек (например, scipy.optimize.minimize_scalar, MATLAB). Он представляет собой гибридный подход, объединяющий сильные стороны сразу двух методов:

  • Метод золотого сечения – гарантирует сходимость.

  • Метод парабол – обеспечивает быстроту на гладких участках.

Принцип работы метода золотого сечения#

  1. Золотое сечение – деление отрезка на две неравные части так, чтобы отношение всего отрезка к большей части равнялось отношению большей части к меньшей (6) :

(6)#\[\frac{b-a}{\tau\cdot(b-a)} = \frac{\tau \cdot(b-a)}{(1-\tau)\cdot(b-a)} \quad \Rightarrow \quad \tau = \frac{\sqrt{5}-1}{2} \approx 0.618\]

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

  1. Алгоритм шага:

  • Задаём начальный интервал [a, b].

  • Вычисляем две внутренние точки:

\(x_1 = a + (1-\tau)(b-a) \quad \text{(левая точка)}\)

\(x_2 = a + \tau(b-a) \quad \text{(правая точка)}\)

  • Вычисляем значения функции \(f(x_1)\) и \(f(x_2)\).

    • Если \(f(x_1) < f(x_2)\), то минимум находится в \([a, x_2]\); новый интервал \([a, x_2]\).

    • Если \(f(x_1) \ge f(x_2)\), то минимум в \([x_1, b]\); новый интервал \([x_1, b]\).

    • На каждом шаге интервал уменьшается в \(\tau \approx 1.618\) раз.

    • Процесс повторяется до достижения требуемой точности по длине интервала.

Как работает метод Брента#

  1. Сначала пытается аппроксимировать функцию параболой через три точки.

  2. Если аппроксимация удачна и приводит к быстрому спуску к минимуму – используется параболический шаг.

  3. Если на каком-то участке парабола начинает вести себя нестабильно, алгоритм автоматически переключается на метод золотого сечения, чтобы гарантированно не потерять минимум и продолжить поиск.

Почему это лучший выбор#

 Комбинация двух подходов делает метод Брента одновременно:

  • Очень быстрым (благодаря параболической интерполяции);

  • Исключительно надёжным (благодаря золотому сечению).

3.3.  Метод Нелдера-Мида#

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

 Суть метода заключается в последовательном перемещении и деформировании симплекса вокруг точки экстремума.

Принцип работы#

 Алгоритм Нелдера-Мида (метод деформируемого многогранника) заключается в построении симплекса и последующем его деформировании в направлении минимума с помощью четырёх операций (см. рис. 4 ниже):

  1. Отражение (reflection)

  2. Растяжение (expansion)

  3. Сжатие (contraction) — бывает внешним и внутренним

  4. Редукция (shrinkage) — также называется сокращением

Рис. 4 Алгоритм Нелдера-Мида

Рис. 4 Алгоритм Нелдера-Мида#

 Симплекс — это геометрическая фигура, являющаяся n-мерным обобщением треугольника. Для одномерного пространства — это отрезок, для двумерного — треугольник. Таким образом, n-мерный симплекс имеет n+1 вершину.

 Алгоритм (одна итерация):

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

  2. К худшей точке последовательно применяются операции деформирования в следующем порядке:

    • Отражение — зеркальное отображение худшей вершины относительно центроида.

    • В зависимости от значения функции в отражённой точке, алгоритм решает:

      • если отражённая точка лучше лучшей → растяжение;

      • если отражённая точка средняя (между лучшей и худшей) → внешнее сжатие;

      • если отражённая точка хуже худшей → внутреннее сжатие;

      • если сжатие не дало улучшения → редукция (сжатие всего симплекса к лучшей вершине).

 На схеме показаны все четыре базовые операции: центроид обозначен серой точкой, худшая вершина — красной.

Примечание

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

 Условия остановки (любое из них):

  1. Выполнено заданное максимальное число итераций.

  2. Характерный размер симплекса (например, максимальное расстояние между вершинами или объём) стал меньше заданной величины xatol.

  3. Значение функции в лучшей вершине изменилось менее чем на fatol за несколько итераций (достигнута необходимая точность по функции).

Примечание

Метод Нелдера-Мида относится к прямым методам оптимизации (нулевого порядка) — он не требует вычисления производных. Сходимость не гарантирована для всех функций, но на практике часто эффективен, особенно для задач малой размерности (обычно \(n \le 20\)).

3.4.  Метод L-BFGS-B#

L-BFGS-B это усовершенствованная версия классического градиентного спуска.

Ограничения метода#

  • Требование к функции: Целевая функция должна быть гладкой и дифференцируемой. L-BFGS-B вычисляет и активно использует градиент. Для негладких или зашумленных функций существуют более устойчивые прямые методы, такие как упомянутый ранее метод Нелдера-Мида.

  • Поиск локального минимума: Как и большинство градиентных методов, L-BFGS-B находит локальный минимум.

Примечание

L-BFGS-B — это золотой стандарт для задач, где функция гладкая, количество параметров — от сотен до миллионов, но при этом есть простые ограничения на их значения.

3.5.  Какой метод выбрать#

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

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

 Для более сложных задач оптимизации, где несколько параметров (2–10) и неизвестная зависимость, рекомендуется использовать метод Нелдера-Мида.

 Если же целевая функция гладкая и дифференцируемая по всем параметрам, что часто встречается в физических моделях, и можно задать границы изменения переменных, то можно применить алгоритм L-BFGS-B, который окажется быстрее метода Нелдера-Мида.

3.6.  Сравнение методов на поставленной задаче#

Рис. 5 Сравнение 1

Рис. 5 Сравнение 1#

Рис. 6 Сравнение 2

Рис. 6 Сравнение 2#

 Как видно все методы дали одинаковое значение оптимальной частоты с точностью до 7 знаков после запятой. Самая быстрая скорость сходимости у метода Брента для поиска корня.

4.  Код для применения методов#

 Сначала представлена общая часть кода для всех методов. В нее входит подключение к модели на платформе REPEAT VISION, инициализация целевого значения и необходимых функций и т.д. Перед запуском необходимо вставить свой api token и id проекта, а также указать ваше время моделирования в [мс] в t_interval.

 1import repeat as rp
 2import pandas as pd
 3import numpy as np
 4import matplotlib.pyplot as plt
 5from scipy.optimize import minimize, minimize_scalar, root_scalar
 6import time
 7
 8user = rp.User(token = "TOKEN") #API
 9app = rp.Application(user)
10project = 'ID'
11t_interval = rp.TimeInterval(start = 0, end = 15000)
12model = app.get_exploration_model(project, t_interval)
13
14TARGET = 326.6   # целевой расход
15bounds = [(140.0, 200.0)]   # границы
16w0 = 130.0   # ваше начальное приближение
17HIGH_PRECISION = {  # параметры оптимизаторов
18    'Nelder-Mead':   {'maxiter': 1000, 'xatol': 1e-8, 'fatol': 1e-8, 'disp': False},
19    'L-BFGS-B':      {'maxiter': 1000, 'gtol': 1e-8, 'disp': False},
20    'Brent':         {'maxiter': 1000, 'xatol': 1e-8, 'fatol': 1e-8},
21    'root_scalar':   {'maxiter': 1000, 'xtol': 1e-8, 'rtol': 1e-8}
22}
23
24
25def result(w, model):
26    v1 = model.existing_variables.copy()
27    v2 = v1.copy()
28    v2['w'] = w
29
30    with model as md:
31        md.run(v2)
32        m = md.get_results('m')
33
34    results = {
35        'm': m
36    }
37
38    return results
39
40
41def get_rashod(w):
42    def run_simulation():
43        with model as md:
44            model_data = result(w, md)
45            rashod_val = model_data["m"].to_numpy()[-1]
46            return rashod_val
47
48
49    return run_simulation()
50
51
52def objective(w):
53    r = get_rashod(w)
54    obj = (r - TARGET) ** 2
55    return obj
56
57
58def root_function(w):
59    r = get_rashod(w)
60    f_val = r - TARGET
61    return f_val
  • Строка 14. Целевой расход.

  • Строка 15. Границы.

  • Строка 16. Ваше начальное приближение.

  • Строка 17. Параметры оптимизаторов.

  • Строки 26-28. Получаем вектор значений расхода от времени (для оптимизации берем последнее значение этого вектора).

  • Строки 42-46. Запускает модель с параметром w и возвращает расход в конце интервала.

  • Строки 52-55. Возвращает целевую функцию(квадрат ошибки).

  • Строки 58-61. Возвращает уравнение для посика корня.

 Теперь перейдем к фрагментам кода для каждого из методов (его необходимо вставить после общего фрагмента кода):

4.1.  Метод Брента для поиска корня#

 1opt_result = root_scalar(
 2            root_function,
 3            method='brentq',
 4            bracket=(bounds[0][0], bounds[0][1]),
 5            options={'maxiter': HIGH_PRECISION['root_scalar']['maxiter'],
 6                     'xtol': HIGH_PRECISION['root_scalar']['xtol'],
 7                     'rtol': HIGH_PRECISION['root_scalar']['rtol']}
 8        )
 9
10optimal_w = opt_result.root
11optimal_rashod = get_rashod(optimal_w)
12print(f"\nОптимум: w = {optimal_w:.4f}, расход = {optimal_rashod:.4f} (цель {TARGET})")

4.2.  Оптимизационный метод Брента#

 1opt_result = minimize_scalar(
 2    objective,
 3    w0,
 4    method='bounded',
 5    bounds=(bounds[0][0], bounds[0][1]),
 6    options=HIGH_PRECISION['Brent']
 7)
 8optimal_w = opt_result.x
 9optimal_rashod = get_rashod(optimal_w)
10print(f"\nОптимум: w = {optimal_w:.4f}, расход = {optimal_rashod:.4f} (цель {TARGET})")

4.3.  Метод Нелдера-Мида#

 1opt_result = minimize(
 2    objective,
 3    w0,
 4    method='Nelder-Mead',
 5    bounds=bounds,
 6    options=HIGH_PRECISION['Nelder-Mead']
 7)
 8optimal_w = opt_result.x[0]
 9optimal_rashod = get_rashod(optimal_w)
10print(f"\nОптимум: w = {optimal_w:.4f}, расход = {optimal_rashod:.4f} (цель {TARGET})")

4.4.  Метод L-BFGS-B#

 1opt_result = minimize(
 2    objective,
 3    w0,
 4    method='L-BFGS-B',
 5    bounds=bounds,
 6    options=HIGH_PRECISION['L-BFGS-B']
 7)
 8optimal_w = opt_result.x[0]
 9optimal_rashod = get_rashod(optimal_w)
10print(f"\nОптимум: w = {optimal_w:.4f}, расход = {optimal_rashod:.4f} (цель {TARGET})")

4.5.  Преимущества реализации оптимизации на Python перед встроенными алгоритмами на платформе REPEAT VISION#

  1. Гибкость в выборе и настройке методов:

    • На платформе REPEAT VISION предлагается фиксированный набор алгоритмов с ограниченным числом параметров.

    • Python + SciPy/NumPy дают доступ к десяткам методов и полному контролю над их гиперпараметрами (точность, максимальное число итераций, допуски, и т.д.).

  2. Визуализация процесса оптимизации:

    • Платформа REPEAT VISION показывает только итоговый результат. Python позволяет строить графики сходимости (значение функции от номера итерации), траектории движения в пространстве параметров, эволюцию симплекса, гистограммы вызовов.

    • Это помогает понять, почему алгоритм остановился, не застрял ли он в локальном минимуме, не слишком ли быстро уменьшился шаг. Визуализация — мощный инструмент для диагностики и тонкой настройки.

  3. Сбор и анализ диагностической информации:

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