应用随机过程建模报告

1、Harbin Institute of Technology 课程设计（论文）课程名称： 应用随机过程 设计题目： 建 模 院 系： 电子与信息工程学院 班 级： 通信1班 设 计 者： 学 号： 指导教师： 设计时间： 2013-11-9 哈尔滨工业大学线性模型电力负荷时间序列建模1 电力系统负荷预测的意义随着我国电力事业的发展，电网的管理日趋现代化，对电力系统负荷预测问题的研究也越来越引起人们的注意。电力负荷预测是电力系统调度、用电、计划、规划等管理部门的重要工作之一。提高负荷预测技术水平，有利于计划用电管理，有利于合理安排电网运行方式和机组检修计划，有利于节煤、节油和降低发电成本，有利于

2、制定合理的电源建设规划，有利于提高电力系统的经济效益和社会效益。电力负荷预测，为编制电力规划提供依据，是电网规划的基础，它规定了电力工业的发展水平、发展速度、源动力资源的需求量，电力工业发展的资金需求量，以及电力工业发展对人力资源的需求量。因此，国内外许多专家和学者开始致力于现代负荷预测方法的研究，而时间序列模型在国际和国内的电力系统短期负荷预测中得到了广泛应用。2 平稳时间序列及其随机线性模型时间序列是指随时间改变而随机的变化的序列。时间序列分析分为时域分析和频域分析，前者是对时间序列在时间域上的各种平均值进行分析研究，后者是进行傅里叶变换以后在频率域进行谱分析。随着计算机技术的飞速发展，时

3、域分析方法为人们所关注。本文所要研究的就是时域分析。平稳时间序列是平稳序列，它满足期望为0，且任意两个时刻的相关函数与时间无关，仅与两个时刻的时间差相关。因为我们所掌握的为平稳时间序列的线性随机模型，而在实际中所遇到的一般都不是平稳时间序列，这就要对其进行相关的处理，使其变化为平稳序列。均值为0且具有有理谱密度的平稳时间序列必可表示为下面三种形式中的一种（其中为白噪声）：（1）自回归模型AR模型AR（p）模型由p2参数来刻画;（2）滑动平均模型MA模型MA(q)模型由q2参数刻画;（3）自回归滑动平均模型或混合模型ARMA模型 ARMA(p,q)混和模型由pq3参数刻画；通过以上介绍可以看出我

4、们可以把AR(p)和MA(q)模型看成APMA(p,q)的两种特例。线性模型中有两个重要的参数：自相关函数和和偏相关函数。其中偏相关函数刻画了平稳序列任意一个长的片段在中间量固定的条件下，两端的线性密切程度，而自相关函数也是刻画两端的线性密切程度，但并不需要中间数值固定。与这两个参数相关的性质为拖尾性和截尾性。拖尾性是指（或）随k无限增长以负指数的速度趋向于0，其图像像拖一条尾巴；截尾性（或）是指参数在k>p或k>q后，其值变为零，其图像像截断了的尾巴一样。总结上面的分析，得出时间序列模型的基本特性，如表1。表1 时间序列模型的基本特性类型 基本方程平稳条件的根全在单位圆外无平稳条

5、件的根全在单位圆外可逆条件无可逆条件根全在单位圆外的根全在单位圆外自相关函数拖尾截尾拖尾偏相关函数截尾拖尾拖尾由模型的基本特性可知，自相关函数具有拖尾性，若算出的偏相关函数具有截尾性，模型为自回归模型（AR）；当算出的自相关函数具有截尾性，偏相关函数具有拖尾性，则能判断模型为滑动平均模型（MA）；若自相关函数和偏相关函数都具有拖尾性，模型为自回归滑动平均模型模型（ARMA）。以此可作为模型识别的标准。3 建模与仿真3.1确定线性模型的步骤（1）对一个时间序列作n次测量得到一个样本 ；（2）数据预处理：作 (为样本数据的算术平均值)，得到n个数据：；（3）计算样本自协方差函数，样本自相关函数，偏

6、相关函数数值，k=0,1,2,K；一般取K<n/4。（4）模型识别：按样本自相关函数和样本偏相关函数的值分别作出点图，按“截尾”，“拖尾”情况，查表确定模型的类别与阶数p，q。（5）参数估计：计算参数估计值（6）模型检验：模型能否正确地描述或比较好的反映所研究的变化过程的特性,还需要进行检验。如模型符合要求，就可以进行预测工作了，如果不符合要求，还需要进行修改或重新识别模型。检验的标准是：由实际观测到的样本序列值与经过模型识别和参数估计得到的随机模型计算出的估计值之差构成的残差序列：是不是白噪声的一个样本序列，若是，则所要检验的随机模型是合适的；否则即为不合适的随机模型。若，在置信度为9

7、5%下，认为序列是白噪声样本序列。（7）预测。3.2 电力负荷时间序列建模与仿真辽宁某电力公司2010年3月1日起连续52日用电量数据如下表： 表2，辽宁某电力公司2010年3月1日起连续52日用电量（单位：万kW·h）时间序号123456789用电量215221752209221822262228223422442253时间序号101112131415161718用电量225722712284228222812274229222932284时间序号192021222324252627用电量229821812223229722652291226722122206时间序号2829303

8电量225722382275229323322336231923352315时间序号373839404142434445用电量235923502355236524342360243624412447时间序号46474849505152用电量2438239723962386236223862378（1） 分析样本数据图1 电力负荷原始数据由上图可知，样本数据整体呈增长趋势，并且两端数据偏离平均值较大，可以认为样本数据是一个非平稳的时间序列样本。（2） 预处理图2 电力负荷一阶差分数据经过一阶差分后的数据如图2，基本可认为是平稳序列。图3 电力负荷一阶差分零均值数据再经过零

9、均值化后的数据如图3。（3） 计算协方差、自相关函数及偏相关函数、图4 自协方差函数图5 自相关函数和偏相关函数（4） 模型识别所谓模型识别就是判定这个要找的预测模型是AR、MA、ARMA中的哪一类。判定的依据是自相关函数及偏相关函数的特性。易见而当时，10个点均满足 ，在95%以上，则可基本认为截尾在处。通过自相关函数和偏相关函数的图像可以基本认为，经过一阶差分后的样本数据自相关函数拖尾，而偏自相关函数是截尾的，由此可以识别模型为AR（l）模型。（5） 参数估计模型AR(1)的参数估计，预测模型为：（6） 模型检验图6 自相关函数样本残差自相关函数计算如下：由计算可知：，可认为是高斯白噪声序

10、列，模型通过检验。（7） 预测结果运行程序预测结果为：，4 对非线性随机模型具体操作办法的思考在现实生活中更多的模型并不能很好的符合线性模型，这就需要引入非线性模型的概念，比如金融时间序列通常经历夹杂着动荡突发的漫长平静期。但由于非线性问题较为复杂，估计难度较大，所以，其本质处理手段是将非线性过程线性化处理。非线性估计领域的经典方法：广义卡尔曼滤波（Extended Kalman Filter，EKF）。EKF是将非线性模型在最优状态估计处进行泰勒展开，相对于其它算法的综合性能较好，算法简单，易于实现，是比较常用的非线性滤波方法。但EKF本质上是一种粗糙的方法，对于强非线性系统，EKF滤波性能

11、极不稳定，甚至发散。自EKF提出的近三十年间，众多学者尝试各种方法对非线性估计问题进行了深入的研究，提出了不少新方法、新思路，比如1997年Julier和Uhlman提出无迹卡尔曼滤波（Unscented Kalman Filter，UKF），其采用非线性变换（UT变换）代替传统的线性变换，体现了一种先进的思想，这种先进思想就是“非线性估计算法应更接近系统的非线性本质”，UKF已逐渐应用于导航、跟踪领域。就像前面所说的非线性模型的本质处理手段是线性化处理，即施加某种数学变换，变换成线性模型。而之前的一些专家学者已经提出了诸多非线性模型，我们建立非线性模型常采用的方法就是：在已提出的非线性模型中

12、挑选进行拟合，如果拟合效果好，则采用，不符合只能另寻出路。这是线性估计理论的基本思想。5 小结电力系统的负荷预测具有很重要的现实意义，学生通过运用在应用随机过程中学到的相关知识，利用MATLAB软件，针对辽宁省某电力公司的电力负荷时间序列建立了随机线性模型并进行了预测，由于掌握的资料广度和学生的知识水平的限制，只能根据本文的数据进行分析，模型还比较粗糙，该模型还可以通过其它的估计方法进行进一步的完善。另外，由于在数据搜索方面的能力不足，只找到了连续52个点的数据，如果数据长度再多一些的话，预测结果将更精确一些。参考文献1 田波平, 吴玉东, 张兴华. 应用随机过程M. 哈尔滨工业大学出版社.

13、2012.2 关薇. 时间序列在电力系统负荷预测中的应用D. 北方工业大学, 2010.3 万志宏. 基于时间序列的电力系统短期负荷预测研究D. 华南理工大学, 2012.4 朱琳. 电力系统负荷预测方法及其应用D. 华北电力大学 (北京), 2011.附录%功能：某电力公司电力负荷时间序列建模%编程时间：2013-11-05clc;clear;close all;data=2152,2175,2209,2218,2226,2228,2234,2244,2253,2257,2271,2284,2282,2281,.2274,2292,2293,2284,2298,2181,2223,2297,

14、2265,2291,2267,2212,2206,2257,.2238,2275,2293,2332,2336,2319,2335,2315,2359,2350,2355,2365,2434,2360,.2436,2441,2447,2438,2397,2396,2386,2362,2386,2378;%辽宁某电力公司2010年3月1日起连续52日用电量% data=1928,1926,1915,1968,1953,1972,1967,1992,2009,2016,2009,2033,2013,2144,.% 2088,2124,2096,2051,2056,2092,2034,2111,21

15、03,2139,2120,2121,2147,2161,.% 2110,2108,2040,2063,2055,2067,2100,2113,2042,2052,2075,2080,2010,2083,.% 2128,2093,2061,2106,2077,2100,2019,2154,2141,2150;%2141,2150%做出原始数据图形x_axis=1:length(data);plot(x_axis,data,'-ob');grid on;title('原始数据');set(gca,'XLim',1 length(data);%X轴的数

16、据显示范围xlabel('时间');ylabel('原始数据');title('电力负荷量');%样本数据整体呈增长趋势,并且两端数据偏离平均值较大%因此，可以认为样本数据是一个非平稳的时间序列样本%计算一阶差分序列并将其零均值化%一阶差分序列for i=1:(length(data)-1) data_D(i)=data(i+1)-data(i);endfigure;plot(1:length(data_D),data_D,'-ob');xlabel('时间');ylabel('一阶差分数据');t

17、itle('对电力负荷量作一阶差分');grid on;%零均值化data_D_average=mean(data_D);data_D_zero=data_D-data_D_average;figure;plot(1:length(data_D_zero),data_D_zero,'-ob');xlabel('时间');ylabel('零均值电力负荷量');title('零均值电力负荷量');grid on;%这样经过处理之后，data_D_zero已经是均值为零的平稳时间序列%求data_D_zero的平均值da

18、ta_D_zero_avedata_D_zero_ave=mean(data_D_zero);%计算零均值化序列协方差，自相关函数及偏相关函数%一般取K<n/4,常用K在n/10左右%我们取K=10%注意：这里以后会调整K=10;%先求自协方差函数rfor i=1:K r_temp=0; for j=1:length(data_D_zero)-i r_temp=r_temp+data_D_zero(j)*data_D_zero(j+i)/(length(data_D_zero); end r(i)=r_temp;endr0=0;for j=1:length(data_D_zero) r0

19、=r0+data_D_zero(j)*data_D_zero(j);endr0=r0/length(data_D_zero);r=r0 r;figure;plot(0:length(r)-1,r,'-ob');title('自协方差函数r(k)');xlabel('k');ylabel('自协方差函数r(k)');grid on;set(gca,'XLim',0 length(r)-1);%X轴的数据显示范围%自相关函数pp=r/r0;figure;subplot(2,1,1);plot(0:length(p)-1

20、,p,'-ob');title('自相关函数p');xlabel('k');ylabel('偏相关函数pk)的值');grid on;set(gca,'XLim',0 length(p)-1);%X轴的数据显示范围%计算偏相关函数m(1,1)=p(2);for k=1:K-1 temp1=0; temp2=0; for j=1:k temp1=temp1+p(k+1-j)*m(k,j); temp2=temp2+p(j)*m(k,j); end m(k+1,k+1)=(p(k+1)-temp1)./(1-temp2

21、); for j=1:k m(k+1,j)=m(k,j)-m(k+1,k+1)*m(k,k-j+1); endendmm=(diag(m)'mm=1 mm;i=1:length(m);subplot(2,1,2);plot(0:length(mm)-1,mm,'-ob');xlabel('k');ylabel('偏相关函数m(k,k)的值');title('偏相关函数m(k,k)');disp('偏相关函数m(k,k)为');mmset(gca,'XLim',0 length(mm)-1);%X轴的数据显示范围grid on;%以上全部算完%转入随机线性模型的参数估计%参数估计%确定为AR,MA,ARMA模型%基本可确定为AR(1)B=p(2);sita2=r0*(1-p(2)2);%模型检验for i=1:50 a(i)=data_D_zero(i+1)-B*data_D_zero(i);end%先求自协方差函数rfor i