Оценка моделей передаточных функций для теплообменника

Открыть скрипт

Этот пример показывает, как оценить передаточную функцию из данных измеренного сигнала.

Теплообменник

В этом примере мы оцениваем передаточную функцию для теплообменника. Теплообменник состоит из температуры хладагента, температуры продукта и температуры окружающего воздуха нарушения порядка. Мы оценим теплоноситель к передаточной функции температуры продукта.

Конфигурирование измеренных данных

Измеренные данные хранятся в файле MATLAB и включают измерения изменений температуры хладагента вокруг номинальной и изменения температуры продукта вокруг номинальной.

load iddemo_heatexchanger_data

Соберите измеренные данные с помощью iddata команда и постройте график.

data = iddata(pt,ct,Ts);
data.InputName  = '\Delta CTemp';
data.InputUnit  = 'C';
data.OutputName = '\Delta PTemp';
data.OutputUnit = 'C';
data.TimeUnit   = 'minutes';
plot(data)

Оценка передаточной функции

Из физики задачи мы знаем, что теплообменник может быть описан системой первого порядка с задержкой. Используйте tfest команда, задающая один полюс, отсутствие нулей и неизвестную задержку ввода/вывода для оценки передаточной функции.

sysTF = tfest(data,1,0,nan)

sysTF =
 
  From input "\Delta CTemp" to output "\Delta PTemp":
                    1.467
  exp(-0.0483*s) * --------
                   s + 1.56
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 1   Number of zeros: 0
   Number of free coefficients: 2
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using TFEST on time domain data "data".
Fit to estimation data: 36%                      
FPE: 0.02625, MSE: 0.02622

The compare и resid команды позволяют нам исследовать, насколько хорошо оцененная модель соответствует измеренным данным.

set(gcf,'DefaultAxesTitleFontSizeMultiplier',1,...
   'DefaultAxesTitleFontWeight','normal',...
   'Position',[100 100 780 520]);
resid(sysTF,data);


clf
compare(data,sysTF)

Рисунок показывает, что невязки сильно коррелируют, что означает, что в измеренных данных есть информация, которая не была адекватно получена оценочной моделью.

Оценка передаточной функции из начальной системы

Ранее мы оценивали передаточную функцию из данных, но отдельно для системного порядка не включали много знаний априори. Из физики задачи мы знаем, что система стабильна и имеет положительный коэффициент усиления. Проверяя измеренные данные, мы также знаем, что задержка ввода/вывода составляет около 1/5 минуты. Мы используем эту информацию, чтобы создать начальную систему и оценить передаточную функцию, используя эту систему как начальное предположение.

sysInit = idtf(NaN,[1 NaN],'ioDelay',NaN);
sysInit.TimeUnit = 'minutes';

Ограничьте числитель передаточной функции и условия знаменателя так, чтобы система была стабильной с положительным усилением.

sysInit.Structure.num.Value   = 1;
sysInit.Structure.num.Minimum = 0;
sysInit.Structure.den.Value   = [1 1];
sysInit.Structure.den.Minimum = [0 0];

Ограничьте задержку ввода/вывода области значений [0 1] минутой и используйте 1/5 минуту в качестве начального значения.

sysInit.Structure.ioDelay.Value   = 0.2;
sysInit.Structure.ioDelay.Minimum = 0;
sysInit.Structure.ioDelay.Maximum = 1;

Используйте систему как начальное предположение для задачи оценки

sysTF_initialized = tfest(data,sysInit)

sysTF_initialized =
 
  From input "\Delta CTemp" to output "\Delta PTemp":
                    1.964
  exp(-0.147*s) * ---------
                  s + 2.115
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 1   Number of zeros: 0
   Number of free coefficients: 2
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using TFEST on time domain data "data".
Fit to estimation data: 54.09%                   
FPE: 0.01351, MSE: 0.01349                       

resid(sysTF_initialized,data);


clf
compare(data,sysTF,sysTF_initialized)

Оценка модели процесса

В вышесказанном мы рассматривали задачу оценки как задачу оценки передаточной функции, но знаем, что существует некоторая дополнительная структура, которую мы можем навязать. В частности, известно, что система теплообменника является процессом 1-го порядка с задержкой или моделью 'P1D' формы:

$e^{-T_d s}\frac{K_p}{T_p s+1}$

Используйте procest команда для дальнейшего навязывания этой структуры задаче.

sysP1D  = procest(data,'P1D')

sysP1D =
Process model with transfer function:
             Kp                      
  G(s) = ---------- * exp(-Td*s)     
          1+Tp1*s                    
                                     
        Kp = 0.90548                 
       Tp1 = 0.32153                 
        Td = 0.25435                 
                                     
Parameterization:
    {'P1D'}
   Number of free coefficients: 3
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                            
Estimated using PROCEST on time domain data "data".
Fit to estimation data: 70.4%                      
FPE: 0.005614, MSE: 0.005607                       

resid(sysP1D,data);


clf
compare(data,sysTF,sysTF_initialized,sysP1D)

Оценка модели процесса с помощью модели возмущения

Остаточные графики всех оценок, выполненных до сих пор, показывают, что остаточная корреляция все еще высока, что означает, что модель недостаточно богата, чтобы объяснить всю информацию в измеренных данных. Ключом к отсутствующей части является температура окружающего воздуха нарушения порядка, которую мы еще не включили в модель.

Создайте модель процесса 'P1D' с ограничением по значениям задержки и постоянной времени и используйте ее как начальное предположение для задачи оценки.

sysInit = idproc('P1D','TimeUnit','minutes');

Ограничьте модель положительным усилением и задержкой в область значений [0 1] минуту.

sysInit.Structure.Kp.Value    = 1;
sysInit.Structure.Kp.Minimum  = 0;
sysInit.Structure.Tp1.Value   = 1;
sysInit.Structure.Tp1.Maximum = 10;
sysInit.Structure.Td.Value    = 0.2;
sysInit.Structure.Td.Minimum  = 0;
sysInit.Structure.Td.Maximum  = 1;

Задайте опцию для использования модели первого порядка ('ARMA1') для компонента нарушения порядка. Используйте модель шаблона sysInit наряду с набором опций для оценки модели.

opt = procestOptions('DisturbanceModel','ARMA1');
sysP1D_noise = procest(data,sysInit,opt)

sysP1D_noise =
Process model with transfer function:                
             Kp                                      
  G(s) = ---------- * exp(-Td*s)                     
          1+Tp1*s                                    
                                                     
        Kp = 0.91001                                 
       Tp1 = 0.3356                                  
        Td = 0.24833                                 
                                                     
An additive ARMA disturbance model has been estimated
      y = G u + (C/D)e                               
                                                     
      C(s) = s + 591.6                               
      D(s) = s + 3.217                               
                                                     
Parameterization:
    {'P1D'}
   Number of free coefficients: 5
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                            
Estimated using PROCEST on time domain data "data".
Fit to estimation data: 96.86% (prediction focus)  
FPE: 6.307e-05, MSE: 6.294e-05                     

resid(sysP1D_noise,data);


clf
compare(data,sysTF,sysTF_initialized,sysP1D,sysP1D_noise)

Остаточный график ясно показывает, что невязки некоррелированы, что означает, что у нас есть модель, которая объясняет измеренные данные. Компонент нарушения порядка 'ARMA1', который мы оценили, сохранен как значения числителя и знаменателя в свойстве «NoiseTF» модели.

sysP1D_noise.NoiseTF

ans = 

  struct with fields:

    num: {[1 591.6038]}
    den: {[1 3.2172]}

Сравнение различных моделей

Хотя мы определили модель, которая объясняет измеренные данные, мы отмечаем, что модель соответствует измеренным данным составляет около 70%. Потеря значения подгонки является следствием сильного эффекта нарушения порядка окружающей температуры, которое показано следующим образом.

Измеренные данные были получены из модели Simulink со следующими точными значениями (d является нарушением порядка Гауссова шума)

   y = (1-pi/100)*exp(-15s)/(21.3s+1)*u + 1/(25s+1)*d

Создайте модель 'P1D' с этими значениями и посмотрите, насколько хорошо эта модель соответствует измеренным данным.

sysReal = idproc('P1D','TimeUnit','minutes');
sysReal.Kp  = 1-pi/100;
sysReal.Td  = 15/60;
sysReal.Tp1 = 21.3/60;
sysReal.NoiseTF = struct('num',{[1 10000]},'den',{[1 0.04]});

compare(data,sysReal,sysTF,sysTF_initialized,sysP1D,sysP1D_noise);

График сравнения показывает, что истинная система, подобранная к измеренным данным, также составляет около 70%, что подтверждает, что наши оценочные модели не хуже реальной модели в подборе кривой измеренных данных - следствие сильного эффекта нарушения порядка окружающей температуры.

Дополнительная информация

Для получения дополнительной информации об идентификации динамических систем с помощью System Identification Toolbox посетите информационную страницу System Identification Toolbox.

Документация по System Identification Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте