时间序列预测和预后预测

打开脚本

这个示例展示了如何创建时间序列模型，并使用该模型进行预测、预测和状态估计。测量数据来自一个感应炉，其槽尺寸随着时间的推移而腐蚀。不能直接测量槽的大小，但可以测量炉电流和消耗功率。众所周知，随着缝隙尺寸的增大，缝隙电阻减小。因此，被测电流的平方与被测功率的比值与槽的大小成比例。您可以使用测量到的电流-功率比(电流和功率测量都是嘈杂的)来创建一个时间序列模型，并使用该模型来估计当前槽的大小并预测未来槽的大小。通过物理检查，及时知道感应炉槽的尺寸。

加载并绘制测量数据

测量的电流-功率比数据存储在iddata_TimeSeriesPredictionMATLAB文件。每小时测量一次的数据显示，随着时间的推移，比率增加，表明炉槽的侵蚀。您可以使用这些数据开发一个时间序列模型。首先将数据分离为标识段和验证段。

负载iddata_TimeSeriesPredictionn =元素个数(y);ns =地板(n / 2);y_id = y (1: ns,:);y_v = y ((ns + 1:结束):);data_id = iddata(y_id， []， Ts，)“TimeUnit”，“小时”）;data_v = iddata(y_v， []， Ts，)“TimeUnit”，“小时”，“Tstart”ns + 1);情节(data_id data_v)传说(的识别数据，验证数据的，“位置”，“东南”）;

模式识别

槽蚀可以建模为输入噪声和输出电流功率比的状态空间系统。测量的电流-功率比与系统状态成正比，或

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

$y_n = Cx_n + e_n$

在哪里 x_n美元状态向量，包含槽大小; 推出美元为测量的电流-功率比; e_n美元噪音和 $ A、C、K美元都有待确认。

使用党卫军()命令从测量数据中识别离散状态空间模型。

sys = ss (data_id 1“t”Ts,“形式”，“规范”）

sys =离散时间识别的状态空间模型: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样本时间:1小时参数化:有指标的规范形式:1。使用“idssdata”、“getpvec”、“getcov”等参数及其不确定性。状态:使用SSEST对时域数据“data_id”进行估计。拟合估计数据:67.38%(预测焦点)FPE: 0.09575, MSE: 0.09348

该模型最大限度地减少了一步预测。验证模型使用10步前预测，即，给定 y_0美元,美元…,推出用这个模型来预测 $美元y_ {n + 10} $$ ．请注意，测量值和预测值之间的误差， $y_0 - \帽子{y_0},…,推出——\帽子{最大}$$ ，用来制作 $美元y_ {n + 10} $$ 预测。

对识别数据和独立验证数据使用10步前预测器。

nstep = 10;比较(sys data_id nstep)% 10步预测与估计数据的比较网格(“上”）;

图;比较(sys data_v nstep)与验证数据的比较网格(“上”）;

以上两组数据表明，预测器与实测数据吻合。

通过预测进一步验证了模型的有效性。预测使用测量数据记录 $y_0 y_1…,推出——\帽子{最大}$ 计算时间步长n时的模型状态。该值作为预测未来时间跨度模型响应的初始条件。我们预测模型响应在验证数据的时间跨度，然后比较两者。我们还可以计算预测的不确定性，并绘制其值的+/- 3sd。

MeasuredData = iddata(y， []， Ts，)“TimeUnit”，“小时”）;% = [data_id; data_v]t0 = MeasuredData.SamplingInstants;地平线=大小(data_v, 1);%预测地平线[yF， ~， ~， yFSD] = forecast(sys, data_id, Horizon);%注意:yF是IDDATA对象，而yFSD是一个双向量t = yF.SamplingInstants;提取时间样本百分比yFData = yF.OutputData;%提取响应作为双向量情节(MeasuredData)在情节(t, yFData' r . - 't yFData + 3 * yFSD“r——”t yFData-3 * yFSD“r——”)举行从标题(“验证数据时间范围内的预测响应”网格)在

图中显示，具有置信区间(红色虚线表示)的模型响应与验证数据的测量值重叠。综合预测和预测结果表明，该模型能较好地反映实测电流功率比。

预测结果还表明，在大的视野下，模型方差很大，为实际目的，未来的预测应限于短期的视野。对于感应炉模型来说，200小时的炉位是合适的。

最后，我们利用模型预测200步后502-701小时的反应。

地平线= 200;%预测地平线[yFuture， ~， ~， yFutureSD] = forecast(sys, MeasuredData, Horizon);t = yFuture.SamplingInstants;提取时间样本百分比yFutureData = yFuture.OutputData;%提取响应作为双向量情节(t0, y,．．.t yFutureData' r . - '，．．.t yFutureData + 3 * yFutureSD“r——”，．．.t yFutureData-3 * yFutureSD“r——”)标题(“预测响应(200步)”网格)在

蓝色曲线显示的是1-501小时内的测量数据。红色曲线是超出测量数据时间范围200小时的预测响应。红色虚线显示了基于识别模型随机抽样的预测响应的3sd不确定性。

状态估计

所建立的模型与实测的电流-功率比相吻合，但我们感兴趣的是模型中的炉槽尺寸。所标识的模型具有任意状态，可以对其进行转换，使状态具有意义，在本例中是槽大小。

为任意状态创建一个预测器。确定的模型协方差需要转换为预测模型使用translatecov ()命令。的createPredictor ()方法的第三个输出参数预测()要使用的函数translatecov ()．

类型createPredictorest = translatecov(@(s) createPredictor(s,data_id)，sys)

function pred = createPredictor(mdl,data) % createPredictor返回提前一步的预测器。% % sys = createPredictor(mdl,data) % %为指定的模型mdl %和测量数据创建一个提前1步的预测模型sys。该函数由% |TimeSeriedPredictionExample|和|translatecov()|命令使用，将识别的模型协方差%转换为预测器。% Copyright 2015 The MathWorks, Inc. [~，~，pred] = predict(mdl,data,1);美国东部时间=离散时间状态空间模型发现:x (t + Ts) = x (t) + B u C (t) y (t) = x (t) + D u (t) = x1 x1 0.9064 B = y1 x1 0.09465 C = (x1, y1 1 D = y₁y1 0样品时间:1小时参数化:规范形式指数:1。参数及其不确定性使用“idssdata”，“getpvec”，“getcov”。状态:由直接构造或转换产生。不估计。

该模型美国东部时间提前一步的预测器是否在与原始模型相同的状态坐标中表示sys．我们如何转换状态坐标，使模型的状态对应于(时间依赖的)槽大小?解决方案是依靠实际的，直接的测量槽大小采取间歇。这在实践中并不少见，因为进行直接测量的成本很高，而且只能定期进行(比如替换组件)。

具体来说，改变预测者的状态， x_n美元 , z_n美元 ,所以美元最大= Cz_n $ 在哪里推出美元测量的电流-功率比，和为炉槽尺寸。在这个例子中，四个直接测量炉槽尺寸，sizeMeasured，和炉电流功率比，ySizeMeasured，用于估计美元加元．在转换预测器时，预测器的协方差也需要进行转换。因此我们使用translatecov ()命令进行状态坐标转换。

Cnew = sizeMeasured \ ySizeMeasured;est = translatecov(@(s) ss2ss(s,s.C/Cnew)，est)

美国东部时间=离散时间状态空间模型发现:x (t + Ts) = x (t) + B u C (t) y (t) = x (t) + D u (t) = x1 x1 0.9064 B = y1 x1 0.9452 C = (x1, y1 0.1001 D = y₁y1 0样品时间:1小时参数化:规范形式指数:1。参数及其不确定性使用“idssdata”，“getpvec”，“getcov”。状态:由直接构造或转换产生。不估计。

预测器现在以所需的状态坐标表示。它有一个输入是测量的系统输出(炉电流功率比)和一个输出是预测的系统输出(炉槽大小)。对预测器进行仿真，估计系统输出和系统状态。

选择= simOptions;选择。InitialCondition = sizeMeasured (1);U = iddata ([], [data_id.Y;data_v。Y], Ts,“TimeUnit”，“小时”）;(你们、ye_sd xe] = sim(美国东部时间U选择);

将估计的输出和槽大小与测量值和已知值进行比较。

yesdp =你们;yesdp。Y = ye.Y + 3 * ye_sd;yesdn =你们;yesdn。Y = ye.Y-3 * ye_sd;n =元素个数(xe);人物,情节(yesdp [data_id; data_v],你们,‘g’yesdn,‘g’)传说(测量输出的，“估计输出”，“99.7%的束缚”，“位置”，“东南”)网格(“上”)图、情节(tSizeMeasured sizeMeasured,的r *1: n, xe 1: n, yesdp.Y / est.C,‘g’1: n, yesdn.Y / est.C,‘g’）;传奇(“测量状态”，“估计状态”，“99.7%的束缚”，“位置”，“东南”)包含(的时间(小时)) ylabel (“振幅”）;网格(“上”）