时间序列预测与预后预测

打开脚本

本例展示了如何创建时间序列模型，并使用该模型进行预测、预测和状态估计。测量数据来自一个感应炉，其槽尺寸随时间侵蚀。槽尺寸无法直接测量，但会测量炉电流和消耗功率。众所周知，随着槽尺寸的增加，温度会降低e插槽电阻减小。因此，测量电流平方与测量功率之比与插槽大小成正比。使用测量电流功率比（电流和功率测量都有噪声）创建一个时间序列模型，并使用该模型估计当前槽尺寸和预测未来槽尺寸。通过物理检查，在某些时间点上感应炉槽尺寸已知。

加载并绘制测量数据

测量的电流功率比数据存储在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）;plot（data_id，data_v）图例（“识别数据”,“验证数据”,'地点',“东南”）;

模型识别

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

$x_ {n + 1} = ax_n + ke_n$

$y\u n=Cx\u n+e\u n$

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

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

sys = ssest（data_id，1，“t”Ts,“形式”,“规范的”)

sys =离散时间识别的状态空间模型：x（t + ts）= x（t）+ k e（t）y（t）y（t）= c x（t）+ e（t）a = x1 x1 1.001 c= x1 y1 1 k = y1 x1 0.09465采样时间：1小时参数化：指数规范形式：1。干扰组件：估计自由系数的数量：2使用“idssdata”，“getpvec”，“getcov”参数的参数及其不确定性．状态：使用SSEST上的时域数据“data_id”估计。适合估算数据：67.38％（预测焦点）FPE：0.09575，MSE：0.09348

确定的模型最小化了1步超前预测。使用10步超前预测器验证模型，即 y_0美元,美元…,推出使用模型预测 $美元y_ {n + 10} $$ ．请注意，测量值与预测值之间的错误， $y_0-\hat{y_0}，…，y_n-\hat{y_n}$ ，用来制作 $美元y_ {n + 10} $$ 预言

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

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

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

上述两个数据集的练习表明，预测值与测量数据相匹配。

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

测量数据= iddata（y，[]，ts，“TimeUnit”,'小时'）;% = [data_id; data_v]t0=测量数据。采样率；水平线=尺寸（数据，1）；%预测期[YF，〜，〜，YFSD] =预测（SYS，DATA_ID，地平线）;％注意：YF是IDDATA对象，而YFSD是一个双矢量t = yf.samplinginstants;提取时间样本百分比yFData=yF.OutputData；%提取响应作为双向量绘图（测量模具）持有在绘图（t，yf数据，'r.-'，t，yFData+3*yFSD，'r--'，t，yFData-3*yFSD，'r--')举行离开标题（“验证数据时间范围内的预测响应”） 网格在

该图表明，具有置信区间的模型响应（由红色彩色虚线曲线表示）与验证数据的测量值重叠。组合预测和预测结果表明该模型代表了测量的电流功率比。

预测结果还表明，在大地平线上，模型方差大，实际目的，未来预测应限于空间。对于感应炉模型，200小时的地平线是合适的。

最后，我们使用该模型预测了502-701小时内未来200步的响应。

地平线= 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）

函数pred = createpredictor（mdl，data）％createpredictor返回1步前预测器。%% sys = createpredictor（mdl，data）%%为指定型号MDL％和测量数据创建一个1步前方预测测量模型SYS。该功能由％| TimesEriedPredictionExample使用和|翻译（）|命令百分比将标识的模型协方差转换为预测器。％Copyright 2015 MathWorks，Inc。[〜，〜，pred] =预测（MDL，数据，1）;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.09465 c = x1 y1 1 d = y1 y1 0采样时间：1小时参数化：指数的规范形式：1。馈通：无干扰组件：没有自由系数：2使用“idsdata”，“getpvec”，“getcov”参数及其不确定性。状态:由直接构造或转换产生。未估计。

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

具体来说，转换预测器状态， $ x_n $ ，到 $ z_n $ 因此 $y\u n=Cz\u n$ 哪里 $ y_n $ 测量的电流功率比，以及 $ z_n $ 是炉槽尺寸。在本例中，炉槽尺寸的四个直接测量值，sizeMeasured和炉电流功率比，ysizemeasured.，用于估算 $ C $ 。在转换预测值时，还需要转换预测值协方差。因此，我们使用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;opts.initialcondition = sizemeasured（1）;u = iddata（[]，[data_id.y; data_v.y]，ts，“TimeUnit”,'小时'）;[ye，ye_sd，xe] = sim（est，你，选择）;

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

yesdp = ye;yesdp.y = ye.y + 3 * ye_sd;是的，=你;yesdn.y = ye.y-3 * ye_sd;n = numel（xe）;图，绘图（[data_id; data_v]，ye，yesdp，'G'，yesdn，'G')传奇(“测量输出”,“估计产量”,'9.7％绑定','地点',“东南”)网格(“开”）图，情节（TsizeMeasured，SizeMeasured，的r *，1：n，xe，1：n，yesdp.y / est.c，'G'，1:n，yesdn.Y/est.C，'G'）;传奇(“测量状态”,“估计状态”,'9.7％绑定','地点',“东南”)包含('时间（小时）') ylabel (“振幅”）;网格(“开”)