主要Content

时间序列预测和预测预测

此示例显示了如何创建时间序列模型并使用该模型进行预测,预测和状态估计。测量的数据来自一个感应炉,其插槽大小随时间流失。不能直接测量插槽尺寸,但是测量了炉电流和消耗的功率。众所周知,随着插槽尺寸的增加,插槽电阻会降低。因此,测得的电流与测量功率的比率与插槽大小成正比。您可以使用测量的电流比(电流和功率测量值嘈杂)来创建时间序列模型,并使用该模型估算当前插槽大小并预测未来的插槽大小。通过物理检查,感应炉插槽大小在某些时间点已知。

加载和绘制测量数据

测得的电流比数据存储在iddata_TimeSeriesPredictionMATLAB file. The data is measured at hourly intervals and shows that over time the ratio increases indicating erosion of the furnace slot. You develop a time series model using this data. Start by separating the data into an identification and a validation segment.

加载iddata_TimeSeriesPredictionn = 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(“识别数据”,,,,“验证数据”,,,,'location',,,,'SouthEast');

模型标识

插槽侵蚀可以建模为具有噪声输入的状态空间系统,并测量的电流比为输出。测得的电流比与系统状态成正比,或

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

$ y_n = cx_n + e_n $

在哪里$x_n$状态向量包含插槽大小;$y_n$is the measured current-power ratio;e_n美元noise and$A, C, K$are to be identified.

使用ssest()command to identify a discrete state-space model from the measured data.

系统= ssest(data_id,1,'Ts',,,,Ts,'形式',,,,'典范'
系统= 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$use the model to predict$y_{n+10}$。Note that the error between the measured and predicted values,$ y_0- \ hat {y_0},...,y_n- \ hat {y_n} $,,,,are used to make the$y_{n+10}$预言。

使用10 identifica提前一步预测tion data and the independent validation data.

nstep = 10; compare(sys,data_id,nstep)百分比比较10步预测与估计数据网格('在');

figure; compare(sys,data_v,nstep)% comparison to validation data网格('在');

上面的练习两个数据集都表明预测变量与测量数据匹配。

Forecasting is used to further verify the model. Forecasting uses the measured data record$y_0,y_1,...,y_n-\hat{y_n}$to compute the model state at time step n. This value is used as initial condition for forecasting the model response for a future time span. We forecast the model response over the time span of the validation data and then compare the two. We can also compute the uncertainty in forecasts and plot +/- 3 sd of their values.

MeasuredData = iddata(y, [], Ts,'TimeUnit',,,,'hours');% = [data_id;data_v]t0 =测量data.smplingInstants;Horizo​​n = size(data_v,1);%预测范围[yF, ~, ~, yFSD] = forecast(sys, data_id, Horizon);% Note: yF is IDDATA object while yFSD is a double vectort = yF.SamplingInstants;% extract time samplesyfdata = yf.outputdata;% extract response as double vectorplot(MeasuredData) holdon情节(t,yfdata,'r.-',t,yfdata+3*yfsd,'r--',t,yfdata-3*yfsd,'r--')holdofftitle('Forecasted response over the validation data''s time span')gridon

The plot shows that the model response with confidence intervals (indicated by the red colored dashed curves) overlap the measured value for the validation data. The combined prediction and forecasting results indicate that the model represents the measured current-power ratio.

The forecasting results also show that over large horizons the model variance is large and for practical purposes future forecasts should be limited to short horizons. For the induction furnace model a horizon of 200 hours is appropriate.

最后,我们使用该模型预测响应在502-701小时的时间跨度中对未来的200个步骤。

Horizon = 200;%预测范围[yFuture, ~, ~, yFutureSD] = forecast(sys, MeasuredData, Horizon); t = yFuture.SamplingInstants;% extract time samplesyfuturedata = yfuture.outputdata;% extract response as double vector情节(t0,y,...t, yFutureData,'r.-',,,,...t, yFutureData+3*yFutureSD,'r--',,,,...t, yFutureData-3*yFutureSD,'r--')title(“预测响应(200个步骤)”)gridon

The blue curve shows the measured data that spans over 1-501 hours. The red curve is the forecasted response for 200 hours beyond the measured data's time range. The red dashed curves shows the 3 sd uncertainty in the forecasted response based on random sampling of the identified model.

State Estimation

确定的模型与测得的电流比率匹配,但我们对模型中的状态的炉插槽大小感兴趣。确定的模型具有可以转换的任意状态,以使状态具有含义,在这种情况下为插槽大小。

Create a predictor for the arbitrary state. The identified model covariances need to be translated to the predictor model using thetranslatecov()command. TheCreatePredictor()function simply extracts the third output argument of the预测()function to be used withtranslatecov()

typeCreatePredictorest = 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.

The modelest是在与原始模型相同的状态坐标中表达的1步预测器系统。How do we transform the state coordinates so that the model's state corresponds to the (time dependent) slot size? The solution is to rely on actual, direct measurements of the slot size taken intermittently. This is not uncommon in practice where the cost of taking direct measurements is high and only be done periodically (such as when the component is being replaced).

具体而言,改变预测指标状态,$x_n$,,,,to$z_n$, 以便$ y_n = cz_n $在哪里$y_n$测得的电流比,并$z_n$是炉子槽的大小。在此示例中,炉插槽大小的四个直接测量值,sizeMeasured,,,,and furnace current-power ratio,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.

现在,预测因子在所需的状态坐标中表示。它具有一个输入,该输入是测量的系统输出(炉电流比率),一个输出是预测的系统输出(炉插槽尺寸)。模拟预测变量以估计系统输出和系统状态。

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') 传奇(“测量输出”,,,,“估计输出”,,,,'99.7% bound',,,,'location',,,,'SouthEast')网格('在')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');网格('在'

使用预测和预测预后

预测模型和预测的组合使我们能够对感应炉进行预后。

预测变量模型使我们能够根据测量数据估算当前的炉插槽大小。如果估计值为或接近临界值,则可以安排检查或维护。预测使我们能够从估计的当前状态预测未来的系统行为,使我们能够预测何时需要检查或维护。

Further the predictor and forecast model can be re-identified as more data becomes available. In this example one data set was used to identify the predictor and forecast models but as more data is accumulated the models can be re-identified.

也可以看看

||

Related Topics