熱交換器向け伝達関数モデルの推定

この例では、測定信号データから伝達関数を推定する方法を示します。

熱交換器

この例では、熱交換器の伝達関数を推定します。熱交換器は、冷却水温度、製品温度、および外乱周囲温度で構成されます。ここでは、冷却水温度から製品温度への伝達関数を推定します。

測定データの設定

測定データは 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 つ、零点なしおよび未知の入出力遅延を指定して 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)

この図は、残差の相関が強いことを示しています。これは、推定モデルによって的確に再現されていない情報が測定データに含まれていることを意味します。

初期システムからの伝達関数の推定

ここまでは、データから伝達関数を推定しましたが、システムの次数以外の事前情報を含んでいませんでした。問題の物理現象からは、システムが安定しており、正のゲインをもっていることが判明しています。また、測定データを調べた結果からは、入出力遅延が約 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% であり、周囲温度による外乱の影響が強いために測定データへの適合において推定モデルが実際のモデルを下回っていないことが確認されました。