Оценка задержки
Частой особенностью систем управления является наличие задержек вследствие времени, требуемого для обработки или передачи по сети, либо другого физического явления. Внутренние задержки системы в условиях дискретного времени могут добавлять большое количество состояний: для представления задержки длительностью секунд в дискретной системе со временем дискретизации требуется переменных состояния. Распространенным особым случаем является ситуация, когда задержка происходит либо на входе, либо на выходе системы, например из-за задержек связи. Зачастую полезно оценивать эту задержку перед оценкой модели, так как это позволяет сократить количество оцениваемых параметров. Ниже генерируется набор данных с большой входной задержкой и рассматривается, как можно оценить .
Оценка задержки
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase
τ = 2.0 # Задержка в секундах
Ts = 0.1 # Время дискретизации
P = c2d(tf(1, [1, 0.5, 1])*delay(τ), Ts) # Динамика описывается с помощью простой системы второго порядка с входной задержкой
u = sin.(0.1 .* (0:Ts:30).^2) # Интересный входной сигнал
res = lsim(P, u', x0=[0.5; 0; zeros(20)])
d = iddata(res)
plot(d)
Если рассмотреть эту (дискретную) систему, можно увидеть, что у нее на самом деле многомерное состояние.
P.nx
22
Простой непараметрический способ определения наличия задержки в системе — по взаимной корреляции между входом и выходом. Для этого можно воспользоваться функцией crosscorplot
:
crosscorplot(d)
График показывает, что между входом и выходом имеется незначительная (ниже пунктирных линий) корреляция для запаздываний короче 2 секунд, что в точности соответствует нашей задержке .
Другой способ определения задержки — оценка импульсной характеристики системы.
impulseestplot(d, 60; λ=0.5, lab="Estimated", title="Impulse response of delay system")
plot!(impulse(P, 6), lab="True", framestyle=:zerolines)
Оценка импульсных характеристик при наличии таких больших задержек является численно сложной задачей, поэтому для получения приемлемого результата потребовалась регуляризация . Таким образом, если набор данных небольшой, а предполагаемая задержка велика, рекомендуется использовать метод на основе взаимной корреляции, ведь на практике, если ожидаемая импульсная характеристика неизвестна, подстроить параметр регуляризации невозможно. Однако при небольших задержках и больших наборах данных метод на основе импульсной характеристики работает достаточно неплохо. Ниже в целях демонстрации показана оцененная импульсная характеристика для набора данных гораздо большего размера.
u2 = sin.(0.01 .* (0:Ts:300).^2) .+ randn.() # Интересный длительный входной сигнал с некоторым шумом
res2 = lsim(P, u2')
d2 = iddata(res2)
impulseestplot(d2, 200; λ=0.0, lab="Estimated", title="Impulse response with large dataset")
plot!(impulse(P, 20), lab="True", framestyle=:zerolines)
Третьим способом, возможно менее элегантным, является использование метода подбора модели для определения задержки. Можно предположить, что модели с порядком, меньшим задержки, будут показывать достаточно плохой результат. Это должно быть очевидно при попытке использовать модели различного порядка. Ниже используется функция find_nanb
, которая пытается определить подходящий порядок для модели arx
.
find_nanb(d, 2:8, 10:30, xrotation=90, size=(800, 300), margin=5Plots.mm, xtickfont=6)
График показывает, что информационный критерий Акаике (AIC) оказывается минимальным при достижении nb
значения 22, что равно количеству выборок задержки + количество параметров в числителе для системы без задержки.
c2d(tf(1, [1, 0.5, 1]), Ts)
TransferFunction{ControlSystemsBase.Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}}
0.004913614996952642z + 0.004832374721150279
-------------------------------------------------
1.0z^2 - 1.9414834347826107z + 0.9512294245007139
Sample Time: 0.1 (seconds)
Discrete-time transfer function model
Учет задержки при оценке модели
Если известно, что имеется задержка, ее необходимо как-то учитывать при оценке модели. Самый примитивный способ — просто выбрать модель с достаточно высоким порядком и на этом все. Однако такой подход сопряжен с численными проблемами и, как правило, не рекомендуется. У некоторых методов есть специальный именованный аргумент inputdelay
. Он позволяет указать известные входные задержки, которые будут учитываться методом автоматически. Если у метода нет такого аргумента, всегда можно устранить задержку, подвергнув данные предварительной обработке. Далее будет показано, как это делается.
τ_samples = 20 # Количество выборок в задержке
d2 = iddata(d.y[:, τ_samples+1:end], d.u[:, 1:end-τ_samples], d.Ts)
crosscorplot(d2)
На графике взаимной корреляции теперь видно, что задержка устранена. Имейте в виду, что задержка в одну выборку обычно является вполне естественной: действительно, реакция на выходе физической системы вряд ли может происходить одновременно с событием на входе. После оценки модели с набором данных, из которого устранена задержка, необходимо вернуть задержку в оцененную модель, например, так:
P̂ = subspaceid(d2, focus=:simulation) # Оценка модели без задержки
tf(d2c(P̂)) # Соответствует ли это исходной непрерывной системе без задержки?
TransferFunction{ControlSystemsBase.Continuous, ControlSystemsBase.SisoRational{Float64}}
1.9105787402143033e-16s^2 + 0.036935897478138344s + 0.9989182628886182
----------------------------------------------------------------------
1.0s^2 + 0.49999999999999306s + 1.000000000000001
Continuous-time transfer function model
Приведенная выше оцененная модель должна быть очень близка к системе tf(1, [1, 0.5, 1])
(в числителе предполагаются небольшие члены более высокого порядка). Чтобы ввести задержку в эту модель, сделаем следующее:
P̂τ = P̂*delay(τ, Ts) # Обратите внимание, что этот метод отличается от метода delay(τ, Ts)*P̂, который добавляет выходную задержку вместо входной.
bodeplot([P, P̂τ], lab=["True system" "" "Estimated system" ""])
Оцененная модель теперь должна весьма точно соответствовать подлинной системе, включая сильный спад по фазе при более высоких частотах вследствие задержки.
!!! warning "Internal delays" Если в системе есть внутренние задержки, на графике взаимной корреляции это может выглядеть так, как будто в системе есть задержка между входом и выходом, но в таком случае выполнять смещение данных, как было показано выше, может быть неуместно. При наличии внутренних задержек используйте вместо этого количество выборок в оцененной задержке в качестве нижней границы для модели. Примером внутренней задержки может служить замыкание контура управления вокруг канала задержки, показанное ниже.
Внутренние задержки
В условиях дискретного времени задержка в одну выборку на входе или выходе выглядит как полюс в начальной точке. Однако, если задержка происходит внутри системы, обнаружить ее не так легко. Ниже производится замыкание контура вокруг приведенной выше системы с помощью ПИ-регулятора и строится новый набор данных, в котором вход связан с ПИ-регулятором, а не является непосредственно входом объекта.
C = pid(0.1, 0.5; Ts)
L = P*C
G = feedback(L)
plot(nyquistplot(L), plot(step(G, 50)))
ref = sign.(sin.(0.02 .* (0:Ts:50).^2)) # Интересный эталонный сигнал
res = lsim(G, ref')
dG = iddata(res)
plot(plot(dG), crosscorplot(dG))
Как видно на графике взаимной корреляции, реакция на выходе по-прежнему занимает более 2 секунд, но, исходя из графика подбора модели, для получения хорошего результата необходимо включить порядки до 23 и выше.
find_nanb(dG, 3:30, 5:30, xrotation=90, size=(800, 300), margin=5Plots.mm, xtickfont=6)
Можно также изучить схемы полюсов и нулей систем с разомкнутым и замкнутым контуром.
plot(
pzmap(P, title="Open-loop system"),
pzmap(G, title="Closed-loop system"),
plot_title="Pole-zero maps",
layout=(1, 2),
ratio=:equal,
size=(800, 400),
xlims=(-1.1, 1.1),
ylims=(-1.1, 1.1),
)
Из-за наличия обратной связи все полюсы задержек в исходной точке устранены (на графике слева имеется несколько полюсов, наложенных друг на друга).
В этом случае для точной регистрации процессов в системе необходимо оценить модель достаточно высокого порядка (23). Если сделать это, можно увидеть, что оцененная модель очень хорошо соответствует реальной системе (так как не было введено никаких возмещений).
model = subspaceid(dG, G.nx)
simplot(model, dG, zeros(model.nx))
На практике оценка модели такого высокого порядка может представлять сложность. Поэтому может быть целесообразно попытаться напрямую оценить систему с открытым контуром более низкого порядка , правильно учтя задержку, которая в этом случае будет происходить на входе, а не внутри системы.
При использовании одного из методов, поддерживающих именованный аргумент inputdelay
, можно определить, возможно ли аппроксимировать систему с замкнутым контуром с помощью модели, имеющей задержку только на входе.
model2 = arx(dG, 4, 4, inputdelay=τ_samples)
plot(
simplot(model2, dG),
pzmap(
model2,
ratio=:equal,
xlims=(-1.1, 1.1),
ylims=(-1.1, 1.1),
title="Estimated ARX model with input delay",
)
)
Мы видим, что модель порядка 3 (4 свободных параметра в знаменателе) в сочетании с заданной входной задержкой в 20 выборок на самом деле дает хорошую аппроксимацию системы с замкнутым контуром.
Если сравнить графики Боде реальной системы с замкнутым контуром с двумя указанными моделями, окажется, что они имеют высокую степень соответствия.
bodeplot([G, model, model2])
Таким образом, систему с внутренними задержками можно аппроксимировать с помощью модели, имеющей задержку только на входе.
Для полноты картины рассмотрим пример, в котором такой возможности нет. Систему в приведенном ниже примере можно представить как эхокамеру, в которой входной сигнал проходит через резонансный канал, прежде чем попасть на выход. 60 % выходной энергии поступает обратно на вход через тот же канал (эхо), из-за чего возникает интересная импульсная характеристика.
ref = sign.(sin.(0.02 .* (0:Ts:100).^2)) # Интересный эталонный сигнал
Pc = feedback(tf(100, [1, 10, 100]), -tf(60, [1, 10, 100])*delay(τ)) # 60 % выходной энергии подается обратно на вход с задержкой в 2 секунды (наподобие эха)
Pd = c2d(Pc, Ts)
res = lsim(Pd, ref')
decho = iddata(res)
plot(bodeplot(Pd, lab=""), pzmap(Pd), plot(impulse(Pd, 10), title="Impulse response"), plot(decho))
На графике подбора модели ниже видно, что для достижения хорошего соответствия требуются порядки модели до 24.
find_nanb(decho, 3:30, 5:30, xrotation=90, size=(800, 300), margin=5Plots.mm, xtickfont=6)
Попытка оценить модель 4-го порядка с входной задержкой в 20 выборок на этот раз совершенно не срабатывает, однако аппроксимация модели 24-го порядка дает результат, а аппроксимация модели 4-го порядка с использованием идентификации подпространства с долгим внутренним горизонтом прогнозирования также оказывается достаточно неплохой.
model3 = arx(decho, 4, 4, inputdelay=τ_samples)
model4 = subspaceid(decho, 24)
model5 = subspaceid(decho, 4, r=50)
figsim = simplot(model3, decho, zeros(ss(model3).nx), sysname="ARX 4")
simplot!(model4, decho, zeros(model4.nx), ploty=false, plotu=false, sysname="Subspace 24")
simplot!(model5, decho, zeros(model5.nx), ploty=false, plotu=false, sysname="Subspace 4")
Имейте в виду, что при таком моделировании не вводится никаких возмущений, а оценка моделей 24-го порядка на практике, скорее всего, будет непростой задачей.