Набор данных ввода-вывода

В этом примере переменные данных, хранящиеся в 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). Мы постараемся использовать эти значения при оценке нелинейных моделей.

Нелинейные модели ARX (IDNLARX)

В модели 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

Нелинейная модель ARX с функцией SIGMOIDNET

Вместо WAVENET в модели могут использоваться другие нелинейные функции. Попробуйте функцию SIGMOIDNET, которая преобразует регрессоры в выходные с помощью суммы сигмоидных модулей (плюс линейный и смещенный термин). Сконфигурируйте сигмоиднет, чтобы использовать 8 модули сигмоида (расширенный и переведенный на различные количества).

ms1 = nlarx(z1,[5 1 3], sigmoidnet(8));
compare(z1,ms1)

Теперь оцените новую модель на наборе данных z2.

compare(z2,ms1);

и на наборе данных z3.

compare(z3,ms1);

Другие нелинейные функции отображения в модели IDNLARX

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

Если у вас есть доступная 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

Модели Гаммерштейна-Винера (IDNLHW)

Модель Гаммерштейна-Винера состоит из до 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 требует дифференцируемых нелинейных функций.