Введение

В модели IDNLARX каждый выход является функцией регрессоров, которые являются преобразованиями прошлых входов и прошлых выходов. Типичные регрессоры являются просто задержанными входными или выходными переменными, представленными объектом спецификации linearRegressor или полиномами задержанных переменных, представленных объектом polynomialRegressor. Однако, кроме возможности навязывать значение абсута (например, abs (y (t-1))), вы не можете создать сложные математические формулы, используя эти типы регрессоров. Например, если ваша выходная функция требует регрессора формы sin (u (t-3)). * exp (-abs (u (t))), вам нужен способ ввести свою пользовательскую формулу в выражение для регрессора. Это обеспечивается объектами customRegressor. Как следует из имени, объект customRegressor используется для включения произвольных, пользовательских формул в качестве регрессоров для вашей нелинейной модели ARX.

Рассмотрим пример электрической системы с напряжением V и током I в качестве входов. Если известно, что электроэнергия является важным количеством системы, то имеет смысл сформировать пользовательский регрессор V * I. Может быть более эффективным использовать соответственно определенные пользовательские регрессоры, чем использовать только линейные и полиномиальные регрессоры.

Пример SISO: Моделирование Engine внутреннего сгорания

Файл icEngine.mat содержит один набор данных с 1500 входно-выходными выборками, собранными со частотой дискретизации 0,04 секунды. Вход u (t) является напряжением [V], управляющим перепускным воздушным клапаном на холостом ходу (BPAV), и выход y (t) является скоростью вращения двигателя [RPM/100]. Данные загружаются и разделяются на набор данных ze для оценки модели и другой набор данных zv для валидации модели.

load icEngine
z = iddata(y,u,0.04);
ze = z(1:1000); 
zv = z(1001:1500);
plot(ze,zv)
legend('Estimation data','Validation data')

Линейные регрессоры как порядки ARX

Матрица порядка [na nb nk], которая также используется в линейной модели ARX, помогает легко определять линейные регрессоры, которые являются просто входными/выходными переменными, задержанными на определенное количество выборок. Выбор порядков модели требует проб и ошибок. В этом примере давайте будем использовать [na nb nk] = [4 2 10], соответствующие линейным регрессорам y (t-1), y (t-2), y (t-3), y (t-4), u (t-10), u (t-11). Выберите линейную выходную функцию, чтобы выход модели был всего лишь взвешенной суммой из его шести регрессоров плюс смещение.

sys0 = nlarx(ze, [4 2 10], linear);

Имя входа, имя выхода и список регрессоров этой модели отображаются ниже. Заметьте, что используются имена по умолчанию 'u1', 'y1'.

sys0.InputName

ans = 1×1 cell array
    {'u1'}

sys0.OutputName

ans = 1×1 cell array
    {'y1'}

disp(getreg(sys0))

    {'y1(t-1)' }
    {'y1(t-2)' }
    {'y1(t-3)' }
    {'y1(t-4)' }
    {'u1(t-10)'}
    {'u1(t-11)'}

Это линейные регрессоры. Они хранятся в свойстве модели Regressors как объект linearRegressor.

sys0.Regressors

ans = 
Linear regressors in variables y1, u1
       Variables: {'y1'  'u1'}
            Lags: {[1 2 3 4]  [10 11]}
     UseAbsolute: [0 0]
    TimeVariable: 't'

  Regressors described by this set

Свойство «Регрессоры» неявно хранит информацию о регрессоре, используя объект спецификации регрессии. Можно представить матрицу порядка [4 2 10] как ярлык для определения линейных регрессоров, где - лаги смежны, и минимальная задержка в выходной переменной фиксирована на 1.

Линейные регрессоры с использованием объекта спецификации linearRegressor

Более гибким способом, который позволяет выбирать переменные с произвольными лагами, является использование объекта linearRegressor. В вышеописанном строении переменная y1 имеет лаги [1 2 3 4], в то время как переменная u1 имеет лаги [10 11]. Используя объект linearRegressor, то же строение может быть достигнута следующим образом:

LinReg = linearRegressor({'y1','u1'}, {1:4, [10, 11]});
sys1 = nlarx(ze, LinReg, linear)

sys1 = 
Nonlinear ARX model with 1 output and 1 input
  Inputs: u1
  Outputs: y1

Regressors:
  Linear regressors in variables y1, u1
  List of all regressors

Model output is linear in regressors.
Sample time: 0.04 seconds

Status:                                          
Estimated using NLARX on time domain data "ze".  
Fit to estimation data: 94.35% (prediction focus)
FPE: 0.001877, MSE: 0.00183

Сравните синтаксисы оценок моделей sys1 и sys0; при создании sys1 матрица порядка была заменена на объект спецификации линейного регрессора. Получившиеся модели sys0 и sys1 идентичны.

[getpvec(sys0), getpvec(sys1)] % estimated parameter vectors of sys0 and sys1

ans = 7×2

    0.7528    0.7528
    0.0527    0.0527
   -0.0621   -0.0621
   -0.0425   -0.0425
   -0.0165   -0.0165
   -0.0289   -0.0289
    0.2723    0.2723

Объект linearRegressor будет использоваться для задания регрессоров вместо матрицы порядка в следующих сценариях:

LinRegWithAbs = linearRegressor({'y1','u1'},{[1 10], [0 2]},[true, false])

LinRegWithAbs = 
Linear regressors in variables y1, u1
       Variables: {'y1'  'u1'}
            Lags: {[1 10]  [0 2]}
     UseAbsolute: [1 0]
    TimeVariable: 't'

  Regressors described by this set

Полиномиальные регрессоры

Часто требуются регрессоры, которые являются полиномами переменных задержки ввода-вывода. Они могут быть добавлены к модели с помощью объекта polynomialRegressor. Предположим, вы хотите добавить $$u1(t-10{)}^2$$, $$y1(t-1{)}^2$$ в качестве регрессоров к модели. Сначала вы создаете объект спецификации с этими настройками, а затем добавляете этот объект к модели.

% Create order 2 regressors in variables 'u1' and 'y1', with lags 10 and 1
% respectively
PolyReg = polynomialRegressor({'u1','y1'},{10, 1},2);

% Use LinReg and PolyReg in the model
sys2 = nlarx(ze, [LinReg; PolyReg], linear)

sys2 = 
Nonlinear ARX model with 1 output and 1 input
  Inputs: u1
  Outputs: y1

Regressors:
  1. Linear regressors in variables y1, u1
  2. Order 2 regressors in variables u1, y1
  List of all regressors

Model output is linear in regressors.
Sample time: 0.04 seconds

Status:                                          
Estimated using NLARX on time domain data "ze".  
Fit to estimation data: 94.47% (prediction focus)
FPE: 0.001804, MSE: 0.001752

Регрессоры на основе пользовательских формул

Хотя наиболее часто используются линейный и полиномиальный регрессоры, иногда нужно использовать другую формулу, которая не описывается полиномом. Примером являются тригнометрические функции, такие как $$sin(),cos(),\dots$$Другой пример - функция насыщения $$x=max(LowerBound,min(x,UpperBound)$$. В этих ситуациях можно использовать массив объектов customRegressor, которые инкапсулируют определенное математическое выражение.

В следующем примере регрессор создается как косинусная функция переменной с именем 'u1' и задерживает 10 выборок, другими словами: $$x=cos(u_1(t-10))$$. Логическое значение в последнем входном параметре указывает, векторизован ли пользовательский регрессор или нет. Векторизованные регрессоры быстрее в расчеты, но требуют ухода в функции, указанной на первом входном параметре.

x = customRegressor('u1', 10, @cos)

x = 
Custom regressor: cos(u1(t-10))
    VariablesToRegressorFcn: @cos
                  Variables: {'u1'}
                       Lags: {[10]}
                 Vectorized: 1
               TimeVariable: 't'

  Regressors described by this set

Заданная формула (@ cos здесь) сохранена в VariablesToRegressorFcn свойство объекта регрессора. По умолчанию вычисление функции принято векторизованным, то есть если вход является матрицей с N строками, то выход функции x.VariablesToRegressorFcn будет вектором-столбцом длины N. Векторизация помогает ускорить оценку регрессора в процессе оценки модели и симуляции. Однако у вас есть опция отключить его при желании путем установки значения Vectorized свойство FALSE.

Можно создать массив регрессоров, все из которых имеют одну и ту же базовую формулу, но используют различные значения задержки. Например, предположим, что формула $$x=u(t-a)*|y(t-b)|$$, где 'u' и 'y' являются двумя переменными, для которых доступны измеренные данные, и лаги (a, b) могут принимать несколько значений в области значений 1:10. Можно создать массив из этих регрессоров с одним вызовом customRegressor конструктор

C = customRegressor({'u','y'},{1:10, 1:10}, @(x,y)x.*abs(y))

C = 
Custom regressor: @(x,y)x.*abs(y)
    VariablesToRegressorFcn: @(x,y)x.*abs(y)
                  Variables: {'u'  'y'}
                       Lags: {[1 2 3 4 5 6 7 8 9 10]  [1 2 3 4 5 6 7 8 9 10]}
                 Vectorized: 1
               TimeVariable: 't'

  Regressors described by this set

C представляет набор из 100 регрессоров, сгенерированных при помощи формулы $$u(t-a).*abs(y(t-b))$$ для всех комбинаций $$a$$ и $$b$$значения в области значений 1:10.

Использование всех трех типов регрессоров в модели для динамики двигателя IC

Предположим, что проба/ошибка или физическое понимание предполагают, что нам нужно использовать набор регрессоров $$R=\{u1(t-10{)}^3,u1(t-11{)}^3,y1(t-1),sin(y1(t-4))\}$$в модели. У нас есть смесь линейных, полиномиальных и тригонометрических формул. Мы действуем следующим образом:

LinReg    = linearRegressor('y1',1);
PolyReg   = polynomialRegressor('u1',[10 11],3);
CustomReg = customRegressor('y1',4,@(x)sin(x));
% for now no nonlinearity; output is a linear function of regressors
sys3 = nlarx(ze,[LinReg; PolyReg; CustomReg], [])

sys3 = 
Nonlinear ARX model with 1 output and 1 input
  Inputs: u1
  Outputs: y1

Regressors:
  1. Linear regressors in variables y1
  2. Order 3 regressors in variables u1
  3. Custom regressor: sin(y1(t-4))
  List of all regressors

Model output is linear in regressors.
Sample time: 0.04 seconds

Status:                                          
Estimated using NLARX on time domain data "ze".  
Fit to estimation data: 92.11% (prediction focus)
FPE: 0.003647, MSE: 0.00357

getreg(sys3)

ans = 4×1 cell
    {'y1(t-1)'     }
    {'u1(t-10)^3'  }
    {'u1(t-11)^3'  }
    {'sin(y1(t-4))'}

Мы можем расширить этот рабочий процесс, включив нелинейные функции отображения, такие как Sigmoid Network в модель, а также обозначить только подмножество набора регрессоров, который будет использоваться в качестве входов для его линейных и нелинейных компонентов (примечание: Сигмоидная сеть является суммой 3 компонентов - линейная функция, член смещения и нелинейная функция, которая является суммой сигмоидных модулей). В следующем примере мы используем рабочий процесс, основанный на модели шаблона, в котором мы серапически готовим шаблонную модель IDNLARX и набор опций оценки перед использованием их в команде NLARX для оценки параметра.

sys4 = idnlarx(ze.OutputName, ze.InputName, [4 2 10], sigmoidnet);
% generate 'u1(t-10)^2', 'y1(t-1)^2', 'u1(t-10)*y1(t-1)'
P = polynomialRegressor({'y1','u1'},{1, 10},2, false, true);
% generate cos(u1(t-10))
C1 = customRegressor('u1',10,@cos);
% generate sin(y1(t-1).*u1(t-10)+u1(t-11))
C2 = customRegressor({'y1','u1','u1'},{1, 10, 11},@(x,y,z)sin(x.*y+z));

% add the polynomial and custom regressors to the model sys2
sys4.Regressors = [sys4.Regressors; P; C1; C2];

% view the regressors and how they are used in the model
disp(sys4.RegressorUsage)

                                       y1:LinearFcn    y1:NonlinearFcn
                                       ____________    _______________

    y1(t-1)                               true              true      
    y1(t-2)                               true              true      
    y1(t-3)                               true              true      
    y1(t-4)                               true              true      
    u1(t-10)                              true              true      
    u1(t-11)                              true              true      
    y1(t-1)^2                             true              true      
    u1(t-10)^2                            true              true      
    cos(u1(t-10))                         true              true      
    sin(y1(t-1).*u1(t-10)+u1(t-11))       true              true      

% designate only the linear regressors to be used in the nonlinear
% component of the sigmoid network
Usage = sys4.RegressorUsage;
Usage{7:end,2} = false;
sys4.RegressorUsage = Usage;
disp(sys4.RegressorUsage)

                                       y1:LinearFcn    y1:NonlinearFcn
                                       ____________    _______________

    y1(t-1)                               true              true      
    y1(t-2)                               true              true      
    y1(t-3)                               true              true      
    y1(t-4)                               true              true      
    u1(t-10)                              true              true      
    u1(t-11)                              true              true      
    y1(t-1)^2                             true              false     
    u1(t-10)^2                            true              false     
    cos(u1(t-10))                         true              false     
    sin(y1(t-1).*u1(t-10)+u1(t-11))       true              false     

% Prepapre estimation options: use Levenberg-Marquardt solver with 30 maximum iterations.
% Turn progress display on and set estimation focus to 'simulation'
opt = nlarxOptions;
opt.Focus = 'simulation';
opt.Display = 'on';
opt.SearchMethod = 'lm';
opt.SearchOptions.MaxIterations = 30;

% estimate parameters of sys4 to fit data ze
sys4 = nlarx(ze, sys4, opt)

sys4 = 
Nonlinear ARX model with 1 output and 1 input
  Inputs: u1
  Outputs: y1

Regressors:
  1. Linear regressors in variables y1, u1
  2. Order 2 regressors in variables y1, u1
  3. Custom regressor: cos(u1(t-10))
  4. Custom regressor: sin(y1(t-1).*u1(t-10)+u1(t-11))
  List of all regressors

Output function: Sigmoid Network with 10 units

Sample time: 0.04 seconds

Status:                                          
Estimated using NLARX on time domain data "ze".  
Fit to estimation data: 87.81% (simulation focus)
FPE: 0.001491, MSE: 0.008517

Теперь мы можем подтвердить модели, сравнив их ответы с выходами набора данных validata zv.

compare(zv, sys1, sys2, sys3, sys4)

Сравнительные графики указывают sys2 как лучшая модель. Обратите внимание, что регрессоры, используемые в этом примере, были выбраны произвольно, главным образом, чтобы показать другой способ создания регрессоров и оценки моделей. Лучшие результаты можно получить, выбирая регрессоры более разумно.

Пример MIMO: моделирование моторизованной камеры

Файл 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');
plot(z)

Использование различных типов регрессоров в случае MIMO не сильно отличается от случая SISO. Все регрессоры используются в динамике для всех выходов системы. The RegressorUsage свойство может использоваться, чтобы назначить определенные регрессоры, которые будут использоваться для определенных выходов.

    nanbnk = [ones(2,2), 2*ones(2,6), ones(2,6)];
    sysMIMO = idnlarx(z.OutputName, z.InputName, nanbnk, linear);
    C = customRegressor({'u5','u6'},{1 1},@(x,y)x.*y);
    P1 = polynomialRegressor('u1',1,2);
    P2 = polynomialRegressor('y2',1,3);
    sysMIMO.Regressors = [sysMIMO.Regressors; C; P1; P2];
    getreg(sysMIMO)

ans = 17×1 cell
    {'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)'         }
    {'u1(t-1)^2'       }
    {'y2(t-1)^3'       }
    {'u5(t-1).*u6(t-1)'}

    sysMIMO = nlarx(z, sysMIMO)

sysMIMO = 
Nonlinear ARX model with 2 outputs and 6 inputs
  Inputs: u1, u2, u3, u4, u5, u6
  Outputs: y1, y2

Regressors:
  1. Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6
  2. Order 2 regressors in variables u1
  3. Order 3 regressors in variables y2
  4. Custom regressor: u5(t-1).*u6(t-1)
  List of all regressors

 All model outputs are linear in their regressors.
Sample time: 0.02 seconds

Status:                                                      
Estimated using NLARX on time domain data "Motorized Camera".
Fit to estimation data: [99.05;98.85]% (prediction focus)    
FPE: 0.2067, MSE: 0.7496

Сравните ответ модели с данными оценки.

    compare(z,sysMIMO)

Анализируйте невязки модели для их автокорреляции (тест белизны) и перекрестной корреляции с входным сигналом (корреляционный тест)

    fig = figure;
    fig.Position(3) = 2*fig.Position(3);
    resid(z,sysMIMO)