MIMO LQR/LQG: линейно квадратичный и линейно гауссовский регуляторы с практическим примером и кодом (Часть 1) / Хабр

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

Цель публикации — демонстрация методов оптимального управления (LQR) и стохастической фильтрации (LQG/Kalman Filter) для решения задачи точного поддержания режима в условиях взаимосвязи физических параметров (температура и давление в замкнутом объеме) и зашумленных измерений.

Проект реализован на языке Python в парадигме Model-Based Design, разделяющей физику процесса, модель управления и среду моделирования.

Он включает постановку задачи, описание решения и его программную реализацию в виде python-пакета, адаптированного для работы в Google Colab. Среда автоматически сохраняет все артефакты каждого цикла моделирования — конфигурацию и результаты.

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

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

Ссылка на полный код проекта в конце статьи.

Оглавление

Зачем от SISO PID переходить к MIMO LQR/LQG
Постановка задачи
Формализация постановки
Решение: синтез регуляторов
Сравнение LQR и LQG по результатам моделирования
Рекомендации по применению модели (симуляция)

Зачем от SISO PID переходить к MIMO LQR/LQG

Классический PID-регулятор - инструмент эмпирического типа с одним входом и одним выходом (SISO), что не вполне соответствует природе обычно многосвязных (MIMO) объектов управления.
Проблема обходится применением комбинаций и каскадов из PID-регуляторов. Однако это лишь обходное решение. Децентрализация не способна явно компенсировать перекрестные связи между каналами, характерные для MIMO-объектов, и лишь косвенно подавляет их как внешние возмущения. Вдобавок, все это требует сложной настройки и часто обладает недостаточной устойчивостью.
Иногда комбинации и каскады оборачиваются внешней моделью, которая координирует работу PID или развязывает перекрестные связи. Но это уже выход за рамки классического PID-управления.
В противовес, линейно квадратичный LQR и линейно гауссовский LQG регуляторы решают проблему архитектурно "внизу" - на уровне регулятора, а не надсистемы.

Типовые области применения регуляторов

PID
SISO медленные или умеренно - быстрые технологические объекты с хорошо разделёнными контурами. Массовая промышленная автоматизация, когда перекрёстные связи слабы или подавлены конструктивно. Преимущественно PI архитектура.
LQR
MIMO системы, где критически важен учет перекрестных связей между каналами управления. Применяется там, где необходимо оптимизированное решение задачи баланса между качеством регулирования и затратами энергии, будь то быстрая робототехника, автопилоты, системы стабилизации или сложные технологические процессы.
LQG
Те же классы задач, что LQR, но при наличии существенного шума измерений и необходимости оптимальной фильтрации: управление полетом, навигация, управление сложными технологическими установками с зашумлёнными датчиками.
MPC (Model Predictive Control)
Системы со сложными ограничениями на управление и на переменные состояния. Позволяет явно учитывать ограничения и прогнозы поведения системы, но требует больших вычислительных ресурсов.
APC (Advanced Process Control)
Надстройка над базовыми PID-регуляторами для улучшения качества и экономических показателей при модернизации существующей архитектуры управления. Могут внутри содержать элементы MPC.

Первые три типа моделей определены достаточно компактно. Два последних - скорее не типы, а классы моделей.

Решаемая задача требует точного и быстрого двухканального управления параметрами, которые в замкнутом объеме автоклава вблизи рабочей точки находятся в термодинамической взаимосвязи, а измерения параметров зашумлены.
Пара PID регуляторов в таких условиях взаимно конкурирует, требует сложной настройки, неустойчива к смене рабочей точки (ее положение зависит от вида продукта). Управление имеет осциллирующий характер что приводит к повышенному износу исполнительных механизмов и перерасходу энергии.

В отличие от них, LQR/LQG регуляторы изначально адаптированы к особенностям задачи и позволяют синтезировать управление:

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

Вариант постановки задачи

Состав исходных данных

Уставка по температуре в диапазоне 118-127 град.
Уставка по давлению с учетом наддува над равновесным давлением, бар. Величина наддува в диапазоне 0.8-1.4 бар.
Начальное отклонение температуры от уставки, град.
Начальное отклонение давления от уставки, бар.
Лимиты по параметрам нагрева и наддува в долях от максимальных значений на инженерной шкале [0, 10].
Параметры производительности приводов по нагреву и наддуву, подобранные под шкалу [0, 10], единица сигнала на (град/с и бар/с).
Шумы измерений температуры и давления в стандартных отклонениях, (град и бар).
Инерция фазового перехода паровоздушной среды tau_phase, с.
Параметры пассивной динамики объекта:
- Инерция тепловых потерь tau_T, с.
- Инерция потерь давления tau_leak, с.

Примечание. Диапазоны в которых задаются уставки по температуре и давлению определяются спецификой прикладной задачи.

Файл конфигурации config.py содержит набор значений параметров для цикла моделирования.

Требования к системе управления (пример)

Требуется синтезировать MIMO-регулятор LQR/LQG, который обеспечивает переход объекта управления из начального состояния к заданным уставкам по температуре и давлению с соблюдением следующих критериев качества:

Статическая точность: в установившемся режиме отклонение от уставки не должно превышать ±0.2°C для температуры и ±0.1 бар для давления.
Качество переходного процесса:
- Перерегулирование: должно отсутствовать или быть минимальным (не более 2%). Целевой характер процесса — апериодический.
- Время установления: вхождение в 2%-зону от уставки не позднее 40 секунд.
Эффективность управления:
- Управляющие сигналы должны находиться в пределах [0, 10] без насыщения.
- Необходимо обеспечить плавность управления для минимизации износа приводов.
Робастность и развязка:
- Эффективное подавление измерительных шумов (фильтрация).
- Минимизация перекрестного влияния каналов управления (развязка управления температурой и давлением).
Оптимизация критерия затрат:
- Алгоритм должен минимизировать квадратичный функционал качества, связывающий точность и затраты на управление , с возможностью явной настройки их баланса.
Сравнительная эффективность:
- Интегральные показатели качества (ISE, IAE, ITAE) синтезированного регулятора должны превосходить показатели базового PID-регулирования.

Формальная постановка задачи

Для последовательного решения задачи разделим процесс синтеза на два этапа:

детерминированный LQR и
стохастический LQG.

Введем векторы состояния и управления для MIMO-системы ():

Вектор состояния (ошибки регулирования):
$x(t) = \begin{bmatrix} T(t) - T_{set} \\ P(t) - P_{set} \end{bmatrix}$
Вектор управления (сигналы на приводы):
$u(t) = \begin{bmatrix} u_{heat}(t) \\ u_{valve}(t) \end{bmatrix}$

Детерминированный LQR

Задача:
Найти закон управления для идеализированной системы, где все параметры состояния доступны для измерения и отсутствуют внешние возмущения.

Модель объекта (детерминированная):

$\dot{x}(t) = A x(t) + B u(t)$

Матрица описывает внутреннюю динамику объекта (термодинамику и перекрестные связи), а матрица — эффективность управляющих воздействий.

Критерий оптимизации и его решение:

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

$J_{LQR} = \int_{0}^{\infty} \left( x^T Q x + u^T R u \right) dt \to \min$

Где:

Матрица (state-cost matrix) определяет "штраф" за отклонение состояния от целевого (нулевого). Её диагональные элементы позволяют задать приоритеты: например, можно сделать поддержание температуры более важным, чем поддержание давления.
Матрица (control-cost matrix) определяет "штраф" за расход управляющего ресурса . Чем выше значения в , тем более "экономным" и плавным будет управление.

Такой постановке соответствует регулятор, оптимальный в идеальных условиях, но не учитывающий реалий зашумленных датчиков и модельных возмущений.

Система с зашумленными измерениями и возмущениями модели LQG

Задача:
Адаптировать управление к реальным условиям, где измерения зашумлены, а в системе присутствуют стохастические возмущения.

Модель объекта (стохастическая):

$\begin{cases}\dot{x}(t) = A x(t) + B u(t) + w(t) \\y(t) = C x(t) + v(t)\end{cases}$ .

Где:

: матрица выхода является единичной, так как мы напрямую измеряем оба состояния системы (температуру и давление ).
$w(t) \sim N(0, Q_w)$ : процессный гауссовский шум с ковариационной матрицей . Он описывает неучтенные флуктуации процесса (неравномерность теплоотдачи, микро-утечки).
$v(t) \sim N(0, R_v)$ : измерительный гауссовский шум датчиков с ковариационной матрицей .
: реальные показания датчиков (доступная информация).

Задача оценивания (Фильтр Калмана):

Так как точное состояние нам неизвестно из-за шумов, необходимо восстановить его оценку $\hat{x}(t)$ . Для этого требуется найти матрицу усиления наблюдателя , минимизирующую ковариацию ошибки оценивания :

$P_{err} = \mathbb{E}[(x - \hat{x})(x - \hat{x})^T] \to \min$

По аналогии с LQR итоговый закон управления для LQG: $u(t) = -K \hat{x}(t)$ с заменой состояния на его оценку.

Такой постановке соответствует регулятор, учитывающий реалии зашумленных датчиков и модельных возмущений.

Решение задачи

Построение модели

(модули проекта thermodynamics.py и model.py)

«В сущности, все модели неправильны, но некоторые полезны».
— Джордж Э. П. Бокс (статистик)

Термодинамические процессы в автоклаве объективно сложны: нагрев осуществляется подачей перегретого пара (~152°C, 5 бар), еще непрогретый продукт орошается через систему форсунок, а вентиляторы перемешивают смесь пара, воздуха и воды. Создание «цифрового двойника», учитывающего гидродинамику струй, пара и локальные зоны перегрева/конденсации, сделало бы модель громоздкой и непригодной для синтеза регулятора.

Мы зададим упрощенную структуру, которая захватывает главную динамику объекта, игнорируя высокочастотные шумы.

Матрица (внутренняя динамика)

Матрица описывает эволюцию состояния системы ( $\dot{x} = Ax$ ) пр�� отсутствии управления. В нашем случае она связывает медленные тепловые процессы и быстрые фазовые переходы.

Для построения матрицы мы используем параметры инерции, заданные в конфигурации проекта (config.py):

tau_T = 3600.0 с : инерция остывания
tau_leak = 900.0 с: инерция утечек
tau_phase = 10.0 с: инерция фазового перехода / восстановления $P_{sat}$

$A = \begin{bmatrix}-1/\tau_T & 0 \\K_{PT}/\tau_{phase} & -(1/\tau_{phase} + 1/\tau_{leak})\end{bmatrix}$

Параметры инерции относятся к оцениваемым параметрам.

Физический смысл элементов:

$a_{11} = -1/\tau_{T}$ :
описывает пассивное остывание автоклава. При $\tau_T = 3600$ система теряет тепло очень медленно. Это "медленная" мода системы.
$a_{12} \approx 0$ :
влиянием давления на температуру (адиабатическое сжатие/расширение) в данном диапазоне пренебрегаем, так как оно несоизмеримо меньше теплоемкости автоклава, загруженного продукцией.
$a_{21} = K_{PT} / \tau_{phase}$ (термодинамическая связь):
самый важный элемент. Он показывает, как температура "тянет" за собой давление.
- Коэффициент $K_{PT}$ (получаемый из уравнения Антуана) определяет, сколько бар добавится при росте на 1°C (производная в рабочей точке).
- Делитель $\tau_{phase}$ делает этот процесс быстрым. Давление "догоняет" температуру насыщения всего за ~10-30 секунд. Именно этот элемент задает перекрестную связь ( $T \to P$ ).
$a_{22} = -(1/\tau_{phase} + 1/\tau_{leak})$ :
скорость изменения давления. Она складывается из двух процессов:
- Медленного падения из-за утечек ( $1/\tau_{leak}$ ).
- Быстрого стремления к равновесию ( $1/\tau_{phase}$ ).

Давление — это "быстрая" мода системы: при фиксированной температуре любые отклонения давления релаксируют к равновесному значению $P_{sat}(T)$ на порядки быстрее, чем происходит изменение температуры.

Матрица (эффективность управления)

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

$B = \begin{bmatrix}k_{heat} & 0 \\0 & k_{valve}\end{bmatrix}$

Физический смысл и значения элементов:

$b_{11} = k_{heat}$ : коэффициент эффективности нагрева. Показывает, что 1 услов��ая единица сигнала управления нагревом приводит к росту температуры на 0.4 °C/с.
$b_{22} = k_{valve}$ : коэффициент эффективности наддува. 1 условная единица сигнала управления давлением приводит к его росту на 0.4 бар/с.
$b_{21}$ и $b_{12}$ равны 0: это модельное допущение означает, что мы считаем привод нагрева и привод наддува полностью развязанными. Такое допущение позволяет в чистом виде оценить способность LQR-регулятора компенсировать внутреннюю термодинамическую связь ( $T \to P$ ), заданную в матрице элементом $a_{21} = K_{PT}/\tau_{phase}.$

Дискретизация

Поскольку модель матриц и непрерывна, а регулятор работает в дискретном времени, требуется дискретизация модели. В модуле model.py используется метод Zero-Order Hold (ZOH) из scipy.signal.cont2discrete.

# Дискретизация Zero-Order Hold из scipy.signal.cont2discrete
    C_dummy = np.eye(2)
    D_dummy = np.zeros((2, 2))

    sys_d = cont2discrete((A, B, C_dummy, D_dummy), dt, method='zoh')
    Ad = sys_d[0]
    Bd = sys_d[1]

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

Для задач управления и фильтрации ZOH (Zero-Order Hold) де-факто является стандартом, так как он математически точно моделирует работу реального цифрового контроллера, где управляющий сигнал удерживается постоянным в течение всего такта (ступенчатый вход).

Использование здесь, например, метода Эйлера опасно тем, что вносимая им ошибка дискретизации может стать соизмеримой с малыми невязками (инновациями), которые фильтр Калмана использует для коррекции. Это приведет к тому, что фильтр будет бороться не с реальными шумами измерений, а с паразитными артефактами модели, что разрушит точность оценки, особенно в жестких системах с разномасштабной динамикой (как наша связь $T \to P$ ).

Синтез оптимального регулятора LQR

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

$J_{LQR} = \int_{0}^{\infty} \left( x^T Q x + u^T R u \right) dt$

Выбор весовых матриц

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

Матрица состояния $Q = \text{diag}([5.0, 2.0])$
- $q_{11} = 5.0$ : высокий штраф за отклонение температуры. Это приоритетный параметр стерилизации.
- $q_{22} = 2.0$ : штраф за отклонение давления - ниже. Стабильность важна, но небольшие колебания допустимы.
Матрица управления $R = \text{diag}([3.0, 8.0])$
- $r_{11} = 3.0$ : стоимость нагрева умеренная. Разрешаем активное использование нагревателя для быстрого выхода на режим.
- $r_{22} = 8.0$ : высокая стоимость работы приводом давления. Мы хотим избежать "дерганой" работы пневматики.

Численное решение

Для нахождения матрицы решается алгебраическое уравнение Риккати (ARE). В проекте используется решатель (сольвер) из библиотеки scipy.linalg (функция solve_discrete_are), который реализует метод Шура для дискретных систем.

В результате расчета для модели с учетом динамики фазового перехода ( $\tau_{phase}$ ) получена следующая матрица коэффициентов обратной связи (/content/autoclave_control/runs/.../json/run_info.json):

json
"K_LQR_gain": [

      [
        1.2574922820043983,
        0.0068017683810473765
      ],
      [
        0.0026608100302188934,
        0.30403542363709946
      ]
    ],

Интерпретация:

Диагональ:
- Высокий коэффициент $k_{11} \approx 1.26$ обеспечивает быструю реакцию на отклонения по температуре.
- Коэффициент $k_{22} \approx 0.30$ значительно ниже, что отражает стратегию "экономии ресурса" при управлении процессом выравнивания давления.
Перекрестные связи:
- Появление ненулевых недиагональных элементов подтверждает MIMO-природу синтезированного регулятора. Он способен вносить превентивные поправки в один канал на основе информации о состоянии другого, компенсируя малые возмущения до их развития.

Синтез LQG (Linear Quadratic Gaussian)

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

С учетом этого необходим переход от LQR к LQG-регулятору.

Оптимальное стохастическое управление базируется на принципе разделения (separation principle). Он гласит, что задачи оценивания состояния и управления им можно решать независимо друг от друга, не теряя оптимальности всей системы. Соответственно, синтез LQG разбивается на две независимые задачи:

1. Детерминированная задача (LQR): Расчет матрицы оптимальных обратных связей , как описано выше.

2. Стохастическая задача (Фильтр Калмана): Построение оптимального наблюдателя состояния, который минимизирует ошибку оценки $\hat{x}(t)$ на фоне шумов.

Итоговый закон управления принимает вид:

$u(t) = -K \cdot \hat{x}(t)$

где $\hat{x}(t)$ — оценка состояния, формируемая фильтром Калмана.

Настройка фильтра Калмана

Для синтеза наблюдателя необходимо задать две ковариационные матрицы, характеризующие интенсивность шумов в системе.

Матрица шумов процесса

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

$Q_w = \begin{bmatrix} 10^{-4} & 0 \\ 0 & 10^{-4} \end{bmatrix}$

Матрица шумов измерений

Эта матрица характеризует точность датчиков. Элементы на диагонали соответствуют дисперсиям ошибок измерений $\sigma^2$ .
Исходя из паспортных характеристик датчиков промышленного уровня (стандартное отклонение для датчика температуры $\sigma_T \approx 0.1^\circ C$ , для датчика давления $\sigma_P \approx 0.05$ бар), матрица имеет вид:

$R_v = \begin{bmatrix} \sigma_T^2 & 0 \\ 0 & \sigma_P^2 \end{bmatrix} = \begin{bmatrix} 0.01 & 0 \\ 0 & 0.0025 \end{bmatrix}$

Фильтр Калмана в реальном времени корректирует предсказание модели на основе поступающих измерений $y_{k+1}$ :

Этап 1: Прогноз (Time Update).

Мы используем модель, чтобы шагнуть в будущее. Получаем прогнозную (априорную) оценку, обозначаемую символом со знаком «минус» ( $\hat{x}^-$ ):

$\hat{x}_{k+1}^- = A_d \hat{x}_k + B_d u_k$

Этап 2: Коррекция (Measurement Update).

В момент времени приходит свежее измерение $y_{k+1}$ . Мы сравниваем его с прогнозом и получаем финальную (апостериорную) оценку:

$\hat{x}_{k+1} = \hat{x}_{k+1}^- + L_k (y_{k+1} - C \hat{x}_{k+1}^-)$

( $y_{k+1} - C \hat{x}_{k+1}^-$ ) — это инновация (ошибка предсказания показания датчика)
К прогнозу (априору) $\hat{x}_{k+1}^-$ добавляем ошибку, взвешенную коэффициентом доверия ,

где — матрица усиления Калмана (в установившемся режиме становится константой: $L_k \to L = \text{const}$ ), а $\hat{x}_{k+1}^-$ — прогнозное значение.

Описанный двухэтапный алгоритм реализован в классе KalmanFilter2D модуля control.py :

class KalmanFilter2D:
    """
    Калмановский фильтр для 2D-состояния:
        x_{k+1} = Ad x_k + Bd u_k + w_k
        y_k     = Cd x_k + v_k
    где:
        x_k ∈ R^2 — состояние [T_dev, P_dev]
        u_k ∈ R^2 — управление [u_heat, u_press]
        y_k ∈ R^2 — измерения температуры и давления
        w_k ~ N(0, Qw), v_k ~ N(0, Rv)
    """

Он соответствует классическому дискретному линейному фильтру Калмана в рекурсивной форме.

Рекурсивность / нестационарность (Recursive / Time-Varying Gain) означает, что матрица ковариации self.P и коэффициент усиления Kk (соответствует в описании алгоритма) пересчитываются на каждом шаге.

Для нашей задачи с относительно низкой динамикой это может быть избыточно и излишне вычислительно затратно. Можно было обойтись стационарным (Steady-State) фильтром Калмана, где коэффициент усиления рассчитывается один раз заранее решением дискретного уравнения Риккати через solve_discrete_are (как для LQR) и замораживается. Но, в самый раз, для моделей с высокой динамикой.

Специфика вычисления коэффициента усиления в рекурсивном фильтре Калмана

В классе KalmanFilter2D нет явного вызова функции решения уравнения Риккати (solve_discrete_are) как это делалось при нахождении решения для LQR.

Уравнение Риккати здесь решается итеративно внутри методов predict и update .

1. В методе predict (это первая часть рекурсии Риккати - прогноз ковариации):

self.P = self.Ad @ self.P @ self.Ad.T + self.Qw

2. В методе update (это вторая часть рекурсии):

    # Вычисление S (инновационная ковариация)
    S = self.Cd @ self.P @ self.Cd.T + self.Rv
    
    # Вычисление K (усиление Калмана)
    Kk = self.P @ self.Cd.T @ np.linalg.inv(S)
    
    # Обновление P (коррекция ковариации)
    self.P = (np.eye(self.n) - Kk @ self.Cd) @ self.P

При многократном вычислении в цикле матрица self.P в итоге сойдется к тому же самому значению, которое выдала бы функция solve_discrete_are (для стационарного случая). Так реализовано "динамическое решение уравнения Риккати" в дискретном времени. Это позволяет фильтру корректно работать в начале (переходный процесс), когда мы еще мало знаем о системе (большая начальная ).

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

kf.reset(x0=np.zeros_like(x0), P0=10.0 * np.eye(len(x0)))

Выбор большого (относительно начальных отклонений $T_{dev} \sim 2.5^\circ C, P_{dev} \sim 1.0$ бар) значения обеспечивает:

Более быструю сх��димость: Фильтр с большим будет вычислять большой коэффициент усиления на первых шагах.
Робастность (устойчивость к ошибкам): Алгоритм становится нечувствительным к тому, какое именно значение задано для $(\hat{x}_0)$ . Фильтр все равно быстро сойдется благодаря большому .

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

Рисунок 2. Сходимость рекурсивного усиления к стационарному за 300 шагов симуляции в условиях различной стартовой неопределенности ().

Графики на рисунке 1 демонстрируют преимущество рекурсивного фильтра перед стационарным (нерекурсивным) в режиме "холодного старта".

Слева (Причина):

Чем выше начальная неопределенность (фиолетовая линия), тем выше стартовый коэффициент усиления $L \to 1$ . Фильтр "понимает", что его начальная оценка ненадежна, и максимально доверяет первым измерениям датчиков. При малых (желтая линия) фильтр излишне "самоуверен" и игнорирует новые данные, стартуя с низким усилением, близким к стационарному $L_{ss}$ .

Справа (Следствие):

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

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

Интерпретация рисунка 2:

Результаты подтверждают корректность рекурсивного алгоритма: независимо от выбора начальной неопределенности (от агрессивного 10.0 до консервативного 0.005), фильтр Калмана неизбежно сходится к одному и тому же стационарному решению $L_{ss}$ , п��лученному через уравнение Риккати (DARE). Ошибка сходимости пренебрежимо мала ( $\approx 10^{-6}$ ), что демонстрирует: влияет только на динамику переходного процесса, но не на финальную точность фильтра в установившемся режиме.

Особенность расчета интегральных метрик

В модуле metrics.py производится вычисление интегральных метрик производительности и метрик стоимости.

Для стохастических систем важно вычислять метрики производительности на основе истинного состояния объекта ( $x_{true}$ ), которое в реальности скрыто, но доступно в симуляции. Распространенная ошибка — использовать для метрик оценку состояния ( $\hat{x}_{est}$ ), которую выдает сам фильтр Калмана.

Это опасно тем, что создает иллюзию идеальной работы: фильтр может "думать", что он управляет прекрасно (ошибка $\hat{x}_{est}\to 0$ ), в то время как реальный объект ( $x_{true}$ ) может находиться далеко от уставки из-за неучтенных смещений или ошибок модели, которые фильтр просто "не видит". Метрика, посчитанная по оценке, оценит лишь "самоуверенность" фильтра, а не реальное качество управления.

Это принципиальный момент: конечной целью системы управления является поддержание фактического режима , а не стабилизация показаний датчиков. Расчет по $y_{meas}$ исказил бы метрику шумом, который невозможно устранить управлением, а расчет по $\hat{x}_{est}$ скрыл бы ошибки оценивания наблюдателя. Использование $x_{true}$ позволяет дать объективную оценку физического результат�� управления.

Сравнение LQR и LQG по итогам моделирования

Диаграммы, приведенные ниже демонстрируют, что LQG-регулятор показал способность эффективно сглаживать высокочастотный шум датчиков, предотвращая "дребезг" управляющих воздействий, который наблюдался бы при использовании «чистого» LQR на зашумленных данных. Это важно для продления ресурса исполнительных механизмов автоклава.

Рисунок 3. Динамика параметров и сигналов управления.

Рисунок 4. Фазовые траектории LQR и LQG.

Рисунок 5. Интегральные метрики производительности.

Рисунок 6. Метрики стоимости для LQR и LQG.

Проведем более детальное сравнение эффективности двух подходов: детерминированного LQR и стохастического LQG регуляторов.

1. Качество переходных процессов

Оба регулятора успешно справляются с задачей выхода на уставки по температуре ( $121^\circ \text{C}$ ) и давлению ( бар).

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

2. Работа с шумами и управляющие воздействия

Наиболее показательным является график управляющих сигналов.

LQR (идеальный): Демонстрирует гладкое экспоненциальное падение управляющего воздействия по мере приближения к цели.
LQG (реальный): Несмотря на сильный шум измерений (красные точки на графиках и ), управляющий сигнал (фиолетовая и бирюзовая линии) остается практически гладким, без высокочастотного "дребезга". Это подтверждает правильную настройку фильтра Калмана: регулятор реагирует на тренд, а не на мгновенный шум, что критически важно для долговечности исполнительных механизмов.

3. Анализ фазовых траекторий

Фазовый портрет в координатах отклонений ( $T_{dev}, P_{dev}$ ) демонстрирует стратегию движения системы к равновесию ().

Траектории LQR и LQG практически совпадают, двигаясь вдоль "хребта" оптимальности.
На графике LQG видно, как оценка (синие точки) сглаживает разброс истинного зашумленного состояния (красный пунктир), удерживая систему в узком коридоре управления.

В контексте фазового графика LQG красный пунктир ("истина") отображает реальное физическое состояние системы с учетом всех внутренних шумов процесса.
Стохастическая модель объекта :
$\dot{x}_{true} = Ax + Bu + w(t)$
Где — это процессный шум.
Истина: красный пунктир — это именно $x_{true}$ .
Это "честная" физика, на каждом шаге состояние возмущает случайный шум.
Измерения (то, что видит датчик): $y = x_{true} + v(t)$ . На фазовом графике они не показаны. Но на графиках динамики они есть — это облака красных точек.
Оценка - синие точки (то, что оценил фильтр): $\hat{x}$ . Это результат работы фильтра Калмана, который пытается угадать $x_{true}$ , глядя на зашумленный .
При дисперсии процессного шума $\sigma^2 =$ $10^{-4}$ дрожание красного пунктира микроскопическое и на общем плане фазового портрета ( $T_{dev} \sim 2.5^\circ C, P_{dev} \sim 1.0$ бар) оно практически неразличимо глазом — линия кажется гладкой.
Вывод. Для LQR и LQG "истина" имеет разную природу хотя и показана одинаковыми линиями.

4. Метрики качества

Количественное сравнение подтверждает ожидаемое небольшое снижение производительности LQG из-за работы в условиях неопределенности:

Функционал стоимости : полная стоимость для LQG () лишь на ~5% выше, че�� для идеального LQR (). Это "цена", которую мы платим за фильтрацию шума. Составляющие стоимости сбалансированы.
Ошибки регулирования (ISE, IAE, ITAE): Метрики ошибки для LQG закономерно выше (на 1-3%), чем для LQR. Например, интегральная квадратичная ошибка (ISE) по давлению составляет для LQG против для LQR. Разница минимальна, что говорит о высокой эффективности синтезированного наблюдателя.

Вывод: Реализованный LQG-регулятор обеспечивает качество управления, сопоставимое с теоретически оптимальным LQR, при этом успешно подавляя измерительный шум. Система демонстрирует робастность. Незначительный рост ошибки регулирования является приемлемой платой за стабильность исполнительных механизмов.

Что дальше: к астатическому LQG-регулятору

Текущая реализация LQG-регулятора демонстрирует высокую эффективность в стабилизации системы и подавлении шумов. Однако, как и любой пропорциональный регулятор по состоянию, он может допускать небольшую статическую ошибку при наличии неучтенных постоянных возмущений, таких как, например, изменение теплопотерь или дрейф датчиков.

Для устранения этого недостатка и достижения нулевой ошибки в установившемся режиме, возможным следующим шагом является модификация контроллера путем введения интегрального действия. Как показано в работе Tran \& Vu (2025) на примере MIMO-системы управления мобильным роботом, этого можно достичь через расширение вектора состояния. Вектор дополняется интегралом ошибки слежения $x_{int} = \int (r - y) dt$ , после чего новый регулятор синтезируется уже для расширенной системы $x_{aug} = [x, x_{int}]^T$ .

Это трансформирует классический LQR из PD-регулятора (учитывающего состояние и скорость изменения) в полноценный оптимальный MIMO PID-регулятор. В отличие от традиционного подхода с настройкой независимых PID-петель для температуры и давления, здесь интегральные составляющие каналов математически связаны друг с другом, а все коэффициенты усиления (P, I и D) настраиваются автоматически и согласованно через решение уравнения Риккати.

Кроме того, подход обладает и теоретическим обоснованием. Авторы статьи демонстрируют, что сама функция стоимости LQR может выступать в роли функции Ляпунова. Это позволяет доказать, что синтезированная система гарантированно сохраняет устойчивость и удерживает состояние в ограниченной области около нуля даже при наличии стохастических возмущений.

Такая модернизация сделает систему управления не только точной, но и доказательно робастной.

Взгляд из будущего (относительно настоящей статьи).
На Хабре опубликована статья с кодом и описанием синтеза астатического регулятора: Augmented LQR: расширяем пространство состояний, чтобы убрать статическую ошибку (Часть 2).

MIMO LQR/LQG: линейно квадратичный и линейно гауссовский регуляторы с практическим примером и кодом (Часть 1)