Идентификация неустойчивых систем

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

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

Моделирование неустойчивой системы

В этом небольшом примере мы строим очень простую неустойчивую линейную систему:

где — управляющий вход, а — возмущение. Сначала мы моделируем систему со стабилизирующей обратной связью и сохраняем данные. Затем мы снова моделируем ту же систему с сохраненными данными в качестве входных, но с немного иным начальным условием. Интуитивно может казаться, что эти две модели должны быть одинаковыми: система одна и та же и входные данные одни и те же, лишь начальное условие немного поменялось.

using ControlSystemsBase, Plots, LinearAlgebra
sys = tf(1.0, [1, -1]) |> ss # Неустойчивая линейная система с одним состоянием
sys = c2d(sys, 0.01) # Дискретизируем

Q = R = I
L = lqr(sys, Q, R) # Проектируем контроллер стабилизации

Ts = 0.01 # Интервал выборки
t = 0:Ts:15 # Вектор времени
inp(x, t) = -L*x .+ sin(t) # Входная функция представляет собой линейную обратную связь + sin(t)

res = lsim(sys, inp, t)
res2 = lsim(sys, res.u, t, x0=[1e-6]) # Повторяем моделирование с теми же входными данными, но начальным условием, отличающимся на 1e-6
plot(res)
plot!(res2)

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

Эта система является экспоненциально устойчивой, пока больше 1. Если мы собираем экспериментальные данные для неустойчивой системы со стабилизирующей обратной связью, мы по сути собираем данные из устойчивой системы.

Во втором случае моделировалась исходная система с разомкнутым контуром:

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

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

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

Как можно исправить эту ситуацию к лучшему? При чистом моделировании особо ничего не поделаешь. Однако при моделировании системы в целях оценки ее параметров есть очень эффективный прием: Решением является введение в модель также и обратной связи по измерению. Иными словами, вместо моделирования неустойчивой системы мы производим оценку ее состояния с использованием доступных данных измерений. Средство оценки состояния, например фильтр Калмана, имеет внутреннюю обратную связь по измерению, которая корректирует отклонения между состоянием системы и данными измерений. Для установившегося фильтра Калмана в дискретном времени эта коррекция выглядит так:

То есть мы производим линейную коррекцию на основе разницы между прогнозируемым выходом и измеренным выходом . Подобная коррекция должна быть вам знакома — сравните ее с уравнением замкнутой системы управления с опорным узлом :

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

Оценка состояния для идентификации системы

Чтобы оценить параметры динамической системы с помощью средства оценки состояния, можно запустить средство оценки состояния, например фильтр Калмана, вперед во времени, сравнить оцененные состояния с измеренными данными и вычислить градиент погрешности относительно параметров. Именно так работает метод на основе погрешности измерения (PEM), реализованный в функции newpem. У этой функции есть именованный аргумент focus со значением по умолчанию focus = :prediction, которое, однако, можно изменить на focus = :simulation, чтобы отключить обратную связь по измерению. Однако для успешной идентификации неустойчивых систем обратная связь по измерению имеет ключевое значение и любой алгоритм моделирования без такой связи имеет тенденцию к отклонению.

Алгоритм PEM применяется для идентификации модели неустойчивой системы шара и балки в руководстве Шар и балка. Ближе к концу этого руководства эффективность оцениваемых моделей сравнивается с использованием все более длительных горизонтов прогнозирования, но уже при горизонте в 20 шагов эффективность прогнозирования существенно снижается вследствие неустойчивости.

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

Особенности метода на основе погрешности прогнозирования

По сути, метод PEM преобразует задачу минимизации потерь в плане производительности моделирования в задачу минимизации потерь в плане длительности прогнозирования ^[1] ^[2]. Это дает ряд преимуществ, два из которых показаны в следующем примере:

Потери зачастую проще оптимизировать.
Помимо точного симулятора, вы также получаете прогноз для системы.
С помощью PEM можно оценивать модели возмущений.

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

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

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

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

Начнем с определения модели маятника. Модель принимает параметр , соответствующий длине маятника.

using Plots, Statistics, DataInterpolations, LowLevelParticleFilters

Ts = 0.01 # Интервал выборки
tsteps = range(0, stop=20, step=Ts)

x0 = [0.0, 3.0] # Начальный угол и угловая скорость

function pendulum(x, u, p, t) # Динамика маятника
    g = 9.82 # Постоянная силы тяжести
    L = p isa Number ? p : p[1] # Длина маятника
    gL = g / L
    θ = x[1]
    dθ = x[2]
    [dθ
     -gL * sin(θ)]
end

Мы предполагаем, что истинная длина маятника равна , и генерируем данные для этой системы.

using LowLevelParticleFilters: rk4, rollout

discrete_pendulum = rk4(pendulum, Ts) # Дискретизируем непрерывную динамику с помощью RK4

function simulate(fun, p)
    x = rollout(fun, x0, tsteps, p; Ts)[1:end-1]
    y = first.(x) # Это имеющиеся данные для оценки параметров (измерение угла)
    x, y
end

x, y = simulate(discrete_pendulum, 1.0) # Выполняем моделирование при L = 1,0
plot(tsteps, y, title = "Pendulum simulation", label = "angle")

Мы также определяем функции, которые моделируют систему и вычисляют потери для заданного параметра p, соответствующего длине.

function simloss(p)
    x,yh = simulate(discrete_pendulum, p)
    yh .= abs2.(y .- yh)
    return mean(yh)
end

Далее характер потерь рассматривается как функция от длины маятника:

Ls = 0.01:0.01:2
simlosses = simloss.(Ls)
fig_loss = plot(Ls, simlosses,
    title  = "Loss landscape",
    xlabel = "Pendulum length",
    ylabel = "MSE loss",
    lab    = "Simulation loss"
)

Этот рисунок вызывает интерес: потери при истинном значении , конечно, равны 0, но для значений наклон в целом неверный! Более того, потери имеют колебательный характер, а это значит, что оптимизировать данную функцию крайне непросто и для схождения локального поиска к истинному значению потребуется очень хорошее начальное предположение. Обратите внимание: мы намеренно выбрали одномерный пример для наглядности, и, хотя одномерные задачи обычно решаются легко, тот же ход рассуждений применим и к более сложным задачам высокой размерности.

Теперь перейдем к определению модели предиктора. Наш предиктор будет очень простым: на каждом временном шаге будет вычисляться погрешность между смоделированным углом и измеренным углом . Доля этой погрешности будет использоваться для коррекции состояния маятника. Используемая коррекция является линейной и выглядит так: . Мы построили то, что обычно называется (линейным) наблюдателем. Фильтр Калмана — это линейный наблюдатель особого рода, где рассчитывается на основе статистической модели возмущений, действующих на систему. Чтобы не усложнять задачу, мы ограничимся простым наблюдателем с постоянным усилением.

Для подачи выборки данных в непрерывную модель мы используем интерполятор. Мы также определим новые функции: функцию predictor, которая содержит динамику маятника с коррекцией наблюдателя, функцию prediction, которая выполняет развертывание (во избежание путаницы мы не используем слово «моделирование»), и функцию потерь.

y_int = LinearInterpolation(y, tsteps; extrapolate=true)

function predictor(x, u, p, t)
    g = 9.82
    L, K, y = p # Длина маятника, коэффициент усиления наблюдателя и измерения
    gL = g / L
    θ = x[1]
    dθ = x[2]
    yt = y(t)
    e = yt - θ
    [dθ + K * e
    -gL * sin(θ)]
end

discrete_predictor = rk4(predictor, Ts)

function predloss(p)
    p_full = (p..., y_int)
    x, yh = simulate(discrete_predictor, p_full)
    yh .= abs2.(y .- yh)
    return mean(yh)
end

predlosses = map(Ls) do L
    K = 1.0 # Коэффициент усиления обратной связи с наблюдателем
    p = (L, K)
    predloss(p)
end

plot!(Ls, predlosses, lab = "Prediction loss")

VYmk4WHh8vl8oSEBIFAYBtlD0CHADVCAEBb4TielpaWl5dnNBrDw8MHDhxIJpM9HRQAjoJECAAAwKvBrVEAAABeDRIhAAAArwaJEAAAgFf7f3YxvFQG8C97AAAAAElFTkSuQmCC

Потери снова рассматриваются как функция от параметра, но на этот раз ситуация выглядит гораздо лучше. Потери не являются выпуклыми, но градиент имеет правильное направление на гораздо большем интервале. Здесь мы произвольно задали коэффициент усиления наблюдателя . Однако при оценке методом PEM оптимизатор обычно определяет этот параметр самостоятельно, что и происходит внутри newpem.

Может возникнуть вопрос, почему мы использовали коррекцию в форме , а не задали угол в модели равным результату измерения. Причин две.

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

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

Ссылки:

1. Ljung, Lennart. System identification---Theory for the user.

2. Larsson, Roger и др. Direct prediction-error identification of unstable nonlinear systems applied to flight test data.