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

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

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

где критерий обычно является взвешенной квадратичной нормой ошибок прогноза. Упорядоченный критерий изменяется как:

Общий особый случай этого когда. Это называется гребенчатой регрессией в статистике, например, смотрите команду ridge в Statistics and Machine Learning Toolbox™.

Полезный образ мыслей о регуляризации, это представляет предварительные знания о неизвестном векторе параметра, и это описывает уверенность в этом знании. (Чем больше, тем более высокая уверенность). Формальная интерпретация в установке Bayesian, это имеет предшествующее распределение, которое является Гауссовым со средней и ковариационной матрицей, где отклонение инноваций.

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

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

Рассмотрите проблему оценки импульсного ответа линейной системы как модель FIR:

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

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

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

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.

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

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

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

Используя ориентировочные стоимости и в 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, заказывают линейную систему с цветным шумом измерения:

где 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. Существует два обычно используемых метода:

Применение этих методов к данным с максимальным порядком, 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: 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*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 к угловой скорости и временной константе:

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

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

От предварительных знаний и опыта мы думаем, что это - приблизительно 4 и является приблизительно 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 и = 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 имеет параметры = 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                               

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

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

Другое использование регуляризации должно численно стабилизировать оценку больших (часто нелинейный) модели. Мы дали записи данных 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, линейных моделей серого поля и моделей Nonlinear ARX. Регуляризация может оказать значительное влияние на качество идентифицированной модели, предоставленной константы регуляризации, Lambda и R выбраны соответственно. Для моделей ARX это может быть сделано очень легко с помощью функции arxRegul. Этот автоматический выбор также подает в специализированный алгоритм оценки пространства состояний ssregest.

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

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