Этот пример показывает идентификацию нелинейных моделей черного ящика входа single-single выхода (SISO) с использованием измеренных входно-выходных данных. Пример использует измеренные данные из системы с двумя емкостями для исследования различных структур модели и вариантов идентификации.
В этом примере переменные данных, хранящиеся в twotankdata.mat будет использоваться файл. Он содержит 3000 выборок входно-выходных данных системы с двумя емкостями, сгенерированных со шаг расчета 0,2 секунды. Вход u (t) является напряжением [V], приложенным к насосу, который генерирует приток к верхнему баку. Довольно небольшое отверстие в нижней части этого верхнего бака приводит к выходу, который поступает в нижний бак, и выход y (t) системы двух баков затем является уровнем жидкости [ m ] нижнего бака. Создадим объект IDDATA z для инкапсуляции загруженных данных:
Этот набор данных также используется в нелинейном сером прямоугольном примере «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}. Более комплексные, произвольные формулы в задержанных переменных также возможны.
В System Identification Toolbox™ наборы регрессоров не должны создаваться явно; большой набор может быть сгенерирована неявно путем простого определения их формы (такой как полиномиальный порядок), способствующих переменных и их лагов. Объекты спецификации 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 = polynomialRegressor ({'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 производит свой выход с помощью 2-ступенчатого преобразования: 1. Сопоставьте обучающие данные (входно-выходные сигналы) с набором регрессоров. Численно набор регрессоров является двойной матрицей R с одним столбцом данных для каждого регрессора. 2. Сопоставьте регрессоры с одним выходом, y = Fcn ®, где Fcn () является скалярной нелинейной функцией.
Чтобы оценить модель IDNLARX, после выбора порядков моделей, нам нужно выбрать тип нелинейной функции отображения Fcn () для использования. Давайте сначала попробуем функцию Wavelet Network (wavenet).
mw1 = nlarx(z1,[5 1 3], wavenet);
Сеть Вейвлета использует комбинацию вейвлета модулей, чтобы сопоставить регрессоры с выходными данными модели. Модель mw1 является объектом @ idnalrx. Он хранит нелинейную функцию отображения (wavenet here) в своей 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 является статической функцией ее регрессоров. Обычно эта функция отображения имеет линейный блок, нелинейный блок плюс термин смещения (выходное смещение). Выходные данные модели являются суммой выходов этих блоков и смещения. В порядок уменьшения сложности модели может быть выбран только подмножество регрессоров, чтобы войти в нелинейный блок и линейный блок. Назначение регрессоров линейным и нелинейным компонентам функции отображения выполняется с помощью свойства RegressorUsage модели.
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
Различные модели, как линейные, так и нелинейные, можно сравнить вместе в одной команде COMPARE.
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 лучше выполняет оба теста (авто- и перекрестная корреляция).
График функций может использоваться, чтобы просмотреть нелинейные отклики различных моделей IDNLARX. График по существу показывает характеристики нелинейной функции отображения Fcn ().
plot(mw2,mw3) % plot nonlinearity response as a function of selected regressors
Вместо WAVENET в модели могут использоваться другие нелинейные функции. Попробуйте функцию 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 пользователи могут задать собственную единичную функцию в файлах 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. Оцените модель, используя функцию treepartition:
mt1 = nlarx(z1,[5 1 3], treepartition);
Сравните отклики моделей mc1 и mt1 к данным z1.
compare(z1, mc1, mt1)
Если у вас есть доступная Deep Learning Toolbox™, можно также использовать нейронную сеть, чтобы задать функцию отображения 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).
где порядки = [nb bf nk] определяет порядки и задержку линейной передаточной функции. InputNonlinearity и OutputNonlinearity задают нелинейные функции для двух нелинейных блоков на входе и выходном конце линейного блока. M является моделью @ idnlhw. Линейная модель Output-Error (OE) соответствует случаю тривиальной (единичного усиления) нелинейности, то есть, если входные и выходные нелинейные функции модели IDNLHW являются единичными усилениями, структура модели эквивалентна линейной модели OE.
Нелинейная функция PWLINEAR (кусочно-линейная) полезна в общих моделях IDNLHW.
Заметьте, что в порядки = [nb nf nk] nf задает количество полюсов линейной передаточной функции, несколько связанных с порядком «na» модели IDNLARX. «nb» связан с числом, если нули.
mhw1 = nlhw(z1, [1 5 3], pwlinear, pwlinear);
Результат может быть оценен как раньше с помощью COMPARE.
compare(z1,mhw1)
Свойства модели могут быть визуализированы командой plot.
Щелкните на блок-схеме, чтобы выбрать представление входа нелинейности (UNL), линейного блока или выхода нелинейности (YNL).
Для линейного блока можно выбрать тип графиков в пределах переходная характеристика, диаграмма Боде, импульсная характеристика и карты полюса-нули.
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
Функции НАСЫЩЕНИЕ и 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 для НАСЫЩЕНИЕ может быть задано при оценке модели.
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 требует дифференцируемых нелинейных функций.