Оценка передаточных функций

На этой странице описывается, как оценивать передаточные функции, которые иногда называют моделями ARX или ARMAX, то есть модели любой из следующих форм:

Данный пакет предназначен для оценки дискретных моделей, но их можно преобразовывать в непрерывные с помощью функции d2c из пакета ControlSystemsBase.jl.

Доступны следующие методы:

Функции

arx: Transfer-function estimation using least-squares fitting with a closed-form solution. Supports multiple datasets and custom estimator functions.
arma: Estimate an ARMA model (no control input).
ar: Estimate an AR model (no input).
arma_ssa Estimate an ARMA model with estimated noise as input (no control input).
plr: Transfer-function estimation (ARMAX model) using pseudo-linear regression. armax is an alias for this function. This method estimates a noise model as well.
arxar: Transfer-function estimation using generalized least-squares method. This method estimates a noise model as well.
getARXregressor/getARregressor: For low-level control over the estimation

Дополнительную справочную информацию см. в документации docstring.

Большинство методов для оценки передаточных функций работает только с системами SISO, SIMO и MISO. Для оценки систем MIMO можно прибегнуть к методам на основе пространства состояний и при необходимости преобразовать результат в передаточную функцию после оценки с помощью tf.

Пример использования:

N  = 2000     # Количество временных шагов
t  = 1:N
Δt = 1        # Интервал выборки
u  = randn(N) # Случайный входной управляющий сигнал
G  = tf(0.8, [1,-0.9], 1)
y  = lsim(G,u,t)[1][:]
yn = y
d  = iddata(y,u,Δt)

na,nb = 1,1   # Количество коэффициентов многочлена

Gls = arx(d,na,nb,stochastic=false) # присваиваем аргументу stochastic значение true, чтобы получить передаточную функцию MonteCarloMeasurements.Particles
@show Gls
# TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     0.8000000000000005
# --------------------------
# 1.0*z - 0.8999999999999997

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

e  = randn(N)
yn = y + e    # Сигнал измерения с шумом
d  = iddata(yn,u,Δt)

na,nb,nc = 1,1,1

Gls      = arx(d,na,nb, stochastic=true)     # Обычная оценка методом наименьших квадратов
Gtls     = arx(d,na,nb, estimator=tls)       # Оценка методом наименьших полных квадратов
Gwtls    = arx(d,na,nb, estimator=wtls_estimator(y,na,nb)) # Взвешенная оценка методом наименьших полных квадратов
Gplr, Gn = plr(d,na,nb,nc, initial_order=20) # Псевдолинейная регрессия
@show Gls; @show  Gtls; @show  Gwtls; @show  Gplr; @show  Gn;
# TransferFunction{ControlSystemsBase.SisoRational{MonteCarloMeasurements.Particles{Float64,500}}}
#     0.824 ± 0.029
# ---------------------
# 1.0*z - 0.713 ± 0.013

# Gtls = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     1.848908051191616
# -------------------------
# 1.0*z - 0.774385918070221

# Gwtls = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#    0.8180228878106678
# -------------------------
# 1.0*z - 0.891939152690534

# Gplr = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     0.8221837077656046
# --------------------------
# 1.0*z - 0.8896345125395438

# Gn = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     0.9347035105826179
# --------------------------
# 1.0*z - 0.8896345125395438

Теперь видно, что оценка с использованием стандартного метода наименьших квадратов имеет сильное смещение и является ложнодостоверной (обратите внимание на знак ± в коэффициентах передаточной функции). Обычный метод наименьших полных квадратов плохо работает в этом примере, так как не все переменные в регрессии содержат одинаковое количество шума. Взвешенным методом наименьших полных квадратов получается достаточно хорошо восстановить истинную модель. Псевдолинейная регрессия также показывает себя неплохо, одновременно оценивая модель шума. Вспомогательная функция wtls_estimator(y,na,nb) возвращает функцию, которая выполняет операцию wtls с использованием ковариационных матриц подходящего размера, зависящего от длины y и порядков модели. Взвешенная оценка методом наименьших полных квадратов реализована в пакете TotalLeastSquares.jl. Дополнительные сведения см. в интерактивных скриптах с примерами.

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

См. также описание функции arma, которая позволяет оценивать модели сигналов без входов.

Видео, посвященное использованию функции arx, доступно здесь:

Моделирование временных рядов

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

ar: Estimate an AR model (no input).
arma_ssa Estimate an ARMA model with estimated noise as input (no control input).

Генерация кода

Для генерации кода на C, например с целью моделирования системы или фильтрации через оцененную передаточную функцию, используйте пакет SymbolicControlSystems.jl.

API передаточных функций

# ControlSystemIdentification.arx — Function

Gtf = arx(d::AbstractIdData, na, nb; inputdelay = ones(Int, size(nb)), λ = 0, estimator=\, stochastic=false)

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

nb и na представляют количество коэффициентов многочленов числителя и знаменателя.

Входную задержку можно добавить с помощью inputdelay = d, что соответствует дополнительной задержке z^-d. inputdelay = 0 Приводит к получению прямой составляющей. Высший порядок многочлена B задается с помощью nb + inputdelay - 1. λ > 0 можно указать для регуляризации L₂. estimator по умолчанию \ (наименьшие квадраты), альтернативами являются estimator = tls для полной оценки по методу наименьших квадратов. arx(Δt,yn,u,na,nb, estimator=wtls_estimator(y,na,nb) потенциально более устойчив в условиях сильного шума измерений. Количество свободных параметров равно na+nb

stochastic: если задано значение true, возвращает передаточную функцию с неопределенными параметрами, представленными MonteCarloMeasurements.Particles.

Поддерживает оценку MISO, предоставляя iddata с матрицей u, с nb = [nb₁, nb₂…] и дополнительным inputdelay = [d₁, d₂…].

Эта функция поддерживает несколько наборов данных, предоставляемых в виде вектора объектов iddata.

# ControlSystemIdentification.ar — Function

ar(d::AbstractIdData, na; λ=0, estimator=\, scaleB=false, stochastic=false)

Вычисляет авторегрессионную передаточную функцию G = 1/A; авторегрессионный процесс определяется как A(z⁻¹)y(t) = e(t)

Аргументы

d: IdData, см. описание iddata.
na: порядок модели
λ: λ > 0 можно указать для регуляризации L₂.
estimator: например, \,tls,irls,rtls
scaleB: нужно ли масштабировать числитель с помощью дисперсии ошибки прогнозирования.
stochastic: если задано значение true, возвращает передаточную функцию с неопределенными параметрами, представленными MonteCarloMeasurements.Particles.

Известно, что оценка авторегрессионных моделей методом наименьших квадратов затруднена при сильном шуме измерений, поэтому использование estimator = tls может улучшить результат в этом случае.

Пример

julia> N = 10000
10000

julia> e = [-0.2; zeros(N-1)] # шум e
10000-element Vector{Float64}:
[...]

julia> G = tf([1, 0], [1, -0.9], 1) # Авторегрессионная передаточная функция
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   1.0z
----------
1.0z - 0.9

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> y = lsim(G, e, 1:N)[1][:] # Получение выходного сигнала авторегрессионной передаточной функции из входного шума e
10000-element Vector{Float64}:
[...]

julia> Gest = ar(iddata(y), 1) # Вычисление авторегрессионной передаточной функции из выходного y
TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}}
          1.0z
-------------------------
1.0z - 0.8999999999999998

Sample Time: 1.0 (seconds)
Discrete-time transfer function model

julia> G ≈ Gest # Проверка правильности вычисления
true

julia> eest = lsim(1/Gest, y, 1:N)[1][:] # восстановление входного шума e из выходного y и вычисленной передаточной функции Gest
10000-element Vector{Float64}:
[...]

julia> isapprox(eest, e, atol = eps()) # входной шум корректно восстановлен
true

# ControlSystemIdentification.arma — Function

model = arma(d::AbstractIdData, na, nc; initial_order=20, method=:ls)

Оценивает модель авторегрессионного скользящего среднего с коэффициентами na в знаменателе и коэффициентами nc в числителе. Возвращает модель и оцененную последовательность шумов, управляющих системой.

Аргументы

d: iddata
initial_order: для оценки остатков используется начальная авторегрессионная модель этого порядка.
estimator: Функция (A,y)->minimizeₓ(Ax-y) по умолчанию равна \, но есть и другой вариант — wtls_estimator(1:length(y)-initial_order,na,nc,ones(nc)).

См. также описание estimate_residuals.

# ControlSystemIdentification.plr — Function

G, Gn = plr(d::AbstractIdData,na,nb,nc; initial_order = 20)

Выполняет псевдолинейную регрессию для оценки модели в форме. Ay = Bu + Cw Остаточная последовательность оценивается сначала путем оценки модели arx высокого порядка, после чего оцененная остаточная последовательность включается во вторую задачу оценки. Возвращаемыми значениями являются оцененная модель системы и оцененная модель шума. G и Gn всегда будут иметь один и тот же многочлен в знаменателе.

armax является псевдонимом для plr. См. также описание pem, ar, arx и arxar.

# ControlSystemIdentification.arxar — Function

G, H, e = arxar(d::InputOutputData, na::Int, nb::Union{Int, Vector{Int}}, nd::Int)

Оценивает модель ARXAR Ay = Bu + v, где v = He и H = 1/D, используя обобщенный метод наименьших квадратов. Дополнительные сведения см. в работе Сёдерстёма (Söderström) Convergence properties of the generalized least squares identification method, 1974.

Аргументы

d: iddata
na: порядок A
nb: количество коэффициентов в B, порядок определяется с помощью nb + inputdelay - 1. При оценке MISO он имеет вид nb = [nb₁, nb₂...].
nd: порядок D

Именованные аргументы

H = nothing: предварительные знания об авторегрессионной модели шума.
inputdelay = ones(Int, size(nb)): необязательная задержка входных данных, inputdelay = 0 приводит к получению прямой составляющей, принимает вид inputdelay = [d₁, d₂…] при оценке MISO.
λ = 0: λ > 0 можно указать для регуляризации L₂.
estimator = \: например, \,tls,irls,rtls, для трех последних требуется using TotalLeastSquares.
δmin = 10e-4: минимальное изменение степени e, указывающее на сходимость.
iterations = 10: максимальное количество итераций.
verbose = false: если задано значение true, выводятся дополнительные сведения.

См. также описание plr и arx.

Пример:

julia> N = 500
500

julia> sim(G, u) = lsim(G, u, 1:N)[1][:]
sim (generic function with 1 method)

julia> A = tf([1, -0.8], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
1.0z - 0.8
----------
   1.0z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> B = tf([0, 1], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Int64}}
1
-
z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> G = minreal(B / A)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   1.0
----------
1.0z - 0.8

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> D = tf([1, 0.7], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
1.0z + 0.7
----------
   1.0z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> H = 1 / D
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   1.0z
----------
1.0z + 0.7

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> u, e = randn(1, N), randn(1, N)
[...]

julia> y, v = sim(G, u), sim(H * (1/A), e) # моделирование процесса
[...]

julia> d = iddata(y .+ v, u, 1)
InputOutput data of length 500 with 1 outputs and 1 inputs

julia> na, nb , nd = 1, 1, 1
(1, 1, 1)

julia> Gest, Hest, res = arxar(d, na, nb, nd)
(G = TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   0.9987917259291642
-------------------------
1.0z - 0.7937837464682017

Sample Time: 1 (seconds)
Discrete-time transfer function model, H = TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
          1.0z
-------------------------
1.0z + 0.7019519225937721

Sample Time: 1 (seconds)
Discrete-time transfer function model, e = [...]

# ControlSystemIdentification.getARXregressor — Function

getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb; inputdelay = ones(Int, size(nb)))

Возвращает сокращенный выходной сигнал y и матрицу регрессии A — такую, что оценка модели ARX методом наименьших квадратов порядка na,nb равна y\A.

Возвращает матрицу регрессии, используемую для подгонки модели ARX, например в виде A(z)y = B(z)u с выходными данными y и входными данными u, где порядок авторегрессии равен na, порядок входного скользящего среднего равен nb, и необязательная входная задержка равна inputdelay. Внимание: при изменении входной задержки порядок изменяется на nb + inputdelay - 1. inputdelay = 0 Приводит к получению прямой составляющей.

Пример

A     = [1,2*0.7*1,1] # Коэффициенты A(z)
B     = [10,5] # Коэффициенты B(z)
u     = randn(100) # Моделирование 100 временных шагов с гауссовым входом
y     = filt(B,A,u)
yr,A  = getARXregressor(y,u,3,2) # Предполагается, что известен системный порядок — 3,2
x     = A\yr # Оценка многочленов модели
plot([yr A*x], lab=["Signal" "Prediction"])

Нелинейные модели ARX см. в описании BasisFunctionExpansions.jl. См. также описание arx.

# ControlSystemIdentification.getARregressor — Function

yt,A = getARregressor(y::AbstractVector, na)

Возвращает такие значения, что x = A\yt. Дополнительные сведения см. в описании getARXregressor.

Видеоруководства

Видеоруководства по оценке передаточных функций доступны здесь: