Создавание нелинейных моделей ARX с нелинейными и пользовательскими регрессорами

В этом примере показано, как использовать полиномиальные и пользовательские регрессоры в моделях Nonlinear ARX (IDNLARX).

Введение

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

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

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

Файл icEngine.mat содержит один набор данных с 1 500 выборками ввода - вывода, собранными в частота дискретизации 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')

Figure contains 2 axes objects. Axes object 1 with title y1 contains 2 objects of type line. These objects represent Estimation data, Validation data. Axes object 2 with title u1 contains 2 objects of type line. These objects represent 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], idLinear);

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

sys0.InputName
ans = 1x1 cell array
    {'u1'}

sys0.OutputName
ans = 1x1 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

Свойство "Regressors" хранит информацию на регрессоре неявно с помощью объекта спецификации регрессии. Можно ли думать о матрице порядка [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, idLinear)
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

Output function: None
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 задать регрессоры вместо матрицы порядка в следующих сценариях:

  • вы хотите использовать произвольные задержки (состоящие из нескольких несмежных участков), в переменных, таких как набор {y1(t-1),y1(t-10),u(t-3),u(t-11)}

  • вы хотите, чтобы минимальная задержка в выходной переменной (переменных) отличалась, чем 1, например, набор {y1(t-4),y1(t-5),...}

  • вы хотите использовать абсолютные значения переменных такой как в наборе {|y1(t-1)|,|y1(t-10)|,u(t),u(t-2)} где только абсолютное значение переменной 'y1' используется. Это может быть достигнуто путем выполнения

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

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

Часто регрессоры, которые являются полиномами задержанных переменных I/O, требуются. Они могут быть добавлены к модели с помощью объекта 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], idLinear)
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

Output function: None
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

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

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

В следующем примере регрессор создается как косинусная функция переменной, названной 'u1', и задержал 10 выборок, другими словами: x=cos(u1(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

Output function: None
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 = 4x1 cell
    {'y1(t-1)'     }
    {'u1(t-10)^3'  }
    {'u1(t-11)^3'  }
    {'sin(y1(t-4))'}

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

sys4 = idnlarx(ze.OutputName, ze.InputName, [4 2 10], idSigmoidNetwork);
% 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

Мы можем теперь подтвердить модели путем сравнения их ответов на выход подтверждать набора данных zv.

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

Figure contains an axes object. The axes object contains 5 objects of type line. These objects represent zv (y1), sys1: 73.84%, sys2: 80.24%, sys3: 42.74%, sys4: 37.97%.

Сравнить графики показывают 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)

Figure contains 8 axes objects. Axes object 1 with title y1 contains an object of type line. This object represents Motorized Camera. Axes object 2 with title y2 contains an object of type line. This object represents Motorized Camera. Axes object 3 with title u1 contains an object of type line. This object represents Motorized Camera. Axes object 4 with title u2 contains an object of type line. This object represents Motorized Camera. Axes object 5 with title u3 contains an object of type line. This object represents Motorized Camera. Axes object 6 with title u4 contains an object of type line. This object represents Motorized Camera. Axes object 7 with title u5 contains an object of type line. This object represents Motorized Camera. Axes object 8 with title u6 contains an object of type line. This object represents Motorized Camera.

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

    nanbnk = [ones(2,2), 2*ones(2,6), ones(2,6)];
    sysMIMO = idnlarx(z.OutputName, z.InputName, nanbnk, idLinear);
    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 = 17x1 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

Output functions:
  Output 1:  None
  Output 2:  None

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)

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Motorized Camera (y1), sysMIMO: 97.11%. Axes object 2 contains 2 objects of type line. These objects represent Motorized Camera (y2), sysMIMO: 98.38%.

    

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

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

Figure contains 14 axes objects. Axes object 1 with title AutoCorr contains 2 objects of type line. This object represents sysMIMO. Axes object 2 contains 2 objects of type line. This object represents sysMIMO. Axes object 3 with title XCorr (u1) contains 2 objects of type line. This object represents sysMIMO. Axes object 4 contains 2 objects of type line. This object represents sysMIMO. Axes object 5 with title XCorr (u2) contains 2 objects of type line. This object represents sysMIMO. Axes object 6 contains 2 objects of type line. This object represents sysMIMO. Axes object 7 with title XCorr (u3) contains 2 objects of type line. This object represents sysMIMO. Axes object 8 contains 2 objects of type line. This object represents sysMIMO. Axes object 9 with title XCorr (u4) contains 2 objects of type line. This object represents sysMIMO. Axes object 10 contains 2 objects of type line. This object represents sysMIMO. Axes object 11 with title XCorr (u5) contains 2 objects of type line. This object represents sysMIMO. Axes object 12 contains 2 objects of type line. This object represents sysMIMO. Axes object 13 with title XCorr (u6) contains 2 objects of type line. This object represents sysMIMO. Axes object 14 contains 2 objects of type line. This object represents sysMIMO.