Измеренный набор данных

Данные, сохраненные в файле 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 (IDNLARX) - регрессоры комплектации

Давайте сначала попробуем нелинейные модели 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.

Нелинейная модель ARX (IDNLARX) - предварительная оценка с использованием Wavenet

Для начала выберем порядки 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)

Нелинейная модель ARX - попытка высших порядков

Давайте проверим, сможем ли мы лучше справляться с более высокими порядками. Обратите внимание, что при идентификации моделей с использованием расширений базиса функций для выражения нелинейности, количество параметров модели может превысить количество выборок данных. В таких случаях некоторые метрики оценки, такие как Отклонение шума и Окончательная Ошибка Предсказания (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

Нелинейная модель ARX - настраивающее число модулей нелинейных функций

Количество модулей измерения в функции 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

Нелинейная модель ARX - Попытка других нелинейных функций отображения

Вместо функции 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 в качестве дополнительного входного параметра.

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

Возможно использовать различные нелинейные функции на разных выходных каналах в одной модели. Предположим, что мы хотим использовать функцию древовидного разбиения, чтобы описать первый выход и использовать вейвлет для второго выхода. Оценка модели показана ниже. Третий входной параметр (нелинейная функция отображения) теперь является массивом двух разных функций.

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, могут использоваться на предполагаемых моделях так же, как и в случае линейных моделей.

Модель Гаммерштейна-Винера (IDNLHW) - Предварительная оценка

Модель Гаммерштейна-Винера состоит из до 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

Начальное предположение для линейного интервала НАСЫЩЕНИЯ и нулевого интервала 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