Этот пример показывает, как оценить нелинейные модели черного ящика с мультивходами с несколькими выходами (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. Необходимо выбрать два важных элемента: регрессоры модели и нелинейные функции отображения.
Регрессоры являются простыми формулами, основанными на задержанных по времени переменных ввода-вывода, самым простым случаем является переменные, отстающие от небольшого набора последовательных значений. Например, если «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)}. В общем случае моделей с выходами NY и входами NU, NN является матрицей с строками NY и строками 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 (wavenet) в качестве функции нелинейного отображения для обоих выходов. Функция NLARX оценивает параметры Нелинейной модели ARX.
NN = [ones(2,2), ones(2,6), ones(2,6)]; % the orders
mw1 = nlarx(z, NN, wavenet);
Исследуйте результат путем сравнения выхода, моделируемой с оценочной моделью и выходом в измеренных данных z:
compare(z,mw1)
Давайте проверим, сможем ли мы лучше справляться с более высокими порядками. Обратите внимание, что при идентификации моделей с использованием расширений базиса функций для выражения нелинейности, количество параметров модели может превысить количество выборок данных. В таких случаях некоторые метрики оценки, такие как Отклонение шума и Окончательная Ошибка Предсказания (FPE), не могут быть определены надежно. В данном примере мы отключаем предупреждение, информирующее нас об этом ограничении.
ws = warning('off','Ident:estimation:NparGTNsamp'); % mw2 = nlarx(z, [ones(2,2), 2*ones(2,6), ones(2,6)], wavenet) 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 =
14×1 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)'}
Числа модулей ("wa вейвлетов) двух функций WAVENET были автоматически выбраны алгоритмом оценки. Эти номера отображаются ниже.
mw2.OutputFcn(1).NonlinearFcn.NumberOfUnits
ans =
27
Количество модулей измерения в функции WAVENET может быть явным образом задано вместо автоматического выбора алгоритмом оценки. Предположим, что мы хотим использовать 10 модули для первой нелинейной функции отображения, и 5 для второй (обратите внимание, что модель для каждого выхода использует свою собственную функцию отображения; массив всех функций отображения хранится в свойстве модели «OutputFcn»).
Fcn1 = wavenet(10); % output function for the first output Fcn2 = wavenet(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
Вместо функции WAVENET могут также использоваться другие нелинейные функции. Давайте попробуем TREEPARTITION для обоих выходов.
mt1 = nlarx(z, nanbnk, 'treepartition');
В вышеописанном вызове мы использовали строку 'treepartition' вместо объекта, чтобы задать функцию отображения. Этот синтаксис облегчает быстрый способ построения модели с ограничением, что все функции отображения должны быть одного типа, и они должны использоваться с их значениями свойств по умолчанию. В следующем примере конструктор объекта treepartition вызывается, чтобы непосредственно создать объект нелинейной функции Tree-partition.
mt2 = nlarx(z, nanbnk, treepartition)
mt2 = <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: Tree Partition with 7 units Output 2: Tree Partition with 7 units Sample time: 0.02 seconds Status: Estimated using NLARX on time domain data "Motorized Camera". Fit to estimation data: [98.56;98.34]% (prediction focus) MSE: 1.615
Точно так же мы можем использовать Sigmoid Network (SIGMOIDNET) в качестве функции отображения. Опции оценки, такие как максимальные итерации (MaxIterations) и отображение итерации, могут быть заданы с помощью nlarxOptions команда.
opt = nlarxOptions('Display','on'); opt.SearchOptions.MaxIterations = 5; ms1 = nlarx(z, nanbnk, 'sigmoidnet', opt);
Этот синтаксис вызова для NLARX очень похож на таковые, используемые ранее - указание данных, порядков и нелинейных функций отображения в качестве первичного входного параметра аргументов. Однако, чтобы изменить опции оценки из их значений по умолчанию, мы создали набор опций, opt, используя команду nlarxOptions и передав ее в команду NLARX в качестве дополнительного входного параметра.
Возможно использовать различные нелинейные функции на разных выходных каналах в одной модели. Предположим, что мы хотим использовать функцию древовидного разбиения, чтобы описать первый выход и использовать вейвлет для второго выхода. Оценка модели показана ниже. Третий входной параметр (нелинейная функция отображения) теперь является массивом двух разных функций.
Fcn1 = treepartition; Fcn2 = wavenet; mtw = nlarx(z, nanbnk, [Fcn1; Fcn2]);
Иногда функция, которая преобразует регрессоры в выходные данные модели, не должна быть нелинейной; это может быть просто взвешенная сумма векторов регрессора плюс необязательное смещение. Это похоже на линейные модели ARX (кроме срока смещения). Отсутствие нелинейности в канале выхода может быть указано путем выбора функции LINEAR. Следующий пример означает, что в model_output (t) = F (y (t-1), u (t-1), u (t-2)), функция F состоит из линейного компонента для первого выхода и нелинейного компонента (SIGMOIDNET) для второго выхода.
Fcn1 = linear; Fcn2 = sigmoidnet(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: Linear Function 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: [1×1 idpack.Channel]
Output: [1×1 idpack.Channel]
LinearFcn: [1×1 nlident.internal.ProjectedLinearFcn]
Offset: [1×1 nlident.internal.Offset]
У нас есть опция сделать его чисто линейным, выбрав не использовать термин смещения. Этот выбор может быть сделан путем фиксации смещения к нулевому значению перед оценкой.
Fcn1.Offset.Value = 0; Fcn1.Offset.Free = false; mlsNoOffset = nlarx(z, nanbnk, [Fcn1; Fcn2], opt);
Различные модели можно сравнить в одной команде COMPARE.
compare(z, mw2, ms1, mls)
График функций может использоваться, чтобы просмотреть нелинейные отклики различных моделей.
plot(mt1,mtw,ms1,mls)
Обратите внимание, что панель управления с правой стороны графика позволяет выбрать и строение регрессор.
Другие функции, такие как RESID, PREDICT и PE, могут использоваться на предполагаемых моделях так же, как и в случае линейных моделей.
Модель Гаммерштейна-Винера состоит из до 3 блоков: блоку линейной передаточной функции предшествует нелинейный статический блок и/или последующий другой нелинейный статический блок. Он называется моделью Винера, если первый нелинейный статический блок отсутствует, и моделью Гаммерштейна, если второй нелинейный статический блок отсутствует.
Модели IDNLHW обычно оцениваются с помощью синтаксиса:
M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).
Where Orders = [nb nf nk] определяет порядки и задержку линейной передаточной функции. InputNonlinearity и OutputNonlinearity задают нелинейные функции для двух нелинейных блоков. Модель линейной выходной ошибки (OE) соответствует случаю тривиальной (UNITGAIN) нелинейности.
Давайте выберем порядки = [nb bf nk] = [one (2,6), one (2,6), one (2,6)]. Это означает, что в линейном блоке каждый выход является суммой 6 передаточных функций первого порядка, управляемых 6 входами. Попробуйте модель Гаммерштейна (без выходной нелинейности) с входной нелинейностью, описанной кусочно-линейной функцией. Заметьте, что введенный один объект PWLINEAR автоматически расширяется на все 6 входных каналов. UNITGAIN указывает на отсутствие нелинейности.
opt = nlhwOptions(); opt.SearchOptions.MaxIterations = 15; NN = [ones(2,6), ones(2,6), ones(2,6)]; InputNL = pwlinear; % input nonlinearity OutputNL = unitgain; % 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 Input 2: Piecewise Linear Input 3: Piecewise Linear Input 4: Piecewise Linear Input 5: Piecewise Linear Input 6: Piecewise Linear 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)];
Давайте попробуем модель Винера. Заметьте, что отсутствие входа нелинейности может быть указано как «[]» вместо объекта UNITGAIN.
opt.SearchOptions.MaxIterations = 10;
mhw2 = nlhw(z, nbnfnk, [], 'pwlinear', 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 Output 2: Piecewise Linear 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
Укажите как входную, так и выходную нелинейности для модели Гаммерштейна-Винера. Как и в случае нелинейных моделей ARX, мы можем использовать вектор символов (а не объект), чтобы задать нелинейные функции отображения.
mhw3 = nlhw(z, nbnfnk,'saturation', 'pwlinear', 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 Input 2: Saturation Input 3: Saturation Input 4: Saturation Input 5: Saturation Input 6: Saturation Output nonlinearities: Output 1: Piecewise Linear Output 2: Piecewise Linear 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
К предельным значениям функции НАСЫЩЕНИЯ можно обращаться следующим образом:
mhw3.InputNonlinearity(1).LinearInterval % view Linear Interval of SATURATION
ans =
0.0103 0.0562
Точно так же к пропускам функции PWLINEAR можно обращаться следующим образом:
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 - нелинейность Sigmoid Network на входном канале # 6 с 5 модулями
Мы можем создать массив объектов нелинейных функций отображения выбранных типов и передать его в функцию оценки NLHW как входную нелинейность.
inputNL = [pwlinear; saturation; deadzone; sigmoidnet; unitgain; sigmoidnet(5)];
inputNL(1).NumberOfUnits = 3;
opt.SearchOptions.MaxIterations = 25;
mhw4 = nlhw(z, nbnfnk, inputNL, [], opt); % "[]" for 4th input denotes no output nonlinearity
Начальное предположение для линейного интервала НАСЫЩЕНИЯ и нулевого интервала DEADZONE может быть непосредственно указано при создании этих объектов; можно также задать ограничения для этих значений, например, привязаны ли они к заданным значениям (путем установки атрибута Free на false) или если их оценки зависят от минимальных/максимальных границ (с помощью атрибутов Minimum и Maximum).
Предположим, что мы хотим задать линейный интервал насыщения [10 200] и нулевой интервал deadzone [11 12] как начальные предположения. Кроме того, мы хотим, чтобы верхний предел насыщения оставался фиксированным. Мы можем достичь этого строения следующим образом.
% Create nonlinear functions with initial guesses for properties. OutputNL1 = saturation([10 200]); OutputNL1.Free(2) = false; % the upper limit is a fixed value OutputNL2 = deadzone([11 12]); mhw5 = idnlhw(nbnfnk, [], [OutputNL1; OutputNL2], 'Ts',z.Ts);
Обратите внимание на использование конструктора объектов модели IDNLHW idnlhw, а не оценщик nlhw. Получившийся объект модели mhw5 еще не оценен из данных. Можно получить доступ к предельным значениям функций насыщения и deadzone. У них все еще есть свои начальные значения, потому что они еще не оценены по данным.
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
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
Output 2: Dead Zone
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
