Идентификация замкнутого контура
В этом примере рассматривается эффективность различных алгоритмов идентификации применительно к данным замкнутого контура (когда входной сигнал системы производится контроллером с использованием обратной связи от выхода).
Мы возьмем очень простую систему с цветным выходным шумом и различными входными сигналами, образуемыми посредством обратной связи от выхода , где будет меняться от эксперимента к эксперименту.
Известно, что в отсутствии и при использовании простого регулятора идентифицируемость системы низкая, ведь
мы получаем замкнутую систему,
в которой и невозможно различить. Введение позволяет разрешить эту проблему.
В самом первом из приведенных ниже экспериментов демонстрируется эта проблема при отсутствии возбуждения через .
Сначала мы определим модель реальной системы, функцию, моделирующую данные и добавляющую цветной выходной шум, а также функцию, оценивающую три разные модели и строящую график их частотных характеристик. Мы рассмотрим три метода оценки.
-
arx, a prediction-error approach based on a least-squares estimate. -
Метод на основе подпространства
subspaceid, который, как известно, подвержен смещению при наличии обратной связи по выходу. -
Метод на основе погрешности прогнозирования (PEM)
newpem.
Методы ARX и PEM теоретически не подвержены смещению при наличии обратной связи по выходу [1], в то время как метод на основе подпространства подвержен. (Обратите внимание: метод на основе подпространства используется с целью получения начального предположения для итеративного алгоритма PEM.)
using ControlSystemsBase, ControlSystemIdentification, Plots
G = tf(1, [1, -0.9], 1) # Истинная система
function generate_data(u; T)
E = c2d(tf(1 / 100, [1, 0.01, 0.1]), 1) # Модель цветного шума
e = lsim(E, randn(1, T)).y # Реализация шума
function u_noise(x, t)
y = x .+ e[t] # Добавляем измеренный шум к состоянию для моделирования обратной связи шума измерения
u(y, t)
end
res = lsim(G, u_noise, 1:T, x0 = [1])
d = iddata(res)
d.y .+= e # Добавляем шум измерения к выходу, хранящемуся в объекте данных
d
end
function estimate_and_plot(d, nx=1; title, focus=:prediciton)
Gh1 = arx(d, 1, 1)
sys0 = subspaceid(d, nx; focus)
tf(sys0)
Gh2, _ = ControlSystemIdentification.newpem(d, nx; sys0, focus)
tf(Gh2)
figb = bodeplot(
[G, Gh1, sys0.sys, Gh2.sys];
ticks = :default,
title,
lab = ["True system" "ARX" "Subspace" "PEM"],
plotphase = false,
)
figd = plot(d)
plot(figb, figd)
endestimate_and_plot (generic function with 2 methods)В первом эксперименте нет эталонного возбуждения и при небольшом объеме данных ( ) оценки получаются крайне плохими.
L = 0.5 # Коэффициент усиления обратной связи u = -L*x
u = (x, t) -> -L * x
title = "-Lx"
estimate_and_plot(generate_data(u, T=80), title=title*", T=80")
При большем объеме данных оценки столь же плохие.
estimate_and_plot(generate_data(u, T=8000), title=title*", T=8000")
Это означает, что, если система определяется только шумом, оценить ее практически невозможно.
Теперь рассмотрим небольшое периодическое возбуждение
L = 0.5 # Коэффициент усиления обратной связи u = -L*x
u = (x, t) -> -L * x .+ 5sin(t)
title = "-Lx + 5sin(t)"
estimate_and_plot(generate_data(u, T=80), title=title*", T=80")
В данном случае все методы, кроме метода на основе подпространства, работают достаточно эффективно.
estimate_and_plot(generate_data(u, T=8000), title=title*", T=8000")
Увеличение объема данных не улучшает работу метода на основе подпространства.
При более сложном возбуждении (случайный шум белого спектра) все методы работают хорошо.
L = 0.5 # Коэффициент усиления обратной связи u = -L*x
u = (x, t) -> -L * x .+ 5randn()
title = "-Lx + 5randn()"
estimate_and_plot(generate_data(u, T=80), title=title*", T=80")
А при увеличении объема данных — еще лучше.
estimate_and_plot(generate_data(u, T=8000), title=title*", T=8000")
При сильной обратной связи, но слабом возбуждении результаты работы всех методов достаточно плохие, поэтому важно, чтобы энергия возбуждения была достаточной в сравнении с каналом обратной связи.
L = 1 # Коэффициент усиления обратной связи u = -L*x
u = (x, t) -> -L * x .+ 0.1randn()
title = "-Lx + 0.1randn()"
estimate_and_plot(generate_data(u, T=80), title=title*", T=80")
В данном случае можно попытаться повысить порядок модели для метода PEM и метода на основе подпространства, чтобы узнать, смогут ли они изучить модель шума (с двумя полюсами).
estimate_and_plot(generate_data(u, T=8000), 3, title=title*", T=8000")
Изучение модели шума иногда может давать достаточно неплохой результат, но требует больше данных. Извлечь изученную модель шума можно с помощью noise_model.
Определение наличия обратной связи
Иногда определить наличие обратной связи в наборе данных можно по взаимной корреляции между входом и выходом. В причинной системе не должно быть корреляции с отрицательными запаздываниями. Выходной сигнал в прямом смысле поступает обратно на вход, что приводит к обратной причинности:
L = 0.5 # Коэффициент усиления обратной связи u = -L*x
u = (x, t) -> -L * x .+ randn.()
title = "-Lx + 5sin(t)"
crosscorplot(generate_data(u, T=500), -5:10, m=:circle)
На этом графике явно видна значительная корреляция как с положительным, так и с отрицательным запаздыванием, что указывает на наличие обратной связи. Здесь используется статический P-регулятор, что дает корреляцию на один шаг назад во времени. С динамическим регулятором (например, ПИ-регулятором) эффект будет существеннее.
Если удалить обратную связь, получится следующее:
L = 0.0 # нет обратной связи
u = (x, t) -> -L * x .+ randn.()
title = "-Lx + 5sin(t)"
crosscorplot(generate_data(u, T=500), -5:10, m=:circle)
теперь корреляция для отрицательных задержек и нулевой задержки практически несущественная (ниже пунктирных линий).