В этом примере показана идентификация нелинейных моделей черного ящика с одним входом-одним выходом (SISO) с использованием измеренных данных ввода-вывода. В примере используются измеренные данные из двухбачковой системы для изучения различных структур моделей и вариантов идентификации.
В этом примере переменные данных, сохраненные в twotankdata.mat будет использован файл. Он содержит 3000 выборок данных ввода-вывода двухбачковой системы, сформированных с использованием времени выборки 0,2 секунды. Вход u (t) представляет собой напряжение [V], приложенное к насосу, который генерирует приток в верхний резервуар. Довольно небольшое отверстие в нижней части этого верхнего резервуара дает выход, который поступает в нижний резервуар, и выход y (t) системы двух резервуаров представляет собой уровень жидкости [m] в нижнем резервуаре. Создание объекта IDDATAz для инкапсуляции загруженных данных:
Этот набор данных также используется в нелинейном сером поле примера «Two Tank System: C MEX-File Modeling of Time-Continuous SISO System», где предполагаются физические уравнения, описывающие систему. Этот пример относится к моделям черного ящика, поэтому не предполагается знание фактических физических законов.
load twotankdata z = iddata(y, u, 0.2, 'Name', 'Two tank system');
Разбейте этот набор данных на 3 подмножества одинакового размера, каждое из которых содержит 1000 образцов.
Мы будем оценивать модели с z1и оценить их на z2 и z3. z1 аналогичен z2, но не z3. Поэтому в некотором смысле z2 поможет нам оценить способность интерполяции моделей, в то время как z3 поможет нам оценить их способность экстраполяции.
z1 = z(1:1000); z2 = z(1001:2000); z3 = z(2001:3000); plot(z1,z2,z3) legend('Estimation','Validation 1', 'Validation 2')
Первым шагом в оценке нелинейных моделей черного ящика является выбор набора регрессоров. Регрессоры представляют собой простые формулы, основанные на запаздывающих по времени входных/выходных переменных, простейший случай - переменные, запаздывающие на небольшой набор последовательных значений. Например, если «u» - имя входной переменной, а «y» - имя выходной переменной, то примерным набором регрессоров может быть {y (t-1), y (t-2), u (t), u (t-1), u (t-2)}, где «t» обозначает переменную времени. Другим примером, включающим полиномиальные члены, могут быть {abs (y (t-2)), y (t-2) * u (t-1), u (t-4) ^ 2}. Более сложные, произвольные формулы в задержанных переменных также возможны.
В Toolbox™ Идентификация системы (System Identification) наборы регрессоров не должны создаваться явно; большая коллекция может быть создана неявно, просто задав их форму (например, порядок многочленов), переменные и их задержки. Объекты спецификации linearRegressor, polynomialRegressor и customRegressor используются для этой цели. Например, R = linearRegressor ({'y', 'u'}, {[1 2 4], [2 10]}) подразумевает набор регрессоров {y (t-1), y (t-2), y (t-4), u (t-2), u (t-10)}. Аналогично, R = полиномный регрессор ({'y'}, {[1 2]}, 3) подразумевает полиномиальные регрессоры 3-го порядка в переменной 'y' с лагами 1 и 2, то есть множество {y (t-1) ^ 3, y (t-2) ^ 3}. Доступны дополнительные варианты конфигурации, такие как использование только абсолютных значений и/или разрешение смешивания лагов в формулах-регрессорах.
Когда в модели используются только линейные регрессоры, основанные на смежных лагах, их часто легче задать с помощью матрицы порядка, а не с помощью linearRegressor объект. Матрица порядка [na nb nk] определяет количество прошлых выходов, прошлых входов и входную задержку, используемую в формулах регрессора. Это та же идея, которая используется при оценке линейных моделей ARX (см. ARX, IDPOLY). Например, NN = [2 3 4] подразумевает, что выходная переменная использует задержки (1,2), а входная переменная использует задержки (4,5,6), ведущие к набору регрессоров {y (t-1), y (t-2), u (t-4), u (t-5), u (t-6)}.
Обычно заказы моделей выбираются методом проб и ошибок. Иногда полезно сначала попробовать линейные модели ARX, чтобы получить подсказки о лучших заказах. Итак, давайте сначала попробуем определить оптимальные порядки для линейной модели ARX, используя функции arxstruc и selstruc.
V = arxstruc(z1,z2,struc(1:5, 1:5,1:5)); % select best order by Akaike's information criterion (AIC) nn = selstruc(V,'aic')
nn =
5 1 3
Критерий AIC имеет выбранный nn = [na nb nk] = [5 1 3], что означает, что в выбранной структуре модели ARX y (t) прогнозируется 6 регрессорами y (t-1), y (t-1),..., y (t-5) и u (t-3). Мы постараемся использовать эти значения при оценке нелинейных моделей.
В модели IDNLARX выходной сигнал модели является нелинейной функцией регрессоров, таких как y (t) = Fcn (y (t-1), y (t-1),..., y (t-5), u (t-3)).
Модели IDNLARX обычно оцениваются с помощью синтаксиса:
M = NLARX(Data, Orders, Fcn).
Она аналогична линейной команде ARX с дополнительным аргументом «Fcn», указывающим тип нелинейной функции. Нелинейная модель ARX создает выходные данные с помощью двухступенчатого преобразования: 1. Сопоставьте обучающие данные (входные-выходные сигналы) с набором регрессоров. Численно набор регрессоров представляет собой двойную матрицу R с одним столбцом данных для каждого регрессора. 2. Сопоставьте регрессоры с одним выходом, y = Fcn ®, где Fcn () является скалярной нелинейной функцией.
Для оценки модели IDNLARX после выбора заказов моделей необходимо выбрать тип используемой функции нелинейного отображения Fcn (). Сначала попробуем использовать функцию вейвлет-сети (wavelet Network).
mw1 = nlarx(z1,[5 1 3], wavenet);
Вейвлет-сеть использует комбинацию вейвлет-единиц для отображения регрессоров на выходные данные модели. Модель mw1 является объектом @ idnalrx. Он сохраняет нелинейную функцию отображения (здесь) в своемOutputFcn собственность. Количество этих единиц для использования может быть указано заранее, или мы можем оставить его алгоритму оценки, чтобы определить оптимальное значение автоматически. Свойство NonlinearFcn.NumberOfUnits хранит выбранное количество единиц в модели.
NLFcn = mw1.OutputFcn; NLFcn.NonlinearFcn.NumberOfUnits
ans =
3
Анализ результата путем сравнения отклика модели mw1 к измеренному выходному сигналу в наборе данных z1.
compare(z1,mw1);
Первая модель была не вполне удовлетворительной, поэтому попробуем изменить некоторые свойства модели. Вместо того, чтобы позволить алгоритму оценки автоматически выбирать количество единиц в функции WAVENET, это число может быть явно указано.
mw2 = nlarx(z1,[5 1 3], wavenet(8));
По умолчанию используются InputName и OutputName.
mw2.InputName mw2.OutputName
ans =
1×1 cell array
{'u1'}
ans =
1×1 cell array
{'y1'}
Регрессоры модели IDNLARX можно просмотреть с помощью команды GETREG. Регрессоры выполнены из mw2.InputName, mw2.OutputName и временные задержки.
getreg(mw2)
ans =
6×1 cell array
{'y1(t-1)'}
{'y1(t-2)'}
{'y1(t-3)'}
{'y1(t-4)'}
{'y1(t-5)'}
{'u1(t-3)'}
Обратите внимание, что матрица порядка эквивалентна линейному набору регрессоров, созданному с помощью объекта спецификации, как в:
R = linearRegressor([z1.OutputName; z1.InputName],{1:5, 3});
mw2_a = nlarx(z1, R, wavenet(8));
getreg(mw2_a)
ans =
6×1 cell array
{'y1(t-1)'}
{'y1(t-2)'}
{'y1(t-3)'}
{'y1(t-4)'}
{'y1(t-5)'}
{'u1(t-3)'}
mw2_a является по существу той же моделью, что и mw2. Использование объекта спецификации обеспечивает большую гибкость; используйте его, если вы не хотите использовать параллельные задержки или если минимальное отставание в выходных переменных больше единицы.
Выходные данные модели IDNLARX являются статической функцией ее регрессоров. Обычно эта функция отображения имеет линейный блок, нелинейный блок, плюс член смещения (выходное смещение). Выходной сигнал модели представляет собой сумму выходных сигналов этих блоков и смещения. Чтобы уменьшить сложность модели, только подмножество регрессоров может быть выбрано для входа в нелинейный блок и линейный блок. Назначение регрессоров линейным и нелинейным компонентам функции отображения выполняется с помощью свойства модели RegingUsage.
RegUseTable = mw2.RegressorUsage
RegUseTable =
6×2 table
y1:LinearFcn y1:NonlinearFcn
____________ _______________
y1(t-1) true true
y1(t-2) true true
y1(t-3) true true
y1(t-4) true true
y1(t-5) true true
u1(t-3) true true
RegUseTable - это таблица, показывающая, какие регрессоры используются в качестве входных данных для линейных и нелинейных компонентов вейвлет-сети. Каждая строка соответствует одному регрессору. Предположим, что мы хотим использовать только те регрессоры, которые являются функциями входной переменной «u1» в нелинейном компоненте. Эта конфигурация может быть достигнута следующим образом.
RegNames = RegUseTable.Properties.RowNames; I = contains(RegNames,'u1'); RegUseTable.("y1:NonlinearFcn")(~I) = false;
В качестве примера рассмотрим модель, которая использует полиномиальные регрессоры порядка 2 в дополнение к линейным, используемым mw2. Сначала создайте набор спецификаций полиномиального регрессора.
P = polynomialRegressor({'y1','u1'},{1:2, 0:3},2)
P =
<strong>Order 2 regressors in variables y1, u1</strong>
Order: 2
Variables: {'y1' 'u1'}
Lags: {[1 2] [0 1 2 3]}
UseAbsolute: [0 0]
AllowVariableMix: 0
AllowLagMix: 0
TimeVariable: 't'
Теперь добавьте спецификацию P к свойству модели Regressors.
mw3 = mw2; % start with the mw2 structure for new estimation mw3.Regressors = [mw3.Regressors; P]; % add the polynomial set to the model structure RegNames = mw3.RegressorUsage.Properties.RowNames; Ind = contains(RegNames,'u1'); mw3.RegressorUsage.("y1:NonlinearFcn")(~Ind) = false;
Наконец, обновите параметры модели в соответствии с данными путем минимизации ошибок прогнозирования на 1 шаг вперед.
mw3 = nlarx(z1, mw3)
mw3 = <strong>Nonlinear ARX model with 1 output and 1 input</strong> Inputs: u1 Outputs: y1 Regressors: 1. Linear regressors in variables y1, u1 2. Order 2 regressors in variables y1, u1 Output function: Wavelet Network with 5 units Sample time: 0.2 seconds Status: Estimated using NLARX on time domain data "Two tank system". Fit to estimation data: 96.72% (prediction focus) FPE: 3.583e-05, MSE: 3.438e-05
Различные модели, как линейные, так и нелинейные, можно сравнивать вместе в одной команде СРАВНИТЬ.
mlin = arx(z1,[5 1 3]); % linear ARX model % compare(z1,mlin,mw2,mw3)
Проверка модели может быть выполнена с помощью команды COMPARE в наборах данных проверки. Во-первых, проверьте модели, проверив, насколько хорошо их смоделированные ответы соответствуют измеренному выходному сигналу в z2.
compare(z2,mlin,mw2,mw3)
Аналогично, проверьте оцененные модели по набору данных z3.
compare(z3,mlin,mw2,mw3)
Анализ остатков модели также может быть выполнен для проверки качества модели. Мы ожидаем, что остатки будут белыми (некоррелированными с самим собой, за исключением запаздывания 0) и некоррелированными с входным сигналом. Остатки - это ошибки прогнозирования на 1 шаг вперед. Остатки можно визуализировать с помощью границ незначительных значений с помощью команды RESID. Например, в наборе данных оценки z2остатки для трех идентифицированных моделей показаны ниже.
resid(z2,mlin,mw2,mw3)
legend show
Все остатки не коррелируются с входным сигналом z2 (z2. InputData). Тем не менее, они демонстрируют некоторую автокорреляцию для моделей mlin и mw3. Модель mw2 лучше работает в обоих тестах (авто- и кросс-корреляция).
Функция PLOT может использоваться для просмотра нелинейных ответов различных моделей IDNLARX. График по существу показывает характеристики функции нелинейного отображения Fcn ().
plot(mw2,mw3) % plot nonlinearity response as a function of selected regressors
Вместо WAVENET в модели могут использоваться другие нелинейные функции. Попробуйте использовать функцию SIGMOIDNET, которая отображает регрессоры в выходные данные, используя сумму сигмоидных единиц (плюс линейный и смещенный член). Настройте sigmoidnet для использования 8 единиц сигмоида (расширенных и переведенных на различные количества).
ms1 = nlarx(z1,[5 1 3], sigmoidnet(8)); compare(z1,ms1)
Теперь проанализируйте новую модель в наборе данных z2.
compare(z2,ms1);
и на наборе данных z3.
compare(z3,ms1);
CUSTOMNET аналогичен SIGMOIDNET, но вместо функции sigmoid unit пользователи могут определить свою собственную функцию unit в файлах MATLAB. В этом примере используется функция gaussunit.
type gaussunit.m
function [f, g, a] = gaussunit(x) %GAUSSUNIT customnet unit function example % %[f, g, a] = GAUSSUNIT(x) % % x: unit function variable % f: unit function value % g: df/dx % a: unit active range (g(x) is significantly non zero in the interval [-a a]) % % The unit function must be "vectorized": for a vector or matrix x, the output % arguments f and g must have the same size as x, computed element-by-element. % Copyright 2005-2006 The MathWorks, Inc. % Author(s): Qinghua Zhang f = exp(-x.*x); if nargout>1 g = - 2*x .* f; a = 0.2; end % FILE END
Для модели на основе customnet используйте только третий, пятый и шестой регрессоры в качестве входных данных для нелинейного компонента функции customnet.
mc1_ini = idnlarx('y1','u1', [5 1 3], customnet(@gaussunit)); Use = false(6,1); Use([3 5 6]) = true; mc1_ini.RegressorUsage{:,2} = Use; mc1 = nlarx(z1,mc1_ini);
Функция секционирования дерева - это еще одна функция отображения, которую можно использовать в модели нелинейного ARX. Оценка модели с использованием функции treepartion:
mt1 = nlarx(z1,[5 1 3], treepartition);
Сравнение ответов моделей mc1 и mt1 к данным z1.
compare(z1, mc1, mt1)
При наличии Toolbox™ Deep Learning можно также использовать нейронную сеть для определения функции отображения Fcn () модели IDNLARX. Эта сеть должна быть прямой сетью (без повторяющихся компонентов).
Здесь мы создадим сеть с одним выходом с неизвестным количеством входов (обозначим по размеру входа нулем в сетевом объекте). Количество входов устанавливается равным количеству регрессоров во время процесса оценки автоматически.
if exist('feedforwardnet','file')==2 % run only if Deep Learning Toolbox is available ff = feedforwardnet([5 20]); % use a feed-forward network to map regressors to output ff.layers{2}.transferFcn = 'radbas'; ff.trainParam.epochs = 50; % Create a neuralnet mapping function that wraps the feed-forward network OutputFcn = neuralnet(ff); mn1 = nlarx(z1,[5 1 3], OutputFcn); % estimate network parameters compare(z1, mn1) % compare fit to estimation data end
Модель Хаммерштейна-Винера состоит из до 3 блоков: блоку линейной передаточной функции предшествует нелинейный статический блок и/или следует другой нелинейный статический блок. Она называется моделью Винера, если первый нелинейный статический блок отсутствует, и моделью Хаммерштейна, если второй нелинейный статический блок отсутствует.
Модели IDNLHW обычно оцениваются с помощью синтаксиса:
M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).
где Orders = [nb bf nk] определяет порядки и задержку линейной передаточной функции. InputNonlinearity и OutputNonlinearity задают нелинейные функции для двух нелинейных блоков на входном и выходном концах линейного блока. M является моделью @ idnlhw. Линейная модель вывода-ошибки (OE) соответствует случаю тривиальных (единичного усиления) нелинейностей, то есть если входная и выходная нелинейные функции модели IDNLHW являются как единичными коэффициентами усиления, то структура модели эквивалентна линейной модели OE.
Нелинейная функция PWLINEAR (кусочно-линейная) полезна в общих моделях IDNLHW.
Обратите внимание, что в Orders = [nb nf nk] nf указывает количество полюсов линейной передаточной функции, несколько связанных с порядком «na» модели IDNLARX. «nb» относится к числу, если нули.
mhw1 = nlhw(z1, [1 5 3], pwlinear, pwlinear);
Результат можно оценить с помощью команды СРАВНИТЬ, как и ранее.
compare(z1,mhw1)
Свойства модели можно визуализировать с помощью команды ПЕЧАТЬ.
Щелкните по блок-схеме, чтобы выбрать вид входной нелинейности (UNL), линейного блока или выходной нелинейности (YNL).
Для линейного блока можно выбрать тип графиков в пределах Step response, Bode plot, Impulse response и карты полюсов-нулей.
plot(mhw1)
Точки разрыва двух кусочно-линейных функций могут быть изучены. Это матрица из двух строк, первая строка соответствует входу, а вторая строка - выходу функции PWLINEAR.
mhw1.InputNonlinearity.BreakPoints
ans =
Columns 1 through 7
2.6179 3.3895 4.1611 4.9327 5.6297 6.3573 7.0850
-2.2681 -1.6345 0.1033 2.2066 3.7696 4.9984 5.8644
Columns 8 through 10
7.8127 8.2741 9.2257
6.4075 6.5725 9.1867
Функции SATURATION и DEADZONE - это нелинейные функции, вдохновленные физикой. Они могут использоваться для описания сценария насыщения привода/датчика и эффектов сухого трения.
mhw2 = nlhw(z1, [1 5 3], deadzone, saturation); compare(z1,mhw2)
Отсутствие нелинейности указывается объектом UNITGAIN или пустым значением «[]». Единичный коэффициент усиления - это только подача входного сигнала на выход без какой-либо модификации.
mhw3 = nlhw(z1, [1 5 3], 'deadzone', unitgain); % no output nonlinearity mhw4 = nlhw(z1, [1 5 3], [],'saturation'); % no input nonlinearity
Предельные значения DEADZONE и SATURATION могут быть изучены после оценки.
mhw2.InputNonlinearity.ZeroInterval mhw2.OutputNonlinearity.LinearInterval
ans =
-0.9974 4.7573
ans =
0.2603 0.5846
При оценке модели можно задать начальное предположение ZeroInterval для DEADZONE или LinearInterval для SATURATION.
mhw5 = nlhw(z1, [1 5 3], deadzone([0.1 0.2]), saturation([-1 1]));
Параметры алгоритма оценки можно задать с помощью команды NLHWOPTIONS.
opt = nlhwOptions();
opt.SearchMethod = 'gna';
opt.SearchOptions.MaxIterations = 7;
mhw6 = nlhw(z1, [1 5 3], deadzone, saturation, opt);
Уже оцененная модель может быть уточнена с помощью большего количества итераций оценки.
Оценка на основе оценочных данных z1, mhw7 должны иметь лучшую форму, чем mhw6.
opt.SearchOptions.MaxIterations = 30; mhw7 = nlhw(z1, mhw6, opt); compare(z1, mhw6, mhw7)
Нелинейные функции PWLINEAR, SATURATION и DEADZONE в основном предназначены для использования в моделях IDNLHW. Они могут оценивать только единичные переменные нелинейности. Другие более общие нелинейные функции, SIGMOIDNET, CUSTOMNET и WAVENET, также могут использоваться в моделях IDNLHW. Однако недифференцированные функции TREEPARTITION и NEURALNET не могут быть использованы, поскольку оценка моделей IDNLHW требует дифференцированных нелинейных функций.