Оценка Моделей Непрерывного времени с помощью Данных Simulink
Этот пример иллюстрирует, как модели, симулированные в Simulink®, могут быть идентифицированы с помощью System Identification Toolbox™. Пример описывает, как иметь дело с системами непрерывного времени и задержками, а также важностью междемонстрационного поведения входа.
if exist('start_simulink','file')~=2 disp('This example requires Simulink.') return end
Получение данных моделирования из модели Simulink
Считайте систему описанной следующей моделью 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 первого порядка дает очень хорошее описание данных. Таким образом мы выберем эту модель как наш конечный результат.
Преобразование дискретной модели к непрерывному времени (LTI)
Преобразуйте эту модель в непрерывное время и представляйте его в форме передаточной функции:
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');