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

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

На этой странице описываются доступные средства для оценки линейных моделей в пространстве состояний, входы которых имеют следующую форму:

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

Существует несколько методов для идентификации моделей в пространстве состояний: subspaceid, n4sid, newpem и era. subspaceid — это самый современный алгоритм для идентификации на основе подпространства, в то время как n4sid — это более старая реализация. newpem решает задачу погрешности прогнозирования с помощью метода итеративной оптимизации (из пакета Optim.jl) и, как правило, немного более точен, хотя и требует больших вычислительных затрат. Если вы не уверены, какой метод использовать, попробуйте сначала subspaceid (если только данные не получены в результате операции с замкнутым контуром; в этом случае используйте newpem).

Идентификация на основе подпространства с помощью n4sid и subspaceid

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

using ControlSystemIdentification, ControlSystemsBase, Plots
Ts = 0.1
G  = c2d(DemoSystems.resonant(), Ts)
u  = randn(1,1000)
y  = lsim(G,u).y
y .+= 0.01 .* randn.() # Добавляем шум измерения
d  = iddata(y,u,Ts)
sys = subspaceid(d, :auto; verbose=false, zeroD=true)

# либо используем устойчивую версию svd, если в y есть выбросы или отсутствующие значения
# using TotalLeastSquares
# sys = n4sid(d, :auto; verbose=false, svd=x->rpca(x)[3])
bodeplot([G, sys.sys], lab=["True" "" "subspace" ""])
julia
AAAAAElFTkSuQmCC

N4SIDStateSpace is a subtype of AbstractPredictionStateSpace, a statespace object that contains an observer gain matrix sys.K (Kalman filter) as well as estimated covariance matrices etc. Если воспользоваться вместо этого функцией n4sid, то получим следующее:

sys2 = n4sid(d, :auto; verbose=false, zeroD=true)
bodeplot!(sys2.sys, lab=["n4sid" ""])
julia
5lS1eLEAAAAASUVORK5CYII=

subspaceid allows you to choose the weighting between :MOESP, :CVA, :N4SID, :IVM and is generally preferred over n4sid. Обе функции позволяют выбрать функции, применяемые для оценки наименьших квадратов и вычисления SVD, например средства оценки, устойчивые к выбросам, и т. д.

Настройка подгонки модели

Алгоритмы оценки на основе подпространства имеют ряд параметров, которые можно настроить, если начальная подгонка модели неудовлетворительная.

  • Параметр focus определяет цель подгонки модели. По умолчанию используется значение :prediction, при котором минимизируется погрешность прогнозирования. Если при этом значении для устойчивой системы получается неустойчивая модель либо если производительность моделирования низкая, может лучше подойти значение focus = :simulation.

  • Настроить можно также несколько параметров горизонта. С помощью именованного аргумента r выбирается горизонт прогнозирования. Значение должно быть больше порядка модели, по умолчанию nx + 10. Аргументы s1 и s2 определяют прошлые горизонты для выхода и входа соответственно. По умолчанию используется s1 = s2 = r. Прошлые горизонты можно настроить только для subspaceid.

  • (Дополнительно.) Метод, используемый для вычисления SVD, а также выполнения подгонки методом наименьших квадратов, можно изменить с помощью именованных аргументов svd, Aestimator, Bestimator.

  • zeroD позволяет сделать оцениваемую матрицу нулевой (строго собственная модель).

  • Можно выбрать тип веса W: :MOESP, :CVA, :N4SID, :IVM. По умолчанию используется :MOESP.

Дополнительные аргументы и более подробные сведения см. в docstring функций subspaceid и n4sid.

ERA и OKID

Алгоритмы реализации собственных значений и идентификации наблюдателя Калмана доступны в виде функций era и okid. При вызове era с объектом данных автоматически используется функция okid с целью получения марковских параметров для алгоритма ERA.

sys3 = era(d, 2) # У функции era много параметров, которые могут требовать настройки.
bodeplot!(sys3, lab=["ERA" ""])
julia
PHHXZ87tLW1lbnLcNKkSdBGBD0HWoQA9KC33nprz549169ft2MMTz75ZEZGRlRU1DfffGN9bwYAgAGJEAAAwIAGt08AAAAY0CARAgAAGND+Pzn8gzOimbXOAAAAAElFTkSuQmCC

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

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

using ControlSystemIdentification, ControlSystemsBase, Plots, Statistics
Ts = 0.1
G = c2d(tf(1, [1,1,1]), Ts) # Истинная система

# Создаем несколько экспериментов
ds = map(1:5) do i
    u = randn(1, 1000)
    y, t, x = lsim(G, u)
    yn = y + 0.2randn(size(y))
    iddata(yn, u, Ts)
end

Ys = okid.(ds, 2, round(Int, 10/Ts), smooth=true, λ=1)    # Оцениваем импульсную характеристику для каждого эксперимента
Y = mean(Ys)                            # Усредняем все импульсные характеристики

imp = impulse(G, 10)
f1 = plot(imp, lab="True", l=5)
plot!(imp.t, vec.(Ys), lab="Individual estimates", title="Impulse-response estimates")
plot!(imp.t, vec(Y), l=(3, :black), lab="Mean")

models = era.(Ys, Ts, 2, 50, 50)    # Оцениваем модели на основе отдельных экспериментов
meanmodel = era(Y, Ts, 2, 50, 50)   # Оцениваем модель на основе средней импульсной характеристики

f2 = bodeplot([G, meanmodel], lab=["True" "" "Combined estimate" ""], l=2)
bodeplot!(models, lab="Individual estimates", c=:black, alpha=0.5, legend=:bottomleft)

plot(f1, f2)
julia
zkb8JMgQIhAIBGJNg8onEAgEArGmQYYQgUAgEGuafwAANhYkoyLFOgAAAABJRU5ErkJggg==

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

era(ds, 2, 50, 50, round(Int, 10/Ts), p=1, λ=1, smooth=true) # Должно быть идентично приведенной выше усредненной модели
julia
StateSpace{ControlSystemsBase.Discrete{Float64}, Float64}
A =
 0.9714275037672085   -0.08556668857687873
 0.08556668857687827   0.9234713852927754
B =
 -0.0655718881394695
  0.06269691994965396
C =
 -0.655718881394695  -0.6269691994965387
D =
 0.011478353122140328

Sample Time: 0.1 (seconds)
Discrete-time state-space model
output

Метод на основе погрешности прогнозирования (PEM)

Функция pem считается устаревшей и не рекомендуется к использованию; вместо нее предпочтительнее использовать функцию newpem, которая является более универсальной и производительной.

Метод на основе погрешности прогнозирования представляет собой простой, но эффективный алгоритм для идентификации дискретных ЛПП-систем в форме пространства состояний:

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

Результатом идентификации с помощью функции newpem является пользовательский тип с дополнительными полями для ковариационных матриц усиления Калмана и шума.

Идентификация серого ящика

Сведения об оценке линейных и нелинейных моделей с фиксированной структурой см. в описании функции ControlSystemIdentification.nonlinear_pem.

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

Ниже мы генерируем систему и моделируем ее вперед во времени. Затем мы пытаемся оценить модель на основе входной и выходной последовательностей с помощью функции newpem.

using ControlSystemIdentification, ControlSystemsBase, Random, LinearAlgebra
using ControlSystemIdentification: newpem
sys = c2d(tf(1, [1, 0.5, 1]) * tf(1, [1, 1]), 0.1)

Random.seed!(1)
T   = 1000                      # Количество временных шагов
nx  = 3                         # Количество полюсов в истинной системе
nu  = 1                         # Количество входов
x0  = randn(nx)                 # Начальное состояние
sim(sys,u,x0=x0) = lsim(ss(sys), u, x0=x0).y # Вспомогательная функция
u   = randn(nu,T)               # Генерируем случайные входные данные
y   = sim(sys, u, x0)           # Моделируем систему
y .+= 0.01 .* randn.()          # Добавляем шум измерения
d   = iddata(y,u,0.1)

sysh,opt = newpem(d, nx, focus=:prediction) # Оцениваем модель

yh = predict(sysh, d) # Выполняем прогнозирование с помощью оцененной модели
predplot(sysh, d)     # Строим график спрогнозированных и истинных выходных данных
julia
AcTfTWNbzBwsAAAAAElFTkSuQmCC

Дополнительные графики можно найти в интерактивных скриптах с примерами, а несколько примеров — в разделе примеров данной документации.

Аргументы

У алгоритма есть несколько параметров:

  • По умолчанию оптимизация начинается с начального предположения, предоставляемого subspaceid, но его можно переопределить, передав начальное предположение в newpem с помощью именованного аргумента sys0.

  • Параметр focus определяет цель подгонки модели. По умолчанию используется значение :prediction, при котором минимизируется погрешность прогнозирования. Если при этом значении для устойчивой системы получается неустойчивая модель либо если производительность моделирования низкая, может лучше подойти значение focus = :simulation.

  • Регуляризатор можно передать с помощью именованного аргумента regularizer.

  • Сделать модель устойчивой можно с помощью аргумента stable = true.

  • Матрицу можно сделать нулевой с помощью аргумента zeroD = true.

  • Добиться необходимого баланса между точностью прогнозирования и производительностью моделирования можно путем оптимизации погрешности прогнозирования на шаге . По умолчанию используется значение , которое соответствует среднеквадратичной погрешности прогнозирования. Его можно изменить с помощью именованного аргумента h. Чем выше значение h, тем больше вычислительных затрат требуется для решения задачи оптимизации.

Дополнительные аргументы и более подробные сведения см. в docstring функции newpem.

Внутренние компоненты

На внутреннем уровне для оптимизации параметров системы используется пакет Optim.jl, который вычисляет градиенты (и при необходимости гессианы) посредством автоматического дифференцирования. Параметры решателя Optim можно контролировать, передавая именованные аргументы в newpem и передавая созданный вручную объект решателя. Решатель по умолчанию — BFGS()

Фильтрация, прогнозирование и моделирование

При оценке моделей иногда можно выбрать «цель» оценки: производительность прогнозирования (:prediciton) или моделирования (:simulation). При выборе моделирования приоритетом является точность низкочастотных свойств, особенно при интеграции систем, а при выборе моделирования — точность модели на более высоких частотах. Если на систему влияют существенные входные возмущения или если система является неустойчивой, в качестве цели предпочтительнее выбирать прогнозирование.

При проверке оцененных моделей их можно имитировать с помощью функции lsim из пакета ControlSystemsBase.jl или с помощью функции simulate. Модель можно также преобразовать в объект KalmanFilter из пакета LowLevelParticleFilters.jl, вызвав KalmanFilter(sys), после чего можно выполнять фильтрацию, сглаживание и другие операции с помощью средств, доступных для KalmanFilter.

Кроме того, имеются перечисленные ниже вспомогательные функции.

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

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

API Statespace

# ControlSystemIdentification.newpem — Function

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,
)
julia

Новая реализация метода на основе ошибки прогнозирования (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
julia

Если начальное предположение 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
julia

Обратите внимание, 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)]
julia

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

# ControlSystemIdentification.subspaceid — Function

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,
)
julia

Оценивает модель пространства состояний с помощью идентификации на основе подпространства. Доступно несколько различных алгоритмов, основанных на подпространстве, которые можно выбрать с помощью ключевого слова 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, и использует эту функцию (по умолчанию) для формирования начального предположения для оптимизации.

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

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

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

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

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

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

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
)
julia

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

Если результаты неудовлетворительны, попробуйте изменить 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.n4sid — Function

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

Оценивает модель пространства состояний с помощью метода 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.era — Function

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

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

Аргументы

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

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

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

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

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

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)
julia

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

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

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

Аргументы

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

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

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

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

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

# ControlSystemIdentification.okid — Function

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

Идентификация с помощью фильтра наблюдателя Калмана. Возвращает размер 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, для полного решения по методу наименьших квадратов (ошибки в переменных)

# ControlSystemsBase.observer_predictor — Function

observer_predictor(sys::N4SIDStateSpace; h=1)
julia

Возвращает систему предикторов

с входным уравнением [B-KD K] * [u; y]

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

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

# ControlSystemsBase.observer_controller — Function

observer_controller(sys::AbstractPredictionStateSpace, L)
julia

Возвращает контроллер обратной связи по измерению, который принимает y и формирует управляющий сигнал u = -Lx̂. См. также описание ff_controller.

# ControlSystemIdentification.prediction_error — Function

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

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

# ControlSystemIdentification.prediction_error_filter — Function

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

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

# ControlSystemIdentification.noise_model — Function

noise_model(sys::AbstractPredictionStateSpace)
julia

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

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

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

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

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