Этот пример показывает, как использовать полиномиальные и пользовательские регрессоры в нелинейных моделях ARX (IDNLARX).
В модели IDNLARX каждый выход является функцией регрессоров, которые являются преобразованиями прошлых входов и прошлых выходов. Типичные регрессоры являются просто задержанными входными или выходными переменными, представленными объектом спецификации linearRegressor или полиномами задержанных переменных, представленных объектом polynomialRegressor. Однако, кроме возможности навязывать значение абсута (например, abs (y (t-1))), вы не можете создать сложные математические формулы, используя эти типы регрессоров. Например, если ваша выходная функция требует регрессора формы sin (u (t-3)). * exp (-abs (u (t))), вам нужен способ ввести свою пользовательскую формулу в выражение для регрессора. Это обеспечивается объектами customRegressor. Как следует из имени, объект customRegressor используется для включения произвольных, пользовательских формул в качестве регрессоров для вашей нелинейной модели ARX.
Рассмотрим пример электрической системы с напряжением V и током I в качестве входов. Если известно, что электроэнергия является важным количеством системы, то имеет смысл сформировать пользовательский регрессор V * I. Может быть более эффективным использовать соответственно определенные пользовательские регрессоры, чем использовать только линейные и полиномиальные регрессоры.
Файл 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')
Матрица порядка [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. В вышеописанном строении переменная 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 sys1ans = 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 будет использоваться для задания регрессоров вместо матрицы порядка в следующих сценариях:
вы хотите использовать произвольные (несмежные) лаги в переменных, таких как набор
вы хотите, чтобы минимальная задержка в выход переменных (переменных ) (ах) отличалась от 1, для примера, набора
вы хотите использовать абсолютные значения переменных, таких как в наборе где используется только абсолютное значение переменной '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
Часто требуются регрессоры, которые являются полиномами переменных задержки ввода-вывода. Они могут быть добавлены к модели с помощью объекта polynomialRegressor. Предположим, вы хотите добавить , в качестве регрессоров к модели. Сначала вы создаете объект спецификации с этими настройками, а затем добавляете этот объект к модели.
% 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
Хотя наиболее часто используются линейный и полиномиальный регрессоры, иногда нужно использовать другую формулу, которая не описывается полиномом. Примером являются тригнометрические функции, такие как Другой пример - функция насыщения . В этих ситуациях можно использовать массив объектов customRegressor, которые инкапсулируют определенное математическое выражение.
В следующем примере регрессор создается как косинусная функция переменной с именем 'u1' и задерживает 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.
Можно создать массив регрессоров, все из которых имеют одну и ту же базовую формулу, но используют различные значения задержки. Например, предположим, что формула , где '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 регрессоров, сгенерированных при помощи формулы для всех комбинаций и значения в области значений 1:10.
Предположим, что проба/ошибка или физическое понимание предполагают, что нам нужно использовать набор регрессоров в модели. У нас есть смесь линейных, полиномиальных и тригонометрических формул. Мы действуем следующим образом:
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 как лучшая модель. Обратите внимание, что регрессоры, используемые в этом примере, были выбраны произвольно, главным образом, чтобы показать другой способ создания регрессоров и оценки моделей. Лучшие результаты можно получить, выбирая регрессоры более разумно.
Файл 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)