Построение и оценка моделей процессов с использованием Toolbox™ идентификации системы

В этом примере используются:

В этом примере показано, как создавать простые модели процессов с помощью 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). Кроме того, на смену реальным полюсам может прийти недостаточно затухающая (сложная) пара полюсов.

Представление моделей процессов с использованием объектов 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 секунды. Как было отмечено, система является функцией непрерывного переноса времени и, следовательно, может быть описана с использованием объектов модели в панели инструментов идентификации системы, таких как 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')

Реакция Боде с доверительной областью, соответствующей 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) м2 (полученные для непрерывного ввода, с предположением 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, когда тестируется на другом наборе данных2:

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)

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