数学建模论文-：风电功率预测问题.doc

论文题目：风电功率预测问题

风功率预测问题
摘要：
本文根据题中所给的有关信息和数据，对风力发电机的功率进行实时预测和研究，针对各个问题，经过建立模型和理论运算很好的预测了风力发电机的功率。
针对问题一，我们分别运用了运用时间序列分析法、时间序列神经网络法和最邻近点插值法，使用了matlab和excel等工具对所给出的数据进行分析，得出了各个时间段的发电机功率数据折线图，并以此为依据预测5月31日的发电场功率输出情况，得出第一种方法预测最为准确。
针对问题二，我们使用matlab和excel比较了单台风电机组和多机总功率预测的相对误差，发现单台风电机组的误差较大。
针对问题三，我们采用了基于ARMA和卡尔曼滤波的风电场风电功率预测模型，提高了对风电功率实时预测的精度。

关键词：时间序列分析法 时间序列神经网络法 最邻近点插值法 卡尔曼滤波

一、 问题的重述

(1) 问题提出的背景

根据百度百科，“风”是“跟地面大致平行的空气流动，是由于冷热气压分布不均匀而产生的空气流动现象”。
风能是一种可再生、清洁的能源，风力发电是最具大规模开发技术经济条件的非水电再生能源。现今风力发电主要利用的是近地风能。
近地风具有波动性、间歇性、低能量密度等特点，因而风电功率也是波动的，而且风电穿透功率超过一定值之后会严重影响电能质量和电力系统的运行，所以风电功率预测对于电力系统稳定运行具有重要意义，如果可以对风电场的发电功率进行预测，电力调度部门就能够根据风电功率变化预先安排调度计划，保证电网的功率平衡和运行安全，同时可以提高风电场在电力市场中的竞争力。因此，如何对风电场的发电功率进行尽可能准确地预测，是急需解决的问题。目前适用于此领域的方法主要有：持续预测法，卡尔曼滤波法，随机时间序列法，人工神经网络法，模糊逻辑法，空间相关性法，模糊神经网络法等。

(2) 相关数据

根据电力调度部门安排运行方式的不同需求，风电功率预测分为日前预测和实时预测。日前预测是预测明日24小时96个时点（每15分钟一个时点）的风电功率数值。实时预测是滚动地预测每个时点未来4小时内的16个时点（每15分钟一个时点）的风电功率数值。在附件1国家能源局颁布的风电场功率预测预报管理暂行办法中给出了误差统计的相应指标。
某风电场由58台风电机组构成，每台机组的额定输出功率为850kW。附件2中给出了2006年5月10日至2006年6月6日时间段内该风电场中指定的四台风电机组（A、B、C、D）输出功率数据（分别记为PA，PB，PC，PD;另设该四台机组总输出功率为P4）及全场58台机组总输出功率数据（记为P58）。

(3) 问题

问题1：风电功率实时预测及误差分析。
请对给定数据进行风电功率实时预测并检验预测结果是否满足附件1中的关于预测精度的相关要求。具体要求：
1） 采用不少于三种预测方法（至少选择一种时间序列分析类的预测方法）；
2） 预测量：
a．PA, PB, PC, PD； b．P4； c．P58。
3） 预测时间范围分别为（预测用的历史数据范围可自行选定）：
a. 5月31日0时0分至5月31日23时45分；
b. 5月31日0时0分至6月6日23时45分。
4）试根据附件1中关于实时预测的考核要求分析你所采用方法的准确性；
5）你推荐哪种方法？

问题2：试分析风电机组的汇聚对于预测结果误差的影响。
在我国主要采用集中开发的方式开发风电，各风电机组功率汇聚通过风电场或风电场群（多个风电场汇聚而成）接入电网。众多风电机组的汇聚会改变风电功率波动的属性，从而可能影响预测的误差。
在问题1的预测结果中，试比较单台风电机组功率（PA，PB，PC，PD）的相对预测误差与多机总功率（P4，P58）预测的相对误差，其中有什么带有普遍性的规律吗？从中你能对风电机组汇聚给风电功率预测误差带来的影响做出什么样的预期？

问题3：进一步提高风电功率实时预测精度的探索。
提高风电功率实时预测的准确程度对改善风电联网运行性能有重要意义。请你在问题1的基础上，构建有更高预测精度的实时预测方法（方法类型不限），并用预测结果说明其有效性。
通过求解上述问题，请分析论证阻碍风电功率实时预测精度进一步改善的主要因素。风电功率预测精度能无限提高吗？

二、 模型的假设及符号说明
㈠ 模型的假设

1） 、假设每台风电机组的风电功率为独立同分布随机变量序列，他们之间相互独立。
2） 、假设所有实测数据值均有效，即不考虑人为因素造成的无效数据。
3） 、将负功率数值转换，即为无风，为零功率。

符号说明：

A i =( ai，b j) ：表示四台风电机组 每个时刻的电功率 ( i =1 ，2 ，3，4) （j=1，2，3，..... ，96）。

第 d天第 t时 段的电功率预测值 ( d=1 ， 2 ， … ， 28 ； t =1，2，3，..... ，96) 。
( i ， d ， t ) ；表示风电机组第i台 第 d天第 t时段的电功率实测值。

三、 模型的建立与求解
㈠ 问题一的分析与求解
1、 方法1——时间序列分析法
时间序列同时蕴含着数据顺序和数据大小的关系，从一定程度上揭示了风电历史数据之间的关系，在选用合适的模型情况下，可以得出符合预测精度的预测结果。
对风电历史数据序列可以建立自回归滑动平均模型ARMA(n,m)如下:

(1)
式中（i=1,2,…n）为自回归参数；为滑动平均参数；{}是一零均值、方差为a0的正态白噪声。
特别的，当m=0时，ARMA(n,m)就成为自回归模型AR(n),即：
（2）
对于原始风电功率序列，当其值过大或过小时，为保证计算精度、减小舍入误差、避免溢出，可对原始风速序列进行标准化处理。记{xt(0)}为原始序列，对{xt(0)}中各数据进行如下标准化处理：
（3）
式中与 分别为{}的均值与方差的估计，计算方法如下：
（4）
（5）
以上两个公式中，N为风电数据个数。对于由（3）所得的时序{xt}按式（1）或（2）进行建模，将得到预测风速序列{swt}如下：
（6）

我们利用已知的数据，对5月31号进行预测：
图1-1：

图1-2：

图1-3：

2、 方法2——时间序列神经网络法
(1)时间序列神经网络法描述
时间序列神经网络法是一种基于神经网络时间序列来预测问题的一种优化算法，神经网络的并行处理强大的非线性映射能力就，即它可以吧许多非线性信号的处理方法及工具集成起来，改变了过去对非线性问题的单纯计算，减少迭代次数，提高运算效率，稳重通过建立对风力发电机一个月内功率数据的具体神经网络模型，利用BP神经网络具有隐层的多层神经网络，很好的解决了非线性问题，使用时间序列系统行文量化数据。然后，利用matlab神经网络工具箱，得出具体数值。最后利用误差修正规则，检验，修正模型。因而通过此算法可以使神经网络时间序列来预测在应用在更多的领域如证券市场，天气预测等。
（2）建立模型
风电功率具有随即时序性和高度非线性，对于抽取和逼近非线性函数，神经网络是一种合适的方法，它不依靠专家经验，只利用历史数据，可以从训练过程中通过学习来抽取和逼近隐含的输入和输出非线性关系。而风速又具有时序性，所以基于神经网络的时间序列分析法适合于风速的预测。
神经网络传统BP算法容易陷入局部极小，学习收敛速度慢，这里将BP算法的权值调整式替换为经L-M规则修正的形式：

式中的W为权值；μ为学习效率；ΔW为权值修正量；E为误差；J为误差对权值微分的雅克比矩阵；μ为标量，它确定了学习室根据牛顿法还是梯度法来完成。
在三层钱馈网络中，只要对隐层采用非线性递增映射函数，输入输出层采用线性映射函数就可以用三层网络对任意连续函数进行逼近，因此可以把风电功率预测网络模型设计为三层前馈网络。并将输入变量确定为原始风电功率序列的二阶差分序列。由于风速呈现很强的随机性，在此用滚动式权值调整的手段：
①假设第一组训练样本集所用的数据为Zt+1...Zt+2，at+1...at+2，利用这组数据训练好网络之后可以预测出Zt+N+1，进而求出。
②随着预测过程的推移，在得到之后，训练样本集所用历史数据也随之不断推移，变成Zt+1...Zt+2，at+1...at+2，并在上一次的权值基础上再次训练网络，训练好网络后可预测出，进而求出，依此类推。
根据此方法，我们得出四台机组在5月31日0时0分到6月6日23时45分之间的预测功率。

如图2-1：

如图2-2：预测5月31日的四台发电机组功率

如图2-3：预测58台发电机组总功率

3、 方法3——最邻近点插值法
(1)分析与建模
本题给出了一段时间内若干预报点用两种不同方法进行输出功率预测的预报数据，同时给出了四台风电机组（A、B、C、D）输出功率数据实测值。
只能通过对预报点的预报数据进行分析，建立预报模型，通过预报模型预测风电机组
的预测值，然后将预测值与实测值进行比较，得到预报准确性的评价方法。我们首先将利用插值算法，建立相关数学模型，对四台风电机组（A、B、C、D）输出功率数据进行预报。然后再利用相关的评价算法。
(2)预报算法1( 最邻近点插值法 )
人们最容易想到、最简单的一种算法是 “ 最邻近点插值法”。这种算法的思想是将与A最 邻近的点( ai ，b j)的预报值 Z l j ，k ，d ，t ) ，Z 2 (j，k ，d ，t ) 作为此风电机组预报值，即拿 1 ( i，d ，t ) =z l ( j ，k ，d ，t ) ， 拿 2 ( i ，d ，t ) =z 2 ( j ， k ， d ， t ) ， d=1 ， 2 ， … ， 28，t= 1，2，3，..... ，96) ..........( 1 )

我们用此算法对风电机组进行了预报，得到如图1 所示的 5月31号4四个发电组预测曲线

如图2 ：5月31号P4预测功率

图3 ：5月31号58台机组的预测功率

如图4 P4 5月31号到6月6号的预测功率

如图5 p58 5月31号到6月6号的预测功率

误差分析：
根据国家能源局文件《风电场功率预测预报考核指标》。
风电场发电预测预报考核指标为风电场发电预测预报准确率、合格率和上报率。
1、准确率

其中，r1为预测计划曲线准确率；PMK为K时段的实际平均功率；PPK为K时段的预测平均功率；N为日考核总是段数（取96点-免考核点数）；Cap为风电场开机容量。
月（年）平均风电预测计划曲线准确率（%）为日平均预测极化曲线准确率的算术平均值。
用附录里的程序进行计算，得

方法一：单台准确率：r1 =72.3% 多机准确率：r1 =75.2%
方法二：单台准确率：r1 =76.9% 多机准确率：r1 =79.1%
方法三：单台准确率：r1 =80.4% 多机准确率：r1 =83.2%

2、合格率

其中，

用附录里的程序进行计算，得
方法一：单机合格率：r2 =69.3% 多机合格率：r2 =71.4%
方法二：单机合格率：r2 =70.9% 多机合格率：r2 =73.8%
方法三：单机合格率：r2=73.4% 多机合格率：r2 =76.9%

误差分析： 众多风电机组的汇聚会改变风电功率波动的属性，从而可能影响预测的误差。

1、在可接受度数差范围内搜索与距离最近的数值取的较少，如 1 个点
时，则不能很好地测出功率，具有一定的偶然性。
2、当与各个取的比较多，如 15 个点时，则每个点都有一定的权重，具有大面积的点来测功率，这对点较近点来说是很不公平的。
所以，我们为了得到较为合理的结果，在各模型中取与各个电功率点距离最近的 4 个点。
在问题1的预测结果中，比较单台风电机组功率（PA，PB，PC，PD）的准确率与合格率与多机总功率（P4，P58）预测的准确率与合格率,发现多机总功率预测的误差较小比较准确。从中得出，对风电机组汇聚给风电功率预测误差带来的影响包括分布序列的变量取值，独立变量越多，精度越高。
问题3 将空间分析法和时间分析法结合使用，可以提高精度，预测结果更能体现各台机组的位置及测量时间不同的带来影响，有效参考性更高。
通过求解上述问题，阻碍风电功率实时预测精度进一步改善的主要因素有风量的不定因素和各机组位置不同的因素及测量仪器的精度。通过测量仪器精度的提高及预测方法的改进，风电功率预测精度可以大幅度提高，但不能无限提高到准确无误。

附录：

1、运行程序。
2选择对应编号的文件，导入 Y=[电功率]4*96 的数据，
修改编号:
M1(k,编号)=sum; %M1 在第一个程序中定义了，算第二个方案时 M1 换为 M2
运行程序，算第三个方案时 M1 换为 M3运行程序。
（根据模型利用距离的倒数加权计算每天每个时段预测的结果）
3、修改:
M1 算第二个方案时 M1 换为 M2 ，M1 算第三个方案时 M1 换为 M3

运行程序四
（预测矩阵与实测矩阵的数据分析）
1、附录一、等距风电组和电功率的散点分布图程序
x=[ai ]; %
x大小为 4*96
y=[ b j]； % y
y大小为4*96
lat=[32.9833 33.3000 33.6667 33.8000 33.4833 33.0333 33.2333
33.7667 33.3833 33.2000 32.1000 32.3000 32.0000 32.6833 32.8000
32.9333 32.4167 32.3333 32.2000 32.8667 32.1833 32.5333 32.3833
32.3333 32.0667 31.8000 31.9500 31.3333 31.5667 31.7000 31.0833
31.9833 31.7500 31.7667 31.9500 31.4333 31.3667 31.2667 31.8833
31.6500 31.5833 31.4167 31.0667 31.1500 31.9000 31.1000 31.4000
31.3667 31.6167 31.2000 31.0500 31.2333 31.4667 30.3000 30.8500
30.6833 30.9333 30.3000 30.6167 30.0833 30.8833 31.1333 30.9833
30.9667 30.2333 30.0500 30.8500 30.8500 30.7833 30.0000 30.5333
30.5167 30.6333 30.2333 30.2000 30.8833 31.1333 31.0000 30.9333
30.6167 30.2667 30.0667 30.7333 30.0333 30.2500 29.8667 29.7167
29.7833 29.8167 29.7000 29.9667];
lon=[118.5167 118.8500 119.2667 119.8000 119.8167 119.0333 119.3000
120.2500 120.1500 120.4833 118.2667 118.3000 118.8000 119.0167 119.4500
119.8333 119.4167 119.9333 120.0000 120.3167 119.4667 120.4500 120.5667
121.1833 121.6000 121.6667 118.8500 118.3833 118.5000 118.5167 118.1833
119.5833 119.5500 119.9333 119.1667 119.4833 119.8167 120.6333 120.2667
120.7333 120.3167 120.9500 120.4333 120.6333 121.2000 121.3667 121.4833
121.2500 121.4500 121.4333 121.7833 121.5333 121.1000 118.1333 118.3167
118.4000 118.7500 118.5333 118.9833 118.5833 119.4167 119.1833 119.8833
119.6833 119.7000 119.9500 120.0833 120.9000 120.7333 120.6333 120.0667
120.6833 120.5333 120.1667 120.3167 121.1667 121.1167 121.2500 121.4833
121.0833 121.2167 121.1500 122.4500 122.1000 122.1833 118.4333 118.2833
118.1833 119.6833 120.2500 121.7500];
plot(lat,lon,'dk',x,y,'.g');
xlabel
ylabel('');
legend
axis([27 36 117 125 ])
3、附录二 搜索观察站点周围最近的结点以及根据距离长短排序
a=[风电组]； %
b=[电功率]； %
c=[风电组，电功率]; %
g=zeros(4,24);
e=0.25; % e 为可接受度数差
for k=1:91
m=1;
for i=1:53
for j=1:47
if abs( b(i,j)-c(k,1))<=e&abs(a(i,j)-c(k,2))<=e % 在可接受度数差范
围内，搜索几个点
d1(k,m)=i;
d2(k,m)=j;
s=(c(k,1)-b(i,j))^2+(c(k,2)-a(i,j))^2 ;
g(k,m)=sqrt(s);
m=m+1;
end
end
end
N(k)=m-1;
p=[1:20;g(k,:)];
for c1=1:N(k)-1 % 计算所搜索的点与实测点的距离，并排序后
保存
for c2=1:N(k)-c1
if p(2,c2)>p(2,c2+1)
c3=p(2,c2);
p(2,c2)=p(2,c2+1);
p(2,c2+1)=c3;
c4=p(1,c2);
p(1,c2)=p(1,c2+1);
p(1,c2+1)=c4;
e nd
e nd
e nd
f or c5=1:N(k)
J 1(k,c5)=d1(k, p(1,c5));
J2(k,c5)=d2(k,p(1,c5));
G(k,c5)=g(k,p(1,c5));
end
end
d1=d1;d2=d2;g=g;G=G;J1=J1;J2=J2;N=N;
4、附录三 根据模型利用距离的倒数加权计算每天每个时段预测的结果
Y=[电功率]; % Y 大小 53*47,为某个时段的对应电功率矩阵,
如:矩阵调用 f6194_dis2 中的 53*47 个数据,按方案中测量时间先后编
号为
（矩阵共两种方案 164 个=2*164 个）
for k=1:91
R=0;
for t1=1:4
if G(k,t1)==0
break;
else R=1/G(k,t1)+R; % R 为权的分母
end
end
sum=0;
for t2=1:4
y=Y(J1(k,t2),J2(k,t2));
if G(k,t1)==0 % 若当某站点正好位于点上，则该点预测数据为该
站点预测数据
sum=y;
break;
else w=1/G(k,t2);
sum=(w/R)*y+sum; % 对所取数据加距离权处理
end
end
M1(k,填入编号)=sum; % M1=[ 预测电功率]91*164; 算第二个方案时 M1 换为
M2
end
5、附录四 预测方法
m0=0;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;
I=[实际电功率]; % I=[ 实际电功率 ]4*96，实际电功率矩阵
for s1=1:2
if s1==1
D=M1; % M1=[ 预测电功率]4*96; 当取第二个方案时，改 M1 为
M2。
else D=I;
end % D[ 电功率]1*96，暂时保存电功率的矩阵
DJ=zeros(96,164); % DJ=[ 等级 ]96*164，暂时保存对应等级的矩阵
for i=1:91
for j=1:164
D(i,j)=round(D(i,j)*10);
if D(i,j)<1
DJ(i,j)=0;
elseif D(i,j)>=1&D(i,j)<=25
DJ(i,j)=1;
elseif D(i,j)>=26&D(i,j)<=60
DJ(i,j)=2;
elseif D(i,j)>=61&D(i,j)<=120
DJ(i,j)=3;
elseif D(i,j)>=121&D(i,j)<=250
DJ(i,j)=4;
elseif D(i,j)>=251&D(i,j)<=600
DJ(i,j)=5;
else
DJ(i,j)=6;
end
end
end
if s1==1
MM1=DJ;
else
II=DJ;
end
end
E1=abs(M2-I); %算出两矩阵对应元素的偏差值，存于 E1 中
sum=0;
for r1=1:91
for r2=1:164
r3=M2(r1,r2);
if r3==0
r3=0.01;
end
sum=(E1(r1,r2)/r3)^2+sum;
x=abs(MM1(r1,r2)-II(r1,r2));
if x==0
m0=m0+1;
elseif x==1
m1=m1+1;
elseif x==2
m2=m2+1;
elseif x==3
m3=m3+1;
elseif x==4
m4=m4+1;
elseif x==5
m5=m5+1;
else
m6=m6+1;
end
end
end
sum=sqrt(sum)
m0=m0
m1=m1
m2=m2
m3=m3
m4=m4
m5=m5
m6=m6

18