Этот пример иллюстрирует, как модели, моделируемые в Simulink®, могут быть идентифицированы с помощью System Identification Toolbox™. Пример описывает, как иметь дело с непрерывно-разовыми системами и задержками, а также важностью междемонстрационного поведения входа.
if exist('start_simulink','file')~=2 disp('This example requires Simulink.') return end
Считайте систему описанной следующей моделью Simulink:
open_system('iddemsl1') set_param('iddemsl1/Random Number','seed','0')
Красная часть является системой, синяя часть является контроллером, и ссылочный сигнал является развернутой синусоидой (сигнал щебета). Шаг расчета данных установлен в 0,5 секунды.
Эта система может быть представлена с помощью структуры idpoly:
m0 = idpoly(1,0.1,1,1,[1 0.5],'Ts',0,'InputDelay',1,'NoiseVariance',0.01)
m0 =
Continuous-time OE model: y(t) = [B(s)/F(s)]u(t) + e(t)
B(s) = 0.1
F(s) = s + 0.5
Input delays (listed by channel): 1
Parameterization:
Polynomial orders: nb=1 nf=1 nk=0
Number of free coefficients: 2
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
Давайте моделируем модель iddemsl1 и давайте сохраним данные в объекте iddata:
sim('iddemsl1') dat1e = iddata(y,u,0.5); % The IDDATA object
Давайте сделаем вторую симуляцию режима в целях валидации.
set_param('iddemsl1/Random Number','seed','13') sim('iddemsl1') dat1v = iddata(y,u,0.5);
Давайте посмотрим на данные об оценке, полученные во время первой симуляции:
plot(dat1e)
Давайте начнем путем оценки порядка по умолчанию дискретная модель, чтобы получить некоторое предварительное понимание характеристик данных:
m1 = n4sid(dat1e, 'best') % A default order model
m1 =
Discrete-time identified state-space model:
x(t+Ts) = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 0.7881 0.1643 0.1116
x2 -0.1214 0.4223 -0.8489
x3 0.155 0.7527 0.2119
B =
u1
x1 -0.0006427
x2 -0.02218
x3 0.07347
C =
x1 x2 x3
y1 -5.591 0.871 1.189
D =
u1
y1 0
K =
y1
x1 -0.001856
x2 0.002363
x3 -0.06805
Sample time: 0.5 seconds
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 18
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using N4SID on time domain data "dat1e".
Fit to estimation data: 86.17% (prediction focus)
FPE: 0.01281, MSE: 0.01251
Проверяйте, как хорошо модель воспроизводит данные о валидации
compare(dat1v,m1)
Как наблюдается, данные о валидации предсказаны хорошо моделью. Чтобы заняться расследованиями больше в характеристики данных, давайте осмотрим непараметрический импульсный ответ, вычисленный с помощью dat1e, где отрицательные задержки для анализа автоматически определяются:
ImpModel = impulseest(dat1e,[],'negative');
clf
h = impulseplot(ImpModel);
showConfidence(h,3)
ImpModel является моделью FIR чей порядок (нет. из коэффициентов), автоматически определяются. Мы также принимаем решение анализировать эффекты обратной связи путем вычисления импульсного ответа для отрицательных задержек. Влияния от отрицательных задержек не все незначительны. Это происходит из-за регулятора (выходная обратная связь). Это означает, что импульсная оценка ответа не может использоваться, чтобы определить задержку. Вместо этого создайте несколько моделей ARX низкоуровневых с различными задержками и узнайте лучшую подгонку:
V = arxstruc(dat1e,dat1v,struc(1:2,1:2,1:10));
nn = selstruc(V,0) %delay is the third element of nn
nn =
2 2 3
Задержка определяется к 3 задержкам. (Это правильно: deadtime 1 секунды дает две задержки задержки и ZOH-блок другой.) Соответствующая модель ARX может также быть вычислена, можно следующим образом:
m2 = arx(dat1e,nn) compare(dat1v,m1,m2);
m2 =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)
A(z) = 1 - 0.2568 z^-1 - 0.3372 z^-2
B(z) = 0.04021 z^-3 + 0.04022 z^-4
Sample time: 0.5 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nk=3
Number of free coefficients: 4
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARX on time domain data "dat1e".
Fit to estimation data: 85.73% (prediction focus)
FPE: 0.01346, MSE: 0.0133
Эти две модели m1 и m2 ведут себя так же в симуляции. Давайте теперь попытаемся подстроить порядки и задержки. Зафиксируйте задержку с 2 (который вместе с отсутствием сквозного соединения дает сетевую задержку 3 выборок), и найдите модель в пространстве состояний порядка по умолчанию с той задержкой:
m3 = n4sid(dat1e,'best','InputDelay',2,'Feedthrough',false); % Refinement for prediction error minimization using pem (could also use % |ssest|) m3 = pem(dat1e, m3);
Позвольте как взгляд на предполагаемую системную матрицу
m3.a % the A-matrix of the resulting model
ans =
0.7332 -0.3784 0.1735
0.4705 0.3137 -0.6955
-0.0267 0.7527 0.6343
Третья динамика порядка автоматически выбрана, который вместе с 2 "дополнительными" задержками дает 5-ю модель в пространстве состояний порядка.
Всегда желательно не вслепую положиться на автоматический выбор порядка. Они под влиянием случайных ошибок. Хороший путь состоит в том, чтобы посмотреть на нули и полюса модели, наряду с областями уверенности:
clf
h = iopzplot(m3);
showConfidence(h,2) % Confidence region corresponding to 2 standard deviations
Очевидно эти два полюса/нуля в модульном кругу, кажется, отменяют, указывая, что динамика первого порядка может быть достаточной. Используя эту информацию, давайте сделаем новую оценку первого порядка:
m4 = ssest(dat1e,1,'Feedthrough',false,'InputDelay',2,'Ts',dat1e.Ts); compare(dat1v,m4,m3,m1)
График compare показывает, что простая модель m4 первого порядка дает очень хорошее описание данных. Таким образом мы выберем эту модель как наш конечный результат.
Преобразуйте эту модель в непрерывное время и представляйте его в форме передаточной функции:
mc = d2c(m4); idtf(mc)
ans =
From input "u1" to output "y1":
0.09828
exp(-1*s) * ----------
s + 0.4903
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
Хорошее описание системы было получено, как отображено выше.
Непрерывная модель времени может также быть оценена непосредственно. Дискретная модель m4 имеет 2 демонстрационных входных задержки, которые представляют 1 вторую задержку. Мы используем команду ssest для этой оценки:
m5 = ssest(dat1e,1,'Feedthrough',false,'InputDelay',1); present(m5)
m5 =
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
x1 -0.4903 +/- 0.008225
B =
u1
x1 0.01345 +/- 1.13e+11
C =
x1
y1 7.307 +/- 6.138e+13
D =
u1
y1 0
K =
y1
x1 -0.02227 +/- 1.871e+11
Input delays (seconds): 1
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 4
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 4, Number of function evaluations: 9
Estimated using SSEST on time domain data "dat1e".
Fit to estimation data: 87.34% (prediction focus)
FPE: 0.01054, MSE: 0.01047
More information in model's "Report" property.
Параметры модели m5 показывают высокий уровень неуверенности даже при том, что модель соответствует данным 87%. Это вызвано тем, что модель использует больше параметров, чем абсолютно необходимое продвижение к потере уникальности в оценках параметра. Чтобы просмотреть истинный эффект неуверенности в модели, существует два возможных подхода:
Просмотрите неуверенность как доверительные границы на ответе модели, а не на параметрах.
Оцените модель в канонической форме.
Позвольте использованию попробовать оба подхода. Сначала мы оцениваем модель в канонической форме.
m5Canon = ssest(dat1e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); present(m5Canon)
m5Canon =
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
x1 -0.4903 +/- 0.007881
B =
u1
x1 0.09828 +/- 0.001559
C =
x1
y1 1
D =
u1
y1 0
K =
y1
x1 -0.1628 +/- 0.03702
Input delays (seconds): 1
Parameterization:
CANONICAL form with indices: 1.
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 3
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 4, Number of function evaluations: 9
Estimated using SSEST on time domain data "dat1e".
Fit to estimation data: 87.34% (prediction focus)
FPE: 0.01054, MSE: 0.01047
More information in model's "Report" property.
m5Canon использует каноническую параметризацию модели. Это соответствует данным об оценке, столь же хорошим как модель m5. Это показывает маленькую неуверенность в значениях его параметров, дающих свидетельские показания его надежности. Однако когда мы видели m5, большая неуверенность не обязательно означает "плохую" модель. Чтобы установить качество этих моделей, позвольте использованию просмотреть их ответы вовремя и частотные диапазоны с областями уверенности, соответствующими 3 стандартным отклонениям. Мы также строим исходную систему m0 для сравнения.
Диаграмма Боде.
clf opt = bodeoptions; opt.FreqScale = 'linear'; h = bodeplot(m0,m5,m5Canon,opt); showConfidence(h,3) legend show
График шага.
clf
showConfidence(stepplot(m0,m5,m5Canon),3)
legend show
Границы неуверенности для этих двух моделей фактически идентичны. Мы можем так же сгенерировать нулевую полюсом карту (iopzplot) и годограф Найквиста (nyquistplot) с областями уверенности для этих моделей.
idtf(m5)
ans =
From input "u1" to output "y1":
0.09828
exp(-1*s) * ----------
s + 0.4903
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by conversion from idss model.
Когда сравнение непрерывных моделей времени вычислило из выборочных данных, важно рассмотреть междемонстрационное поведение входного сигнала. В примере до сих пор, вход к системе был кусочной константой, из-за "Нулевого порядка Содержат" (zoh) схему в контроллере. Теперь демонтируйте эту схему и рассмотрите действительно непрерывную систему. Сигналы ввода и вывода все еще выбираются 2 Гц, и все остальное - то же самое:
open_system('iddemsl3') sim('iddemsl3') dat2e = iddata(y,u,0.5);
Модели дискретного времени все еще преуспеют на этих данных, с тех пор когда они будут настроены к измерениям, они включат свойства выборки и междемонстрационное входное поведение (для текущего входа). Однако, когда создавание непрерывных моделей времени, знание междемонстрационных свойств важны. Сначала создайте модель так же, как для случая ZOH:
m6 = ssest(dat2e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); idtf(m6)
ans =
From input "u1" to output "y1":
0.1119
exp(-1*s) * ----------
s + 0.5607
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by conversion from idss model.
Давайте сравним предполагаемую модель (m6) с истинной моделью (m0):
step(m6,m0) % Compare with true system
Соглашение теперь не так хорошо. Мы можем, однако, включать в информацию об объекте данных о входе. Как приближение, позвольте ему быть описанным как кусочные линейный (Содержите первый порядок, FOH) между моментами выборки. Эта информация затем используется средством оценки для соответствующей выборки:
dat2e.Intersample = 'foh'; m7 = ssest(dat2e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); % new estimation with correct intersample behavior idtf(m7)
ans =
From input "u1" to output "y1":
0.09937
exp(-1*s) * ----------
s + 0.4957
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by conversion from idss model.
Давайте посмотрим на сравнение переходного процесса снова:
step(m7,m0) % Compare with true system
Эта модель (m7) дает намного лучший результат, чем m6. Это завершает этот пример.
bdclose('iddemsl1'); bdclose('iddemsl3');