Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

今後の見通しについての時系列による予測および予想

この例では、時系列モデルを作成し、そのモデルを予測、予想および状態推定に使用する方法を示します。測定データは誘導炉のものであり、そのスロット サイズは時間の経過とともに摩損します。スロット サイズは直接には測定できませんが、炉の電流と消費電力が測定されます。スロット サイズが大きくなるほどスロットの抵抗は減ることがわかっています。したがって、測定電力に対する測定電流の二乗の比率はスロット サイズに比例します。測定電流と測定電力の比率 (電流と電力のいずれの測定値もノイズを含む) を使用して時系列モデルを作成し、そのモデルを使用して現在のスロット サイズの推定と今後のスロット サイズの予想を行います。物理的検査により、誘導炉のスロット サイズはいくつかの時点で把握されます。

測定データの読み込みとプロット

測定電流と測定電力の比率データは、MATLAB ファイル iddata_TimeSeriesPrediction に格納されます。データは 1 時間間隔で測定され、時間の経過とともに比率が大きくなって、炉のスロットが摩損していることが示されます。このデータを使用して時系列モデルを作成します。まず、データを同定と検証のセグメントに分けます。

load iddata_TimeSeriesPrediction
n = numel(y);
ns = floor(n/2);
y_id = y(1:ns,:);
y_v = y((ns+1:end),:);
data_id = iddata(y_id, [], Ts, 'TimeUnit', 'hours');
data_v  = iddata(y_v, [], Ts, 'TimeUnit', 'hours', 'Tstart', ns+1);

plot(data_id,data_v)
legend('Identification data','Validation data','location','SouthEast');

モデルの同定

スロットの摩損は、ノイズ入力があり測定電流と測定電力の比率を出力とする状態空間システムとしてモデル化できます。測定電流と測定電力の比率はシステム状態に比例し、次のように表されます。

$x_{n+1} = Ax_n + Ke_n$

$y_n = Cx_n + e_n$

ここで、$x_n$ は状態ベクトルであり、スロット サイズを含みます。$y_n$ は測定電流と測定電力の比率です。$e_n$ のノイズと $A, C, K$ は特定されます。

ssest() コマンドを使用して、測定データから離散状態空間モデルを同定します。

sys = ssest(data_id,1,'Ts',Ts,'form','canonical')
sys =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + K e(t)
       y(t) = C x(t) + e(t)
 
  A = 
          x1
   x1  1.001
 
  C = 
       x1
   y1   1
 
  K = 
            y1
   x1  0.09465
 
Sample time: 1 hours

Parameterization:
   CANONICAL form with indices: 1.
   Disturbance component: estimate
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                             
Estimated using SSEST on time domain data "data_id".
Fit to estimation data: 67.38% (prediction focus)   
FPE: 0.09575, MSE: 0.09348                          
 

同定されたモデルでは 1 ステップ先の予測が最小化されます。10 ステップ先の予測子を使用してモデルを検証します。すなわち、$y_0,...,y_n$ が与えられた場合に、モデルを使用して $y_{n+10}$ を予測します。測定値と予測値との間の誤差 $y_0-\hat{y_0},...,y_n-\hat{y_n}$$y_{n+10}$ の予測に使用されることに注意してください。

同定データおよび別個の検証データに、10 ステップ先の予測子を使用します。

nstep = 10;
compare(sys,data_id,nstep)  % comparison of 10-step prediction to estimation data
grid('on');

figure; compare(sys,data_v,nstep)  % comparison to validation data
grid('on');

上記においては、両方のデータセットで予測子が測定データと一致していることが示されています。

さらにモデルを検証するために予想を使用します。予想では測定データの記録 $y_0,y_1,...,y_n-\hat{y_n}$ を使用して、タイム ステップ n でのモデルの状態を計算します。この値は、未来の時間範囲でのモデル応答を予想するための初期条件として使用されます。検証データの時間範囲にわたるモデル応答を予想したうえで、この 2 つを比較します。予想に伴う不確かさを計算して、その値の +/- 3 sd をプロットすることもできます。

MeasuredData = iddata(y, [], Ts, 'TimeUnit', 'hours'); % = [data_id;data_v]
t0 = MeasuredData.SamplingInstants;

Horizon = size(data_v,1); % forecasting horizon
[yF, ~, ~, yFSD]  = forecast(sys, data_id, Horizon);
% Note: yF is IDDATA object while yFSD is a double vector
t = yF.SamplingInstants; % extract time samples
yFData = yF.OutputData;      % extract response as double vector
plot(MeasuredData)
hold on
plot(t, yFData, 'r.-', t, yFData+3*yFSD, 'r--', t, yFData-3*yFSD, 'r--')
hold off
title('Forecasted response over the validation data''s time span')
grid on

プロットでは、信頼区間 (赤の破線による曲線) を伴うモデル応答が、検証データに対する測定値とオーバーラップすることが示されています。予測と予想との結果を組み合わせると、モデルが測定電流と測定電力の比率を表現していることが示されます。

また、予想結果には、対象期間が長いとモデルの分散が大きく、実用目的としては将来の予想期間を短期に限定すべきであることも示されています。誘導炉モデルでは、対象期間を 200 時間にするのが妥当です。

最後に、モデルを使用して、将来の 200 ステップでの応答を 502 ~ 701 時間の時間範囲について予想します。

Horizon = 200; % forecasting horizon
[yFuture, ~, ~, yFutureSD] = forecast(sys, MeasuredData, Horizon);
t = yFuture.SamplingInstants; % extract time samples
yFutureData = yFuture.OutputData;      % extract response as double vector
plot(t0, y,...
   t, yFutureData, 'r.-', ...
   t, yFutureData+3*yFutureSD, 'r--', ...
   t, yFutureData-3*yFutureSD, 'r--')
title('Forecasted response (200 steps)')
grid on

青の曲線は、1 ~ 501 時間の範囲にわたる測定データを示しています。赤の曲線は、測定データの時間範囲を超えた 200 時間で予想される応答です。赤の破線による曲線は、同定されたモデルの無作為抽出に基づいた、予想される応答における 3 sd の不確かさを示しています。

状態の推定

同定されたモデルは測定電流と測定電力の比率と一致しますが、問題となっているのはモデルの状態である炉のスロット サイズです。同定されたモデルには任意の状態があり、これを変換して、今の場合であればスロット サイズとして意味をもたせることができます。

その任意の状態についての予測子を作成します。同定されたモデルの共分散は、translatecov() コマンドを使用して予測子モデルに変換する必要があります。関数 createPredictor() は、translatecov() で使用するために、関数 predict() の 3 番目の出力引数を単純に抽出します。

type createPredictor
est = translatecov(@(s) createPredictor(s,data_id),sys)
function pred = createPredictor(mdl,data)
%CREATEPREDICTOR Return 1-step ahead predictor.
%
%   sys = createPredictor(mdl,data)
%
%   Create a 1-step ahead predictor model sys for the specified model mdl
%   and measured data. The function is used by
%   |TimeSeriedPredictionExample| and the |translatecov()| command to
%   translate the identified model covariance to the predictor.

% Copyright 2015 The MathWorks, Inc.
[~,~,pred] = predict(mdl,data,1);

est =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t)
       y(t) = C x(t) + D u(t)
 
  A = 
           x1
   x1  0.9064
 
  B = 
            y1
   x1  0.09465
 
  C = 
       x1
   y1   1
 
  D = 
       y1
   y1   0
 
Sample time: 1 hours

Parameterization:
   CANONICAL form with indices: 1.
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.
 

モデル est は 1 ステップ先の予測子であり、元のモデル sys と同じ状態座標で表現されます。モデルの状態を (時間に依存する) スロット サイズに対応付けるには、状態座標をどのように変換すればよいでしょうか。その解決法は、断続的に取得される、スロット サイズの実際の直接測定値に基づきます。これは、測定値を直接取得するためのコストが高く、(コンポーネントを交換する際など) 定期的にのみ行われるような場合、実際問題として珍しいことではありません。

具体的には、$y_n = Cz_n$ となるように、予測子の状態 $x_n$$z_n$ に変換します。ここで、$y_n$ は測定電流と測定電力の比率、$z_n$ は炉のスロット サイズです。この例では、炉のスロット サイズの直接測定値 sizeMeasured を 4 つと、炉の電流と電力の比率 ySizeMeasured を使用して、$C$ を推定します。予測子の変換では、予測子の共分散も変換する必要があります。このため、translatecov() コマンドを使用して状態座標の変換を実行します。

Cnew = sizeMeasured\ySizeMeasured;
est  = translatecov(@(s) ss2ss(s,s.C/Cnew),est)
est =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t)
       y(t) = C x(t) + D u(t)
 
  A = 
           x1
   x1  0.9064
 
  B = 
           y1
   x1  0.9452
 
  C = 
           x1
   y1  0.1001
 
  D = 
       y1
   y1   0
 
Sample time: 1 hours

Parameterization:
   CANONICAL form with indices: 1.
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.
 

これで、予測子が望ましい状態座標で表現されるようになりました。測定されたシステム出力 (炉の電流と電力の比率) である入力が 1 つ、予測されたシステム出力 (炉のスロット サイズ) である出力が 1 つ設けられています。予測子がシミュレートされて、システム出力とシステム状態が推定されます。

opts = simOptions;
opts.InitialCondition = sizeMeasured(1);
U = iddata([],[data_id.Y; data_v.Y],Ts,'TimeUnit','hours');
[ye,ye_sd,xe] = sim(est,U,opts);

推定した出力およびスロット サイズを、測定値および既知の値と比較します。

yesdp   = ye;
yesdp.Y = ye.Y+3*ye_sd;
yesdn   = ye;
yesdn.Y = ye.Y-3*ye_sd;
n = numel(xe);
figure, plot([data_id;data_v],ye,yesdp,'g',yesdn,'g')
legend('Measured output','Estimated output','99.7% bound','location','SouthEast')
grid('on')
figure, plot(tSizeMeasured,sizeMeasured,'r*',1:n,xe,1:n,yesdp.Y/est.C,'g',1:n,yesdn.Y/est.C,'g');
legend('Measured state','Estimated state','99.7% bound','location','SouthEast')
xlabel('Time (hours)')
ylabel('Amplitude');
grid('on')

予測と予想を使用した今後の見通し

予測子モデルと予想との組み合わせにより、誘導炉について今後の見通しを立てることができます。

予測子モデルでは、測定データを基に現在の炉のスロット サイズを推定できます。推定値が棄却限界値かそれに近い値であれば、検査または保守作業をスケジュールできます。予想によって、推定された現在の状態から今後のシステムの動作を予測することができ、検査や保守がいつ必要になるか予測できるようになります。

さらに、予測子と予想モデルは、より多くのデータが入手可能になると同定し直すことができます。この例では 1 つのデータセットが予測子と予想モデルの同定に使用されましたが、より多くのデータが蓄積されるとモデルを同定し直すことができます。

参考

| |

関連するトピック