数学建模实验五：风电功率预测问题

一、实验目的

（1）学习和掌握马尔科夫链预测方法； （2）掌握三次指数平滑预测方法； （3）掌握BP神经网络预测方法；

（4）掌握NAR时间序列的动态神经网络进行预测等

二、实验步骤

2.1 问题提出：

某风电场由58台风电机组构成，每台机组的额定输出功率为85kw。附件2中给出了2006年5月10日~2006年6月6日时间段内该风电场中指定的四台风电机组（A，B，C和D）每隔15分钟的输出功率数据（分别记为PA、PB、PC和PD；另设该四台机组总输出功率为P4）及全场58台机组总输出功率数据（记为P58）。 问题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的基础上，构建有更高预测精度的实时预测方法，并用预测结果说明其有效性。

通过求解上述问题，请分析论证阻碍风电功率实时预测精度进一步改善的主要因素。 风电功率预测精度能无限提高吗？

2.2 问题1的分析、建模与求解

（1）数据分析：绘制各发电机输出功率曲线（参考程序ysw.m）,分析数据的波动规律。 程序代码：

clc,clear,close all load('data.mat') a = 23; % 某一天

subplot(321),plot(PA(a,:));grid on xlabel('时点');ylabel('PA输出功率'); title('PA发动机输出功率曲线')

subplot(322),plot(PB(a,:));grid on xlabel('时点');ylabel('PB输出功率'); title('P5发动机输出功率曲线')

subplot(323),plot(PC(a,:));grid on xlabel('时点');ylabel('PC输出功率'); title('PC发动机输出功率曲线')

subplot(324),plot(PD(a,:));grid on xlabel('时点');ylabel('PD输出功率'); title('PD发动机输出功率曲线')

subplot(325),plot(P4(a,:));grid on xlabel('时点');ylabel('P4输出功率'); title('P4发动机输出功率曲线')

subplot(326),plot(P58(a,:));grid on xlabel('时点');ylabel('P58输出功率'); title('P58发动机输出功率曲线')

输出：

图1

分析：如图1所示，各发电机组的数据波动规律基本一致，具有周期性，并没有较强的规律性，具有随机波动性，在0-100的时点内，一开始从一个特定的值逐步上升，到达一个特定时间点（大概是12时点）开始迅速下降，一直降到0，在很长一段时间内都在0的附近微弱波动，从大概40时点开始，在动荡中上升，大概在70时点是上升到最大值，之后又迅速下降，再次将到0，从90时点开始又有小幅度上升。

（2）建立三次指数平滑法模型并求解（参考问题1程序） 建立模型：

指数平滑预测方法是移动平均预测方法加以发展的一种特殊加权移动平均预测方法。它可分为一次指数平滑法和多次指数平滑法，一般常用于时间序列数据资料既有长期趋势变动又有季节波动的场合。根据平滑次数的不同，又分为一次指数平滑法、二次指数平滑法和三次指数平滑法等。

一次指数平滑法是以最后一次指数平滑值为基础，确定市场预测值的一种特殊的加权平均法，其基本预测方程为：

1)St(1)xt(1)St(1

Ft1St(1)xt(1)Ft式中

xt——实际观察值，t1，2，，nSt(1)——时间t观察值的一次指数平滑值

——时间序列的平滑指数，01Ft1——t11)St(1)xt(1)St(1Ft1St(1)xt(1)Ftxt——实际观察值，t1，2，，nSt(1)——时间t观察值的一次指数平滑值——时间序列的平滑指数，01

由上式可知，一次指数平滑法是一种加权预测。指数平滑法是以首项系数为α，公比为（1一α）的等比数列作为权数的加权平均法。体现了“近重远轻”的赋权原则它既不需要存储全部历史数据，也不需要存储一组数据，从而可以大大减少数据存储问题，甚至有时只需一个最新观察值、最新预测值和α值，就可以进行预测。它提供的预测值是前一期预测值加上前期预测值中产生的误差的修正值。由预测模型可见，起到一个调节器的作用。如果值选取得越大，则越加大当前数据的比重，预测值受近期影响越大；如果值选取得越小，则越加大过去数据的比重，预测值受远期影响越大。因此，值大小的选取对预测的结果关系很大，一般根据最小均方差选取的值。

二次指数平滑法是在一次指数平滑的基础上再进行一次指数平滑。并根据一次、二次的最后一项的指数平滑值，建立直线趋势预测模型，并用之进行预测的方法，称之为二次指数平滑预测法。当时间序列的变动呈线性趋势时，可采用二次指数平滑法。二次指数平滑法的计算方法为：

其中

Ft1——t1期预测值二次指数平滑预测模型为： 其中

三次指数平滑法是二次指数平滑法的进一步推广，仿照二次指数平滑法的推导方法，可推得估计值公式：

ˆt3St(1)3St(2)St(3)aˆ265St(1)2(54)St43St(3)bt2(1)2ˆtS(1)2St(2)St(3)c2t21

所以，最终可得三次指数平滑预测模型为：

因此，三次指数平滑预测法的求解流程为：

ˆcyˆtaˆtbˆt2,1,2,...tˆt3St(1)3St(2)St(3)aˆ65St(1)2(54)St243St(3)bt2(1)2ˆtS(1)2St(2)St(3)c2t21(1)(1)Styt(1)St1S(2)S(1)(1)S(2)tt1t(3)(2)StSt(1)St(3)1

程序代码：

针对PA，设置平滑指数为0.9，对5月31日00:00-23:45时段内进行预测，MATLAB程序如下：

%PA 5月31:0:0 计算的MATLAB 程序如下：三次指数平滑预测 clc,clear

load('data.mat') %原始数据以列向量的方式存放在workplace文件中 PA=PA(2:29,:);PA=PA';

yt=PA(:,22); n=length(yt);

alpha=0.9; %平滑系数如果时间序列具有快速明显的变化时，则α宜选用较大的值 st1_0=mean(yt(1:3)); st2_0=st1_0;st3_0=st1_0; st1(1)=alpha*yt(1)+(1-alpha)*st1_0; st2(1)=alpha*st1(1)+(1-alpha)*st2_0; st3(1)=alpha*st2(1)+(1-alpha)*st3_0; for i=2:n

st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1); st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1); st3(i)=alpha*st2(i)+(1-alpha)*st3(i-1); end

xlswrite('PA531.xls',[st1',st2',st3'])

st1=[st1_0,st1];st2=[st2_0,st2];st3=[st3_0,st3]; a=3*st1-3*st2+st3;

b=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1-2*(5-4*alpha)*st2+(4-3*alpha)*st3); c=0.5*alpha^2/(1-alpha)^2*(st1-2*st2+st3); yhat=a+b+c; %预测yt+1的值

xlswrite('PA531.xls',yhat','Sheet1','D1') %求解三次指数平滑预测方程c、b、a系数 xishu=[c(n+1),b(n+1),a(n+1)]; %Sheet1.E1列存放真实值

xlswrite('PA531.xls',yt,'Sheet1','E1') %误差分析，预测值与真实值差之差 st3=st3(2:n+1)'; delta=abs(st3-yt);

xlswrite('PA531.xls',delta,'Sheet1','F1')%绝对误差 deltaxd=delta./yt;

xlswrite('PA531.xls',deltaxd,'Sheet1','G1')%计算相对误差 deltajdminmax=minmax(delta');

xlswrite('PA531.xls',deltajdminmax,'Sheet1','H1')%计算绝对误差的最小值，最大值 deltaxdminmax=minmax(deltaxd');

xlswrite('PA531.xls',deltaxdminmax,'Sheet1','J1')%计算相对误差的最小值，最大值 deltasum=sum(sum(delta));%总误差值

xlswrite('PA531.xls',deltasum,'Sheet1','L1') %函数绘图

plot(1:n,yt,1:n,st3(1:n),'r') legend('实际值','预测值',2) grid on

xlabel('时点x'),ylabel('发电功率y');

title('PA5.31.0.0-5.31.23.45发电功率随时点变化图像')

输出：

图2

对5月31日00:00-6月6日23：45整体上进行预测，MATLAB程序如下：

%PA 5月31:0:0 计算的MATLAB 程序如下：BP神经网络预测 %5月31位于第23行 clc,clear

load PA %原始数据以列向量的方式存放在workplace文件中 PA=PA(2:29,:); %数据的标准化 N=size(PA); for j=1:N(1,2)

PAHminmax=minmax(PA(:,j)'); for i=1:N(1,1)

PA(i,j)=(PA(i,j)-PAHminmax(1,1))/(PAHminmax(1,2)-PAHminmax(1,1)); end end

%以每天的从0时计数起，每隔十五分钟作为输入 P=PA(1:21,:);

%以5月31的间隔十五分钟的发电量作为目标向量 T=PA(22:N(1,1),:);

%创建一个BP神经网络，每一个输入向量的取值范围为[0 ,1]，隐含层有22个神经元， %输出层有一个神经元，隐含层的激活函数为tansig，输出层的激活函数为%logsig， %训练函数为梯度下降函数，即标准学习算法 for i=1:21 a(i,1)=0; a(i,2)=1; end

net=newff(a,[21,7],{'tansig','logsig'},'traingd'); net.trainParam.epochs=25000; net.trainParam.goal=0.01; %设置学习速率为0.1 LP.lr=0.1; %训练网络

net=train(net,P,T);

%预测5月31的发电量数据 T1=sim(net,P);%预测值

%PA.5月31日发电量真实值 T0=PA(22:N(1,1),:);

%预测值与实际值的误差 n=size(T1); k=1;

for i=1:n(1,1) for j=1:n(1,2)

a(k,1)=T1(i,j); b(k,1)=T0(i,j); k=k+1; end end

N=size(a); %绘制误差图

plot(1:N(1,1),b(:,1),1:N(1,1),a(:,1),'r') grid on

legend('实际值','预测值',2)

xlabel('时点x'),ylabel('发电功率y');

title('PA5.31.0.0-6.6.23.45发电功率实时函数图像')

输出：

图3

对于其他机组及总机组，仿照上述三次指数平滑模型的计算方法依次计算，利用MATALB编写程序，可得类似的预测结果。

分析：从图3、图4可看出，三次指数平滑法对风电机功率的非线性、非规律性的适应性较强，可以较好的预测功率的非线性和非规律性特点。但该方法权数的确定具有很强的主观性，当数据特征发生变化时，指数平滑法不能自动调整权数，以适应新数据的要求；同时，当预测对象保持较长时间的稳定后，出现突然上升或下降的趋势时，指数平滑法就难以适应。因此，指数平滑法应用于中短期预测时误差较小，效果较好。在本题中，采用一个时点预测下一个时点，不断的滚动向前，长远预测误差将很大。

（3）建立BP神经网络模型并求解（参考问题1程序） 建立模型：

ＢＰ神经网络是一种典型的多层前向型神经网络，具有一个输入层、一个或多个隐含层和一个输出层。层与层之间采用权连接的方式，同一层的神经元之间不存在相互连接。理论上已经证明，具有一个隐含层的三层网络可以逼近任意非线性函数。

BP网络的学习过程主要由以下四部分组成： 1）输入样本顺传播

输入样本传播也就是样本由输入层经中间层向输出层传播计算。这一过程主要是 输入样本求出它所对应的实际输出。

① 隐含层中第i个神经元的输出为

Ra1if1w1ijpjb1ij1② 输出层中第k个神经元的输出为：

i1,2,,s1

a2kS1f2wab2k,2ki1ii1i1,2,s2

其中f1(·), f2 (·)分别为隐含层和输出层的传递函数。

2）输出误差逆传播

在第一步的样本顺传播计算中我们得到了网络的实际输出值，当这些实际的输出值与期望输出值不一样时，或者说其误差大于所限定的数值时，就要对网络进行校正。

首先，定义误差函数

1s22E(w,b)=(tka2k)

2k1其次，给出权值的变化 ① 输出层的权值变化

从第i个输入到第k个输出的权值为：

w2ki其中：

Ekia1i w2kikiekf2'， eklka2k

② 隐含层的权值变化

从第j个输入到第i个输出的权值为：

wij其中：

Eijpjw1ij01 (η为学习系数)

ijeif1'

eikiw2ki

k1s2由此可以看出：①调整是与误差成正比，即误差越大调整的幅度就越大。②调整量与输入值大小成比例，在这次学习过程中就显得越活跃，所以与其相连的权值的调整幅度就应该越大，③调整是与学习系数成正比。通常学习系数在0.1~0.8之间，为使整个学习过程加快，又不会引起振荡，可采用变学习率的方法，即在学习初期取较大的学习系数随着学习过程的进行逐渐减小其值。最后，将输出误差由输出层经中间层传向输入层，逐层进行校正。

由此可知，BP神经网络的求解流程为：

程序代码：

根据上述建立的输入层和输出层数学模型，对PA利用MATLAB编程求解（以5月31日为例）：

%PA 5月31:0:0 计算的MATLAB 程序如下：BP神经网络预测 %5月31位于第23行 clc,clear

load('data.mat') %原始数据以列向量的方式存放在workplace文件中 PA=PA(2:29,:); %数据的标准化 N=size(PA); for j=1:N(1,2)

PAHminmax=minmax(PA(:,j)'); for i=1:N(1,1)

PA(i,j)=(PA(i,j)-PAHminmax(1,1))/(PAHminmax(1,2)-PAHminmax(1,1)); end end

%以每天的从0时计数起，每隔十五分钟作为输入 P=PA(1:21,:);

%以5月31的间隔十五分钟的发电量作为目标向量 T=PA(22,:);

%创建一个BP神经网络，每一个输入向量的取值范围为[0 ,1]，隐含层有22个神经元，输出层有一个神经元，隐含层的激活函数为tansig，输出层的激活函数为%logsig，训练函数为梯度下降函数，即标准学习算法 for i=1:21 a(i,1)=0; a(i,2)=1;

end

net=newff(a,[21,1],{'tansig','logsig'},'traingd'); net.trainParam.epochs=30000; net.trainParam.goal=0.01; %设置学习速率为0.1 LP.lr=0.1; %训练网络

net=train(net,P,T); %预测5月31的发电量数据 T1=sim(net,P);%预测值 %PA.5月31日发电量真实值 T0=PA(22,:);

%预测值与实际值的误差 for i=1:N(1,2)

error(1,i)=T1(1,i)-T0(1,i); end

%绘制误差图 figure(1)

plot(1:N(1,2),error(1:N(1,2)),'-*') grid on

xlabel('时点x'),ylabel('发电功率误差y');

title('PA5.31.0.0-5.31.23.45发电功率随时点误差变化图像') %函数绘图 figure(2)

plot(1:N(1,2),T0,1:N(1,2),T1,'r') legend('实际值','预测值',2) grid on

xlabel('时点x'),ylabel('发电功率y');

title('PA5.31.0.0-5.31.23.45发电功率随时点变化图像')

输出：

图4 图5

同理，对于PA、PB、PC、PD、P4和P58等机组在5月31日和其他日期，可以利用MATLAB得到形如图4、图5的风电功率随时间变化图和误差图。

分析：由图4、图5可知，预测曲线可以较好地拟合原始数据曲线，在每个时点的值逼近于实际值，预测曲线与实际曲线的走势大概相同，，故该预测模型较合理。

（4）建立马尔科夫链模型并求解 建立模型和求解：

由于风电功率变化的不确定性，而未来一段时间内的风电功率对历史数据有一定的依赖性，又有一定的随机变化性，因此我们利用马尔科夫链模型进行数据预测。以PA为例，我们选用5月29日和5月30日的数据预测5月31日的风电机功率，选用5月30日和5月31日的数据预测6月1日的风电机功率，以此类推。

由附件2中已给出的各天每隔15分钟的各机组风电功率值，我们定义A机组的增长率

为：

以A机组为例，可以利用MATLAB绘出其增长率变化趋势图：

图6

于是对于任意一个状态Ei有状态转移方程：

其中Pi为状态Ei的状态概率向量，P为进一步转移概率矩阵。

以A机组为例，取5月29-5月30之间的194组增长率数据参加计算，其余机组求解类似。由A机组5月29-5月30之间的数据可计算各个状态Ei的频数：

M=Mi=（31,65,0,74,20）

然后计算样本空间中从状态Ei一步转移到状态Ej的样本个数Mij：

计算一步转移概率矩阵P：

计算任一状态Ei的状态概率Pi（n）：

通过MATLAB编程，可得5月31日内任一时点任一状态的概率，MATLAB代码如下：

%PA 5月31:0:0 计算的MATLAB 程序如下：马尔克夫链预测 clc,clear;

%构造马尔科夫链模型

k=1;

for j=20:21 for l=1:n

zsj(1,k)=PA(j,l); k=k+1; end end

n1=size(zsj);

b=zeros(1,n1(1,2)-1);%存放状态值 %2 1 0 -1 -2

%'2表示快速下降、1表示缓慢下降、0表示基本保持不变、1表示缓慢上升、2表示快速上升' c=zeros(1,25);%表示相邻状态发生的概率 j=1;

for i=1:n1(1,2)-1

a(j,i)=(zsj(j,i+1)-zsj(j,i))/zsj(j,i); if a(j,i)>=0.3 b(j,i)=2;

elseif (a(j,i)<0.3)&&(a(j,i)>0) b(j,i)=1;

elseif a(j,i)==0 b(j,i)=0;

elseif a(j,i)>-0.3 && a(j,i)<0 b(j,i)=-1;

elseif a(j,i)<-0.3 b(j,i)=-2; end end

%c（1,25）表示相邻状态发生的概率统计 for i=1:n1(1,2)-2

if (b(j,i)==2&&b(j,i+1)==2) c(j,1)=c(j,1)+1;

elseif(b(j,i)==2&&b(j,i+1)==1) c(j,2)=c(j,2)+1;

elseif(b(j,i)==2&&b(j,i+1)==0) c(j,3)=c(j,3)+1;

elseif(b(j,i)==2&&b(j,i+1)==-1) c(j,4)=c(j,4)+1;

elseif(b(j,i)==2&&b(j,i+1)==-2) c(j,5)=c(j,5)+1;

elseif(b(j,i)==1&&b(j,i+1)==2) c(j,6)=c(j,6)+1;

elseif(b(j,i)==1&&b(j,i+1)==1) c(j,7)=c(j,7)+1;

elseif(b(j,i)==1&&b(j,i+1)==0) c(j,8)=c(j,8)+1;

elseif(b(j,i)==1&&b(j,i+1)==-1) c(j,9)=c(j,9)+1;

elseif(b(j,i)==1&&b(j,i+1)==-2) c(j,10)=c(j,10)+1;

elseif(b(j,i)==0&&b(j,i+1)==2) c(j,11)=c(j,11)+1; elseif(b(j,i)==0&&b(j,i+1)==1) c(j,12)=c(j,12)+1; elseif(b(j,i)==0&&b(j,i+1)==0) c(j,13)=c(j,13)+1;

elseif(b(j,i)==0&&b(j,i+1)==-1) c(j,14)=c(j,14)+1;

elseif(b(j,i)==0&&b(j,i+1)==-2) c(j,15)=c(j,15)+1;

elseif(b(j,i)==-1&&b(j,i+1)==2) c(j,16)=c(j,16)+1;

elseif(b(j,i)==-1&&b(j,i+1)==1) c(j,17)=c(j,17)+1;

elseif(b(j,i)==-1&&b(j,i+1)==0) c(j,18)=c(j,18)+1;

elseif(b(j,i)==-1&&b(j,i+1)==-1) c(j,19)=c(j,19)+1;

elseif(b(j,i)==-1&&b(j,i+1)==-2) c(j,20)=c(j,20)+1;

elseif(b(j,i)==-2&&b(j,i+1)==2) c(j,21)=c(j,21)+1;

elseif(b(j,i)==-2&&b(j,i+1)==1) c(j,22)=c(j,22)+1;

elseif(b(j,i)==-2&&b(j,i+1)==0) c(j,23)=c(j,23)+1;

elseif(b(j,i)==-2&&b(j,i+1)==-1) c(j,24)=c(j,24)+1;

elseif(b(j,i)==-2&&b(j,i+1)==-2) c(j,25)=c(j,25)+1; end end

d=zeros(1,5); i=1;

%统计某一个状态到另一个状态的次数和 for j=1:25 if(j<6)

d(i,1)=d(i,1)+c(i,j); elseif(j>5&&j<11)

d(i,2)=d(i,2)+c(i,j); elseif(j>10&&j<16)

d(i,3)=d(i,3)+c(i,j); elseif(j>15&&j<21)

d(i,4)=d(i,4)+c(i,j); else

d(i,5)=d(i,5)+c(i,j); end end

e= zeros(5,5);%状态概率矩阵 i=1;

for j=1:25 if(j<6)

if(d(i,1)==0) e(5,j)=0; else

e(1,j)=c(i,j)/d(i,1); %该2状态转移到另一个状态的概率 end

elseif(j>5&&j<11) if(d(i,2)==0) e(5,j-5)=0; else

e(2,j-5)=c(i,j)/d(i,2); %该1状态转移到另一个状态的概率 end

elseif(j>10&&j<16) if(d(i,3)==0) e(5,j-10)=0; else

e(3,j-10)=c(i,j)/d(i,3); %该0状态转移到另一个状态的概率 end

elseif(j>15&&j<21) if(d(i,4)==0) e(5,j-15)=0; else

e(4,j-15)=c(i,j)/d(i,4); %该-1状态转移到另一个状态的概率 end else

if(d(i,5)==0) e(5,j-20)=0; else

e(5,j-20)=c(i,j)/d(i,5); %该-2状态转移到另一个状态的概率 end end end

f=b(:,n1(1,2)-1); %最后的增长模式 g=zeros(1,5);

m=1;

if(f(1,1)==2) %预测 h=[1 0 0 0 0 ]*e; for k=1:96 h=h*e;

yuce(m,:)=h; m=m+1; end

elseif(f(1,1)==1)

g(1,:)=[0 1 0 0 0 ]; h=g(1,:)*e; for k=1:96 h=h*e;

yuce(m,:)=h; m=m+1; end

elseif(f(1,1)==0)

g(1,:)=[0 0 1 0 0 ]; h=g(1,:)*e; for k=1:96 h=h*e;

yuce(m,:)=h; m=m+1; end

elseif(f(1,1)==-1)

g(1,:)=[0 0 0 1 0 ]; h=g(1,:)*e; for k=1:96 h=h*e;

yuce(m,:)=h; m=m+1; end

elseif(f(1,1)==-2)

g(1,:)=[0 0 0 0 1 ]; h=g(1,:)*e; for k=1:96 h=h*e;

yuce(m,:)=h; m=m+1; end end

由此可得5月31日各时点预测结果为： 5.31.23.40 5.31.23.45 5.31.00：00 5.31.00：15 ... E4 E4 E4 E4 预测 ... 由于5月31日内A机组风电功率的增长状态为E4，即缓慢下降对每个状态向量，均取其中最大的那个概率值，因此可知，在未来一天内，A机组风电功率趋于缓慢下降的概率较其他状态的概率大得多，因此可以说A机组风电功率5月31日内趋于缓慢下降。 对于PA、PB、PC、PD、P4和P58等机组在5月31日和其他日期，可依照上述方法得出预测结果。

（5）模型比较

进行模型比较的两个考核指标如下： （1）准确率为：

r11（2）合格率为：

1N2PmkPpk100%capk1 N1r2NBk1Nk100%

为了保证数据的合理性，避免偶然误差，三种模型中均省去6月6日当天时段准确率和合格率的计算，取5月31日到6月5日的数据进行分析。计算可知，三个模型的考核指标如下表：

3次指数平滑预测模型的考核指标

表1

BP神经网络预测模型的考核指标

表2

马尔科夫链模型的考核指标

表3

分析：

（1） 三次指数平滑模型预测出来的数据的准确率大于0.95，合格率邓毅1，与实际值很

接近，这个结果与α的取值有关，由于α值很大，实际值的波动特征会被削减，故预测出来的值与实际值很接近，而在实际中，实际值是不会提前知道的，所以三次只是平滑模型不能应用于实际。

（2） BP神经网络模型的准确率在0.6到0.85之间，合格率在0.61到0.86之间，是一个

合理的范围，由于BP神经网络输入值是实际值，根据前期的实际值训练，寻找一个最优的模拟网络来预测下一时段的值，从而BP神经网络的预测结果不受外界约束，主要取决于输入值，因此具有可行性。

（3） 马尔科夫链模型的准确率在0.5-0.8之间，预测还算可以，然而马尔科夫链模型只能

预测某一个状态发生的概率，不能预测出实际值，因此此模型不可行。

2.3 问题2的分析、建模与求解 （1）建立BP神经网络模型 建立模型：

首先根据问题1中选出的较优预测方法即BP神经网络，经过网络的学习训练计算出单台风电机组功率以及多机总功率预测值和他们的相对误差。然后再绘出单机、多机相对误差随各时点变化的相对误差图。最后经过观察分析相对误差图以及机理分析可找出普遍规律。从而可对风电机组汇聚后给风电功率预测误差带来影响做出预期。

在问题1中，我们已经通过分别建立三次指数平滑、BP神经网络预测模型，再对预测结果进行误差分析，并比较准确率，合格率等指标，判断哪一个预测效果最佳。最终，得出BP神经网络较优，从而推荐BP神经网络算法。因此对于本问题，可建立BP神经网络模型。

对于输出层，有：

对于隐层，有：

其中，

模型的求解类似问题1中的模型2，在问题1的基础上将BP神经网络模型算出来的各个机组预测值与实测值间相对误差进行分析。

（2）分别绘制PA、PB、PC、PD与P4、P58的相对误差变化图（参考问题2程序） 程序代码（以PA与P4、P58为例）：

clc,clear,close all figure(1)

%计算PA的预测误差与P4的预测误差的规律 load('data.mat')

%PA的神经网络模拟后的预测相对误差值 %P4的神经网络模拟后的预测相对误差值 n1=size(WA); k=1;

for i=1:n1(1,2) for j=1:n1(1,1)

A(k,1)=WA(j,i);%5.31-6.5合并于一列 k=k+1; end end

l=1;

for i=1:n1(1,2) for j=1:n1(1,1)

w4(l,1)=W4(j,i);%5.31-6.5合并于一列 l=l+1; end end

t=n1(1,1)*n1(1,2);

plot(1:t,A(:,1),'r',1:t,w4(:,1)) grid on

title('PA的预测误差与P4的预测误差随时间变化图'); legend('PA的预测误差','P4的预测误差',2)

xlabel('时点x'),ylabel('发电功率预测误差y');

figure(2)

%计算PA的预测误差与P58的预测误差的规律 load('data.mat')

%PA的神经网络模拟后的预测相对误差值 %P4的神经网络模拟后的预测相对误差值 n1=size(WA); k=1;

for i=1:n1(1,2) for j=1:n1(1,1)

A(k,1)=WA(j,i);%5.31-6.5合并于一列 k=k+1; end end

l=1;

for i=1:n1(1,2) for j=1:n1(1,1)

w58(l,1)=W58(j,i);%5.31-6.5合并于一列 l=l+1; end end

t=n1(1,1)*n1(1,2);

plot(1:t,A(:,1),'r',1:t,w58(:,1)) grid on

title('PA的预测误差与P58的预测误差随时间变化图'); legend('PA的预测误差','P58的预测误差',2)

xlabel('时点x'),ylabel('发电功率预测误差y');

figure(3)

%计算PA的预测误差与P58的预测误差的规律 subplot(2,1,1) load('data.mat')

%PA的神经网络模拟后的预测相对误差值 %P4的神经网络模拟后的预测相对误差值 n1=size(WA); k=1;

for i=1:n1(1,2) for j=1:n1(1,1)

A(k,1)=WA(j,i);%5.31-6.5合并于一列 k=k+1; end end

l=1;

for i=1:n1(1,2) for j=1:n1(1,1)

w58(l,1)=W58(j,i);%5.31-6.5合并于一列 l=l+1; end end

t=n1(1,1)*n1(1,2);

plot(1:t,A(:,1),'r',1:t,w58(:,1)) grid on

title('PA的预测误差与P58的预测误差随时间变化图'); legend('PA的预测误差','P58的预测误差',2)

xlabel('时点x'),ylabel('发电功率预测误差y'); subplot(2,1,2) hold on

%计算PA的预测误差与P4的预测误差的规律 load('data.mat')

%PA的神经网络模拟后的预测相对误差值 %P4的神经网络模拟后的预测相对误差值 n1=size(WA); k=1;

for i=1:n1(1,2) for j=1:n1(1,1)

A(k,1)=WA(j,i);%5.31-6.5合并于一列 k=k+1;

end end

l=1;

for i=1:n1(1,2) for j=1:n1(1,1)

w4(l,1)=W4(j,i);%5.31-6.5合并于一列 l=l+1; end end

t=n1(1,1)*n1(1,2);

plot(1:t,A(:,1),'r',1:t,w4(:,1)) grid on

title('PA的预测误差与P4的预测误差随时间变化图'); legend('PA的预测误差','P4的预测误差',2)

xlabel('时点x'),ylabel('发电功率预测误差y');

输出：

图7：PA和P4、P58相对误差变化图 图8：PB和P4、P58相对误差变化图

图9：PC和P4、P58相对误差变化图 图10：PD和P4、P58相对误差变化图

（3）结果分析 由图7-10可知PA、PB、PC、PD的预测误差都很小，但是P4误差较大，这说明多机组接入时，由于各机组之间的相互影响，使实际输出功率变小，导致预测不准确。多机组产生对的误差几乎都是出现在同一时点，而其他时点不受影响。

因此，在预测多机组同时工作时的总输出功率时必须把这些特殊时点剔除。而这些“异常”时点可通过历史数据分析得到。得出的结果在风电机组原理上也能说通：风电机组即普通的发电机，只不过此时的原动力来源于风，将风能转化成机械能，再通过转换装置拖动发电机旋转产生电能，由于发电机内部有很多线圈，因此发电的同时线圈也会耗能，这就解释了为什么附件中所给的数据会出现负值的现象，即发电机并未正常工作，其内部的线圈反而消耗能量的结果。正因为这样，当多台发电机汇聚运行时，由于线圈间互感的作用以及线圈

本身耗能的原因，输出地总功率会小于单个电机工作时的功率之和。至于为何只出现在某些特殊时点，则需要更深入的研究了。

2.4问题3的分析、建模与求解

（1）建立动态神经网络时间序列模型（NAR模型） 对于问题3，由问题1可知，BP神经网络模型预测结果较好，但从表3可知，BP神经网络虽然在一定程度上可以合理地预测时点的大致趋势，但是预测值与实际值的误差仍较大，为了弥补BP神经网络的弊端，可以改用动态神经网络时间序列模型，该方法的记忆功能对时间序列的滞后性给予了一定的弥补，并且精度较高，其结构原理图如图11：

图11

动态神经网络时间序列模型的算法流程如下：

1） 对PA数据进行标准化，使各值处于0-1之间；

2） 选取训练神经网络的数据（预测某天，则选取该天之前的所有天数据）； 3） 对训练动态神经网络的数据进行矩阵变换，使其变为一列； 4） 启动MATLAB的Neural Network Start，如图12所示；

图12

5） 点击Time Series Tool，进入动态神经网络序列预测工具箱；

6） 使用非线性自回归模型，选择训练该网络的数据作为输入行向量；

7） 对数据进行分割，采用系统默认值，输入数据的70%作为训练数据，15%作为验证数据，

其余15%作为测试数据；

8） 反复调节隐层神经元个数和时间滞后的yi个数，并反复进行训练，直到达到规定要求；

执行上述步骤，对PA在5月1日，6月1日、6月2日、6月3日、6月4日、6月5日、6月6日进行预测，并不断地训练网络，对PB、PC、PD、P4、P58亦如此。

（2）短期预测（参考问题3程序）

建立动态神经网络时间序列模型如图13：

图13

设置隐层神经元个数为15，利用MATLAB求解（以PA在5月31日为例）：

程序代码：

%P4 5月31:0:0 计算的MATLAB 程序如下：BP神经网络预测 %5月31位于第22行 clc,clear

load('data.mat') %原始数据以列向量的方式存放在workplace文件中 PA=PA(2:29,:); %数据的标准化 N=size(PA); for j=1:N(1,2)

PAHminmax=minmax(PA(:,j)'); for i=1:N(1,1)

PA(i,j)=(PA(i,j)-PAHminmax(1,1))/(PAHminmax(1,2)-PAHminmax(1,1)); end end

%以每天的从0时计数起，每隔十五分钟作为输入 P=PA(1:21,:);

%以5月31的间隔十五分钟的发电量作为目标向量 T=PA(22,:);%实际值，应预测的值 T=T';

n=size(P); k=1;

for i=1:n(1,1) for j=1:n(1,2)

P1(k,1)=P(i,j); k=k+1; end end

% Solve an Autoregression Time-Series Problem with a NAR Neural Network % P1 - feedback time series.

targetSeries = tonndata(T,false,false);

% Create a Nonlinear Autoregressive Network feedbackDelays = 1:6; hiddenLayerSize = 15;

net = narnet(feedbackDelays,hiddenLayerSize);

% open loop or closed loop feedback modes.

[inputs,inputStates,layerStates,targets] = preparets(net,{},{},targetSeries);

% Setup Division of Data for Training, Validation, Testing net.divideParam.trainRatio = 70/100; net.divideParam.valRatio = 15/100; net.divideParam.testRatio = 15/100;

% Train the Network

[net,tr] = train(net,inputs,targets,inputStates,layerStates);

% Test the Network

outputs = net(inputs,inputStates,layerStates); errors = gsubtract(targets,outputs);

performance = perform(net,targets,outputs)

% View the Network view(net)

netc = closeloop(net);

[xc,xic,aic,tc] = preparets(netc,{},{},targetSeries); yc = netc(xc,xic,aic);

perfc = perform(net,tc,yc)

nets = removedelay(net);

[xs,xis,ais,ts] = preparets(nets,{},{},targetSeries); ys = nets(xs,xis,ais);

closedLoopPerformance = perform(net,tc,yc)

得到PA的实时预测图如图14：

图14

其均方差图如图15：

图15

分析：由图14可知，经过动态神经网络训练后，预测值与真实值的误差很小，各个时点的预测值与真实值相差无几，大部分的点都能正确预测，因此该网络训练合理。由图15可知，均方误差较小，其中训练用的数据误差最小，大部分在0.1-0.001之间。 同样对6月1日，6月2日……6月6日进行预测，利用今天对明天的预测增加前一天的真实值进行再训练，逐步向前预测，对于PB、PC、PD机组同样如此，即可得出每个机组的预测结果。

（3）长期预测（参考问题3程序） 跟短期预测类似，长期预测只是改变了预测的目标量，由于短期训练每次训练误差都不一样，因此采用重新训练新的动态神经网络进行预测，其MATLAB程序如下（以P4为例）：

%P4 5月31:0:0 计算的MATLAB 程序如下：BP神经网络预测 %5月31位于第22行 clc,clear

load P4 %原始数据以列向量的方式存放在workplace文件中 P4=P4(2:29,:); %数据的标准化 N=size(P4); for j=1:N(1,2)

P4Hminmax=minmax(P4(:,j)'); for i=1:N(1,1)

P4(i,j)=(P4(i,j)-P4Hminmax(1,1))/(P4Hminmax(1,2)-P4Hminmax(1,1)); end end

%以每天的从0时计数起，每隔十五分钟作为输入 P=P4(1:21,:);

%以5月31的间隔十五分钟的发电量作为目标向量 T=P4(22:28,:);%实际值，应预测的值 n=size(P); k=1;

for i=1:n(1,1) for j=1:n(1,2)

P1(k,1)=P(i,j); k=k+1; end end m=1;

n1=size(T);

for i=1:n1(1,1) for j=1:n1(1,2)

T1(m,1)=T(i,j); m=m+1; end end

% Solve an Autoregression Time-Series Problem with a NAR Neural Network % Script generated by NTSTOOL

% Created Sun Nov 27 17:30:25 CST 2011 %

% This script assumes this variable is defined: %

% P1 - feedback time series.

targetSeries = tonndata(T1,false,false);

% Create a Nonlinear Autoregressive Network feedbackDelays = 1:6; hiddenLayerSize = 15;

net = narnet(feedbackDelays,hiddenLayerSize);

% Prepare the Data for Training and Simulation

% The function PREPARETS prepares timeseries data for a particular network, % shifting time by the minimum amount to fill input states and layer states. % Using PREPARETS allows you to keep your original time series data unchanged, while

% easily customizing it for networks with differing numbers of delays, with % open loop or closed loop feedback modes.

[inputs,inputStates,layerStates,targets] = preparets(net,{},{},targetSeries);

% Setup Division of Data for Training, Validation, Testing net.divideParam.trainRatio = 70/100; net.divideParam.valRatio = 15/100;

net.divideParam.testRatio = 15/100;

% Train the Network

[net,tr] = train(net,inputs,targets,inputStates,layerStates);

% Test the Network

outputs = net(inputs,inputStates,layerStates); errors = gsubtract(targets,outputs);

performance = perform(net,targets,outputs)

% View the Network view(net)

% Plots

% Uncomment these lines to enable various plots. %figure, plotperform(tr) %figure, plottrainstate(tr)

%figure, plotresponse(targets,outputs) %figure, ploterrcorr(errors)

%figure, plotinerrcorr(inputs,errors)

% Closed Loop Network

% Use this network to do multi-step prediction.

% The function CLOSELOOP replaces the feedback input with a direct % connection from the outout layer. netc = closeloop(net);

[xc,xic,aic,tc] = preparets(netc,{},{},targetSeries); yc = netc(xc,xic,aic);

perfc = perform(net,tc,yc)

nets = removedelay(net);

[xs,xis,ais,ts] = preparets(nets,{},{},targetSeries); ys = nets(xs,xis,ais);

closedLoopPerformance = perform(net,tc,yc)

图16为P4风电机组长期预测图：

图16

在问题二中，我们讨论过“异常点”对预测风电功率精度有重要影响，避免在那些“异常点”下多机组同时工作；在这些“异常点”的影响下，导致风电机组的发电功率异常。 当然，风电功率预测精度不能无限的提高，风电功率受风力、风速、气温、气压等因素的影响，变化没规律；而现行的预测精度一方面受外界因素的影响，一方面受模型本身因素影响，即存在偶然误差以及系统误差，使得风电功率预测精度不能无限的提高。