В этом примере показано, как оценить нелинейные модели черных ящиков с множеством входов и множеством выходов (MIMO) на основе данных. В Toolbox™ системной идентификации предлагаются два типа нелинейных моделей черных ящиков - модели Nonlinear ARX и Hammerstein-Wiener.
В этом примере будут использованы данные, сохраненные в файле motorizedcamera.mat. Он содержит 188 образцов данных, собранных с моторной камеры со временем выборки 0,02 секунды. Входной вектор u (t) состоит из 6 переменных: 3 компоненты скорости перемещения в ортогональной системе координат X-Y-Z, прикрепленной к камере [м/с], и 3 компоненты скорости вращения вокруг осей X-Y-Z [рад/с]. Выходной вектор y (t) содержит 2 переменные: положение (в пикселях) точки, которая является изображением, снятым камерой фиксированной точки в 3D пространстве. Создание объекта IDDATAz для хранения загруженных данных:
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. Давайте выберем вейвлетную сеть (волновую сеть) в качестве нелинейной функции отображения для обоих выходов. Функция 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)'}
Количество единиц («вейвлетов») двух функций 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');
В приведенном выше вызове мы использовали строку «treepartion» вместо объекта, чтобы указать функцию отображения. Этот синтаксис облегчает быстрый способ построения модели с тем ограничением, что все функции сопоставления должны быть одного типа и должны использоваться со значениями свойств по умолчанию. В следующем примере конструктор объекта treepartion вызывается для непосредственного создания объекта нелинейной функции 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 может использоваться для просмотра нелинейных откликов различных моделей.
plot(mt1,mtw,ms1,mls)
Обратите внимание, что панель управления в правой части графика позволяет выполнять выбор и настройку регрессора.
Другие функции, такие как RESID, PREDICT и PE, могут использоваться в оценочных моделях так же, как в случае линейных моделей.
Модель Хаммерштейна-Винера состоит из до 3 блоков: блоку линейной передаточной функции предшествует нелинейный статический блок и/или следует другой нелинейный статический блок. Она называется моделью Винера, если первый нелинейный статический блок отсутствует, и моделью Хаммерштейна, если второй нелинейный статический блок отсутствует.
Модели IDNLHW обычно оцениваются с использованием синтаксиса:
M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).
где Orders = [nb nf nk] определяет порядки и задержку линейной передаточной функции. InputNonlinearity и OutputNonlinearity задают нелинейные функции для двух нелинейных блоков. Модель линейной выходной ошибки (OE) соответствует случаю тривиальной нелинейности (UNITGAIN).
Выберем Порядки = [nb bf nk] = [единицы (2,6), единицы (2,6), единицы (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(z, mhw1);
Свойства модели можно визуализировать с помощью команды ПЕЧАТЬ.
Щелкните по блок-схеме, чтобы выбрать вид входной нелинейности (UNL), линейного блока или выходной нелинейности (YNL).
Если выбран вид линейного блока, по умолчанию все 12 каналов отображаются в очень уменьшенных размерах. Выберите один из каналов для лучшего просмотра. Можно выбрать тип графиков в пределах Step response, Bode plot, Impulse response и карты полюсов-нулей.
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
Начальное предположение для линейного интервала SATURATION и нулевого интервала DEADZONE может быть указано непосредственно при создании этих объектов; можно также задать ограничения для этих значений, например, фиксированы ли они для указанных значений (путем установки для атрибута Free значения false) или если их оценки подчинены минимальным/максимальным границам (с использованием атрибутов Минимум (Minimum) и Максимум (Maximum)).
Предположим, мы хотим установить линейный интервал насыщения в [10 200], а нулевой интервал мертвой зоны в [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 пока не оценивается по данным. Можно получить доступ к предельным значениям функций насыщения и мертвой зоны. У них все еще есть свои исходные значения, потому что они еще не оценены по данным.
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
