Упорядоченная идентификация динамических систем

Этот пример показывает преимущества регуляризации для идентификации линейных и нелинейных моделей.

Что такое Регуляризация

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

$\hat{\theta} = arg \min_\theta V_N(\theta)$

где критерий обычно является взвешенной квадратичной нормой ошибок предсказания $\varepsilon(t,\theta)$ . $L_2$ Упорядоченный критерий изменяется как:

$\hat\theta = arg \min_\theta V_N(\theta) + \lambda (\theta-\theta^*)^TR(\theta-\theta^*)$

Общий особый случай этого когда $\theta^*=0, R=I$ . Это называется гребенчатой регрессией в статистике, например, смотрите ridge команда в Statistics and Machine Learning Toolbox™.

Полезный образ мыслей о регуляризации, это $\theta^*$ представляет предварительные знания о неизвестном векторе параметра, и это $\lambda*R$ описывает доверие к этому знанию. (Чем больше $\lambda*R$ , тем более высокое доверие). Формальная интерпретация в установке Bayesian, это $\theta$ имеет предшествующее распределение, которое является Гауссовым со средней $\theta^*$ и ковариационной матрицей $\sigma^2/\lambda R^{-1}$ , где $\sigma^2$ отклонение инноваций.

Использование регуляризации может поэтому быть соединено с некоторой предшествующей информацией о системе. Это могло быть довольно мягко, как этот, система устойчива. Переменные регуляризации $\lambda$ и $R$ могут рассматриваться как инструменты, чтобы найти хорошую сложность модели для лучшего компромисса между смещением и отклонением. Константы регуляризации $\lambda$ и $R$ представлены опциями под названием Lambda и R соответственно в System Identification Toolbox™. Выбором $\theta^*$ управляет Nominal опция регуляризации.

Смещение - Компромисс Отклонения в КИХ-моделировании

Рассмотрите задачу оценки импульсной характеристики линейной системы как модель FIR:

$y(t)=\sum^{nb}_{k=0} g(k)u(t-k)$

Они оцениваются командой: m = arx(z,[0 nb 0]). Выбор порядка nb компромисс между смещением (большой nb требуется, чтобы получать медленно затухающие импульсные характеристики без слишком большой ошибки) и отклонение (большой nb дает много параметров, чтобы оценить, который дает большое отклонение).

Давайте проиллюстрируем его с симулированным примером. Мы выбираем простой фильтр Баттерворта второго порядка как систему:

$G(z) = \frac{0.02008 + 0.04017 z^-1 + 0.02008 z^-2}{ 1 - 1.561 z^-1 + 0.6414 z^-2}$

Его импульсную характеристику показывают в рисунке 1:

trueSys = idtf([0.02008 0.04017 0.02008],[1 -1.561 0.6414],1);
[y0,t] = impulse(trueSys);
plot(t,y0)
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')

Рисунок 1: истинная импульсная характеристика.

Импульсная характеристика затухла, чтобы обнулить меньше чем после 50 выборок. Давайте оценим его из данных, сгенерированных системой. Мы симулируем систему с отфильтрованным белым шумом lowpass, как введено и давайте добавим маленький белый шум выходное воздействие с отклонением 0.0025 к выходу. Собраны 1 000 выборок. Эти данные сохранены в regularizationExampleData.mat файле и показанные в рисунке 2.

load regularizationExampleData.mat eData
plot(eData)

Рисунок 2: данные используются для оценки.

Определить хорошее значение для nb мы в основном должны попробовать несколько значений, и некоторой процедурой валидации оценивают, который является лучшим. Это может быть сделано несколькими способами, но поскольку мы знаем истинную систему в этом случае, мы можем определить теоретически самое лучшее значение путем испытания всех моделей с nb=1,...,50 и найдите, какой имеет лучшую подгонку к истинной импульсной характеристике. Такой тест показывает тот nb = 13 дает лучшую ошибочную норму (mse = 0.2522) к импульсной характеристике. Эту предполагаемую импульсную характеристику показывают вместе с истинной в рисунке 3.

nb = 13;
m13 = arx(eData,[0 nb 0]);

[y13,~,~,y13sd] = impulse(m13,t);
plot(t,y0,t,y13)
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True system','13:th order FIR model')

Рисунок 3: истинная импульсная характеристика вместе с оценкой для порядка nb = 13.

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

plot(t,y0,t,y13,t,y13+y13sd,'r:',t,y13-y13sd,'r:')
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True system','13:th order FIR model','Bounds')

Рисунок 4: Предполагаемый ответ с доверительными границами, соответствующими 1 s.d.

Давайте поэтому попытаемся достигнуть хорошего компромисса отклонения смещения гребенчатой регрессией для модели FIR порядка 50. Используйте arxOptions сконфигурировать константы регуляризации. Для этого осуществления мы применяем простой штраф $||\theta||^2$ .

aopt = arxOptions;
aopt.Regularization.Lambda = 1;
m50r = arx(eData, [0 50 0], aopt);

Получившаяся оценка имеет ошибочную норму 0,1171 к истинной импульсной характеристике и показана в рисунке 5 наряду с доверительными границами.

[y50r,~,~,y50rsd] = impulse(m50r,t);
plot(t,y0,t,y50r,t,y50r+y50rsd,'r:',t,y50r-y50rsd,'r:')
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True system','50:th order regularized estimate')

Рисунок 5: истинная импульсная характеристика вместе с упорядоченной гребнем оценкой для порядка nb = 50.

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

Автоматическое определение констант регуляризации для КИХ-моделей

Мы можем сделать еще лучше. При помощи понимания, которое истинная импульсная характеристика затухает, чтобы обнулить и является гладкой, мы можем адаптировать выбор $R, \lambda$ в соответствии с данными. Это достигается arxRegul функция.

[L,R] = arxRegul(eData,[0 50 0],arxRegulOptions('RegularizationKernel','TC'));
aopt.Regularization.Lambda = L;
aopt.Regularization.R = R;
mrtc = arx(eData, [0 50 0], aopt);
[ytc,~,~,ytcsd] = impulse(mrtc,t);

arxRegul использование fmincon от Optimization Toolbox™, чтобы вычислить гиперпараметры, сопоставленные с ядром регуляризации ("TC" здесь). Если Optimization Toolbox не доступен, простая схема поиска Ньютона Гаусса используется вместо этого; используйте "Усовершенствованный. SearchMethod" опция arxRegulOptions выбрать метод поиска явным образом. Предполагаемые гиперпараметры затем используются, чтобы получить значения $R$ и $\lambda$ .

Используя ориентировочные стоимости $R$ и $\lambda$ в ARX приводит к ошибочной норме 0,0461, и ответ показывают в рисунке 6. Этот вид настроенной регуляризации - то, что достигается также impulseest команда. Когда рисунок показывает, подгонка к импульсной характеристике, а также отклонению значительно уменьшается по сравнению с неупорядоченными оценками. Цена является смещением в оценке ответа, которая, кажется, незначительна для этого примера.

plot(t,y0,t,ytc,t,ytc+ytcsd,'r:',t,ytc-ytcsd,'r:')
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True system','50:th order tuned regularized estimate')

Рисунок 6: истинная импульсная характеристика вместе с настроенной упорядоченной оценкой для порядка nb = 50.

Используя упорядоченные модели ARX для оценки моделей в пространстве состояний

Считайте систему m0, который является 30:th, заказывают линейную систему с цветным шумом измерения:

$$y(t) = G(q)u(t) + H(q)e(t)$$

где G(q) передаточная функция входа к выходу и H(q) передаточная функция воздействия. Эта система хранится в regularizationExampleData.mat файле данных. Импульсные характеристики G(q) и H(q) показаны в рисунке 7.

load regularizationExampleData.mat m0
m0H = noise2meas(m0); % the extracted noise component of the model
[yG,t] = impulse(m0);
yH = impulse(m0H,t);

clf
subplot(211)
plot(t, yG)
title('Impulse Response of G(q)'), ylabel('Amplitude')

subplot(212)
plot(t, yH)
title('Impulse Response of H(q)'), ylabel('Amplitude')
xlabel('Time (seconds)')

Рисунок 7: импульсные характеристики G(q) (верхняя часть) и H(q) (нижняя часть).

Мы собрали 210 точек данных путем симуляции m0 с входом u белого шума с отклонением 1, и уровень шума e с отклонением 0.1. Эти данные сохранены в regularizationExampleData.mat и отображены на графике ниже.

load regularizationExampleData.mat m0simdata
clf
plot(m0simdata)

Рисунок 8: данные, которые будут использоваться для оценки.

Оценить импульсные характеристики m0 из этих данных мы можем естественно использовать модели в пространстве состояний в инновационной форме (или эквивалентно модели ARMAX) и вычислить импульсную характеристику с помощью impulse команда как прежде. Для вычисления модели в пространстве состояний мы можем использовать синтаксис, такой как:

mk = ssest(m0simdata, k, 'Ts', 1);

Выгода должна определить хорошее состояние k. Существует два обычно используемых метода:

CV перекрестной проверки: Оцените mk для k = 1,...,maxo использование первой половины данных ze = m0simdata(1:150) и оцените подгонку к второй половине данных zv = m0simdata(151:end) использование compare команда: [~,fitk] = compare(zv, mk, compareOptions('InitialCondition', 'z')). Определите порядок k это максимизирует подгонку. Затем повторно оцените модель с помощью целой записи данных.

Используйте AIC критерия Akaike: Оцените модели для порядков k = 1,...,maxo использование целого набора данных, и затем выбирает ту модель, которая минимизирует aic(mk).

Применение этих методов к данным с максимальным порядком maxo = 30 показывает, что тот CV выбирает k = 15 и AIC выбирает k = 3.

Тест "Oracle": В дополнение к CV и тестам AIC, можно также проверять на какой порядок k подгонка между истинной импульсной характеристикой G(q) (или H(q)) и предполагаемая модель максимизируется. Это, конечно, требует знания истинной системы m0 который непрактичен. Однако, если мы действительно продолжаем это сравнение для нашего примера где m0 известен, мы находим тот k = 12 дает лучший припадок импульсной характеристики предполагаемой модели тому из m0 (= |G (q) |). Точно так же мы находим тот k = 3 дает лучший припадок импульсной характеристики шумового компонента предполагаемой модели тому из шумового компонента m0 (= |H (q) |). Наборы тестов Oracle контрольная точка для сравнения качества моделей, сгенерированных при помощи различных порядков и параметров регуляризации.

Давайте сравним импульсные характеристики, вычисленные для различных критериев выбора порядка:

m3 = ssest(m0simdata, 3, 'Ts', 1);
m12 = ssest(m0simdata, 12, 'Ts', 1);
m15 = ssest(m0simdata, 15, 'Ts', 1);

y3 = impulse(m3, t);
y12 = impulse(m12, t);
y15 = impulse(m15, t);

plot(t,yG, t,y12, t,y15, t,y3)
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True G(q)',...
   sprintf('Oracle choice: %2.4g%%',100*(1-goodnessOfFit(y12,yG,'NRMSE'))),...
   sprintf('CV choice: %2.4g%%',100*(1-goodnessOfFit(y15,yG,'NRMSE'))),...
   sprintf('AIC choice: %2.4g%%',100*(1-goodnessOfFit(y3,yG,'NRMSE'))))

Рисунок 9: истинная импульсная характеристика G(q) по сравнению с предполагаемыми моделями различных порядков.

yH3 = impulse(noise2meas(m3), t);
yH15 = impulse(noise2meas(m15), t);

plot(t,yH, t,yH3, t,yH15, t,yH3)
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True H(q)',...
   sprintf('Oracle choice: %2.4g%%',100*(1-goodnessOfFit(yH3,yH,'NRMSE'))),...
   sprintf('CV choice: %2.4g%%',100*(1-goodnessOfFit(yH15,yH,'NRMSE'))),...
   sprintf('AIC choice: %2.4g%%',100*(1-goodnessOfFit(yH3,yH,'NRMSE'))))

Рисунок 10: истинная импульсная характеристика H(q) по сравнению с предполагаемыми шумовыми моделями различных порядков.

Мы видим, что подгонки, столь же хорошей как 83%, возможно достигнуть для G(q) среди моделей в пространстве состояний, но процедуры отбора порядка может не найти тот лучший порядок.

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

aopt = arxOptions;
[Lambda, R] = arxRegul(m0simdata, [5 60 0], arxRegulOptions('RegularizationKernel','TC'));
aopt.Regularization.R = R;
aopt.Regularization.Lambda = Lambda;
mr = arx(m0simdata, [5 60 0], aopt);
nmr = noise2meas(mr);
ymr = impulse(mr, t);
yHmr = impulse(nmr, t);
fprintf('Goodness of fit for ARX model is: %2.4g%%\n',100*(1-goodnessOfFit(ymr,yG,'NRMSE')))
fprintf('Goodness of fit for noise component of ARX model is: %2.4g%%\n',100*(1-goodnessOfFit(yHmr,yH,'NRMSE')))

Goodness of fit for ARX model is: 83.12%
Goodness of fit for noise component of ARX model is: 78.71%

Оказывается, что эта упорядоченная модель ARX показывает подгонку к истинному G(q) это является четным лучше, чем выбор Oracle. Подгонка к H(q) больше чем 80%, который также лучше что выбор Oracle порядка для лучшей шумовой модели. Это могло быть обсуждено тот mr старший разряд (60 состояний) модель, и несправедливо сравнить его с моделями в пространстве состояний более низкоуровневыми. Но эта высокого уровня модель может уменьшаться до, скажем, порядка 7 при помощи balred команда (требует Control System Toolbox™):

mred7 = balred(idss(mr),7);
nmred7 = noise2meas(mred7);
y7mr = impulse(mred7, t);
y7Hmr = impulse(nmred7, t);

Рисунки 11 и 12 показывают, как упорядоченные модели упорядоченного и уменьшаемого порядка соответствуют выбору Oracle порядка пространства состояний для ssest без любой потери точности.

plot(t,yG, t,y12, t,ymr, t,y7mr)
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True G(q)',...
   sprintf('Oracle choice: %2.4g%%',100*(1-goodnessOfFit(y12,yG,'NRMSE'))),...
   sprintf('High order regularized: %2.4g%%',100*(1-goodnessOfFit(ymr,yG,'NRMSE'))),...
   sprintf('Reduced order: %2.4g%%',100*(1-goodnessOfFit(y7mr,yG,'NRMSE'))))

Рисунок 11: упорядоченные модели по сравнению с выбором Oracle для G(q).

plot(t,yH, t,yH3, t,yHmr, t,y7Hmr)
xlabel('Time (seconds)'), ylabel('Amplitude'), title('Impulse Response')
legend('True H(q)',...
   sprintf('Oracle choice: %2.4g%%',100*(1-goodnessOfFit(yH3,yH,'NRMSE'))),...
   sprintf('High order regularized: %2.4g%%',100*(1-goodnessOfFit(yHmr,yH,'NRMSE'))),...
   sprintf('Reduced order: %2.4g%%',100*(1-goodnessOfFit(y7Hmr,yH,'NRMSE'))))

Рисунок 12: упорядоченные модели по сравнению с выбором Oracle для H(q).

Естественный вопрос спросить состоит в том, является ли выбором порядков в модели ARX столь же чувствительное решение как порядок модели в пространстве состояний в ssest. Простой тест, использование, например, arx(z,[10 50 0], aopt), показывает только незначительные изменения в припадке G(q).

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

Перечисленные выше шаги оценки старшей модели ARX, сопровождаемой преобразованием в пространство состояний и сокращением к желаемому порядку, могут быть автоматизированы с помощью ssregest команда. ssregest значительно упрощает эту процедуру, также упрощая другие полезные опции, такие как поиск оптимального порядка и точно настраивая структуры модели спецификацией значений задержки и сквозного соединения. Здесь мы просто повторно оцениваем упрощенную модель, похожую на mred7 использование ssregest:

opt = ssregestOptions('ARXOrder',[5 60 0]);
mred7_direct = ssregest(m0simdata, 7, 'Feedthrough', true, opt);
compare(m0simdata, mred7, mred7_direct)

Рисунок 13: Сравнение ответов моделей в пространстве состояний к данным об оценке.

h = impulseplot(mred7, mred7_direct, 40);
showConfidence(h,1) %  1 s.d. "zero interval"
hold on
s = stem(t,yG,'r');
s.DisplayName = 'True G(q)';
legend('show')

Рисунок 14: Сравнение импульсных характеристик моделей в пространстве состояний.

В рисунке 14 доверительную границу только показывают для модели mred7_direct поскольку это не было вычислено для модели mred7. Можно использовать translatecov команда для генерации доверительных границ для произвольных преобразований (здесь balred) из идентифицированных моделей. Обратите внимание также что ssregest команда не требует, чтобы вы ввели значение опции "ARXOrder". Это делает автоматический выбор на основе длины данных, когда никакое значение явным образом не установлено.

Основное смещение - компромисс отклонения в серых моделях поля

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

Рассмотрите двигатель постоянного тока (см., например, iddemo7) со статическим усилением G к скорости вращения и постоянной времени $\tau$ :

$G(s) = \frac{G}{s(1+s\tau)}$

В форме пространства состояний мы имеем:

$\dot{x_1} = x_2$

$\dot{x_2} = -1/\tau \cdot x_2$

$$ y = x + e $$

где $x = [x_1; x_2]$ вектор состояния состоит из угла $x_1$ и скорости $x_2$ . Мы наблюдаем оба состояния в шуме, как предложено выходным уравнением.

От предварительных знаний и опыта мы думаем, что это $G$ - приблизительно 4 и $\tau$ является приблизительно 1. Мы собираемся в motorData 400 точек данных от системы, со значительным количеством шума (стандартное отклонение e 50 в каждом компоненте. Мы также сохраняем бесшумные данные моделирования для той же модели в целях сравнения. Данные показывают в рисунке 15.

load regularizationExampleData.mat motorData motorData_NoiseFree
t = motorData.SamplingInstants;
subplot(311)
plot(t,[motorData_NoiseFree.y(:,1),motorData.y(:,1)])
ylabel('Output 1')
subplot(312)
plot(t,[motorData_NoiseFree.y(:,2),motorData.y(:,2)])
ylabel('Output 2')
subplot(313)
plot(t,motorData_NoiseFree.u) % input is the same for both datasets
ylabel('Input')
xlabel('Time (seconds)')

Рисунок 15: зашумленные данные, которые будут использоваться для серой оценки поля, наложенной на бесшумные данные моделирования, которые будут использоваться для проверок. Сверху донизу: Угол, Скорость вращения, Входное напряжение.

Истинными значениями параметров в этой симуляции является G = 2.2 и $\tau$ = 0.8. Чтобы оценить модель, мы создаем idgrey файл модели DCMotorODE.m.

type('DCMotorODE')

function [A,B,C,D] = DCMotorODE(G,Tau,Ts)
%DCMOTORODE ODE file representing the dynamics of a DC motor parameterized
%by gain G and time constant Tau.
%
%   [A,B,C,D,K,X0] = DCMOTORODE(G,Tau,Ts) returns the state space matrices
%   of the DC-motor with time-constant Tau and static gain G. The sample
%   time is Ts.
%
%   This file returns continuous-time representation if input argument Ts
%   is zero. If Ts>0, a discrete-time representation is returned.
%
% See also IDGREY, GREYEST.

%   Copyright 2013 The MathWorks, Inc.

A = [0 1;0 -1/Tau];
B = [0; G/Tau];
C = eye(2);
D = [0;0];
if Ts>0 % Sample the model with sample time Ts
   s = expm([[A B]*Ts; zeros(1,3)]);
   A = s(1:2,1:2);
   B = s(1:2,3);
end

idgrey объект затем создается как:

mi = idgrey(@DCMotorODE,{'G', 4; 'Tau', 1},'cd',{}, 0);

где мы вставили предполагаемое значение параметров как начальные значения. Эта модель настроена к информации в наблюдаемых данных при помощи greyest команда:

m = greyest(motorData, mi)

m =
  Continuous-time linear grey box model defined by @DCMotorODE function:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
           x1      x2
   x1       0       1
   x2       0  -1.741
 
  B = 
          u1
   x1      0
   x2  3.721
 
  C = 
       x1  x2
   y1   1   0
   y2   0   1
 
  D = 
       u1
   y1   0
   y2   0
 
  K = 
       y1  y2
   x1   0   0
   x2   0   0
 
  Model parameters:
   G = 2.138
   Tau = 0.5745
 
Parameterization:
   ODE Function: @DCMotorODE
   (parametrizes both continuous- and discrete-time equations)
   Disturbance component: none
   Initial state: 'auto'
   Number of free coefficients: 2
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                 
Estimated using GREYEST on time domain data "motorData".
Fit to estimation data: [29.46;4.167]%                  
FPE: 6.074e+06, MSE: 4908

Модель m имеет параметры $\tau$ = 0.57 и G = 2.14 и воспроизводит данные, показан в рисунке 16.

copt = compareOptions('InitialCondition', 'z');
[ymi, fiti] = compare(motorData, mi, copt);
[ym, fit] = compare(motorData, m, copt);
t = motorData.SamplingInstants;
subplot(211)
plot(t, [motorData.y(:,1), ymi.y(:,1), ym.y(:,1)])
axis tight
ylabel('Output 1')
legend({'Measured output',...
   sprintf('Initial: %2.4g%%',fiti(1)),...
   sprintf('Estimated: %2.4g%%',fit(1))},...
   'Location','BestOutside')
subplot(212)
plot(t, [motorData.y(:,2), ymi.y(:,2), ym.y(:,2)])
ylabel('Output 2')
axis tight
legend({'Measured output',...
   sprintf('Initial: %2.4g%%',fiti(2)),...
   sprintf('Estimated: %2.4g%%',fit(2))},...
   'Location','BestOutside')

Рисунок 16: Измеренный выход и выходные параметры модели для первоначальных и предполагаемых моделей.

В этом симулированном случае у нас есть также доступ к бесшумным данным (motorData_NoiseFree) и изобразите подгонку к бесшумным данным в рисунке 17.

[ymi, fiti] = compare(motorData_NoiseFree, mi, copt);
[ym, fit] = compare(motorData_NoiseFree, m, copt);
subplot(211)
plot(t, [motorData_NoiseFree.y(:,1), ymi.y(:,1), ym.y(:,1)])
axis tight
ylabel('Output 1')
legend({'Noise-free output',...
   sprintf('Initial: %2.4g%%',fiti(1)),...
   sprintf('Estimated: %2.4g%%',fit(1))},...
   'Location','BestOutside')
subplot(212)
plot(t, [motorData_NoiseFree.y(:,2), ymi.y(:,2), ym.y(:,2)])
ylabel('Output 2')
axis tight
legend({'Noise-free output',...
   sprintf('Initial: %2.4g%%',fiti(2)),...
   sprintf('Estimated: %2.4g%%',fit(2))},...
   'Location','BestOutside')

Рисунок 17: бесшумный выход и выходные параметры модели для первоначальных и предполагаемых моделей.

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

opt = greyestOptions;
opt.Regularization.Lambda = 100;
opt.Regularization.R = [1, 1000]; % second parameter better known than first
opt.Regularization.Nominal = 'model';
mr = greyest(motorData, mi, opt)

mr =
  Continuous-time linear grey box model defined by @DCMotorODE function:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
           x1      x2
   x1       0       1
   x2       0  -1.119
 
  B = 
          u1
   x1      0
   x2  2.447
 
  C = 
       x1  x2
   y1   1   0
   y2   0   1
 
  D = 
       u1
   y1   0
   y2   0
 
  K = 
       y1  y2
   x1   0   0
   x2   0   0
 
  Model parameters:
   G = 2.187
   Tau = 0.8938
 
Parameterization:
   ODE Function: @DCMotorODE
   (parametrizes both continuous- and discrete-time equations)
   Disturbance component: none
   Initial state: 'auto'
   Number of free coefficients: 2
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                 
Estimated using GREYEST on time domain data "motorData".
Fit to estimation data: [29.34;3.848]%                  
FPE: 6.135e+06, MSE: 4933

Мы здесь сказали процессу оценки, что имеем некоторое доверие к начальным значениям параметров и верим больше в наше предположение $\tau$ , чем в нашем предположении G. Получившаяся упорядоченная оценка mr рассматривает эту информацию вместе с информацией в результатах измерений. Они взвешены вместе со справкой Lambda и R. В рисунке 18 показано, как получившаяся модель может воспроизвести выход. Безусловно, упорядоченная модель делает лучшее задание, чем оба первоначальная модель (к которому параметры "привлечены"), и неупорядоченная модель.

[ymr, fitr] = compare(motorData_NoiseFree, mr, copt);
subplot(211)
plot(t, [motorData_NoiseFree.y(:,1), ymi.y(:,1), ym.y(:,1), ymr.y(:,1)])
axis tight
ylabel('Output 1')
legend({'Noise-free output',...
   sprintf('Initial: %2.4g%%',fiti(1)),...
   sprintf('Estimated: %2.4g%%',fit(1)),...
   sprintf('Regularized: %2.4g%%',fitr(1))},...
   'Location','BestOutside')
subplot(212)
plot(t, [motorData_NoiseFree.y(:,2), ymi.y(:,2), ym.y(:,2), ymr.y(:,2)])
ylabel('Output 2')
axis tight
legend({'Noise-free output',...
   sprintf('Initial: %2.4g%%',fiti(2)),...
   sprintf('Estimated: %2.4g%%',fit(2)),...
   sprintf('Regularized: %2.4g%%',fitr(2))},...
   'Location','BestOutside')

Рисунок 18: бесшумный измеренный выход и выходные параметры модели для начальной буквы, оцененных и упорядоченных моделей.

Упорядоченная оценка также уменьшала отклонение параметра по сравнению с неупорядоченными оценками. Это показывают более трудные доверительные границы на Диаграмме Боде mr сравните с тем из m:

clf
showConfidence(bodeplot(m,mr,logspace(-1,1.4,100)),3) % 3 s.d. region
legend('show')

Рисунок 19: Диаграмма Боде m и mr с доверительными границами

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

Разделите данные в две части - оценка и данные о валидации
Вычислите упорядоченную модель с помощью данных об оценке для различных значений Lambda
Оцените, как хорошо эти модели могут воспроизвести данные о валидации: сведите в таблицу значения подгонки NRMSE, обеспеченные compare команда или goodnessOfFit команда.
Выберите тот Lambda который дает модель с лучшей подгонкой к данным о валидации.

Использование регуляризации к большим нелинейным моделям Robustify

Другое использование регуляризации должно численно стабилизировать оценку больших (часто нелинейный) модели. Мы дали записи данных nldata это имеет нелинейную динамику. Мы пробуем нелинейную модель ARX символа нейронной сети со все большим количеством модулей:

load regularizationExampleData.mat nldata
opt = nlarxOptions('SearchMethod','lm');
m10 = nlarx(nldata, [1 2 1], idSigmoidNetwork(10), opt);
m20 = nlarx(nldata, [1 2 1], idSigmoidNetwork(20), opt);
m30 = nlarx(nldata, [1 2 1], idSigmoidNetwork(30), opt);

compare(nldata, m10, m20) % compare responses of m10, m20 to measured response

Рисунок 20: график Сравнения для моделей m10 и m20.

fprintf('Number of parameters (m10, m20, m30): %s\n',...
   mat2str([nparams(m10),nparams(m20),nparams(m30)]))
compare(nldata, m30, m10, m20) % compare all three models
axis([1 800 -57 45])

Number of parameters (m10, m20, m30): [54 104 154]

Рисунок 21: график Сравнения для моделей m10, m20 и m30.

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

opt.Regularization.Lambda = 1e-8;
m30r = nlarx(nldata, [1 2 1], idSigmoidNetwork(30), opt);
compare(nldata, m30r, m10, m20)

Рисунок 22: график Сравнения для моделей m10, m20 и упорядоченная модель m30r.

Подгонка к данным об оценке значительно улучшилась для модели с 30 нейронами. Как обсуждено прежде, систематический поиск Lambda значение, чтобы использовать потребовало бы тестов перекрестной проверки.

Заключения

Мы обсудили преимущество регуляризации для оценки моделей FIR, линейных моделей серого ящика и моделей Nonlinear ARX. Регуляризация может оказать значительное влияние на качество идентифицированной модели, предоставленной константы регуляризации Lambda и R выбраны соответственно. Для моделей ARX это может быть сделано очень легко с помощью arxRegul функция. Этот автоматический выбор также подает в специализированный алгоритм оценки пространства состояний ssregest.

Для других типов оценок необходимо использовать основанный на перекрестной проверке поиск, чтобы определить Lambda. Для структурированных моделей, таких как серые модели поля, R может использоваться, чтобы указать на надежность соответствующего начального значения параметра. Затем использование Nominal опция регуляризации, можно объединить предварительные знания значений параметров с информацией в данных.

Опции регуляризации доступны для всех линейных и нелинейных моделей включая передаточные функции и модели процессов, пространство состояний и полиномиальные модели, Нелинейный ARX, Хаммерстайна-Винера и линейные/нелинейные серые модели поля.

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