Main Content

열 교환기에 대한 전달 함수 모델 추정하기

이 예제에서는 측정된 신호 데이터에서 전달 함수를 추정하는 방법을 보여줍니다.

열 교환기

이 예제에서는 열 교환기에 대한 전달 함수를 추정합니다. 열 교환기는 냉각수 온도, 제품 온도, 외란 주위 온도로 구성됩니다. 냉각수와 제품 간 온도 전달 함수를 추정해 보겠습니다.

측정된 데이터 구성하기

측정된 데이터는 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)

전달 함수 추정

이 문제의 물리학에 따르면, 열 교환기는 지연이 있는 1차 시스템으로 설명할 수 있습니다. 전달 함수를 추정하려면 1개의 극점, 0개의 영점, 미지의 입력/출력을 지정하는 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                       
 

compareresid 명령을 사용하면 추정된 모델이 측정된 데이터와 얼마나 잘 일치하는지 조사할 수 있습니다.

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

clf
compare(data,sysTF)

위 그림은 잔차의 상관관계가 긴밀함을 나타내고 있으며, 이는 추정된 모델에서 적절히 포착되지 않은 정보가 측정된 데이터에 있음을 암시합니다.

초기 시스템으로부터 전달 함수 추정

이전 예제에서 데이터로부터 전달 함수를 추정했지만, 이때 시스템 차수를 제외하고는 선험적 정보를 충분히 포함하지 않았습니다. 문제의 물리학에 따르면, 이 시스템은 안정적이고 양의 이득을 가진다는 사실이 알려져 있습니다. 측정된 데이터를 조사하면 입력/출력 지연이 1분의 약 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;

외란 성분에 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%임을 보여줍니다. 즉, 측정된 데이터의 피팅에서 추정된 모델이 실제 모델보다 더 나쁘지 않음을 확인할 수 있으며, 이러한 손실은 주위 온도 외란의 영향을 강하게 받은 결과입니다.