Создание и оценка моделей процессов Используя System Identification Toolbox™

Этот пример использует:

В этом примере показано, как создать простые модели процессов с помощью System Identification Toolbox™. Методы для создания этих моделей и оценки их параметров с помощью экспериментальных данных описаны. Этот пример требует Simulink®.

Введение

Этот пример иллюстрирует, как создать простые модели процессов, часто используемые в перерабатывающей промышленности. Простые, передаточные функции непрерывного времени младшего разряда обычно используются, чтобы описать поведение процесса. Такие модели описаны объектами IDPROC, которые представляют передаточную функцию в форме нулевого усиления полюса.

Модели процессов имеют основной тип 'Статическое Усиление + Постоянная времени + Задержка'. Они могут быть представлены как:

$P(s) = K.e^{-T_d*s}.\frac{1 + T_z*s}{(1 + T_{p1}*s)(1 + T_{p2}*s)}$

или как процесс интеграции:

$P(s) = K.e^{-T_d*s}.\frac{1 + T_z*s}{s(1 + T_{p1}*s)(1 + T_{p2}*s)}$

где пользователь может определить количество действительных полюсов (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)

Документация