В этом примере показано, как оценить мультивход мультивыводится (MIMO) нелинейные модели черного квадрата из данных. Два типа нелинейных моделей черного квадрата предлагаются в System Identification Toolbox™ - Нелинейный ARX и модели Хаммерстайна-Винера.
Данные, сохраненные в файле motorizedcamera.mat, будут использоваться в этом примере. Это содержит 188 выборок данных, собранных из моторизованной камеры с шагом расчета 0,02 секунд. Входной вектор u (t) состоит из 6 переменных: 3 скоростных компонента перевода в ортогональной системе координат X-Y-Z, зафиксированной к камере [m/s] и 3 скоростным компонентам вращения вокруг X-Y-Z [rad/s]. Выходной вектор y (t) содержит 2 переменные: положение (в пикселях) точки, которая является изображением, взятым камерой фиксированной точки на 3D пробеле. Мы создаем объект IDDATA z содержать загруженные данные:
load motorizedcamera z = iddata(y, u, 0.02, 'Name', 'Motorized Camera', 'TimeUnit', 's');
Давайте сначала попробуем нелинейные модели ARX. Должны быть выбраны два важных элемента: регрессоры модели и нелинейные функции отображения.
Регрессоры являются простыми формулами на основе задержанных временем переменных I/O, самый простой случай, являющийся переменными, изолированными маленьким набором последовательных значений. Например, если "u" от имени входной переменной и "y" имя выходных переменных, то набор регрессора в качестве примера может быть {y (t-1), y (t-2), u (t), u (t-1), u (t-2)}, где "t" обозначает переменную времени. Другой пример включающие полиномиальные термины мог быть {y (t-2) ^4, y (t-2) *u (t-1), u (t-4) ^2}. Более комплексные, произвольные формулы в задержанных переменных также возможны.
Самый легкий способ задать линейные регрессоры при помощи "матрицы порядков". Эта матрица принимает форму NN = [na nb nk], и это задает тем, сколько задержек каждый выход (na) и каждая входная переменная (nb, nk) задерживаются. Это - та же идея, которая используется при оценке линейных моделей ARX (см. ARX, IDPOLY). Например, NN = [2 3 4] подразумевает набор регрессора {y (t-2), u (t-4), u (t-5), u (t-6)}. В общем случае моделей с Нью-Йорком выходные параметры и входные параметры NU, NN является матрицей с нью-йоркскими строками и строками NY+2*NU.
Чтобы запуститься, давайте выберем порядки NN = [na nb nk] = [единицы (2,2), единицы (2,6), единицы (2,6)]. Это означает, что каждая выходная переменная предсказана всеми выходными и входными переменными, каждый задержанный 1 выборкой. Уравнение модели может быть записано как y_i (t) = F_i (y1 (t-1), y2 (t-1), u1 (t-1), u2 (t-1), u3 (t-1)), i = 1, 2. Давайте выберем Wavelet Network (idWaveletNetwork объект) как нелинейная функция отображения для обоих выходных параметров. Функциональный NLARX оценивает параметры модели Nonlinear ARX.
NN = [ones(2,2), ones(2,6), ones(2,6)]; % the orders
mw1 = nlarx(z, NN, idWaveletNetwork);
Исследуйте результат путем сравнения выхода, симулированного с предполагаемой моделью и выходом в результатах измерений z:
compare(z,mw1)
Давайте проверять, можем ли мы добиться большего успеха с высшими порядками. Обратите внимание на то, что при идентификации моделей с помощью расширений основной функции, чтобы описать нелинейность, количество параметров модели может превысить количество выборок данных. В таких случаях некоторые метрики оценки, такие как Шумовое Отклонение и Итоговая ошибка предсказания (FPE) не могут быть определены надежно. Для текущего примера мы выключаем предупреждение, сообщающее нам об этом ограничении.
ws = warning('off','Ident:estimation:NparGTNsamp'); % mw2 = nlarx(z, [ones(2,2), 2*ones(2,6), ones(2,6)], idWaveletNetwork) compare(z,mw2)
mw2 = <strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong> Inputs: u1, u2, u3, u4, u5, u6 Outputs: y1, y2 Regressors: Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6 Output functions: Output 1: Wavelet network with 27 units Output 2: Wavelet network with 25 units Sample time: 0.02 seconds Status: Estimated using NLARX on time domain data "Motorized Camera". Fit to estimation data: [99.22;99.15]% (prediction focus) FPE: 0.1262, MSE: 0.4453
Вторая модель mw2 довольно хороша. Таким образом давайте сохраним этот выбор порядков модели в следующих примерах.
nanbnk = [ones(2,2), 2*ones(2,6), ones(2,6)]; % final orders
Чтобы просмотреть набор регрессора, подразумеваемый этим выбором порядков, используйте команду GETREG:
getreg(mw2)
ans =
14x1 cell array
{'y1(t-1)'}
{'y2(t-1)'}
{'u1(t-1)'}
{'u1(t-2)'}
{'u2(t-1)'}
{'u2(t-2)'}
{'u3(t-1)'}
{'u3(t-2)'}
{'u4(t-1)'}
{'u4(t-2)'}
{'u5(t-1)'}
{'u5(t-2)'}
{'u6(t-1)'}
{'u6(t-2)'}
Количества модулей ("вейвлеты") двух функций idWaveletNetwork были автоматически выбраны алгоритмом оценки. Эти числа отображены ниже.
mw2.OutputFcn(1).NonlinearFcn.NumberOfUnits mw2.OutputFcn(2).NonlinearFcn.NumberOfUnits
ans =
27
ans =
25
Количество модулей в функции idWaveletNetwork может быть явным образом задано вместо того, чтобы быть автоматически выбранным алгоритмом оценки. Предположим, что мы хотим использовать 10 модулей для первой нелинейной функции отображения, и 5 для второго (обратите внимание, что модель для каждого выхода использует свою собственную функцию отображения; массив всех функций отображения хранится в свойстве "OutputFcn" модели).
Fcn1 = idWaveletNetwork(10); % output function for the first output Fcn2 = idWaveletNetwork(5); % output function for the second output mw3 = nlarx(z, nanbnk, [Fcn1; Fcn2])
mw3 = <strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong> Inputs: u1, u2, u3, u4, u5, u6 Outputs: y1, y2 Regressors: Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6 Output functions: Output 1: Wavelet network with 10 units Output 2: Wavelet network with 5 units Sample time: 0.02 seconds Status: Estimated using NLARX on time domain data "Motorized Camera". Fit to estimation data: [99.01;98.89]% (prediction focus) FPE: 0.2273, MSE: 0.7443
Вместо функции idWaveletNetwork могут также использоваться другие нелинейные функции. Давайте попробуем idTreePartition за оба выходных параметров.
mt1 = nlarx(z, nanbnk, idTreePartition);
В вышеупомянутом вызове древовидная функция отображения раздела с настройкой по умолчанию используется для обоих выходные параметры.
Точно так же мы можем использовать Сигмоидальную Сеть (idSigmoidNetwork объект) как функция отображения. Опции оценки, такие как максимальные итерации (MaxIterations) и отображение итерации могут быть заданы с помощью nlarxOptions команда.
opt = nlarxOptions('Display','on'); opt.SearchOptions.MaxIterations = 5; ms1 = nlarx(z, nanbnk, idSigmoidNetwork, opt);
Этот синтаксис вызова для NLARX очень похож на тех используемых прежде - определение данных, порядков и нелинейных функций отображения в качестве их аргументов первичного входного параметра. Однако для того, чтобы изменить опции оценки от их значений по умолчанию, мы создали набор опции, opt, использование nlarxOptions команды и передало его команде NLARX как дополнительный входной параметр.
Возможно использовать различные нелинейные функции на различных выходных каналах в той же модели. Предположим, что мы хотим использовать древовидную функцию раздела, чтобы описать первый выход и использовать сеть вейвлета для второго выхода. Оценку модели показывают ниже. Третий входной параметр (нелинейная функция отображения) является теперь массивом двух различных функций.
Fcn1 = idTreePartition; Fcn2 = idWaveletNetwork; mtw = nlarx(z, nanbnk, [Fcn1; Fcn2]);
Иногда функция, которая сопоставляет регрессоры с выходом модели, не должна быть нелинейной; это могла просто быть взвешенная сумма векторов регрессора плюс дополнительное смещение. Это похоже на линейные модели ARX (за исключением термина смещения). Отсутствие нелинейности в выходном канале может быть обозначено путем выбора функции idLinear. Следующий пример означает, что, в model_output (t) = F (y (t-1), u (t-1), u (t-2)), функция F состоит из линейного компонента для первого выхода и нелинейного компонента (idSigmoidNetwork) для второго выхода.
Fcn1 = idLinear; Fcn2 = idSigmoidNetwork(2); opt.Display = 'off'; % do not show estimation progress anymore mls = nlarx(z, nanbnk, [Fcn1; Fcn2], opt)
mls = <strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong> Inputs: u1, u2, u3, u4, u5, u6 Outputs: y1, y2 Regressors: Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6 Output functions: Output 1: None Output 2: Sigmoid network with 2 units Sample time: 0.02 seconds Status: Estimated using NLARX on time domain data "Motorized Camera". Fit to estimation data: [98.72;98.79]% (prediction focus) FPE: 0.5594, MSE: 1.05
Нет никакой нелинейности в модели для первого выхода. Это не, однако, чисто линейно из-за присутствия термина смещения.
disp(mls.OutputFcn(1))
<strong>Linear Function</strong>
Inputs: y1(t-1), y2(t-1), u1(t-1), u1(t-2), u2(t-1), u2(t-2), u3(t-1), u3(t-2), u4(t-1), u4(t-2), u5(t-1), u5(t-2), u6(t-1), u6(t-2)
Output: y1
Nonlinear Function: None
Linear Function: initialized to [48.3 -32.2 -0.229 -0.0705 -0.113 -0.0516 -0.182 0.297 0.199 -0.133 -0.337 0.583 -0.448 0.167]
Output Offset: initialized to 109
Input: '<Function inputs>'
Output: '<Function output>'
LinearFcn: '<Linear function parameters>'
Offset: '<Offset parameters>'
У нас есть опция создания его чисто линейный путем желания не использовать термин смещения. Этот выбор может быть сделан путем фиксации смещения к нулевому значению перед оценкой.
Fcn1.Offset.Value = 0; Fcn1.Offset.Free = false; mlsNoOffset = nlarx(z, nanbnk, [Fcn1; Fcn2], opt);
Если у вас есть доступ к Statistics and Machine Learning Toolbox™, вы можете также использовать функции гауссова процесса (GP) и сложенный в мешок или повысили древовидные ансамбли, чтобы создать ваши модели Nonlinear ARX. GP представлены объектом idGaussianProcess, в то время как ансамбли дерева регрессии представлены объектом idTreeEnsemble. GP используют ядро в квадрате экспоненциальное по умолчанию. Другие типы ядра могут быть заданы при создании объекта. В примере здесь, мы используем основанный на ядре GP 'ARDMatern52' для первого выхода и уволенный древовидный ансамбль для второго выхода.
if exist('fitrgp','file')==2 warning('off','stats:lasso:MaxIterReached'); NL1 = idGaussianProcess('ardmatern52'); NL2 = idTreeEnsemble; % estimate model mML = nlarx(z, nanbnk, [NL1; NL2]) end
mML = <strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong> Inputs: u1, u2, u3, u4, u5, u6 Outputs: y1, y2 Regressors: Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6 Output functions: Output 1: Gaussian process function using a ARDMatern52 kernel. Output 2: Bagged Regression Tree Ensemble Sample time: 0.02 seconds Status: Estimated using NLARX on time domain data "Motorized Camera". Fit to estimation data: [99.43;98.57]% (prediction focus) MSE: 0.7846
Различные модели могут быть сравнены в той же команде COMPARE. Например, сравните симулированные выходные параметры моделей mw2, ms1, mls, и mlsNoOffset к измеренному выходу в наборе данных "z".
compare(z, mw2, ms1, mls, mlsNoOffset)
График функций может использоваться, чтобы просмотреть нелинейную функцию отображения "форма" для различных моделей.
plot(mt1,mtw,ms1,mls)
Обратите внимание на то, что панель управления на правой стороне графика допускает выбор регрессора и настройку.
Другие функции, такие как RESID, PREDICT и PE могут использоваться на предполагаемых моделях таким же образом как в случае линейных моделей.
Модель Хаммерстайна-Винера состоит максимум из 3 блоков: линейному блоку передаточной функции предшествует нелинейный статический блок и/или сопровождает другой нелинейный статический блок. Это называется Винеровской моделью, если первый нелинейный статический блок отсутствует, и модель Хаммерстайна, если второй нелинейный статический блок отсутствует.
Модели IDNLHW обычно оцениваются с помощью синтаксиса:
M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).
где Порядки = [nb nf nk] задают порядки и задержку линейной передаточной функции. InputNonlinearity и OutputNonlinearity задают нелинейные функции для двух нелинейных блоков. Линейная модель ошибки на выходе (OE) соответствует случаю тривиальных (модульное отображение усиления) нелинейность.
Давайте выберем Orders = [nb bf nk] = [единицы (2,6), единицы (2,6), единицы (2,6)]. Это означает, что в линейном блоке каждый выход является суммой 6 передаточных функций первого порядка, управляемых 6 входными параметрами. Попробуйте модель Хаммерстайна (никакая выходная нелинейность) с входной нелинейностью, описанной кусочной линейной функцией. Заметьте, что вводимый один объект idPiecewiseLinear автоматически расширен до всех 6 входных каналов. idUnitGain указывает на отсутствие нелинейности.
opt = nlhwOptions(); opt.SearchOptions.MaxIterations = 15; NN = [ones(2,6), ones(2,6), ones(2,6)]; InputNL = idPiecewiseLinear; % input nonlinearity OutputNL = idUnitGain; % output nonlinearity mhw1 = nlhw(z, NN, InputNL, OutputNL, opt)
mhw1 = Hammerstein-Wiener model with 2 outputs and 6 inputs Linear transfer function matrix corresponding to the orders: nb = [1 1 1 1 1 1; 1 1 1 1 1 1] nf = [1 1 1 1 1 1; 1 1 1 1 1 1] nk = [1 1 1 1 1 1; 1 1 1 1 1 1] Input nonlinearities: Input 1: Piecewise linear with 10 break-points Input 2: Piecewise linear with 10 break-points Input 3: Piecewise linear with 10 break-points Input 4: Piecewise linear with 10 break-points Input 5: Piecewise linear with 10 break-points Input 6: Piecewise linear with 10 break-points Output nonlinearities: Output 1: absent Output 2: absent Sample time: 0.02 seconds Status: Estimated using NLHW on time domain data "Motorized Camera". Fit to estimation data: [98.46;97.93]% FPE: 7.928, MSE: 2.216
Исследуйте результат с COMPARE.
compare(z, mhw1);
Свойства модели могут визуализироваться Командой plot.
Нажмите на блок-схему, чтобы выбрать представление входной нелинейности (UNL), линейного блока или выходной нелинейности (YNL).
Когда линейное представление блока выбрано, по умолчанию все 12 каналов показывают в очень уменьшаемых размерах. Выберите один из каналов, чтобы иметь лучшее представление. Возможно выбрать тип графиков в Переходном процессе, Диаграммы Боде, Импульсной характеристики и нулевой полюсом карты.
plot(mhw1)
Результат, показанный COMPARE, был довольно хорош, таким образом давайте сохраним те же порядки модели в следующих испытаниях.
nbnfnk = [ones(2,6), ones(2,6), ones(2,6)];
Давайте попробуем Винеровскую модель. Заметьте, что отсутствие входной нелинейности может быть обозначено" []" вместо объекта idUnitGain.
opt.SearchOptions.MaxIterations = 10; mhw2 = nlhw(z, nbnfnk, [], idPiecewiseLinear, opt)
mhw2 = Hammerstein-Wiener model with 2 outputs and 6 inputs Linear transfer function matrix corresponding to the orders: nb = [1 1 1 1 1 1; 1 1 1 1 1 1] nf = [1 1 1 1 1 1; 1 1 1 1 1 1] nk = [1 1 1 1 1 1; 1 1 1 1 1 1] Input nonlinearities: Input 1: absent Input 2: absent Input 3: absent Input 4: absent Input 5: absent Input 6: absent Output nonlinearities: Output 1: Piecewise linear with 10 break-points Output 2: Piecewise linear with 10 break-points Sample time: 0.02 seconds Status: Estimated using NLHW on time domain data "Motorized Camera". Fit to estimation data: [73.85;71.36]% FPE: 1.314e+05, MSE: 503.8
Укажите на обе нелинейности ввода и вывода для модели Хаммерстайна-Винера.
mhw3 = nlhw(z, nbnfnk,idSaturation, idPiecewiseLinear, opt)
mhw3 = Hammerstein-Wiener model with 2 outputs and 6 inputs Linear transfer function matrix corresponding to the orders: nb = [1 1 1 1 1 1; 1 1 1 1 1 1] nf = [1 1 1 1 1 1; 1 1 1 1 1 1] nk = [1 1 1 1 1 1; 1 1 1 1 1 1] Input nonlinearities: Input 1: Saturation with linear interval: [0.0103 0.0562] Input 2: Saturation with linear interval: [-0.00143 0.000909] Input 3: Saturation with linear interval: [-0.0947 -0.0185] Input 4: Saturation with linear interval: [-0.00385 0.00527] Input 5: Saturation with linear interval: [0.0195 0.13] Input 6: Saturation with linear interval: [-0.00302 0.000387] Output nonlinearities: Output 1: Piecewise linear with 10 break-points Output 2: Piecewise linear with 10 break-points Sample time: 0.02 seconds Status: Estimated using NLHW on time domain data "Motorized Camera". Fit to estimation data: [86.88;84.55]% FPE: 1.111e+04, MSE: 137.3
К предельным значениям функции idSaturation можно получить доступ можно следующим образом:
mhw3.InputNonlinearity(1).LinearInterval % view Linear Interval of saturation function
ans =
0.0103 0.0562
Точно так же к точкам останова функции idPiecewiseLinear можно получить доступ можно следующим образом:
mhw3.OutputNonlinearity(1).BreakPoints
ans =
Columns 1 through 7
17.1233 34.2491 51.3726 68.4968 85.6230 102.7478 119.8742
2.6184 16.0645 45.5178 41.9621 62.3246 84.9038 112.2970
Columns 8 through 10
136.9991 154.1238 171.2472
135.4543 156.1016 173.2701
Различные нелинейные функции могут быть смешаны в той же модели. Предположим, что мы хотим модель с: - никакая нелинейность на любых выходных каналах ("Модель Хаммерстайна") - Кусочная линейная нелинейность на входном канале № 1 с 3 модулями - нелинейность Насыщения на входном канале № 2 - Мертвая Зональная нелинейность на входном канале № 3 - Сигмоидальная Сетевая нелинейность на входном канале № 4 - Никакая нелинейность (заданный объектом unitgain) на входном канале № 5 - Сигмоидальная Сетевая нелинейность на входном канале № 6 с 5 модулями
Мы можем создать массив нелинейных функциональных объектов отображения выбранных типов и передать его функции оценки NLHW как входная нелинейность.
inputNL = [idPiecewiseLinear; ... idSaturation; ... idDeadZone; ... idSigmoidNetwork; ... idUnitGain; ... idSigmoidNetwork(5)]; inputNL(1).NumberOfUnits = 3; opt.SearchOptions.MaxIterations = 25; mhw4 = nlhw(z, nbnfnk, inputNL, [], opt); % "[]" as the 4th input argument denotes no output nonlinearity
Исходное предположение для линейного интервала насыщения и нулевого интервала мертвой зоны может быть непосредственно обозначено, когда эти объекты создаются; можно также задать ограничения на эти значения такой как, фиксируются ли они к их заданным значениям (атрибутом установки Free ко лжи), или если их оценки подвергаются минимальным/максимальным границам (использующий Минимальные и Максимальные атрибуты).
Предположим, что мы хотим установить линейный интервал насыщения на [10 200] и нулевой интервал deadzone к [11 12] как исходные предположения. Кроме того, мы хотим, чтобы верхний предел насыщения остался зафиксированным. Мы можем достигнуть этой настройки можно следующим образом.
% Create nonlinear functions with initial guesses for properties. OutputNL1 = idSaturation([10 200]); OutputNL1.Free(2) = false; % the upper limit is a fixed value OutputNL2 = idDeadZone([11 12]); mhw5 = idnlhw(nbnfnk, [], [OutputNL1; OutputNL2], 'Ts',z.Ts);
Заметьте использование конструктора объекта модели IDNLHW idnlhw, не средство оценки nlhw. Получившийся объект модели mhw5 еще не оценивается из данных. К предельным значениям насыщения и функций deadzone можно получить доступ. У них все еще есть свои начальные значения, потому что они еще не оцениваются из данных.
mhw5.OutputNonlinearity(1).LinearInterval % view the linear interval on saturation at first output channel after estimation mhw5.OutputNonlinearity(2).ZeroInterval % view the zero interval on dead zone at second output channel after estimation
ans =
10 200
ans =
11 12
Каковы значения этих пределов после оценки?
opt.SearchOptions.MaxIterations = 15; mhw5 = nlhw(z, mhw5, opt) % estimate the model from data mhw5.OutputNonlinearity(1).LinearInterval % show linear interval on saturation at first output channel after estimation mhw5.OutputNonlinearity(2).ZeroInterval % show zero interval on dead zone at second output channel after estimation
mhw5 =
Hammerstein-Wiener model with 2 outputs and 6 inputs
Linear transfer function matrix corresponding to the orders:
nb = [1 1 1 1 1 1; 1 1 1 1 1 1]
nf = [1 1 1 1 1 1; 1 1 1 1 1 1]
nk = [1 1 1 1 1 1; 1 1 1 1 1 1]
Input nonlinearities:
Input 1: absent
Input 2: absent
Input 3: absent
Input 4: absent
Input 5: absent
Input 6: absent
Output nonlinearities:
Output 1: Saturation with linear interval: [10 200] Output 2: Dead Zone with zero interval: [11 12]
Sample time: 0.02 seconds
Status:
Estimated using NLHW on time domain data "Motorized Camera".
Fit to estimation data: [27.12;6.857]%
FPE: 3.373e+06, MSE: 4666
ans =
9.9974 200.0000
ans =
11.0020 12.0011
Модели различной природы (IDNLARX и IDNLHW) могут быть сравнены в той же команде COMPARE.
compare(z,mw2,mhw1)
warning(ws) % reset the warning state
