Создание и оценка моделей процессов Используя System Identification Toolbox™
Этот пример показывает, как создать простые модели процессов с помощью System Identification Toolbox™. Методы для создания этих моделей и оценки их параметров с помощью экспериментальных данных описаны. Этот пример требует Simulink®.
Введение
Этот пример иллюстрирует, как создать простые модели процессов, часто используемые в перерабатывающей промышленности. Простые, непрерывно-разовые передаточные функции младшего разряда обычно используются, чтобы описать поведение процесса. Такие модели описаны объектами IDPROC, которые представляют передаточную функцию в форме нулевого усиления полюса.
Модели процессов имеют основной тип 'Статическое Усиление + Временная константа + Задержка'. Они могут быть представлены как:
или как процесс интеграции:
где пользователь может определить количество действительных полюсов (0, 1, 2 или 3), а также присутствие нуля в числителе, присутствие термина интегратора (1/s) и присутствие задержки (Td). Кроме того, underdamped (комплексная) пара полюсов может заменить действительные полюса.
Представление Моделей процессов с помощью Объектов IDPROC
Объекты IDPROC задают модели процессов при помощи букв P (для модели процесса), D (для задержки), Z (для нуля) и I (для интегратора). Целое число обозначит количество полюсов. Модели сгенерированы путем вызова idproc с вектором символов, созданным с помощью этих букв.
Например:
idproc('P1') % transfer function with only one pole (no zeros or delay) idproc('P2DIZ') % model with 2 poles, delay integrator and delay idproc('P0ID') % model with no poles, but an integrator and a delay
ans =
Process model with transfer function:
Kp
G(s) = ----------
1+Tp1*s
Kp = NaN
Tp1 = NaN
Parameterization:
'P1'
Number of free coefficients: 2
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
ans =
Process model with transfer function:
1+Tz*s
G(s) = Kp * ------------------- * exp(-Td*s)
s(1+Tp1*s)(1+Tp2*s)
Kp = NaN
Tp1 = NaN
Tp2 = NaN
Td = NaN
Tz = NaN
Parameterization:
'P2DIZ'
Number of free coefficients: 5
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
ans =
Process model with transfer function:
Kp
G(s) = --- * exp(-Td*s)
s
Kp = NaN
Td = NaN
Parameterization:
'P0DI'
Number of free coefficients: 2
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
Создание объекта IDPROC (использование модели Simulink® как пример)
Считайте систему описанной следующей моделью SIMULINK:
open_system('iddempr1') set_param('iddempr1/Random Number','seed','0')
Красная часть является системой, синяя часть является контроллером, и ссылочный сигнал является развернутой синусоидой (сигнал щебета). Шаг расчета данных установлен в 0,5 секунды. Как наблюдается, система является непрерывно-разовой передаточной функцией и может следовательно быть описана с помощью объектов модели в System Identification Toolbox, таких как idss, idpoly или idproc.
Давайте опишем систему с помощью объектов idproc и idpoly. Используя объект idpoly, система может быть описана как:
m0 = idpoly(1,0.1,1,1,[1 0.5],'Ts',0,'InputDelay',1.57,'NoiseVariance',0.01);
Форма IDPOLY, используемая выше, полезна для описания передаточных функций произвольных порядков. Поскольку система, которую мы рассматриваем здесь, является довольно простой (один полюс и никакие нули), и является непрерывно-разовой, мы можем использовать более простой объект IDPROC получить его динамику:
m0p = idproc('p1d','Kp',0.2,'Tp1',2,'Td',1.57) % one pole+delay, with initial values % for gain, pole and delay specified.
m0p =
Process model with transfer function:
Kp
G(s) = ---------- * exp(-Td*s)
1+Tp1*s
Kp = 0.2
Tp1 = 2
Td = 1.57
Parameterization:
'P1D'
Number of free coefficients: 3
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
Оценка параметров моделей IDPROC
Если система описана объектом модели, таким как IDPROC, она может использоваться для оценки его параметров с помощью данных об измерении. Как пример, мы считаем проблему оценки параметров системы модели Simulink (красный фрагмент) использованием данных моделирования. Мы начинаем путем получения данных для оценки:
sim('iddempr1') dat1e = iddata(y,u,0.5); % The IDDATA object for storing measurement data
Давайте посмотрим на данные:
plot(dat1e)
Мы можем идентифицировать модель процесса с помощью команды procest путем предоставления той же информации структуры, указанной, чтобы создать модели IDPROC. Например, 1-pole+delay модель может быть оценена путем вызова procest можно следующим образом:
m1 = procest(dat1e,'p1d'); % estimation of idproc model using data 'dat1e'. % Check the result of estimation: m1
m1 =
Process model with transfer function:
Kp
G(s) = ---------- * exp(-Td*s)
1+Tp1*s
Kp = 0.20045
Tp1 = 2.0431
Td = 1.499
Parameterization:
'P1D'
Number of free coefficients: 3
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "dat1e".
Fit to estimation data: 87.34%
FPE: 0.01069, MSE: 0.01062
To get information about uncertainties, use
present(m1)
m1 =
Process model with transfer function:
Kp
G(s) = ---------- * exp(-Td*s)
1+Tp1*s
Kp = 0.20045 +/- 0.00077275
Tp1 = 2.0431 +/- 0.061216
Td = 1.499 +/- 0.040854
Parameterization:
'P1D'
Number of free coefficients: 3
Use "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 PROCEST on time domain data "dat1e".
Fit to estimation data: 87.34%
FPE: 0.01069, MSE: 0.01062
More information in model's "Report" property.
Параметры модели, K, Tp1 и Td теперь показывают с одной областью значений неуверенности стандартного отклонения.
Вычисление времени и частотной характеристики моделей IDPROC
Модель m1, оцененная выше, является объектом модели IDPROC, к которому могут быть применены все образцовые команды тулбокса:
step(m1,m0) %step response of models m1 (estimated) and m0 (actual) legend('m1 (estimated parameters)','m0 (known parameters)','location','northwest')
Предвещайте ответ с областью уверенности, соответствующей 3 стандартным отклонениям, может быть вычислен путем выполнения:
h = bodeplot(m1,m0); showConfidence(h,3)
Точно так же данные об измерении могут сравниться с моделями выходные параметры с помощью compare можно следующим образом:
compare(dat1e,m0,m1)
Другие операции, такие как sim, impulse, c2d также доступен, как они для других объектов модели.
bdclose('iddempr1')
Размещение эффекта междемонстрационного поведения по оценке
Для медленной выборки) может быть важно (по крайней мере, рассмотреть междемонстрационное поведение входных данных. Чтобы проиллюстрировать это, давайте изучим ту же систему как прежде, но без схемы выборки и хранения:
open_system('iddempr5')
Моделируйте эту систему с тем же шагом расчета:
sim('iddempr5') dat1f = iddata(y,u,0.5); % The IDDATA object for the simulated data
Мы оцениваем модель IDPROC с помощью data1f, также наложение верхней границы на допустимом значении задерживается. Мы будем использовать 'lm' в качестве метода поиска и также принимать решение просмотреть прогресс оценки.
m2_init = idproc('P1D'); m2_init.Structure.Td.Maximum = 2; opt = procestOptions('SearchMethod','lm','Display','on'); m2 = procest(dat1f,m2_init,opt); m2
m2 =
Process model with transfer function:
Kp
G(s) = ---------- * exp(-Td*s)
1+Tp1*s
Kp = 0.20038
Tp1 = 2.01
Td = 1.31
Parameterization:
'P1D'
Number of free coefficients: 3
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "dat1f".
Fit to estimation data: 87.26%
FPE: 0.01067, MSE: 0.01061
Эта модель имеет немного менее точную оценку задержки, чем предыдущая, m1:
[m0p.Td, m1.Td, m2.Td] step(m0,m1,m2) legend('m0 (actual)','m1 (estimated with ZOH)','m2 (estimated without ZOH)','location','southeast')
ans =
1.5700 1.4990 1.3100
Однако путем сообщения процесса оценки, что междемонстрационное поведение является хранением первый порядок (приближение к непрерывной истине) вход, мы добиваемся большего успеха:
dat1f.InterSample = 'foh';
m3 = procest(dat1f,m2_init,opt);
Сравните эти четыре модели m0 (TRUE) m1 (полученный из входа zoh) m2 (полученный для непрерывного входа, с zoh предположением) и m3 (полученный для того же входа, но с foh предположением)
[m0p.Td, m1.Td, m2.Td, m3.Td] compare(dat1e,m0,m1,m2,m3)
ans =
1.5700 1.4990 1.3100 1.5570
step(m0,m1,m2,m3) legend('m0','m1','m2','m3') bdclose('iddempr5')
Моделирование системы, действующей в замкнутом цикле
Давайте теперь рассмотрим более комплексный процесс с интегрированием, которое управляется в замкнутом цикле:
open_system('iddempr2')
Истинная система может быть представлена:
m0 = idproc('P2ZDI','Kp',1,'Tp1',1,'Tp2',5,'Tz',3,'Td',2.2);
Процессом управляет регулятор PD с ограниченной входной амплитудой, и нулевой порядок содержат устройство. Шаг расчета составляет 1 секунду.
set_param('iddempr2/Random Number','seed','0') sim('iddempr2') dat2 = iddata(y,u,1); % IDDATA object for estimation
Две различных симуляции сделаны, первое для оценки и второе в целях валидации.
set_param('iddempr2/Random Number','seed','13') sim('iddempr2') dat2v = iddata(y,u,1); % IDDATA object for validation purpose
Давайте посмотрим на данные (оценка и валидация).
plot(dat2,dat2v) legend('dat2 (estimation)','dat2v (validation)')
Давайте теперь выполним оценку с помощью dat2.
Warn = warning('off','Ident:estimation:underdampedIDPROC'); m2_init = idproc('P2ZDI'); m2_init.Structure.Td.Maximum = 5; m2_init.Structure.Tp1.Maximum = 2; opt = procestOptions('SearchMethod','lsqnonlin','Display','on'); opt.SearchOptions.MaxIterations = 100; m2 = procest(dat2, m2_init, opt)
m2 =
Process model with transfer function:
1+Tz*s
G(s) = Kp * ------------------- * exp(-Td*s)
s(1+Tp1*s)(1+Tp2*s)
Kp = 0.98412
Tp1 = 2
Tp2 = 1.4838
Td = 1.713
Tz = 0.027244
Parameterization:
'P2DIZ'
Number of free coefficients: 5
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "dat2".
Fit to estimation data: 91.51%
FPE: 0.1128, MSE: 0.1092
compare(dat2v,m2,m0) % Gives very good agreement with data
bode(m2,m0)
legend({'m2 (est)','m0 (actual)'},'location','west')
impulse(m2,m0)
legend({'m2 (est)','m0 (actual)'})
Сравните также с параметрами истинной системы:
present(m2) [getpvec(m0), getpvec(m2)]
m2 =
Process model with transfer function:
1+Tz*s
G(s) = Kp * ------------------- * exp(-Td*s)
s(1+Tp1*s)(1+Tp2*s)
Kp = 0.98412 +/- 0.013672
Tp1 = 2 +/- 8.2231
Tp2 = 1.4838 +/- 10.193
Td = 1.713 +/- 63.703
Tz = 0.027244 +/- 65.516
Parameterization:
'P2DIZ'
Number of free coefficients: 5
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Change in cost was less than the specified tolerance..
Number of iterations: 3, Number of function evaluations: 4
Estimated using PROCEST on time domain data "dat2".
Fit to estimation data: 91.51%
FPE: 0.1128, MSE: 0.1092
More information in model's "Report" property.
ans =
1.0000 0.9841
1.0000 2.0000
5.0000 1.4838
2.2000 1.7130
3.0000 0.0272
Предостережение. Идентификация нескольких оперативных констант может иногда быть плохо обусловленной проблемой, особенно если данные собраны в замкнутом цикле.
Чтобы проиллюстрировать это, давайте оценим основанное на модели на данных о валидации:
m2v = procest(dat2v, m2_init, opt) [getpvec(m0), getpvec(m2), getpvec(m2v)]
m2v =
Process model with transfer function:
1+Tz*s
G(s) = Kp * ------------------- * exp(-Td*s)
s(1+Tp1*s)(1+Tp2*s)
Kp = 0.95747
Tp1 = 1.999
Tp2 = 0.60819
Td = 2.314
Tz = 0.0010561
Parameterization:
'P2DIZ'
Number of free coefficients: 5
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "dat2v".
Fit to estimation data: 90.65%
FPE: 0.1397, MSE: 0.1353
ans =
1.0000 0.9841 0.9575
1.0000 2.0000 1.9990
5.0000 1.4838 0.6082
2.2000 1.7130 2.3140
3.0000 0.0272 0.0011
Эта модель имеет намного худшие значения параметров. С другой стороны, это выполняет почти тождественно к истинной системе m0, когда протестировано на другом наборе данных dat2:
compare(dat2,m0,m2,m2v)
Фиксация известных параметров во время оценки
Предположим, что мы знаем из других источников, что одна временная константа равняется 1:
m2v.Structure.Tp1.Value = 1; m2v.Structure.Tp1.Free = false;
Мы можем зафиксировать это значение при оценке других параметров:
m2v = procest(dat2v,m2v)
%
m2v =
Process model with transfer function:
1+Tz*s
G(s) = Kp * ------------------- * exp(-Td*s)
s(1+Tp1*s)(1+Tp2*s)
Kp = 1.0111
Tp1 = 1
Tp2 = 5.3014
Td = 2.195
Tz = 3.231
Parameterization:
'P2DIZ'
Number of free coefficients: 4
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "dat2v".
Fit to estimation data: 92.05%
FPE: 0.09952, MSE: 0.09794
Как наблюдается, фиксирование Tp1 к его известному значению существенно улучшает оценки остающихся параметров в модели m2v.
Это также указывает, что простое приближение должно преуспеть на данных:
m1x_init = idproc('P2D'); % simpler structure (no zero, no integrator) m1x_init.Structure.Td.Maximum = 2; m1x = procest(dat2v, m1x_init) compare(dat2,m0,m2,m2v,m1x)
m1x =
Process model with transfer function:
Kp
G(s) = ----------------- * exp(-Td*s)
(1+Tp1*s)(1+Tp2*s)
Kp = -1.2555
Tp1 = 1.0249e-06
Tp2 = 0.078021
Td = 1.959
Parameterization:
'P2D'
Number of free coefficients: 4
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "dat2v".
Fit to estimation data: -23.87%
FPE: 24.15, MSE: 23.77
Таким образом более простая модель может оценить систему вывод вполне прилично. Однако m1x не содержит интегрирования, таким образом, долговременное поведение области значений разомкнутого цикла будет очень отличаться:
step(m0,m2,m2v,m1x) legend('m0','m2','m2v','m1x') bdclose('iddempr2') warning(Warn)
