Солвер, который решает всё, что может: принцип HYSYS-класса на примере Sim42

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

В литературе выделяют два основных подхода к решению схемы^{[1, §1.1.1]}. В уравнение-ориентированном (equation-oriented) подходе уравнения всех аппаратов собираются в одну систему и решаются одновременно методами типа Ньютона. В модульном (modular) подходе схема решается аппарат за аппаратом; модульные симуляторы делятся на последовательные (sequential modular), где расчёт идёт строго от сырья к продукту, и непоследовательные (non-sequential modular), способные распространять информацию вперёд и назад, так что пользователь может задать данные в любой точке схемы. Определение последнего типа диссертация, на которую опирается эта статья, связывает с работой Svrcek и соавторов 1984 года^{[1, §1.1.1]}.

Коммерческие солверы закрыты, и разобраться, как именно «думает» непоследовательно-модульное ядро, по документации обычно нельзя. Этим объясняется выбор объекта изучения: Sim42 — открытый симулятор, исходный код которого доступен^[7], а принцип работы солвера описан в диссертации его автора^[1] и подробно разобран в аналитических материалах по коду^{[12, 13]}. Это позволяет показать механизм «изнутри» и на конкретных числах.

2.1. Калгарийская школа и HYSYS

Компания Hyprotech создана исследователями Университета Калгари; профессор В. Свэрчек (W. Y. Svrcek) значится партнёром Hyprotech с момента её учреждения в 1976 году^{[9, 10]}. Продукты компании — HYSIM, затем HYSYS (справочник по HYSYS версии 1.1 издан в 1996 году). В мае 2002 года Hyprotech приобрела AspenTech; впоследствии, по решению Федеральной торговой комиссии США, права на часть активов, включая исходный код HYSYS, были переданы Honeywell, на базе чего появился UniSim^[9]. Цены и даты сделок в источниках расходятся, подробности для сути статьи несущественны.

Для истории солвера важнее другое: в диссертации по Sim42 HYSYS перечислен среди «непоследовательно-модульных» симуляторов (Канада), а само понятие со ссылкой на Svrcek et al. (1984)^{[1, §1.1.1–1.1.2]}. Оригинал этой работы найти в открытых источниках не удалось, поэтому ссылка приводится «через диссертацию».

2.2. Рождение Sim42 (2001–2003)

Хронологию проекта дают две группы первоисточников: страница «The History of Sim42» на сайте проекта (архивная копия от 9.12.2002)^[2] и благодарности в диссертации^[1]. В мае 2001 года Raúl Cota (Рауль Кота) начал в Университете Калгари магистерскую работу — создание симулятора с открытым кодом на Python; руководителями были Marco Satyro, Bill Svrcek и Brent Young^{[1, 2]}.

Участники проекта и их вклад по тем же источникам^{[1, 2, 3, 4]}: Raúl Cota — исходная архитектура (аппараты, порты, интерфейс к термодинамике), графическая версия командной строки, испанские сообщения; Craig Morris — солвер, колонна, контроллеры и скриптовые аппараты, первый командный интерфейс, хостинг репозитория; Marco Satyro — термодинамика (через VMG), руководство пользователя, научное руководство; Bill Svrcek и Brent Young — научное руководство; Yau Kun Li — насосы с кривыми, конверсионные реакторы, гипотетические компоненты; Andrew Liu и Rajeev Agarwal — код и спецификации аппаратов; Patrick O’Brien — интерпретатор PyCrust; Jose-Henrique Pinto — переводы сообщений и FAQ. Систему единиц передала компания VMG^[2].

2.3. Откуда пришли идеи солвера

Автор солвера Craig Morris в профиле на сайте проекта пишет, что «когда-то написал симулятор Hysim» и покинул Hyprotech за десять с лишним лет до 2002 года; об устройстве будущего солвера у него были «довольно необычные идеи», и он вызвался написать солвер для проекта Cota^[3]. Satyro, в свою очередь, отмечает, что Sim42 «унаследовал» от Morris двунаправленную передачу информации, вложенность аппаратов и командный язык^[4]. Эти свидетельства самоописательны (профили участников, написанные ими самими), но совпадают с описанием солвера в диссертации и по коду^{[1, 13]}.

Следовательно, правомерно сказать: принцип непоследовательно-модульного солвера, реализованный в Sim42, — того же класса, что и в HYSYS, и создан в среде и руками людей, связанных с Hyprotech. Утверждать идентичность алгоритмов или происхождение кода HYSYS в Sim42 открытые источники оснований не дают.

2.4. Лицензия и сообщество

Код Sim42 распространялся по «SIM42 Software License, Version 1.0» (© 2002 SIM42 Software Foundation)^[6]. Сайт и FAQ называли её BSD-подобной: разрешены использование и изменение без обязательной передачи изменений, что важно для компаний, добавляющих собственные модели^{[1, 5]}. Текст лицензии содержит условия атрибуции: в документации распространяемого продукта должно быть указание на разработку SIM42 Software Foundation; производные продукты нельзя называть SIM42 без письменного разрешения^[6]. В 2008 году участники проекта OpSim сочли эти условия несовместимыми с GPL и разработали собственный симулятор «с нуля», используя Sim42 как «источник идей»^[11] — это мнение участников, а не юридическое заключение.

2.5. После 2004 года

Дальнейшую судьбу проекта в открытых источниках удаётся проследить лишь фрагментарно. Домен sim42.org с 2013 года в веб-архиве представлен страницами-«парковками». Дистрибутив версии 2.0.0.0 сохранился на SourceForge в проекте OllinTS (зарегистрирован в 2007 году, последнее обновление — 2015; лицензия GPLv2) — это термодинамический сервер, который создавался «для завершения Sim42»^[8]. В 2011 году этот дистрибутив был импортирован в репозиторий GitHub, после чего автором импорта начата конвертация с устаревшего расширения Numeric на NumPy^[7]; последний коммит — ноябрь 2011 года. Остатки обращений к Numeric в отдельных модулях, отмеченные в аналитических материалах^{[12, 13]}, — следы этой незавершённой миграции. Признаков активного сопровождения в открытых источниках не обнаружено.

3.1. Главная идея

В классическом последовательном симуляторе пользователь подаёт сырьё, затем считается первый аппарат, потом второй и так далее. Если нужно «задать то, что на выходе», схему приходится перекраивать вручную. Солвер класса Sim42 работает иначе. Представьте большую таблицу судоку: каждая клетка — параметр потока (температура, давление, расход, состав), а каждый аппарат — набор правил («сумма на входе равна сумме на выходе»). Как только в какой-то клетке появляется число, все правила, которым хватает данных, срабатывают и заполняют соседние клетки — в любую сторону^{[1, 12]}.

• Информация идёт в обе стороны. Если известны состав и расход смеси на выходе смесителя и один из двух входов, второй вход вычисляется как разность.
• Порядок ввода не важен. Можно сначала задать давление на выходе, потом температуру на входе — результат тот же.
• Расчёт идёт сразу после каждого ввода (если не включён режим паузы), поэтому схему можно собирать и «чувствовать» в интерактивном режиме.
• Часто рецикл считать не нужно. Авторы отмечают, что за счёт распространения частичной информации во все стороны многие сложные задачи решаются без итераций по рециклам^[5].
• Защита от переопределения. Если одно значение получено из двух источников и они не совпали, программа не молчит, а сообщает о противоречии.

3.2. Из каких «кирпичиков» состоит модель

Аппарат (unit operation) — любой объект схемы: смеситель, нагреватель, колонна, насос, контроллер. В Sim42 потоки — тоже аппараты, а не отдельная сущность; аппараты можно вкладывать друг в друга (теплообменник собран из нагревателя и охладителя)^[1]. Порт (port) — «розетка» аппарата, через которую течёт информация; бывают материальные порты (T, P, расход, состав), энергетические и сигнальные (одно число — перепад давления, управляющий сигнал). Направление порта — это направление потока вещества, а не порядок расчёта: информация через порт ходит в обе стороны.

Каждое значение хранится вместе со статусом происхождения — откуда оно взялось. Именно статус позволяет программе решать, что можно перезаписывать, а что нельзя.

3.3. Степени свободы и контроль избыточности

У каждого аппарата есть физически обоснованное число параметров, которые нужно задать: например, у нагревателя — состав и расход на входе, T и P на входе, плюс одно из трёх: температура на выходе, тепловая нагрузка или перепад давления. Когда задано ровно столько, сколько нужно, аппарат считает остальное. Если задано меньше — часть значений останется None. Если больше — возможны противоречия. Руководство Sim42 называет это «Interactive Simulation and Degrees of Freedom Monitoring»^[12].

Правило, на котором это построено, универсально: если в балансе не хватает ровно одного значения — оно вычисляется; если не хватает двух и более — программа ждёт новых данных; если известно всё — проверяется, что баланс сошёлся. Расчётное значение никогда не перезаписывает уже известное: оно только сверяется с ним. Расхождение сверх допуска считается противоречием (ConsistencyError). Сверка идёт не в абсолютных единицах, а в долях масштаба переменной, чтобы 0,1 К и 0,1 Па не считались равнозначными^{[12, 13]}.

3.4. Порядок ввода не важен, история — тоже

Расчёт запускается не кнопкой «Рассчитать», а любым изменением данных. Команда hold приостанавливает расчёт — удобно, когда нужно ввести много данных; go возобновляет его. После go результат должен совпасть с вводом тех же данных без паузы. Самый сильный тест на корректность: «задал 80, потом изменил на 90» и «сразу задал 90» обязаны дать один и тот же результат^[12]. Как солвер этого добивается, не храня граф зависимостей, разобрано в разделе 5.5.

3.5. Как читать результат

Песочница использует модель из раздела 4. Попробуйте набор «Только с выхода»: по продукту и параметрам аппаратов солвер полностью восстановит вход; «Недоопределён» — получите None без ошибки; «Лишняя, конфликтная» — получите ConsistencyError.

4.1. Постановка

Схема из трёх аппаратов: нагреватель (Heater), насос (Pump) и холодильник (Cooler). Вместо перепадов температур заданы температуры на выходе (Heater — 80 °C, Cooler — 40 °C), у насоса задан КПД. Рабочее тело — вода; расход и прочие условия приведены в таблице 6. Реализованные в лаборатории правила соответствуют описанию в [12] и [13]: флэш при двух известных из (T, P, H) и составе, мольный и энергетический балансы по принципу «ровно одно неизвестное», расчёт перепада давления, изоэнтропа насоса.

Насос — составной аппарат: Ideal (идеальное сжатие по изоэнтропе), Waste (потери превращаются в тепло, ΔP = 0), Set (связывает КПД с энергиями: Q_ид = η·Q_вал) и Mix (энергобаланс Q_вал = Q_ид + Q_пот). Каждая часть решает только свою маленькую задачу; вместе они дают полный расчёт. Это пример вложенности аппаратов, о которой говорится в описании Sim42^[1].

4.2. Что происходит при вводе

1. T и P питания ничего не дают: для флэша нужен ещё состав. Аппарат — Feed — ставится в очередь, расчёт стоит.
2. Состав питания завершает тройку (T, P, z): сработал флэш, появилась энтальпия H(S1). Heater получает «новое знание» и встаёт на вершину стека.
3. ΔP нагревателя и T на выходе: давление на выходе Heater = 100 − 20 = 80 кПа; по (T, P) выхода — энтальпия H(S2). Расхода пока нет, поэтому тепловая нагрузка ещё неизвестна — в энергобалансе два неизвестных (F и Q).
4. Насос и холодильник заполняют то, что можно вычислить без расхода: давления (580 и 550 кПа) и энтальпии потоков. Температура после насоса пока неизвестна: она зависит от потерь в тепло, а потери — от расхода.
5. Последним вводится расход продукта F = 100 кмоль/ч — на выходе схемы. Мольный баланс отправляет его назад по всей схеме до питания; как только расход известен, в каждом энергобалансе остаётся ровно одно неизвестное, и один за другим определяются Q нагревателя, энергии насоса и нагрузка холодильника.

Именно пятый шаг показывает отличие от последовательного расчёта: расход, заданный в конце схемы, вычисляет всё в её начале. Модель можно пройти по шагам: фазы расчёта подсвечены, стек решения виден, значения вспыхивают в момент вычисления.

4.3. Порядок ввода меняет путь, но не итог

Если вводить те же данные от выхода к входу, число крупных шагов (один шаг — один решённый аппарат или ввод) меняется, а итоговые значения остаются теми же. В режиме hold → go очередь копится и решается одним проходом, поэтому шагов заметно меньше.

4.4. Смена спецификации: забывание и решение

Что произойдёт, если после полного расчёта изменить температуру выхода Heater с 80 до 90 °C? Все значения, полученные из 80 °C, устарели — но солвер не хранит карту зависимостей. Сначала выполняется фаза «забывания»: новая спецификация на время притворяется неизвестной, старые рассчитанные значения скрыты, аппараты решаются «вхолостую» и подтверждают только то, что от изменённой величины не зависит; остальное стирается. В лаборатории (сценарий «+ смена T выхода 80→90») при этом стирается 23 величины в шести проходах; затем новая спецификация записывается, и обычный расчёт заново заполняет схему. Результат — тепловая нагрузка 489 450 кДж/ч — в точности совпадает с расчётом «с чистого листа» при 90 °C. Подробнее — в разделе 5.5.

4.5. Рецикл: оценка на потоке и итерации

Если выход схемы возвращается на вход, возникает замкнутый круг: чтобы посчитать смеситель, нужно знать рецикл, а чтобы знать рецикл — посчитать смеситель. Волна распространения упирается в неизвестное. В Sim42 специального аппарата «Recycle» нет: пользователь ставит оценку (оператор ~=) в любой порт — он становится точкой разрыва^[1]. Оценочные значения не передаются дальше как факт; направление вычислений определяется тем, во входном или в выходном порту поставлена оценка (рисунок 4).

Далее схема считается целиком, результат сравнивается с оценкой, ошибка нормируется масштабом переменной, и оценка обновляется: первый шаг — простая подстановка, затем метод Бройдена (квази-ньютоновский, без расчёта производных)^{[12, 13]}. Автоматического выбора разрывного потока в коде не обнаружено — это решает пользователь^[12].

Модель рецикла — учебная (смеситель, реактор, сепаратор, линия возврата), а не схема из Sim42. Параметры солвера соответствуют значениям по умолчанию: MaxNumIterations = 20, MaxRecycleStep = 0,05, MaxError = 0,0001^{[12, 13]}. При них простая подстановка не укладывается в 20 итераций и завершается сообщением MaxSolverIterExceeded, а метод Бройдена сходится за 9.

Раздел адресован разработчикам. Описание построено по аналитическому разбору исходного кода Sim42 v2.0.0.0 (каталоги sim/solver, sim/unitop; репозиторий^[7]) и тексту диссертации^[1]; анализ статический, программа не запускалась^[13]. Имена классов и методов (BasicProperty, Flowsheet.InnerSolve и т. п.) даны по коду; номера строк в источнике могут смещаться. В режиме «Специалист» подразделы свёрнуты; режим «Разработчик» раскрывает их.

InnerSolve():
    consistencyErrorStack.clear()
    iter = 0; jacobian = None
    while iter < maxIter:
        iter += 1
        SolverForget()
        if hold: return
        while op = PopSolveOp():                  # LIFO
            if op is not self and op.Ignored is None:
                op.BlockPush(1)
                try:
                    op.Solve()
                    for port in op.GetPorts(): port.UpdateConnection()
                finally:
                    op.BlockPush(0)
        ResetNewCalc()
        n = len(iterationStack)
        if n == 0: break                          # рециклов нет — готово
        maxError = max(prop.CalculateError(newValue))   # итерации + consistency
        if maxError < tol and no consistencyErrors: break
        if iter + 1 < maxIter: BroydenStep()      # см. 5.6
    if iter >= maxIter: raise SimError("MaxSolverIterExceeded")
    if consistencyErrorStack: raise ConsistencyError(...)

SolverForget():
    _isForgetting = 1
    while op = PopForgetOp():                     # LIFO
        if op is not self:
            op.BlockPush(1)
            try:
                if op.Ignored is None: op.Solve()  # "холостой" прогон
                for port in op.GetPorts():
                    if port.GetParentOp() == op: port.UpdateConnection()
                op.Forget()                        # стереть CALCULATED без NEW
                while port = PopResetCalcPort(): port.ResetNewCalc()
            finally:
                op.BlockPush(0)
    while port = PopResetFixedPort(): port.ResetNewFixed()
    _isForgetting = 0

def Solve(self):
    self.FlashAllPorts()                 # что можно — сфлэшить
    self.balance.DoBalance()             # мольный + энергетический баланс
    hadPressure = self.CalcDP()          # перепад давления
    didBalance = 0
    while self.FlashAllPorts():          # новый флэш дал новое знание?
        self.balance.DoBalance()         # → снова баланс
        didBalance = 1
    if not didBalance: self.balance.DoBalance()
    if not hadPressure: self.CalcDP()

6.1. Что отличает непоследовательно-модульный солвер

6.2. Ограничения и технический долг Sim42

По результатам статического анализа кода^[13] отмечены наблюдения, существенные и для самого алгоритма, и для возможного переноса.

6.3. Что стоит сохранить в любом наследнике

• Идея статусов значений и принцип «расчёт не перезаписывает известное, а проверяет» — простая модель защиты от переопределения.
• Инвалидация без графа зависимостей: экономит память и код, хорошо работает при интерактивной работе.
• Нормировка ошибок шкалами типов — единый механизм допусков.
• Слабая связность ядра с термодинамикой, интерфейсами и единицами.
• Протокол Solve(), позволяющий подключать аппараты любой сложности — от формулы до колонны^[13].

7.1. Для специалиста по моделированию

1. None — не ошибка. Это сигнал о недоопределённости: сосчитайте, каких спецификаций не хватает (например, давление на выходе холодильника). Ошибка — сообщение о противоречии или несходимости.
2. Противоречие — это лишняя спецификация. Если солвер сообщил о конфликте по порту, ищите значение, которое задано дважды разными путями (например, P на выходе и ΔP одновременно).
3. Порядок ввода не должен менять результат. Если меняет — подозревайте оценки рецикла, ограничение числа итераций или дефект, а не «особенность» программы. Сравнивайте «изменил A→B» с «сразу B».
4. Для рецикла задавайте оценку близко к ожидаемой и проверяйте масштаб переменных: слишком жёсткий допуск или далёкая оценка приводят к MaxSolverIterExceeded. Перечень несошедшихся портов показывает, где искать.
5. Отличайте ошибки термодинамики от ошибок солвера. Если одна и та же термодинамическая задача решается стабильно при ручном вводе, но «плывёт» в схеме с рециклом — проблема в солвере; если не решается вовсе — в термо-пакете^[12].
6. Нулевой расход — «почти ноль». Внутри он заменяется на 10⁻⁴⁰ ради избежания деления на ноль, это нужно учитывать при проверке нулевых потоков^{[12, 13]}.

7.2. Для разработчика: свойства, которые стоит проверять

Для солвера с динамическим порядком расчёта главная гарантия — детерминизм результата. Поэтому тестировать нужны метаморфические свойства, а не только значения^[13]:

• независимость от порядка ввода и от режима hold/go;
• «задал A → изменил на B» равно «сразу B» (проверка фазы «забывания»);
• идемпотентность Solve(): повторный вызов не меняет состояния;
• недоопределённость даёт None без ошибки, избыточная несовместимая спецификация — ConsistencyError;
• сходимость рецикла к одному решению при разных начальных оценках;
• модульные тесты ветвей SetValue/GetValue, UpdateConnection (оценочный порт, режим «забывания»), UpdateJacobian.

В репозитории Sim42 уже есть регрессионные сценарии (.tst — текстовые команды CLI: рециклы, колонны, реакторы, компрессоры); при миграции их следует прогнать как «золотой» набор — сначала на оригинале, затем на переносе — сравнивая результаты с допуском шкал^[13]. Именно эти свойства (кроме последнего) проверяют тесты лаборатории: набор из 17 автоматических тестов к двум моделям запускается встроенным тест-раннером Node.js.

Принцип непоследовательно-модульного солвера укладывается в несколько идей: значение хранится вместе со статусом происхождения; расчёт распространяет частичную информацию в обе стороны и никогда не перезаписывает известное, а лишь сверяет; порядок вычислений определяет стек аппаратов, получивших новое знание; устаревшее находится не обходом графа, а «холостым» прогоном; итерации нужны только там, где информация упирается в замкнутый контур, — и тогда применяется подстановка и метод Бройдена. Эти идеи делают поведение программы предсказуемым для специалиста: результат определяется набором спецификаций, а не историей ввода.

История проекта показывает, что открытый Sim42 создан в среде калгарийской школы моделирования и при участии автора симулятора Hysim; это объясняет, почему описанный принцип совпадает по классу с принципом HYSYS. Одновременно важно не выдавать гипотезу за факт: прямой связи кода Sim42 и HYSYS открытые источники не подтверждают. Интерактивная лаборатория, сопровождающая статью, — самостоятельная реализация идей, проверенная автоматическими тестами; она позволяет увидеть механизм на конкретных числах и проверить свойства, которые обычно принимают на веру.

• Анализ кода — статический. Исходный код Sim42 (Python 2, бинарные зависимости) не запускался; описание поведения построено по тексту кода и встроенному руководству^{[12, 13]}. Подробно разобраны Variables.py, Ports.py (ключевые методы), Flowsheet.py, EquationSolver.py, Controller.py, UnitOperations.py, часть Balance.py, Mixer.py, Stream.py, Equation.py, Heater.py (HeaterCooler) и структура Tower.Solve. Не разбирались построчно: тела DoMoleBalance/DoVolBalance, SolveFlowMatrix, Heater.HeatExchanger, KineticReactor, PipeSegment, GUI и веб-интерфейс.
• Лаборатория — упрощённая модель. Рабочее тело — вода; потоки показаны одной карточкой; рецикл — учебная схема, а не пример из Sim42. Корректность проверена тестами на свойства алгоритма, но не сверкой с оригинальным кодом (он не запускался).
• История — по открытым источникам. Первоисточники — страницы sim42.org из веб-архива (2002–2004) и диссертация; профили участников — их самоописания. Не найдены: оригинал работы Svrcek et al. (1984); независимое подтверждение роли Morris в Hyprotech; судьба проекта после 2004 года (репозиторий создан энтузиастом в 2011 году, активное сопровождение не обнаружено).
• Лицензионные замечания носят справочный характер и не являются юридическим заключением. Рисунки из диссертации не воспроизводятся (на титульном листе указан запрет воспроизведения), схемы выполнены заново; два снимка веб-архива приведены уменьшенными как историческая иллюстрация со ссылкой на источник.
• Связь с HYSYS — гипотеза на уровне людей и идей (раздел 2.3), а не утверждение о коде.