Создание и оценка моделей процессов Используя 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.
Давайте опишем систему с помощью idpoly и idproc объекты. Используя 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)
