Создание и оценка моделей процесса с помощью System Identification Toolbox™

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

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

Введение

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

Модели процесса имеют основной тип 'Static Gain + Time Constant + Time Delay'. Они могут быть представлены в виде:

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

Представление моделей процессов с использованием объектов 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-полюсной задержки + может быть оценена путем вызова 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')

Ответ Bode с доверительной областью, соответствующей 3 стандартным отклонениям, может быть вычислен путем:

h = bodeplot(m1,m0);
showConfidence(h,3)

Точно так же данные измерения могут сравниваться с выходами моделей, используя compare следующим образом:

compare(dat1e,m0,m1)

Другие операции, такие как sim, impulse, c2d также доступны, так же как и для других объектов модели.

bdclose('iddempr1')

Учет эффекта поведения интерсampла в оценке

Может быть важным (по крайней мере, для медленной дискретизации) рассмотреть поведение интерсампов входных данных. Чтобы проиллюстрировать это, давайте изучим ту же систему, что и прежде, но без схемы sample-and-hold:

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)

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