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

Открыть скрипт

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

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

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

$$\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$тем выше доверие). Формальная интерпретация в байесовской обстановке является тем, что имеет$\theta$ предшествующее распределение, которое является Гауссовым со средней и$\theta^*$ ковариационной матрицей, $\sigma^2/\lambda R^{-1}$где является$\sigma^2$ отклонением нововведений.

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

Смещение - компромисс отклонений в конечной импульсной характеристике моделирования

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

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

Их оценивает команда: m = arx(z,[0 nb 0]). Выбор порядка <reservedrangesplaceholder0> является компромиссом между смещением (большим 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 к выходу. Собрано 1000 выборки. Эти данные сохраняются в файле regularizationExampleData.mat и показаны на фигуре 2.

load regularizationExampleData.mat eData
plot(eData)

Фигура. 2. Данные, используемые для оценки.

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

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: Истинная импульсная характеристика вместе с оценкой для порядок <reservedrangesplaceholder0>.

Несмотря на 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 с.д.

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

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

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

[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: Истинная импульсная характеристика вместе с упорядоченной по гребню оценкой для порядок <reservedrangesplaceholder0> точности.

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

Автоматическое определение констант регуляризации для конечной импульсной характеристики Моделей

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

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

Рассмотрим системную m0, которая является линейной системой 30-го порядка с цветным измерительным шумом:

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

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

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

Улов состоит в том, чтобы определить хороший порядок <reservedrangesplaceholder0>. Существует два широко используемых метода:

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

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

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

Тест «Oracle»: В дополнение к тестам CV и AIC можно также проверить, какой порядок <reservedrangesplaceholder0> подгонка между истинной импульсной характеристикой 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) это даже лучше, чем выбор 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*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 в каждом компоненте. Мы также сохраняем данные моделирования без шума для той же модели в целях сравнения. Данные показаны на фигура.

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) и изображают подгонку к данным без шума на фигура.

[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: Диаграмма Боде of 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 значение для использования потребует перекрестных валидаций тестов.

Заключения

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

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

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

Документация по System Identification Toolbox

Поддержка

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