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

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

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

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

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

Конфигурирование результатов измерений

Результаты измерений хранятся в файле 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

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

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


clf
compare(data,sysTF)

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

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

Ранее мы оценили, что передаточная функция из данных, но независимо для системного порядка не включала много apriori знания. От физики проблемы мы знаем, что система устойчива и имеет положительное усиление. Осмотр результатов измерений, мы также знаем, что задержка ввода/вывода вокруг 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

Поддержка

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