Документация Engee

Экспортируемые функции и типы

Указатель

# ControlSystemIdentification.FRD — Type

FRD(w, r)

Представляет данные частотной характеристики. w содержит вектор частоты, а r — характеристику. К методам, определенным для этого типа, относятся

  • +-*

  • length, vec, sqrt

  • plot

  • feedback

  • freqvec

  • tfest для оценки рациональной модели

  • Индексирование в частотной области с использованием, например G[1Hz : 5Hz], G[1rad : 5rad]

Если r представляет собой частотную характеристику MIMO, ее измерения равны ny × nu × nω.

График объекта frd::FRD можно построить с помощью plot(frd, hz=false), если была вызвана команда using Plots.

# ControlSystemIdentification.FRD — Method

FRD(w, sys::LTISystem)

Создает объект данных частотной характеристики, оценив частотную характеристику sys на частотах w.

# ControlSystemIdentification.Hz — Type

Представляет частоты в герцах для индексирования объектов FRD: frd[2Hz:10Hz]

# ControlSystemIdentification.InputOutputData — Type

См. описание iddata.

# ControlSystemIdentification.InputOutputFreqData — Type

См. описание iddata.

# ControlSystemIdentification.N4SIDStateSpace — Type

N4SIDStateSpace <: AbstractPredictionStateSpace

Результат оценки модели пространства состояний с помощью метода n4sid.

Поля:

  • sys: оцениваемая модель в виде объекта StateSpace

  • Q: оцениваемая ковариационная матрица состояний

  • R: оцениваемая ковариационная матрица измерений

  • S: оцениваемая матрица взаимной ковариации между состояниями и измерениями

  • K: коэффициент усиления наблюдателя Калмана

  • P: решение уравнения Риккати

  • x: оцениваемая траектория состояния (n4sid) или начальное условие (subspaceid)

  • s: сингулярное разложение

  • fve: доля дисперсии, объясненная сингулярными значениями

# ControlSystemIdentification.PredictionStateSpace — Type

PredictionStateSpace{T, ST <: AbstractStateSpace{T}, KT, QT, RT, ST2} <: AbstractPredictionStateSpace{T}
PredictionStateSpace(sys, K, Q=nothing, R=nothing, S=nothing)

Тип пространства состояний, содержащий дополнительную модель фильтра Калмана для целей прогнозирования.

Аргументы

  • sys: ОПИСАНИЕ

  • K: коэффициент усиления Калмана в бесконечном интервале

  • Q = nothing: ковариация динамики

  • R = nothing: ковариация измерений

  • S = nothing: взаимная ковариация

# ControlSystemIdentification.rad — Type

Представляет частоты в рад/с для индексирования объектов FRD: frd[2rad:10rad]

# ControlSystemIdentification.apply_fun — Function

apply_fun(fun, d::InputOutputData)

Применяет fun(y) ко всем временным рядам y[,u,[x]] ∈ d и возвращает новый iddata с преобразованными рядами.

# ControlSystemIdentification.ar — Method

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 — Method

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.arma_ssa — Method

arma_ssa(d::AbstractIdData, na, nc; L=nothing, estimator=\, robust=false)

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

Аргументы

  • d: iddata

  • na: количество параметров знаменателя

  • nc: количество параметров числителя

  • L: длина внедрения запаздывания, используемого для разделения сигнала и шума. nothing соответствует автоматическому выбору.

  • estimator: функция для решения задачи наименьших квадратов. Примеры \,tls,irls,rtls.

  • robust: использование надежного PCA для обеспечения устойчивости к выбросам.

# ControlSystemIdentification.armax — Function

armax является псевдонимом для plr

# ControlSystemIdentification.arx — Method

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.arxar — Method

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.autocorplot — Function

autocorplot(y, Ts, [lags])

Строит график автокорреляции y для lags, которая по умолчанию равна 1:size(y, 2)÷2.

# ControlSystemIdentification.coherence — Method

κ² = coherence(d; n = length(d)÷10, noverlap = n÷2, window=hamming, method=:welch)

Вычисляет и строит график функции когерентности с квадратом величины. κ² с величиной, близкой к 1, указывает на хорошую объяснимость энергии в выходном сигнале энергией входного сигнала. κ² << 1 указывает на то, что либо система нелинейна, либо сильный шум занимает значительную часть в выходной энергии.

  • κ: Квадратичная функция когерентности в виде FRD.

  • method: :welch или :corr. :welch использует метод Велча для оценки спектральной плотности мощности, а :corr использует метод коррелограммы. Для method = :corr дополнительный именованный аргумент σ определяет ширину гауссова окна, применяемого к оцененным корреляционным функциям перед БПФ. Большее значение σ подразумевает меньшее сглаживание.

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

Расширенная справка

Для модели сигнала y = Gu + v κ² определяется следующим образом

откуда следует, что 0 ≤ κ² ≤ 1 и что κ² близка к 1, если энергия шума S*{vv} мала по сравнению с выходной энергией, обусловленной входными данными S*{uu}\|G(iω)\|\^2.

# ControlSystemIdentification.coherenceplot — Function

coherenceplot(d, [(;n=..., noverlap=...); hz=false)

Вычисляет и строит график функции (квадрата) когерентности κ. κ с величиной, близкой к 1, указывает на хорошую объяснимость энергии в выходном сигнале энергией входного сигнала. κ << 1 указывает на то, что либо система нелинейна, либо сильный шум занимает значительную часть в выходной энергии.

hz указывает на герцы вместо рад/с.

Именованные аргументы для coherence предоставляются в виде именованного кортежа в качестве второго позиционного аргумента.

# ControlSystemIdentification.crosscorplot — Function

crosscorplot(data, [lags])
crosscorplot(u, y, Ts, [lags])

Строит график взаимной корреляции между входом и выходом для lags, который по умолчанию составляет 10 % от длины набора данных с отрицательной стороны и 50 % с положительной стороны, но не более 100 с каждой стороны.

# ControlSystemIdentification.detrend — Method

detrend(d::AbstractArray)
detrend(d::AbstractIdData)

Удаляет среднее значение из d.

# ControlSystemIdentification.era — Function

era(YY::AbstractArray{<:Any, 3}, Ts, nx::Int, m::Int = 2nx, n::Int = 2nx)

Алгоритм реализации собственных значений. Алгоритм возвращает модель пространства состояний.

Аргументы

  • YY: Размер марковских параметров (импульсного отклика) ny × nu × n_time

  • Ts: Интервал выборки

  • nx: порядок модели

  • m: количество строк в матрице Ганкеля

  • n: количество столбцов в матрице Ганкеля

# ControlSystemIdentification.era — Function

era(d::AbstractIdData, nx; m = 2nx, n = 2nx, l = 5nx, p = l, λ=0, smooth=false)
era(ds::Vector{IdData}, nx; m = 2nx, n = 2nx, l = 5nx, p = l, λ=0, smooth=false)

Алгоритм реализации собственных значений. Использует okid для нахождения марковских параметров в качестве начального шага.

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

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

Аргументы

  • nx: порядок модели

  • l: количество оцениваемых марковских параметров.

  • λ: параметр регуляризации (им не стоит злоупотреблять, лучше провести больше экспериментов)

  • smooth: Если задано значение true, регуляризация, заданная с помощью λ, штрафует за кривизну в оцененном импульсном отклике.

  • p: По желанию удалите первые p столбцов во внутренних матрицах Ганкеля, чтобы учесть начальные условия != 0. Если x0 != 0, то для era p по умолчанию будет l, а при прямом вызове okid p по умолчанию будет 0.

# ControlSystemIdentification.estimate_residuals — Method

estimate_residuals(model, y)

Оценивает остатки, определяющие динамику модели ARMA.

# ControlSystemIdentification.estimate_x0 — Function

estimate_x0(sys, d, n = min(length(d), 3 * slowest_time_constant(sys)); fixed = fill(NaN, sys.nx)

Оценивает начальное состояние системы.

Аргументы

  • d: iddata

  • n: количество используемых выборок.

  • fixed: Если задан вектор той же длины, что и x0, конечные значения указывают на фиксированные значения, которые не подлежат оценке, а неограниченные значения являются свободными.

Пример

sys   = ssrand(2,3,4, Ts=1)
x0    = randn(sys.nx)
u     = randn(sys.nu, 100)
y,t,x = lsim(sys, u; x0)
d     = iddata(y, u, 1)
x0h   = estimate_x0(sys, d, 8, fixed=[Inf, x0[2], Inf, Inf])
x0h[2] == x0[2] # Должно быть точное равенство
norm(x0-x0h)    # Должен быть небольшим

# ControlSystemIdentification.filter_bank — Method

filter_bank(basis::AbstractStateSpace{<:Discrete}, signal::AbstractMatrix)

Фильтрует signal во всех системах в basis.

# ControlSystemIdentification.find_na — Function

find_na(y::AbstractVector,n::Int)

Строит графики RMSE и AIC. Для порядков моделей, не превышающих n. Полезна для выбора моделей.

# ControlSystemIdentification.find_nanb — Function

find_nanb(d::InputOutputData, na, nb; logrms = false, method = :aic)

Строит графики RMSE и AIC. Для порядков моделей, не превышающих na, nb. Полезна для выбора моделей. na может быть целым числом или диапазоном. То же самое относится и к nb.

  • logrms определяет, строить ли график логарифма среднеквадратичной ошибки по основанию 10.

  • method: определяет, использовать ли для определения порядка модели информационный критерий Акаике (:aic) или итоговую ошибку прогнозирования (:fpe).

Если цветовая шкала трудночитаема из-за нескольких плиток, представляющих очень большие ошибки, старайтесь не рисовать эти плитки, указывая диапазоны na и nb вместо целых чисел. Например, не показывайте порядок модели меньше 2, используя find_nanb(d, 3:na, 3:nb)

# ControlSystemIdentification.find_similarity_transform — Function

find_similarity_transform(sys1, sys2, method = :obsv)

Находит T так, что ControlSystemsBase.similarity_transform(sys1, T) == sys2

Источник: Minimal state-space realization in linear system theory: an overview, B. De Schutter

Если method == :obsv, для нахождения T используются матрицы наблюдаемости sys1 и sys2, а method == :ctrb использует матрицы управляемости.

julia> T = randn(3,3);

julia> sys1 = ssrand(1,1,3);

julia> sys2 = ControlSystemsBase.similarity_transform(sys1, T);

julia> T2 = find_similarity_transform(sys1, sys2);

julia> T2 ≈ T
true

# ControlSystemIdentification.freqvec — Method

freqvec(h, k)

Возвращает вектор частот длиной k для систем с интервалом выборки h.

# ControlSystemIdentification.getARXregressor — Method

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 — Method

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

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

# ControlSystemIdentification.iddata — Function

iddata(y,       Ts = nothing)
iddata(y, u,    Ts = nothing)
iddata(y, u, x, Ts = nothing)

Создает объект идентификационных данных временной области.

Аргументы

  • y::AbstractArray: выходные данные (обязательно)

  • u::AbstractArray: входные данные (если доступны)

  • x::AbstractArray: данные состояния (если доступны)

  • Ts::Union{Real,Nothing} = nothing: дополнительный интервал выборки

Если временные ряды являются многомерными, время находится в последнем измерении, то есть размеры массивов равны (num_variables, num_timepoints) (см. примеры ниже).

Операции с iddata

  • detrend

  • prefilter

  • resample

  • Добавление двух к временному измерению [d1 d2] (это делается только в том случае, если состояние системы в конце d1 близко к состоянию в начале d2)

  • Индексирование временных рядов d[output_index, input_index]

  • Индексирование оси времени с помощью индексов d[time_indices]

  • Индексирование оси времени с помощью секунд d[3Sec:12Sec] (using ControlSystemIdentification: Sec)

  • Доступ к количеству входных данных, выходных данных и времени выборки: d.nu, d.ny, d.Ts

  • Доступ к вектору времени d.t

  • Умножение в обычном порядке для масштабирования выходных данных C * d. Перед оценкой модели обычно рекомендуется масштабировать выходные данные системы с несколькими выводами так, чтобы они имели примерно одинаковый размер, в случае если они имеют разную величину.

  • Умножение в обычном порядке для масштабирования входных данных d * B

  • writedlm

  • ramp_in, ramp_out

  • plot

  • specplot

  • crosscorplot

Примеры

julia> iddata(randn(10))
Output data of length 10 with 1 outputs, Ts = nothing

julia> iddata(randn(10), randn(10), 1)
InputOutput data of length 10, 1 outputs, 1 inputs, Ts = 1

julia> d = iddata(randn(2, 10), randn(3, 10), 0.1)
InputOutput data of length 10, 2 outputs, 3 inputs, Ts = 0.1

julia> [d d] # Конкатенация по времени
InputOutput data of length 20, 2 outputs, 3 inputs, Ts = 0.1

julia> d[1:3]
InputOutput data of length 3, 2 outputs, 3 inputs, Ts = 0.1

julia> d.nu
3

julia> d.t # доступ к вектору времени
0.0:0.1:0.9

Использование нескольких наборов данных

Некоторые методы оценки поддерживают использование нескольких наборов данных для оценки модели. В этом случае наборы данных предоставляются в виде вектора объектов iddata. Вот эти методы:

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

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

# ControlSystemIdentification.iddata — Method

iddata(y::AbstractArray, u::AbstractArray, w::AbstractVector)

Создает объект входных-выходных данных частотной области. w должен быть указан в рад/с.

# ControlSystemIdentification.iddata — Method

iddata(res::ControlSystemsBase.SimResult)

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

# ControlSystemIdentification.impulseest — Method

ir, t, Σ = impulseest(d::AbstractIdData, n; λ=0, estimator=ls)

Оценивает импульсный отклик системы путем подгонки FIR-модели n-го порядка. Возвращает оценку импульсного отклика, вектор времени и ковариационную матрицу.

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

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

# ControlSystemIdentification.impulseestplot — Function

impulseestplot(data,n; σ = 2)

Оценивает импульсный отклик системы путем подгонки FIR-модели n-го порядка и строит график результата с доверительной полосой 95 % (2σ). Обратите внимание, что доверительная полоса строится вокруг нуля, то есть она строится таким образом, чтобы можно было определить, существенно ли отличается импульсный отклик от нуля.

Этот метод поддерживает только данные с одним выходным значением. Для данных с несколькими выходными значениями используйте okid.

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

# ControlSystemIdentification.kautz — Method

kautz(a::Vector, h)

Создает базис Каутца с дискретным временем с полюсами в точке a и временем выборки h.

# ControlSystemIdentification.laguerre — Method

laguerre(a::Number, Nq)

Создает базис Лагерра длиной Nq с полюсами в точке -a.

# ControlSystemIdentification.laguerre_id — Method

laguerre_id(a::Number, Nq, Ts)

Создает базис Лагерра с дискретным временем длиной Nq с полюсами в точке -a для идентификации системы.

ПРИМЕЧАНИЕ. При больших значениях Nq этот базис может быть численно плохо обусловлен. Рекомендуется применить balance_statespace к полученному базису.

# ControlSystemIdentification.laguerre_oo — Method

laguerre_oo(a::Number, Nq)

Создает выходной ортогонализированный базис Лагерра длиной Nq с полюсами в точке -a.

# ControlSystemIdentification.minimum_phase — Method

minimum_phase(G)

Перемещает нули и полюса G из неустойчивой полуплоскости в устойчивую. Если G является системой с пространством состояний, сначала она преобразуется в передаточную функцию. Это может привести к потере точности.

# ControlSystemIdentification.model_spectrum — Method

model_spectrum(f, h, args...; kwargs...)

Аргументы

  • f: функция оценки модели, например ar,arma

  • h: интервал выборки

  • args: аргументы для f

  • kwargs: именованные аргументы для f

Пример:

using ControlSystemIdentification, DSP
T = 1000
s = sin.((1:T) .* 2pi/10)
S1 = spectrogram(s,window=hanning)
estimator = model_spectrum(ar,1,2)
S2 = spectrogram(s,estimator,window=rect)
plot(plot(S1),plot(S2)) # Требуется пакет LPVSpectral.jl

# ControlSystemIdentification.modelfit — Method

modelfit(y, yh)

Вычисляет соответствие модели yh для y в процентах, иногда называемое нормализованной среднеквадратичной ошибкой (NRMSE).

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

См. также описание rms, sse, mse, fpe, aic.

# ControlSystemIdentification.n4sid — Method

res = n4sid(data, r=:auto; verbose=false)

Оценивает модель пространства состояний с помощью метода n4sid. Возвращает объект типа N4SIDStateSpace, обеспечивающий доступ к модели в виде res.sys.

Реализует упрощенный алгоритм (alg 2), приведенный в работе N4SID: Subspace Algorithms for the Identification of Combined Deterministic Stochastic Systems, авторы: PETER VAN OVERSCHEE и BART DE MOOR

Идеи для частотного взвешивания взяты из работы Frequency Weighted Subspace Based System Identification in the Frequency Domain, Tomas McKelvey (1996). В частности, мы применяем весовую матрицу выходных частот (Fy) так, как она представлена в уравнениях. (16)--(18).

Аргументы

  • data: Идентификационные данные data = iddata(y,u)

  • r: ранг модели (порядок модели)

  • verbose: выводить что-либо?

  • Wf: модель возмущений измерений в частотной области. Чтобы сфокусировать внимание модели на узкой полосе частот, попробуйте использовать что-то вроде Wf = Bandstop(lower, upper, fs=1/Ts), чтобы указать, что существуют возмущения за пределами этой полосы.

  • i: параметр алгоритма, обычно настраивать его не нужно.

  • γ: задайте для него значение в диапазоне (0,1) для стабилизации неустойчивых моделей таким образом, чтобы наибольшее собственное значение имело величину γ.

  • zeroD: значение по умолчанию — false

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

# ControlSystemIdentification.newpem — Method

sys, x0, res = newpem(
    d,
    nx;
    zeroD  = true,
    focus  = :prediction,
    stable = true,
    sys0   = subspaceid(d, nx; zeroD, focus, stable),
    metric = abs2,
    regularizer = (p, P) -> 0,
    output_nonlinearity = nothing,
    input_nonlinearity = nothing,
    nlp = nothing,
    optimizer = BFGS(
        linesearch = LineSearches.BackTracking(),
    ),
    store_trace = true,
    show_trace  = true,
    show_every  = 50,
    iterations  = 10000,
    time_limit  = 100,
    x_tol       = 0,
    f_abstol    = 0,
    g_tol       = 1e-12,
    f_calls_limit = 0,
    g_calls_limit = 0,
    allow_f_increases = false,
)

Новая реализация метода на основе ошибки прогнозирования (PEM). Обратите внимание, что это экспериментальная реализация и в нее могут быть внесены изменения, не соответствующие semver.

Метод на основе ошибки прогнозирования является итеративной задачей оптимизации на основе градиента, поэтому он может быть очень чувствителен к масштабированию сигнала. Перед оценкой рекомендуется выполнить масштабирование к d, например путем умножения слева и справа на диагональные матрицы d̃ = Dy*d*Du, и применить обратное масштабирование к полученной системе В этом случае у нас есть

отсюда G = Dy \ G̃ * Du, где  — это объект, оцениваемый для масштабированных iddata. Пример:

Dy = Diagonal(1 ./ vec(std(d.y, dims=2))) # Нормализация дисперсии
Du = Diagonal(1 ./ vec(std(d.u, dims=2))) # Нормализация дисперсии
d̃ = Dy * d * Du

Если начальное предположение sys0 задано вручную, оно также должно быть соответствующим образом масштабировано.

Аргументы

  • d: iddata

  • nx: порядок модели

  • zeroD: Делает матрицу D нулевой

  • stable: если true, устойчивость оцениваемой системы будет обеспечиваться отражением собственных значений с помощью schur_stab при ϵ=1/100 (по умолчанию). Если stable является вещественным значением, оно используется вместо заданного по умолчанию ϵ.

  • sys0: Начальное предположение. Если оно не указано, то в качестве начального предположения используется subspaceid.

  • focus: prediction или :simulation. Если :simulation, матрица K будет нулевой.

  • optimizer: один из оптимизаторов Optim

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

Остальные аргументы относятся к Optim.Options.

  • regularizer: Функция вектора параметров и соответствующей системы PredictionStateSpace/StateSpace, которая может быть использована для регуляризации оценки.

  • output_nonlinearity: Функция (y::Vector, p), которая работает с выходным сигналом в одной временной точке yₜ и изменяет его на месте. Подробные сведения см. ниже. p является вектором оцениваемых параметров, которые можно оптимизировать.

  • input_nonlinearity: Функция (u::Matrix, p), которая работает со всем входным сигналом u за раз и изменяет его на месте. Подробные сведения см. ниже. p представляет собой вектор оцениваемых параметров, который является общим и для output_nonlinearity.

  • nlp: вектор начального предположения для нелинейных параметров. Если указан output_nonlinearity, он может быть указан дополнительно.

Нелинейная оценка

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

  1. Если существует известная входная нелинейность, вручную примените входную нелинейность к входному сигналу u перед оценкой, т. е. используйте нелинейно преобразованный входной сигнал в объекте iddata d. Если входная нелинейность имеет неизвестные параметры, укажите входную нелинейность в виде функции с именованным аргументом input_nonlinearity для newpem. Предполагается, что эта функция будет работать со всем (матричным) выходным сигналом u и изменять его на месте.

  2. Если выходная нелинейность является инвертируемой, примените инверсию к выходному сигналу y перед оценкой аналогично описанному выше.

  3. Если выходная нелинейность является неинвертируемой, укажите выходное нелинейное преобразование в виде функции с именованным аргументом output_nonlinearity для newpem. Предполагается, что эта функция будет работать с (векторным) выходным сигналом y и изменять его на месте. Пример:

function output_nonlinearity(y, p)
    y[1] = y[1] + p[1]*y[1]^2       # Обратите внимание, как входящий вектор изменяется на месте.
    y[2] = abs(y[2])
end

Обратите внимание, y = f(y) не изменяет y на месте, а создает новый вектор y и присваивает его переменной y. Это не то, что нам нужно.

Второй аргумент для input_nonlinearity и output_nonlinearity — это (необязательный) вектор параметров, которые могут быть оптимизированы. Чтобы использовать эту возможность, передайте именованный аргумент nlp в newpem с вектором начальных предположений для нелинейных параметров. Нелинейные параметры являются общими для выходных и входных нелинейностей, т. е. эти две функции получат один и тот же вектор параметров.

Результатом данной оценки является линейная система без нелинейностей.

Пример

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

using ControlSystemIdentification, ControlSystemsBase Plots
G = DemoSystems.doylesat()
T = 1000  # Количество временных шагов
Ts = 0.01 # Интервал выборки
sys = c2d(G, Ts)
nx = sys.nx
nu = sys.nu
ny = sys.ny
x0 = zeros(nx) # фактическое начальное состояние
sim(sys, u, x0 = x0) = lsim(sys, u; x0)[1]

σy = 1e-1 # Ковариация шума

u  = randn(nu, T)
y  = sim(sys, u, x0)
yn = y .+ σy .* randn.() # Добавляем шум измерения
d  = iddata(yn, u, Ts)

sysh, x0h, opt = ControlSystemIdentification.newpem(d, nx, show_every=10)

plot(
    bodeplot([sys, sysh]),
    predplot(sysh, d, x0h), # Включает оцененное начальное состояние в прогнозирование
)

Возвращаемая модель имеет тип PredictionStateSpace и содержит поле sys с моделью системы, а также ковариационные матрицы и вычисленный коэффициент усиления Калмана для фильтра Калмана.

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

Расширенная справка

В данной реализации используется трехдиагональная параметризация матрицы A, которая, как было показано, является предпочтительной с точки зрения оптимизации¹. Начальное предположение sys0 автоматически преобразуется в специальную трехдиагональную модальную форму. [1] Mckelvey, Tomas & Helmersson, Anders. (1997). Параметризация пространств состояний многомерных линейных систем с помощью трехдиагональных матричных форм.

Вектор параметров, используемый при оптимизации, имеет следующий вид:

p = [trivec(A); vec(B); vec(C); vec(D); vec(K); vec(x0)]

Где ControlSystemIdentification.trivec векторизует диагонали -1,0,1 для A. Если focus = :simulation, K опускается, а если zeroD = true, опускается D.

# ControlSystemIdentification.noise_model — Method

noise_model(sys::AbstractPredictionStateSpace)

Возвращает модель шума, управляющего системой v, в

Модель игнорирует u и задается следующим образом:

Также называется «инновационная форма». Эта функция вызывает ControlSystemsBase.innovation_form.

# ControlSystemIdentification.okid — Function

H = okid(d::AbstractIdData, nx, l = 5nx; p = 1, λ=0, estimator = /)

Идентификация с помощью фильтра наблюдателя Калмана. Возвращает размер H ny × nu × (l+1) марковских параметров (импульсного отклика).

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

Аргументы

  • nx: порядок модели

  • l: количество оцениваемых марковских параметров (длина импульсного отклика).

  • λ: параметр регуляризации

  • smooth: Если задано значение true, регуляризация, заданная с помощью λ, штрафует за кривизну в оцененном импульсном отклике. Если era будет использоваться после okid, следует отдать предпочтение небольшому λ с smooth=true, но, если импульсный отклик будет проверяться на глаз, большее сглаживание может дать визуально более точную оценку импульсного отклика.

  • p: По желанию удалите первые p столбцов во внутренних матрицах Ганкеля, чтобы учесть начальные условия != 0. Если x0 != 0, попробуйте задать для p примерно то же значение, что и для l.

  • estimator: функция, используемая для оценки марковских параметров. По умолчанию используется / (наименьшие квадраты), но может быть и надежным вариантом, таким как TotalLeastSquares.irls / flts или TotalLeastSquares.tls, для полного решения по методу наименьших квадратов (ошибки в переменных)

# ControlSystemIdentification.pem — Method

Эта функция является нерекомендуемой, см. описание newpem.

# ControlSystemIdentification.plr — Method

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

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

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

# ControlSystemIdentification.prediction_error — Method

e = prediction_error(sys::AbstractStateSpace, d::AbstractIdData, args...; kwargs...)

Возвращает ошибки прогнозирования `d.y - predict(sys, d, …​)

# ControlSystemIdentification.prediction_error_filter — Method

prediction_error_filter(sys::AbstractPredictionStateSpace; h=1)
prediction_error_filter(sys::AbstractStateSpace, R1, R2; h=1)

Возвращает фильтр, который принимает [u; y] в качестве входных данных и выдает ошибку прогнозирования e = y - ŷ. См. также описание innovation_form и noise_model. h ≥ 1 является горизонтом прогнозирования. См. описание функции predictiondata для создания iddata с [u; y] в качестве входных данных.

# ControlSystemIdentification.predictiondata — Method

predictiondata(d::AbstractIdData)

Добавляет выходное значение y к входному u_new = [u; y]

# ControlSystemIdentification.predplot — Function

predplot(sys, data, x0=:estimate; ploty=true, plote=false, h=1, sysname="")

Строит графики h-шагового прогнозирования системы и измеренных выходных данных для их сравнения.

По умолчанию начальное условие x0 оценивается с использованием данных. Чтобы начать моделирование из источника, укажите x0 = :zero или x0 = zeros(sys.nx).

  • ploty определяет, строить ли график измеренного сигнала.

  • plote определяет, строить ли график остатка.

  • h является горизонтом прогнозирования.

# ControlSystemIdentification.prefilter — Method

prefilter(f, d::InputOutputData)

Применяет коэффициенты фильтрации к идентификационным данным

# ControlSystemIdentification.prefilter — Method

prefilter(d::AbstractIdData, responsetype::FilterType)

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

# ControlSystemIdentification.prefilter — Method

prefilter(d::AbstractIdData, l::Number, u::Number)

Фильтрует входные и выходные данные с помощью полосового пропускающего фильтра в диапазоне от l до u Гц. Если l = 0, будет использоваться низкочастотный фильтр, а если u = Inf, — высокочастотный.

# ControlSystemIdentification.ramp_in — Method

ramp_in(d::InputOutputData, h::Int; rev = false)

Умножает начальные h выборок входного и выходного сигналов с линейно возрастающим темпом.

# ControlSystemIdentification.ramp_out — Method

ramp_out(d::InputOutputData, h::Int)

Умножает конечные h выборок входного и выходного сигналов с линейно убывающим темпом.

# ControlSystemIdentification.residualplot — Function

residualplot(model, data)

Строит графики остаточной автокорреляции и корреляции между входным сигналом и остатком.

# ControlSystemIdentification.schur_stab — Method

schur_stab(A::AbstractMatrix{T}, ϵ = 0.01)

Стабилизирует собственные значения матрицы A с дискретным временем путем преобразования матрицы A в комплексную. Форма Шура и проецирование неустойчивых собственных значений 1-ϵ < λ ≤ 2 на единичный диск. Для собственных значений больше двух задается нуль.

# ControlSystemIdentification.simplot — Function

simplot(sys, data, x0=:estimate; ploty=true, plote=false, sysname="")

Строит графики моделирования системы и измеренных выходных данных для их сравнения.

По умолчанию начальное условие x0 оценивается с использованием данных. Чтобы начать моделирование из источника, укажите x0 = :zero или x0 = zeros(sys.nx).

  • ploty определяет, строить ли график измеренного сигнала.

  • plote определяет, строить ли график остатка.

# ControlSystemIdentification.specplot — Function

specplot(d::IdData, args...; kwargs...)

Строит спектрограмму входных и выходных временных рядов. См. также описание welchplot.

Дополнительные аргументы передаются в DSP.spectrogram.

# ControlSystemIdentification.subspaceid — Method

model, x0 = subspaceid(frd::FRD, Ts, args...; estimate_x0 = false, bilinear_transform = false, kwargs...)

Если предоставлен объект данных частотной характеристики

  • FRD автоматически преобразуется в InputOutputFreqData.

  • Для estimate_x0 по умолчанию задан 0.

  • bilinear_transform преобразует вектор частот в дискретное время, см. примечание ниже.

Примечание. Если данные частотных характеристик получены в результате анализа частотных характеристик, перед оценкой необходимо провести билинейное преобразование данных. Это преобразование будет применено, если bilinear_transform = true.

# ControlSystemIdentification.subspaceid — Method

subspaceid(
    data::InputOutputData,
    nx = :auto;
    verbose = false,
    r = nx === :auto ? min(length(data) ÷ 20, 50) : nx + 10, # Используемый максимальный горизонт прогнозирования
    s1 = r, # Количество прошлых выходных значений
    s2 = r, # Количество прошлых входных значений
    W = :MOESP,
    zeroD = false,
    stable = true,
    focus = :prediction,
    svd::F1 = svd!,
    scaleU = true,
    Aestimator::F2 = \,
    Bestimator::F3 = \,
    weights = nothing,
)

Оценивает модель пространства состояний с помощью идентификации на основе подпространства. Доступно несколько различных алгоритмов, основанных на подпространстве, которые можно выбрать с помощью ключевого слова W. Возможные варианты: :MOESP, :CVA, :N4SID, :IVM

Источник: Ljung, Theory for the user.

Устойчивость к выбросам можно повысить, предоставив собственный алгоритм разложения и заменив внутренние средства оценки по методу наименьших квадратов. См. документацию по именованным аргументам svd, Aestimator и Bestimator ниже.

Возвращаемая модель имеет тип N4SIDStateSpace и содержит поле sys с моделью системы, а также ковариационные матрицы для фильтра Калмана.

Аргументы

  • data: Идентификационные данные iddata

  • nx: ранг модели (порядок модели)

  • verbose: выводить что-либо?

  • r: горизонт прогнозирования. Модель может более эффективно работать в моделировании, если сделать его длиннее, но это потребует больше времени на вычисления.

  • s1: прошлый горизонт выходных данных

  • s2: прошлый горизонт входных данных

  • W: Тип веса, можно выбрать среди :MOESP, :CVA, :N4SID, :IVM

  • zeroD: Позволяет сделать матрицу D нулевой.

  • stable: стабилизирует неустойчивую систему с помощью отражения собственных значений.

  • focus: :prediction или simulation

  • svd: Функция, используемая для svd. Устойчивость к выбросам можно повысить, используя TotalLeastSquares.rpca для предварительной обработки матрицы данных перед применением svd, как svd = A->svd!(rpca(A)[1]).

  • scaleU: изменяет масштаб входных каналов так, чтобы они имели одинаковую энергию.

  • Aestimator: Функция средства оценки, используемая для оценки A,C. По умолчанию используется `, т. е. наименьшие квадраты, но для повышения устойчивости к выбросам можно использовать робастные средства оценки, такие как `irls, flts, rtls из TotalLeastSquares.jl.

  • Bestimator: Функция средства оценки, используемая для оценки B,D. Взвешенную оценку можно получить, передав wls из TotalLeastSquares.jl вместе с именованным аргументом weights.

  • weights: Вектор весов может быть предоставлен, если Bestimator имеет значение wls.

Расширенная справка

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

# ControlSystemIdentification.subspaceid — Method

model, x0 = subspaceid(data::InputOutputFreqData,
    Ts = data.Ts,
    nx = :auto;
    cont = false,
    verbose = false,
    r = nx === :auto ? min(length(data) ÷ 20, 20) : 2nx, # Внутренний порядок модели
    zeroD = false,
    estimate_x0 = true,
    stable = true,
    svd = svd!,
    Aestimator = \,
    Bestimator = \,
    weights = nothing
)

Оценивает модель пространства состояний с помощью идентификации на основе подпространства в частотной области.

Если результаты неудовлетворительны, попробуйте изменить r, особенно если объем данных невелик.

Пример см. в документации.

Аргументы

  • data: объект идентификационных данных частотной области.

  • Ts: интервал выборки, в который были собраны данные.

  • nx: Желаемый порядок модели, целое число или :auto.

  • cont: возвращать модель с непрерывным временем? Для преобразования оцененной модели с дискретным временем используется билинейное преобразование, см. функцию d2c.

  • verbose: выводить что-либо?

  • r: Внутренний порядок модели, значение должно быть больше или равно nx.

  • zeroD: Позволяет сделать матрицу D нулевой.

  • estimate_x0: оценка дополнительных параметров для учета начальных условий. Это может потребоваться, если данные получены в результате БПФ данных временной области, но может и не потребоваться, если данные собраны с помощью анализа частотных характеристик с точно периодическим входом и надлежащей обработкой переходных процессов.

  • stable: Для обеспечения стабильности модели (используется schur_stab).

  • svd: Используемая функция svd.

  • Aestimator: Средство оценки матрицы A (и исходной C-матрицы).

  • Bestimator: средство оценки матриц B/D и C/D.

  • weights: Необязательный вектор частотных весов, длина которого совпадает с количеством частот в `data.

# ControlSystemIdentification.tfest — Function

tfest(
    data::FRD,
    p0,
    link = log ∘ abs;
    freq_weight = sqrt(data.w[1]*data.w[end]),
    refine = true,
    opt = BFGS(),
    opts = Optim.Options(
        store_trace       = true,
        show_trace        = true,
        show_every        = 1,
        iterations        = 100,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 0,
        f_tol             = 0,
        g_tol             = 1e-8,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
    ),
)

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

На начальном этапе оптимизации решается следующее:

а на втором этапе (если refine=true) решается следующее:

(abs2(link(B/A) - link(l)))

Аргументы

  • data: Объект FRD с данными частотной области.

  • p0: начальное предположение параметров. Может быть NamedTuple или ComponentVector с полями b,a, задающими числитель и знаменатель в том виде, в каком они появляются при вызове tf, т. е. (b = [1.0], a = [1.0,1.0,1.0]). Также может быть экземпляром TransferFunction.

  • link: по умолчанию информация об этапе отбрасывается при подгонке. Чтобы включить этап, измените на link = log.

  • freq_weight: применяет взвешивание с обратной частотой. Это значение определяет частоту среза, до которой вес постоянен и после которой вес линейно уменьшается. По умолчанию используется среднее геометрическое значение наименьшей и наибольшей частоты.

  • refine: указывает, выполняется ли второй этап оптимизации для уточнения результатов первого.

  • opt: используемый оптимизатор Optim.

  • opts: Optim.Options — управление параметрами решателя.

См. также описание minimum_phase для преобразования возможно неминимально фазовой системы в минимально фазовую.

# ControlSystemIdentification.tfest — Function

H, N = tfest(data, σ = 0.05, method = :corr)

Оценивает модель передаточной функции с помощью коррелограммы (задано умолчанию), используя модель сигнала y = H(iω)u + n.

И H, и N имеют тип FRD (данные частотной характеристики).

  • σ определяет ширину гауссова окна, применяемого к оцененным корреляционным функциям перед БПФ. Большее значение σ подразумевает меньшее сглаживание.

  • H = Syu/Suu Передаточная функция процесса

  • N = Sy - |Syu|²/Suu Оцениваемая PSD шума (также оценка дисперсии H). Обратите внимание, PSD связана с моделью шума N*m, используемой в литературе по идентификации систем, как N*{psd}} = N*m\^* N*m. Кривая величины модели шума может быть визуализирована путем построения графика √(N).

  • method: :welch или :corr. :welch использует метод Велча для оценки спектральной плотности мощности, а :corr (по умолчанию) — метод коррелограммы. Если method = :welch, дополнительные именованные аргументы n, noverlap и window определяют количество образцов на сегмент (по умолчанию 10 % данных), количество образцов, перекрываемых между сегментами (по умолчанию 50 %), и используемую функцию окна (по умолчанию hamming) соответственно.

Расширенная справка

Этот метод оценки является несмещенным, если входной сигнал u нескоррелирован с шумом n, но в противном случае является смещенным (например, при идентификации в замкнутом контуре).

# ControlSystemIdentification.tfest — Method

tfest(data::FRD, basis::AbstractStateSpace;
    freq_weight = 1 ./ (data.w .+ data.w[2]),
    opt = BFGS(),
    metric::M = abs2,
    opts = Optim.Options(
        store_trace       = true,
        show_trace        = true,
        show_every        = 50,
        iterations        = 1000000,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 1e-5,
        f_tol             = 0,
        g_tol             = 1e-8,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
)

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

Аргументы

  • data: Объект FRD с данными частотной области.

функция kautz(a::AbstractVector)

  • basis: базис для оценки. См., например, описание laguerre, laguerre_oo, kautz.

  • freq_weight: вектор весов на частоту. По умолчанию это значение равно примерно 1/f.

  • opt: используемый оптимизатор Optim.

  • opts: Optim.Options — управление параметрами решателя.

# ControlSystemIdentification.welchplot — Function

welchplot(d::IdData, args...; kwargs...)

Строит периодограмму Велча входных и выходных временных рядов. См. также описание specplot.

Дополнительные аргументы передаются в DSP.welch_pgram.

# ControlSystemIdentification.wtls_estimator — Function

wtls_estimator(y,na,nb, σu=0)

Создает функцию средства оценки для оценки моделей arx при наличии шума измерений. Если известна дисперсия шума на входе σu (ошибки модели), ее можно задать для повышения точности.

# DSP.Filters.resample — Method

DSP.resample(sys::AbstractStateSpace{<:Discrete}, Qd::AbstractMatrix, newh::Real)

Изменяет интервал выборки ковариационной матрицы Qd, принадлежащей sys, на newh. Эта функция не обрабатывает ковариацию измерений. Способ выполнения этой задачи зависит от контекста. Если сигнал из более быстрой выборки имеет такой же уровень шума при измерении, ничего менять не нужно. Если сигнал из медленной выборки был подвергнут пониженной дискретизации с помощью фильтрации, ковариация измерений должна увеличиться, если система перейдет на более высокую частоту дискретизации. Для поддержки частотной характеристики системы необходимо соответствующим образом изменить ковариацию измерений.

Аргументы

  • sys: Система с дискретным временем, имеющая ковариационную матрицу шума динамики Qd.

  • Qd: ковариационная матрица шума динамики.

  • newh: новый интервал выборки.

# DSP.Filters.resample — Method

resample(sys::AbstractStateSpace{<:Discrete}, newh::Real)

Изменяет интервал выборки sys на newh.

# DSP.Filters.resample — Method

dr = resample(d::InputOutputData, f)

Выполняет повторную выборку iddata d с долей f, например f = fs_new / fs_original.

# LowLevelParticleFilters.simulate — Function

simulate(sys, u, x0 = nothing)
simulate(sys, d, x0 = nothing)

См. также описание simplot, predict.

# StatsAPI.predict — Method

predict(sys, d::AbstractIdData, args...)
predict(sys, y, u, x0 = nothing)

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

# StatsAPI.predict — Method

yh = predict(ar::TransferFunction, y)

Прогнозирование авторегрессионной модели

# StatsAPI.predict — Method

predict(ARX::TransferFunction, d::InputOutputData)

Прогнозирование на шаг вперед для процесса ARX. Длина возвращаемого прогнозирования равна length(d) - max(na, nb).

Пример:

julia> predict(tf(1, [1, -1], 1), iddata(1:10, 1:10))
9-element Vector{Int64}:
  2
  4
  6
  8
 10
 12
 14
 16
 18

# StatsAPI.residuals — Method

residuals(ARX::TransferFunction, d::InputOutputData)

Вычисляет остатки v = Ay - Bu процесса ARX и InputOutputData d. Длина возвращаемых остатков равна length(d) - max(na, nb).

Пример:

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

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

julia> u = 1:5
1:5

julia> y = lsim(ARX, u, 1:5)[1][:]
5-element Vector{Float64}:
  0.0
  1.0
  3.0
  6.0
 10.0

julia> d = iddata(y, u)
InputOutput data of length 5 with 1 outputs and 1 inputs

julia> residuals(ARX, d)
4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0

# DelimitedFiles.writedlm — Function

DelimitedFiles.writedlm(io::IO, d::AbstractIdData, args...; kwargs...)

Записывает идентификационные данные на диск.

# ControlSystemsBase.c2d — Function

c2d(w::AbstractVector{<:Real}, Ts; w_prewarp = 0)
c2d(frd::FRD, Ts; w_prewarp = 0)

Преобразует вектор частот с непрерывным временем w или данные частотных характеристик frd из непрерывного в дискретное время с помощью билинейного преобразования (метод Тастина). Это полезно в тех случаях, когда частотная характеристика получена с помощью анализа частотной характеристики и необходимо использовать функцию subspaceid.