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

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

Непараметрическая оценка

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

using ControlSystemIdentification, ControlSystemsBase, Plots
T          = 100000
h          = 1
sim(sys,u) = lsim(sys, u, 1:T)[1][:]
σy         = 0.5
sys        = tf(1,[1,2*0.1,0.1])
ωn         = sqrt(0.3)
sysn       = tf(σy*ωn,[1,2*0.1*ωn,ωn^2])

u  = randn(1, T)
y  = sim(sys, u)
yn = y + sim(sysn, randn(size(u))) # Добавляем шум, отфильтрованный через sysn
d  = iddata(y,u,h)
dn = iddata(yn,u,h)

InputOutput data of length 100000, 1 outputs, 1 inputs, Ts = 1

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

k = coherence(d)  # Должно быть близко к 1, если система линейная и лишена шума
k = coherence(dn) # Если в системе присутствует шум измерения, получаются немного более низкие значения

Передаточную функцию можно также оценить спектральными методами. Главной отправной точкой здесь является функция tfest, которая возвращает оценку передаточной функции и оценку энергетической спектральной плотности шума (имейте в виду, что единица измерения PSD возводится в квадрат относительно передаточной функции, поэтому в приведенном ниже коде при построении графика берется √N):

G,N = tfest(dn)
bodeplot([sys,sysn], exp10.(range(-3, stop=log10(pi), length=200)), layout=(1,3), plotphase=false, subplot=[1,2,2], ylims=(0.1,300), linecolor=:blue)

coherenceplot!(dn, subplot=3)
plot!(G, subplot=1, lab="G Est", alpha=0.3, title="Process model")
plot!(√N, subplot=2, lab="N Est", alpha=0.3, title="Noise model")

fv36nTt3fsQ2BAUFVVRUODk55ebm4rkZbxl+0opwFrOYxSxmMYu3Kn1iFrOYxSxmMYsXxawinMUsZjGLWfyk8X8wduygrrNuSwAAAABJRU5ErkJggg==

На рисунке слева показана величина ЛАЧХ истинной системы вместе с оценкой (зашумленной), а на рисунке в центре — оцененная модель шума. На рисунке справа представлена функция когерентности (coherenceplot), значение которой близко к 1 везде, кроме резонансного пика шума log10(sqrt(0.3)) = -0.26.

Дополнительные сведения см. в интерактивных скриптах с примерами.

Параметрическая оценка

Передаточные функции

Для оценки параметрической рациональной передаточной функции на основе данных частотной области вызовите tfest с объектом FRD и начальным предположением для модели системы. Начальное предположение определяет количество коэффициентов в числителе и знаменателе оцениваемой модели.

G0 = tf(1.0, [1,1,1]) # Начальное предположение
G = tfest(d::FRD, G0)

На внутреннем уровне Optim использует градиентный оптимизатор для нахождения оптимального соответствия кривой Боде системы. Используемый по умолчанию оптимизатор BFGS можно изменить; см. docstring функции ?tfest.

Сравнение оценки во временной и частотной областях см. в этом интерактивном скрипте.

Если приведенную выше задачу трудно решить, можно параметризовать модель, например путем расширения базиса Лагерра:

basis = laguerre_oo(1, 50) # Используем функции с базисом 50; итоговый порядок модели можно понизить с помощью baltrunc
Gest,p = tfest(d::FRD, basis)

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

Пространство состояний

Функция subspaceid работает с данными частотной области (а также с данными временной области). При передаче объекта InputOutputFreqData (который можно создать с помощью функции iddata) автоматически применяется метод для частотной области. Кроме того, можно передать объект частотной характеристики FRD, который автоматически преобразуется в InputOutputFreqData. Если данные частотной характеристики получены в результате анализа частотной характеристики, может потребоваться произвести билинейное преобразование по оси частот объекта данных, чтобы преобразовать непрерывную ось частот в дискретную. Пример:

Ts    = 0.01 # Интервал выборки
frd_d = c2d(frd_c::FRD, Ts) # Выполняем билинейное преобразование в дискретный вектор частот
Ph, _ = subspaceid(frd_d, Ts, nx)

Это можно сделать автоматически, передав аргумент bilinear_transform = true в функцию subspaceid.

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

using ControlSystemIdentification, ControlSystemsBase, Plots
ny,nu,nx = 2,3,5                        # Количество выходов, входов и состояний
Ts = 1                                  # Интервал выборки
G = ssrand(ny,nu,nx; Ts, proper=true)   # Создаем случайную систему

N = 200             # Количество частотных точек
w = freqvec(Ts, N)
frd = FRD(w, G)     # Создаем объект с данными частотной характеристики (на практике он должен получаться в результате эксперимента)

Gh, x0 = subspaceid(frd, G.Ts, nx; zeroD=true) # Подгоняем частотную характеристику
sigmaplot([G, Gh], w[2:end], lab=["True system" "Estimated model"])

Спектральная оценка на основе модели

Процедуры оценки моделей можно использовать для оценки спектрограмм. В данном пакете некоторые методы из пакета DSP.jl расширяются так, чтобы вторым аргументом они могли принимать функцию оценки. Для создания такой функции предоставляется функция model_spectrum. Ее использование проиллюстрировано ниже.

using ControlSystemIdentification, DSP
T  = 1000
fs = 1
s = sin.((1:1/fs:T) .* 2pi/10) + 0.5randn(T)
S1 = spectrogram(s,window=hanning, fs=fs)            # Стандартная спектрограмма
estimator = model_spectrum(ar,1/fs,6)
S2 = spectrogram(s,estimator,window=rect, fs=fs)     # Спектрограмма на основе модели
plot(plot(S1,title="Standard Spectrogram"),plot(S2,title="AR Spectrogram")) # Требуется пакет LPVSpectral.jl

# 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, но в противном случае является смещенным (например, при идентификации в замкнутом контуре).

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 для преобразования возможно неминимально фазовой системы в минимально фазовую.

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

κ² = 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.laguerre_oo — Function

laguerre_oo(a::Number, Nq)

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

# ControlSystemIdentification.model_spectrum — Function

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