Разворачиваем систему уравнений в граф24.01.2024 16:00

Как-то во время учебы на инженера-теплоэнергетика я наткнулся на одну книгу [Попырин Л.С. Математическое моделирование и оптимизация теплоэнергетических установок. М. Энергия 1978 г.], в которой был описан алгоритм построения расчётных схем энергетических установок, разработанный в Сибирском энергетическом институте (ныне — ИСЭМ СО РАН). Этот алгоритм заложен в основу СМПП (система машинного построения программ) — кодогенератора, который используется в исследованиях в ИСЭМ и по сей день. Собственно алгоритм предназначен для решения систем нелинейных уравнений, и, условно говоря, обобщает метод подстановки, знакомый всем из школьной алгебры.

Зачем это нужно?

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

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

В библиотеке TESPy такие системы уравнений решаются методом Ньютона-Рафсона. Для схем электростанций, состоящих из нескольких агрегатов (котлов, турбин, бойлерных), время единичного расчёта целевой функции может превышать 1 секунду, что для задачи многокритериальной оптимизации по 30–40 переменным в большинстве случаев будет неудовлетворительно.

Особенность таких схем состоит в том, что большинство уравнений связывает между собой небольшое количество переменных. Собственно алгоритм, разработанный в ИСЭМ, предварительно ориентирует каждое уравнение на вычисление одной переменной, и упорядочивает вычисления этих переменных, обеспечивая поочередное решение уравнений на основании результатов предыдущих вычислений. Таким путём, как утверждается в первоисточнике, размерность задачи снижается в 10–20 раз.

Разберём на примере

Возьмём систему:

Предположим, что хотим её решить при заданном значении $x_{3}$ . Понятно, что тогда здесь по сути 2 уравнения и 2 переменные, однако такие избыточные (но не переопредёленные) системы как раз довольно часто встречаются на практике (расход среды на входе = расход на выходе). Представим ее в виде двудольного графа: в одной доле — уравнения, в другой — переменные, а ребра показывают наличие переменной в уравнении.

Задаётся желаемое множество исходных данных, в данном случае $X_{0}=\left\{ x_{3} \right\}$ . Задача — определить, какая переменная должна определяться из каждого уравнения, при этом рёбра между уравнениями и переменными в множестве X0 не должны использоваться. Это является вариантом задачи поиска наибольшего паросочетания. Способов решения существует множество — алгоритмы пополняющего пути (augmenting path), сведение к задаче о назначениях и применение венгерского алгоритма, программирование в ограничениях (Constraint Programming), целочисленное программирование и пр.

Во время учебы я написал в Матлабе венгерский алгоритм, который очень примитивно подходил к требованию, связанным с исходными данными, — на одном из шагов случайно выполнял назначения, когда было несколько вариантов, и потом проверял, не получилось ли так, что вершины из множества X0 оказались задействованы. Решалось это чудо 10 секунд для системы из 22 уравнений.

В какой-то мере это не давало мне покоя, и недавно дошли руки интереса ради нормально реализовать решение такой задачи. Выбор пал на солвер Max-flow из Google OR Tools. Прежде всего, необходимо преобразовать нашу задачу к задаче поиска максимального потока. Для большего разнообразия можно добавить требование к некоторым переменным быть вычисляемыми, например, $Y_{0}=\left\{ x_{1} \right\}$ .

Шаг 1. Максимальный поток

Ориентируем рёбра (теперь уже дуги) от уравнений к переменным, добавляем стартовую вершину s, дуги от нее к уравнениям с пропускной способностью (capacity) = 1, и конечную вершину t.

Для переменных 3 варианта:

Для переменных в Y0: дуги (зеленые) с capacity = 1 направляем сразу к t;
Для переменных (красных) в X0: дуги не выходят, и по свойству сохранения потока условный поток через все вершины в X0 будет в любом решении равен нулю;
Для всех остальных: дуги (синие) выходят с capacity = 1 в промежуточную вершину t0, и одна дуга идет из t0 в t с capacity, равной количеству входящих в t0 дуг.

Выше по факту описано приведение нашей задачи к задаче «Circulation with demands with lower bounds», а затем к Max-flow. Получается такая картина, числами обозначены пропускные способности:

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

По сравнению с моей наивной реализацией в Матлабе, солвер max-flow в OR Tools справляется с задачей из 6000 уравнений менее, чем за 0,01 с.

После решения самой задачи Max-flow отбрасываем левую и правую части графа, выделяем ненулевые потоки и разворачиваем дуги с нулевыми потоками:

Шаг 2. Компоненты сильной связности

Далее объединяем каждую пару уравнение-переменная в одну вершину (слева) и применяем алгоритм Тарьяна для поиска компонент сильной связности и одновременно топологической сортировки полученного графа (справа).

Шаг 3. Разрыв циклов

Далее для каждой компоненты необходимо разорвать все циклы, причем сделать это, удалив желательно минимальное количество дуг. Это задача поиска минимального разрезающего циклы набора рёбер (Minimum feedback arc set), является NP-сложной, существуют точные и аппроксимационные алгоритмы. После удаления этих дуг опять же производится топологическая сортировка полученного ациклического подграфа.

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

Собственно код

Можно генерировать код. Исходные данные — список из одной переменной.

def solve(x):
  # Присваиваем x3 значение из исходных данных
  x3 = x[0]
  # Задаем начальные значения x2 и x1
  x2 = 0.5
  x1 = 0.5
  while True:
    temp_x2 = x2
    temp_x1 = x1
    # Выражаем каждую переменную из соответствующего ей уравнения
    x0 = x2                            #eq2
    x2 = ((x0 - x1**2)/x3)**(1/3)      #eq0
    x1 = x2**0.5 - x0                  #eq1
    # Сравниваем с предыдущей итерацией с относительной погрешностью 1e-4
    # +1e-8 в знаменателе, чтобы избежать деления на ноль
    if abs((x2 - temp_x2) / (temp_x2 + 1e-8)) > 1e-4: continue
    if abs((x1 - temp_x1) / (temp_x1 + 1e-8)) > 1e-4: continue
    break
  # Возвращаем вычисленные переменные в виде массива
  y = [0.0] * 3
  y[0] = x2
  y[1] = x1
  y[2] = x0
  return y

%timeit res = solve([2])

8.76 µs ± 59.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Для сравнения, решаем с помощью fsolve из scipy:

def func(x,x3):
    out = np.empty(3)
    out[0] = x[1]**2 + x3*x[2]**2 - x[0]
    out[1] = (x[0] + x[1])**2 - x[2]
    out[2] = x[0] - x[2]
    return out

x = np.array([0.5,0.5,0.5])
x3 = 2
%timeit res = fsolve(func, x, args=(x3))

43.8 µs ± 1.27 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Сюда выложил код такого генератора, в котором можно задать систему уравнений следующим образом:

from eqsystem import EqSystem, EqSystemGraph, PyScriptGenerator


S = EqSystem()
# Начальное приближение можно задать для всех переменных
S.add_var('x0', 0.5)
S.add_var('x1', 0.5)
S.add_var('x2', 0.5)
S.add_var('x3', 0.5)

eq0 = {'x0': (True, f"x0 = x1**2 + x3 * x2**3"),
       'x1': (True, f"x1 = (x0 - x3 * x2**3)**0.5"),
       'x2': (True, f"x2 = ((x0 - x1**2)/x3)**(1/3)"),
       'x3': (False, f"x1**2 + x * x2**3 - x0")}
eq1 = {'x0': (True, f"x0 = x2**0.5 - x1"),
       'x1': (True, f"x1 = x2**0.5 - x0"),
       'x2': (True, f"x2 = (x0 + x1)**2")}
eq2 = {'x0': (True, f"x0 = x2"),
       'x2': (True, f"x2 = x0")}
S.add_std_eq('eq0', eq0)
S.add_std_eq('eq1', eq1)
S.add_std_eq('eq2', eq2)

S.set_desired_x({'x3'})
S.set_desired_y({'x1'})

G = EqSystemGraph(S)
G.solve_matching_maxflow_ortools()
G.scc()
G.mfas()
gen = PyScriptGenerator("out.py", "solve")
G.generate_code(gen)

Если невозможно (или не хочется) из уравнения выразить какую-либо переменную, можно эту пару решать численно, например, с помощью fsolve. Так, для x3 False означает, что уравнение нужно решить численно относительно x.

Основные тезисы

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