时间序列预测和预后预测

打开脚本

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

加载并绘制测量数据

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

负载iddata_timeseriesprediction.n = numel（y）;ns =地板（n / 2）;y_id = y（1：ns，:);y_v = y（（ns + 1：结束），:);data_id = iddata（y_id，[]，ts，'不时'，“小时”）;data_v = iddata(y_v， []， Ts，)'不时'，“小时”，“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 $ 将被识别。

使用SSEST（）命令从测量数据中识别离散状态空间模型。

sys = ssest（data_id，1，'ts'，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 $ 用这个模型来预测 $y_ {n + 10}$ ．请注意，测量值和预测值之间的误差， $y_0 - \帽子{y_0},…,推出——\帽子{最大}$$ ，用于制造 $y_ {n + 10}$ 预测。

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

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

数字;比较（sys，data_v，nstep）与验证数据的比较网格(“上”）;

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

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

MeasuredData = iddata(y， []， Ts，)'不时'，“小时”）;％= [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] =预测（Sys，Itsertsdata，地平线）;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小时超过测量数据的时间范围。基于所识别模型的随机抽样，红色虚线曲线显示了预测响应中的3个SD不确定性。

状态估计

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

为任意状态创建一个预测器。确定的模型协方差需要转换为预测模型使用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 $ 在哪里推出美元测量的电流-功率比，和为炉槽尺寸。在这个例子中，四个直接测量炉槽尺寸，sem，和炉电流功率比，ySizeMeasured，用于估计美元加元．在转换预测器时，预测器的协方差也需要进行转换。因此我们使用translatecov ()命令进行状态坐标转换。

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

EST =离散时间识别的状态空间模型：x（t + ts）= a x（t）+ b u（t）y（t）= c x（t）+ d u（t）a = x1 x1 0.9064B = y1 x1 0.9452 c = x1 y1 0.1001 d = y1 y1 0采样时间：1小时参数化：指数的规范形式：1。馈通：无干扰组件：没有自由系数：2使用“idsdata”，“getpvec”，“getcov”参数及其不确定性。状态:由直接构造或转换产生。不估计。

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

选择= simOptions;选择。InitialCondition = sizeMeasured (1);U = iddata ([], [data_id.Y;data_v。Y], Ts,'不时'，“小时”）;(你们、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%的束缚”，“位置”，'东南'）xlabel（的时间(小时)）ylabel（'振幅'）;网格(“上”）