Создание структурированных и пользовательских моделей Используя System Identification Toolbox™
Этот пример показывает, как оценить параметры в пользовательских образцовых структурах. Такие структуры заданы IDGREY (линейное пространство состояний) или IDNLGREY (нелинейное пространство состояний) модели. Мы займемся расследованиями, как присвоить структуру, зафиксировать параметры и создать зависимости среди них.
Данные об эксперименте
Мы исследуем данные, произведенные (моделируемым) двигателем постоянного тока. Мы сначала загружаем данные:
load dcmdata
who
Your variables are: text u y
Матричный y содержит эти два выходных параметров: y1 является угловым положением вала двигателя, и y2 является угловой скоростью. Существует 400 выборок данных, и шаг расчета составляет 0,1 секунды. Вход содержится в векторном u. Это - входное напряжение к двигателю.
z = iddata(y,u,0.1); % The IDDATA object z.InputName = 'Voltage'; z.OutputName = {'Angle';'AngVel'}; plot(z(:,1,:))
Фигура: данные об измерении: напряжение, чтобы повернуть
plot(z(:,2,:))
Фигура: данные об измерении: напряжение, чтобы повернуть
Образцовый выбор структуры
d/dt x = A x + B u + K e
y = C x + D u + eМы создадим модель двигателя постоянного тока. Динамика двигателя известна. Если мы выбираем x1 в качестве углового положения и x2 как угловая скорость, легко настроить модель в пространстве состояний следующих воздействий пренебрежения символа: (см. Пример 4.1 в Ljung (1999):
| 0 1 | | 0 |
d/dt x = | | x + | | u
| 0 -th1 | | th2 | | 1 0 |
y = | | x
| 0 1 |Параметр th1 является здесь инверсией, постоянной во времени из двигателя и th2, таков, что th2/th1 является статическим усилением от входа до угловой скорости. (См. Ljung (1987) для того, как th1 и th2 относятся к физическим параметрам двигателя). Мы оценим эти два параметра от наблюдаемых данных. Образцовая структура (параметризованное пространство состояний) описанный выше может быть представлена в MATLAB® с помощью моделей IDSS и IDGREY. Эти модели позволяют вам выполнить оценку параметров с помощью экспериментальных данных.
Спецификация номинальной (первоначальной) модели
Если мы предполагаем, что th1=1 и th2 = 0.28 получаем номинальную или первоначальную модель
A = [0 1; 0 -1]; % initial guess for A(2,2) is -1 B = [0; 0.28]; % initial guess for B(2) is 0.28 C = eye(2); D = zeros(2,1);
и мы группируем это в объект модели IDSS:
ms = idss(A,B,C,D);
Модель характеризуется ее матрицами, их значениями, какие элементы свободны (быть оцененными) и верхние и нижние пределы тех:
ms.Structure.a
ans =
Name: 'A'
Value: [2x2 double]
Minimum: [2x2 double]
Maximum: [2x2 double]
Free: [2x2 logical]
Scale: [2x2 double]
Info: [2x2 struct]
1x1 param.Continuous
ms.Structure.a.Value ms.Structure.a.Free
ans =
0 1
0 -1
ans =
2x2 logical array
1 1
1 1
Спецификация свободных (независимых) параметров Используя модели IDSS
Таким образом, мы должны теперь отметить это, это только (2,2) и B (2,1), которые являются свободными параметрами, которые будут оценены.
ms.Structure.a.Free = [0 0; 0 1]; ms.Structure.b.Free = [0; 1]; ms.Structure.c.Free = 0; % scalar expansion used ms.Structure.d.Free = 0; ms.Ts = 0; % This defines the model to be continuous
Первоначальная модель
ms % Initial model
ms =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2
x1 0 1
x2 0 -1
B =
u1
x1 0
x2 0.28
C =
x1 x2
y1 1 0
y2 0 1
D =
u1
y1 0
y2 0
K =
y1 y2
x1 0 0
x2 0 0
Parameterization:
STRUCTURED form (some fixed coefficients in A, B, C).
Feedthrough: none
Disturbance component: none
Number of free coefficients: 2
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
Оценка свободных параметров модели IDSS
Ошибка прогноза (наибольшее правдоподобие) оценка параметров теперь вычисляется:
dcmodel = ssest(z,ms,ssestOptions('Display','on')); dcmodel
dcmodel =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2
x1 0 1
x2 0 -4.013
B =
Voltage
x1 0
x2 1.002
C =
x1 x2
Angle 1 0
AngVel 0 1
D =
Voltage
Angle 0
AngVel 0
K =
Angle AngVel
x1 0 0
x2 0 0
Parameterization:
STRUCTURED form (some fixed coefficients in A, B, C).
Feedthrough: none
Disturbance component: none
Number of free coefficients: 2
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "z".
Fit to estimation data: [98.35;84.42]%
FPE: 0.001071, MSE: 0.1192
Ориентировочные стоимости параметров вполне близко к используемым, когда данные моделировались (-4 и 1). Чтобы оценить качество модели, мы можем моделировать модель с фактическим входом и сравнить его с фактическим выводом.
compare(z,dcmodel);
Мы можем теперь, например, построить нули и полюса и их области неуверенности. Мы будем чертить области, соответствующие 3 стандартным отклонениям, поскольку модель довольно точна. Обратите внимание на то, что полюс в начале координат является абсолютно бесспорным, поскольку это - часть образцовой структуры; интегратор от угловой скорости до положения.
clf showConfidence(iopzplot(dcmodel),3)
Теперь, мы можем сделать различные модификации. С 1,2 элементами из A-матрицы (зафиксированный к 1) говорит нам, что x2 является производной x1. Предположим, что датчики не калибруются, так, чтобы мог быть неизвестный коэффициент пропорциональности. Чтобы включать оценку такой константы, мы только "выпускаем" A(1,2) и переоценку:
dcmodel2 = dcmodel; dcmodel2.Structure.a.Free(1,2) = 1; dcmodel2 = pem(z,dcmodel2,ssestOptions('Display','on'));
Получившаяся модель
dcmodel2
dcmodel2 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2
x1 0 0.9975
x2 0 -4.011
B =
Voltage
x1 0
x2 1.004
C =
x1 x2
Angle 1 0
AngVel 0 1
D =
Voltage
Angle 0
AngVel 0
K =
Angle AngVel
x1 0 0
x2 0 0
Parameterization:
STRUCTURED form (some fixed coefficients in A, B, C).
Feedthrough: none
Disturbance component: none
Number of free coefficients: 3
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PEM on time domain data "z".
Fit to estimation data: [98.35;84.42]%
FPE: 0.001077, MSE: 0.1192
Мы находим, что предполагаемый A(1,2) близко к 1. Чтобы сравнить две модели, мы используем команду compare:
compare(z,dcmodel,dcmodel2)
Спецификация моделей с двойными параметрами Используя объекты IDGREY
Предположим, что мы точно знаем статическое усиление двигателя постоянного тока (от входного напряжения до угловой скорости, например, из предыдущего эксперимента переходного процесса. Если статическим усилением является G, и временная константа двигателя является t, то модель в пространстве состояний становится
|0 1| | 0 |
d/dt x = | |x + | | u
|0 -1/t| | G/t | |1 0|
y = | | x
|0 1|С известным G существует зависимость между записями в различных матрицах. В порядке описать, что, ранее используемый путь со "Свободными" параметрами не будет достаточен. Мы таким образом должны записать файл MATLAB, который производит A, B, C, и D, и опционально также K и матрицы X0 как выходные параметры, для каждого данного вектора параметра как вход. Это также берет вспомогательные аргументы в качестве входных параметров, так, чтобы пользователь мог изменить определенные вещи в образцовой структуре, не имея необходимость редактировать файл. В этом случае мы позволяем известному статическому усилению G вводиться как таковое аргумент. Файл, который был записан, имеет имя 'motorDynamics.m'.
type motorDynamics
function [A,B,C,D,K,X0] = motorDynamics(par,ts,aux) %MOTORDYNAMICS ODE file representing the dynamics of a motor. % % [A,B,C,D,K,X0] = motorDynamics(Tau,Ts,G) % returns the State Space matrices of the DC-motor with % time-constant Tau (Tau = par) and known static gain G. The sample % time is Ts. % % This file returns continuous-time representation if input argument Ts % is zero. If Ts>0, a discrete-time representation is returned. To make % the IDGREY model that uses this file aware of this flexibility, set the % value of Structure.FcnType property to 'cd'. This flexibility is useful % for conversion between continuous and discrete domains required for % estimation and simulation. % % See also IDGREY, IDDEMO7. % L. Ljung % Copyright 1986-2015 The MathWorks, Inc. t = par(1); G = aux(1); A = [0 1;0 -1/t]; B = [0;G/t]; C = eye(2); D = [0;0]; K = zeros(2); X0 = [0;0]; if ts>0 % Sample the model with sample time Ts s = expm([[A B]*ts; zeros(1,3)]); A = s(1:2,1:2); B = s(1:2,3); end
Мы теперь создаем объект модели IDGREY, соответствующий этой образцовой структуре: принятая временная константа будет
par_guess = 1;
Мы также даем значение 0.25 вспомогательной переменной G (усиление) и шаг расчета.
aux = 0.25; dcmm = idgrey('motorDynamics',par_guess,'cd',aux,0);
Временная константа теперь оценивается
dcmm = greyest(z,dcmm,greyestOptions('Display','on'));
Мы таким образом теперь оценили временную константу двигателя непосредственно. Его значение находится в хорошем соглашении с предыдущей оценкой.
dcmm
dcmm =
Continuous-time linear grey box model defined by "motorDynamics" function:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2
x1 0 1
x2 0 -4.006
B =
Voltage
x1 0
x2 1.001
C =
x1 x2
Angle 1 0
AngVel 0 1
D =
Voltage
Angle 0
AngVel 0
K =
Angle AngVel
x1 0 0
x2 0 0
Model parameters:
Par1 = 0.2496
Parameterization:
ODE Function: motorDynamics
(parametrizes both continuous- and discrete-time equations)
Disturbance component: parameterized by the ODE function
Initial state: parameterized by the ODE function
Number of free coefficients: 1
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using GREYEST on time domain data "z".
Fit to estimation data: [98.35;84.42]%
FPE: 0.00107, MSE: 0.1193
С этой моделью мы можем теперь продолжить тестировать различные аспекты как прежде. Синтаксис всех команд идентичен предыдущему случаю. Например, мы можем сравнить idgrey модель с другой моделью в пространстве состояний:
compare(z,dcmm,dcmodel)
Они ясно очень близки.
Оценка многомерных моделей ARX
Часть пространства состояний тулбокса также обрабатывает многомерный (несколько выходных параметров) модели ARX. Многомерной моделью ARX мы имеем в виду следующее:
A(q) y(t) = B(q) u(t) + e(t)
Здесь (q) ny | нью-йоркская матрица, записи которой являются полиномами в операторе задержки 1/q. k-l элемент обозначается:
где:
Это - таким образом полином в 1/q степени nakl.
Так же B (q) является ny | матрица ню, kj-элемент которой:
Существует таким образом задержка nkkj от входного номера j, чтобы вывести номер k. Наиболее распространенный способ создать тех состоял бы в том, чтобы использовать команду ARX. Порядки заданы как: nn = [na nb nk] с na, являющимся матрицей ny-by-ny, kj которой - запись является nakj; nb и nk заданы так же.
Давайте протестируем некоторые модели ARX на dc-данных. Сначала мы могли просто создать общую модель второго порядка:
dcarx1 = arx(z,'na',[2,2;2,2],'nb',[2;2],'nk',[1;1])
dcarx1 =
Discrete-time ARX model:
Model for output "Angle": A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + e_1(t)
A(z) = 1 - 0.5545 z^-1 - 0.4454 z^-2
A_2(z) = -0.03548 z^-1 - 0.06405 z^-2
B(z) = 0.004243 z^-1 + 0.006589 z^-2
Model for output "AngVel": A(z)y_2(t) = - A_i(z)y_i(t) + B(z)u(t) + e_2(t)
A(z) = 1 - 0.2005 z^-1 - 0.2924 z^-2
A_1(z) = 0.01849 z^-1 - 0.01937 z^-2
B(z) = 0.08642 z^-1 + 0.03877 z^-2
Sample time: 0.1 seconds
Parameterization:
Polynomial orders: na=[2 2;2 2] nb=[2;2] nk=[1;1]
Number of free coefficients: 12
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARX on time domain data "z".
Fit to estimation data: [97.87;83.44]% (prediction focus)
FPE: 0.002157, MSE: 0.1398
Результат, dcarx1, хранится как модель IDPOLY, и все предыдущие команды применяются. Мы могли, например, явным образом перечислить полиномы ARX:
dcarx1.a
ans =
2x2 cell array
{1x3 double} {1x3 double}
{1x3 double} {1x3 double}
как массивы ячеек, где, например, {1,2} элемент dcarx1.a является полиномом (1,2) описанный ранее, имея отношение y2 к y1.
Мы могли также протестировать структуру, где мы знаем, что y1 получен путем пропущения y2 через фильтр первого порядка. (Угол является интегралом угловой скорости). Мы могли затем также постулировать динамику первого порядка от входа, чтобы вывести номер 2:
na = [1 1; 0 1]; nb = [0 ; 1]; nk = [1 ; 1]; dcarx2 = arx(z,[na nb nk])
dcarx2 =
Discrete-time ARX model:
Model for output "Angle": A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + e_1(t)
A(z) = 1 - 0.9992 z^-1
A_2(z) = -0.09595 z^-1
B(z) = 0
Model for output "AngVel": A(z)y_2(t) = B(z)u(t) + e_2(t)
A(z) = 1 - 0.6254 z^-1
B(z) = 0.08973 z^-1
Sample time: 0.1 seconds
Parameterization:
Polynomial orders: na=[1 1;0 1] nb=[0;1] nk=[1;1]
Number of free coefficients: 4
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARX on time domain data "z".
Fit to estimation data: [97.52;81.46]% (prediction focus)
FPE: 0.003452, MSE: 0.177
Сравнить различные модели получило, мы используем
compare(z,dcmodel,dcmm,dcarx2)
Наконец, мы могли сравнить bodeplots, полученный из входа, чтобы вывести один для различных моделей при помощи bode: Сначала вывод:
dcmm2 = idss(dcmm); % convert to IDSS for subreferencing bode(dcmodel(1,1),'r',dcmm2(1,1),'b',dcarx2(1,1),'g')
Second output:
bode(dcmodel(2,1),'r',dcmm2(2,1),'b',dcarx2(2,1),'g')
Две первых модели находятся более или менее в точном соглашении. Модели ARX не так хороши, из-за смещения, вызванного цветным ошибочным шумом уравнения. (У нас был белый шум измерения в симуляциях).
Заключения
Оценка моделей с предварительно выобранными структурами может быть выполнена с помощью тулбокса System Identification. В форме пространства состояний параметры могут быть зафиксированы к их известным значениям или ограничены лечь в предписанной области значений. Если отношение между параметрами или другими ограничениями должно быть задано, объекты IDGREY могут использоваться. Модели IDGREY оценивают заданный пользователями файл MATLAB для оценки системных параметров пространства состояний. Многомерные модели ARX предлагают другую опцию для того, чтобы быстро оценить мультивыходные модели с заданной пользователями структурой.