Exponenta Banner
exponenta event banner

Регуляризованная идентификация динамических систем

Открыть сценарий

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

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

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

$$\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 в Toolbox™ статистики и машинного обучения.

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

Поэтому использование регуляризации может быть связано с некоторой предшествующей информацией о системе. Это может быть довольно мягким, как если бы система была стабильной. Переменные регуляризации$\lambda$ и$R$ могут рассматриваться как инструменты для поиска хорошей сложности модели для наилучшего компромисса между смещением и дисперсией. Константы регуляризации$\lambda$$R$ и представлены опциями, называемыми Lambda и R соответственно в 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 образцов. Оценим его по данным, сгенерированным системой. Мы моделируем систему с фильтруемым белым шумом нижних частот в качестве входного сигнала и добавляем небольшое выходное возмущение белого шума с дисперсией 0,0025 к выходному сигналу. Собрано 1000 образцов. Эти данные сохраняются в файле regulateExampleData.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.

Несмотря на 1000 точек данных с очень хорошим отношением сигнал/шум, оценка не впечатляет. Неопределенность в ответе также достаточно велика, как показано значениями 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 с.д.

Поэтому давайте попытаемся достичь хорошего компромисса между отклонениями смещения и отклонениями путем регрессии хребта для модели 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.

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

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

Мы можем сделать еще лучше. Используя понимание того, что истинная импульсная характеристика распадается до нуля и является гладкой, мы можем адаптировать выбор$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» здесь). Если панель инструментов оптимизации недоступна, вместо нее используется простая схема поиска Гаусса-Ньютона; используйте опцию «Advanced.SeeyMethod» для 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-го порядка с цветным шумом измерения:

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

где G(q) - передаточная функция «вход-выход» и H(q) - функция передачи возмущений. Эта система хранится в файле данных reguledExampleData.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. Эти данные сохраняются в reguledExampleData.mat и отображаются ниже.

load regularizationExampleData.mat m0simdata
clf
plot(m0simdata)

Рис. 8: Данные, используемые для оценки.

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

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: Модели оценки для заказов 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*goodnessOfFit(y12,yG,'NRMSE')),...
   sprintf('CV choice: %2.4g%%',100*goodnessOfFit(y15,yG,'NRMSE')),...
   sprintf('AIC choice: %2.4g%%',100*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*goodnessOfFit(yH3,yH,'NRMSE')),...
   sprintf('CV choice: %2.4g%%',100*goodnessOfFit(yH15,yH,'NRMSE')),...
   sprintf('AIC choice: %2.4g%%',100*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*goodnessOfFit(ymr,yG,'NRMSE'))
fprintf('Goodness of fit for noise component of ARX model is: %2.4g%%\n',100*goodnessOfFit(yHmr,yH,'NRMSE'))

Goodness of fit for ARX model is: 16.88%
Goodness of fit for noise component of ARX model is: 21.29%

Получается, что эта регуляризованная модель ARX показывает соответствие истинному G(q) это даже лучше, чем выбор Оракула. Соответствие H(q) более 80%, что также лучше, чем выбор Oracle для лучшей модели шума. Можно утверждать, что mr является моделью высокого порядка (60 состояний), и несправедливо сравнивать ее с моделями пространства состояния нижнего порядка. Но эта модель высокого порядка может быть уменьшена до, скажем, порядка 7 с помощью balred (требуется 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*goodnessOfFit(y12,yG,'NRMSE')),...
   sprintf('High order regularized: %2.4g%%',100*goodnessOfFit(ymr,yG,'NRMSE')),...
   sprintf('Reduced order: %2.4g%%',100*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*goodnessOfFit(yH3,yH,'NRMSE')),...
   sprintf('High order regularized: %2.4g%%',100*goodnessOfFit(yHmr,yH,'NRMSE')),...
   sprintf('Reduced order: %2.4g%%',100*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 что дает модели наилучшее соответствие с данными проверки.

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

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

load regularizationExampleData.mat nldata
opt = nlarxOptions('SearchMethod','lm');
m10 = nlarx(nldata, [1 2 1], sigmoidnet('NumberOfUnits',10),opt);
m20 = nlarx(nldata, [1 2 1], sigmoidnet('NumberOfUnits',20),opt);
m30 = nlarx(nldata, [1 2 1], sigmoidnet('NumberOfUnits',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], sigmoidnet('num',30), opt);
compare(nldata, m30r, m10, m20)

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

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

Заключения

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

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

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

Документация по инструментам идентификации системы

Поддержка