Harbin Institute of Technology课程设计（论文）课程名称：应用随机过程设计题目：建模院系：电子与信息工程学院班级：通信1班设计者：学号：指导教师：设计时间： 2013-11-9哈尔滨工业大学线性模型——电力负荷时间序列建模1 电力系统负荷预测的意义随着我国电力事业的发展，电网的管理日趋现代化，对电力系统负荷预测问题的研究也越来越引起人们的注意。

电力负荷预测是电力系统调度、用电、计划、规划等管理部门的重要工作之一。

提高负荷预测技术水平，有利于计划用电管理，有利于合理安排电网运行方式和机组检修计划，有利于节煤、节油和降低发电成本，有利于制定合理的电源建设规划，有利于提高电力系统的经济效益和社会效益。

电力负荷预测，为编制电力规划提供依据，是电网规划的基础，它规定了电力工业的发展水平、发展速度、源动力资源的需求量，电力工业发展的资金需求量，以及电力工业发展对人力资源的需求量。

因此，国内外许多专家和学者开始致力于现代负荷预测方法的研究，而时间序列模型在国际和国内的电力系统短期负荷预测中得到了广泛应用。

2 平稳时间序列及其随机线性模型时间序列是指随时间改变而随机的变化的序列。

时间序列分析分为时域分析和频域分析，前者是对时间序列在时间域上的各种平均值进行分析研究，后者是进行傅里叶变换以后在频率域进行谱分析。

随着计算机技术的飞速发展，时域分析方法为人们所关注。

本文所要研究的就是时域分析。

平稳时间序列是平稳序列，它满足期望为0，且任意两个时刻的相关函数与时间t 无关，仅与两个时刻的时间差相关。

因为我们所掌握的为平稳时间序列的线性随机模型，而在实际中所遇到的一般都不是平稳时间序列，这就要对其进行相关的处理，使其变化为平稳序列。

均值为0且具有有理谱密度的平稳时间序列必可表示为下面三种形式中的一种（其中{,0,1,2,}t a t =±±L 为白噪声）： （1）自回归模型——AR 模型1122,0,1,2,t t t p t p t a t ωφωφωφω-------==±±L L AR （p ）模型由p ＋2参数来刻画; （2）滑动平均模型——MA 模型1122,0,1,2,t t t t q t q a a a a t ωθθθ---=---=±±L L MA(q)模型由q ＋2参数刻画;（3）自回归滑动平均模型或混合模型——ARMA 模型11221122,0,1,2,,0,1,2,t t t p t p t t t q t q a a a a t t ωφωφωφωθθθ----------=---=±±=±±L L L LARMA(p,q)混和模型由p ＋q ＋3参数刻画；通过以上介绍可以看出我们可以把AR(p)和MA(q)模型看成APMA(p,q)的两种特例。

线性模型中有两个重要的参数：自相关函数k ρ和和偏相关函数kk φ。

其中偏相关函数kk φ刻画了平稳序列任意一个长1k +的片段在中间量固定的条件下，两端的线性密切程度，而自相关函数k ρ也是刻画两端的线性密切程度，但并不需要中间数值固定。

与这两个参数相关的性质为拖尾性和截尾性。

拖尾性是指k ρ（或kk φ）随k 无限增长以负指数的速度趋向于0，其图像像拖一条尾巴；截尾性k ρ（或kk φ）是指参数在k>p 或k>q 后，其值变为零，其图像像截断了的尾巴一样。

总结上面的分析，得出时间序列模型的基本特性，如表1。

表1 时间序列模型的基本特性类型 ()AR p ()MA q (),ARMA p q 基本方程 ()t t B y a φ=()t ty B a θ=()()t t B y B a φθ=平稳条件 ()0B φ=的根全在单位圆外 无平稳条件()0B φ=的根全在单位圆外可逆条件 无可逆条件()0B θ=根全在单位圆外()0B θ=的根全在单位圆外自相关函数 拖尾 截尾 拖尾 偏相关函数截尾拖尾拖尾由模型的基本特性可知，自相关函数具有拖尾性，若算出的偏相关函数具有截尾性，模型为自回归模型（AR ）；当算出的自相关函数具有截尾性，偏相关函数具有拖尾性，则能判断模型为滑动平均模型（MA ）；若自相关函数和偏相关函数都具有拖尾性，模型为自回归滑动平均模型模型（ARMA ）。

以此可作为模型识别的标准。

3 建模与仿真3.1确定线性模型的步骤（1）对一个时间序列作n 次测量得到一个样本 123,,,,n Z Z Z Z L ；（2）数据预处理：作 t t Z Z ω=- (11ni i Z Z n ==∑为样本数据的算术平均值)，得到n 个数据：123,,,,n ωωωωL ；（3）计算样本自协方差函数ˆk r，样本自相关函数ˆk ρ，偏相关函数ˆkkφ数值，k=0,1,2,…,K ；一般取K<n/4。

（4）模型识别：按样本自相关函数和样本偏相关函数的值分别作出点图，按“截尾”，“拖尾”情况，查表确定模型的类别与阶数p ，q 。

（5）参数估计：计算参数估计值（6）模型检验：模型能否正确地描述或比较好的反映所研究的变化过程的特性,还需要进行检验。

如模型符合要求，就可以进行预测工作了，如果不符合要求，还需要进行修改或重新识别模型。

检验的标准是：由实际观测到的样本序列值t y 与经过模型识别和参数估计得到的随机模型计算出的估计值ˆt y之差构成的残差序列： ˆˆt t t ay y =- 是不是白噪声的一个样本序列，若是，则所要检验的随机模型是合适的；否则即为不合适的随机模型。

若ˆkρ≤（a ），在置信度为95%下，认为序列是白噪声样本序列。

（7）预测。

3.2 电力负荷时间序列建模与仿真辽宁某电力公司2010年3月1日起连续52日用电量数据如下表：（1） 分析样本数据图1 电力负荷原始数据由上图可知，样本数据整体呈增长趋势，并且两端数据偏离平均值较大，可以认为样本数据是一个非平稳的时间序列样本。

（2） 预处理电力负荷量时间原始数据图2 电力负荷一阶差分数据经过一阶差分后的数据如图2，基本可认为是平稳序列。

图3 电力负荷一阶差分零均值数据再经过零均值化后的数据如图3。

（3） 计算协方差、自相关函数及偏相关函数时间一阶差分数据时间零均值电力负荷量零均值电力负荷量、图4 自协方差函数图5 自相关函数和偏相关函数（4） 模型识别所谓模型识别就是判定这个要找的预测模型是AR 、MA 、ARMA 中的哪一类。

判定的依据是自相关函数及偏相关函数的特性。

易见而当1k >时，10个点均满足$0.2773kk φ<==，在95%以上，自协方差函数r(k)k自协方差函数r (k)k偏相关函数p k )的值k偏相关函数m (k ,k )的值则可基本认为kk φ截尾在1k p ==处。

通过自相关函数和偏相关函数的图像可以基本认为，经过一阶差分后的样本数据自相关函数拖尾，而偏自相关函数是截尾的，由此可以识别模型为AR （l ）模型。

（5） 参数估计模型AR(1)的参数估计，预测模型为：11,0,1,2,t t t a t ωφω--==±±L11ˆˆ0.2951φρ==-2201ˆˆˆ(1)1047.9a σγρ=-=（6） 模型检验图6 自相关函数样本残差自相关函数计算如下：由计算可知：ˆk ρ≤（a ），可认为t a 是高斯白噪声序列，模型通过检验。

（7） 预测结果运行程序预测结果为：53ˆ2386.2Z =，53ˆ2394.5Z =4 对非线性随机模型具体操作办法的思考在现实生活中更多的模型并不能很好的符合线性模型，这就需要引入非线性模型的概念，比如金融时间序列通常经历夹杂着动荡突发的漫长平静期。

但由于非线性问题较为复杂，估计难度较大，所以，其本质处理手段是将非线性过程线自相关函数p a自相关函数性化处理。

非线性估计领域的经典方法：广义卡尔曼滤波（Extended Kalman Filter，EKF）。

EKF是将非线性模型在最优状态估计处进行泰勒展开，相对于其它算法的综合性能较好，算法简单，易于实现，是比较常用的非线性滤波方法。

但EKF本质上是一种粗糙的方法，对于强非线性系统，EKF滤波性能极不稳定，甚至发散。

自EKF提出的近三十年间，众多学者尝试各种方法对非线性估计问题进行了深入的研究，提出了不少新方法、新思路，比如1997年Julier和Uhlman提出无迹卡尔曼滤波（Unscented Kalman Filter，UKF），其采用非线性变换（UT变换）代替传统的线性变换，体现了一种先进的思想，这种先进思想就是“非线性估计算法应更接近系统的非线性本质”，UKF已逐渐应用于导航、跟踪领域。

就像前面所说的非线性模型的本质处理手段是线性化处理，即施加某种数学变换，变换成线性模型。

而之前的一些专家学者已经提出了诸多非线性模型，我们建立非线性模型常采用的方法就是：在已提出的非线性模型中挑选进行拟合，如果拟合效果好，则采用，不符合只能另寻出路。

这是线性估计理论的基本思想。

5 小结电力系统的负荷预测具有很重要的现实意义，学生通过运用在应用随机过程中学到的相关知识，利用MATLAB软件，针对辽宁省某电力公司的电力负荷时间序列建立了随机线性模型并进行了预测，由于掌握的资料广度和学生的知识水平的限制，只能根据本文的数据进行分析，模型还比较粗糙，该模型还可以通过其它的估计方法进行进一步的完善。

另外，由于在数据搜索方面的能力不足，只找到了连续52个点的数据，如果数据长度再多一些的话，预测结果将更精确一些。

参考文献1田波平, 吴玉东, 张兴华. 应用随机过程[M]. 哈尔滨工业大学出版社. 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,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,2103,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轴的数据显示范围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('一阶差分数据');title('对电力负荷量作一阶差分');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的平均值data_D_zero_avedata_D_zero_ave=mean(data_D_zero);%%计算零均值化序列协方差，自相关函数及偏相关函数%%一般取K<n/4,常用K在n/10左右%%我们取K=10%%注意：这里以后会调整K=10;%%先求自协方差函数rfor i=1:Kr_temp=0;for j=1:length(data_D_zero)-ir_temp=r_temp+data_D_zero(j)*data_D_zero(j+i)/(length(data_D_zero));endr(i)=r_temp;endr0=0;for j=1:length(data_D_zero)r0=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],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-1temp1=0;temp2=0;for j=1:ktemp1=temp1+p(k+1-j)*m(k,j);temp2=temp2+p(j)*m(k,j);endm(k+1,k+1)=(p(k+1)-temp1)./(1-temp2);for j=1:km(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:50a(i)=data_D_zero(i+1)-B*data_D_zero(i);end%%先求自协方差函数rfor i=1:Kr_a_temp=0;for j=1:length(a)-ir_a_temp=r_a_temp+a(j)*a(j+i)/(length(a));endr_a(i)=r_a_temp;endr0=0;for j=1:length(a)r0=r0+a(j)*a(j);endr0=r0/length(a);figure;plot([1:length(r_a)],r_a,'-ob');title('自协方差函数r_a');ylabel('自协方差');grid on;set(gca,'XLim',[1 length(r_a)]);%X轴的数据显示范围%%%自相关函数pp_a=r_a/r0;figure;plot([1:length(p_a)],p_a,'-ob');title('自相关函数p_a');ylabel('自相关函数');grid on;set(gca,'XLim',[1 length(p_a)]);%X轴的数据显示范围%%%得到的p_a均小于1.95/sqrt(50),这样就可以证明预测模型通过检验disp('自相关函数p_a为');p_a%%%%模型检验完成%%%预测第53个点，第54个点L=length(data_D_zero);for i=1:2data_D_zero(L+i)=B*data_D_zero(L);data_D(L+i)=data_D_zero(L+i)+data_D_average;data(52+i)=data(52+i-1)+data_D(L+i);end。