MATLAB数学建模与仿真（第2版·微课视频版）程序代码 第13章-风电功率预测问题代码

绘制各发电机输出功率曲线，分析数据的波动规律。实现代码如下：clc,clear,closeallload('data.mat')a=23;%某一天subplot(321),plot(PA(a,:));gridonxlabel('时点');ylabel('PA输出功率');title('PA发动机输出功率曲线')subplot(322),plot(PB(a,:));gridonxlabel('时点');ylabel('PB输出功率');title('P5发动机输出功率曲线')subplot(323),plot(PC(a,:));gridonxlabel('时点');ylabel('PC输出功率');title('PC发动机输出功率曲线')subplot(324),plot(PD(a,:));gridonxlabel('时点');ylabel('PD输出功率');title('PD发动机输出功率曲线')subplot(325),plot(P4(a,:));gridonxlabel('时点');ylabel('P4输出功率');title('P4发动机输出功率曲线')subplot(326),plot(P58(a,:));gridonxlabel('时点');ylabel('P58输出功率');title('P58发动机输出功率曲线')图13-2各发电机组输出功率曲线针对PA，设置平滑指数为0.9，对5月31日00:00-23:45时段内进行预测，MATLAB程序如下：clc,clearload('data.mat')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;fori=2:nst1(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);endxlswrite('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('实际值','预测值')gridonxlabel('时点x'),ylabel('发电功率y');title('PA-5发电功率随时点变化图像')图13-4PA风电机功率5月31日预测值对5月31日00:00-6月6日23：45整体上进行预测，MATLAB程序如下：clc,clearloadPAPA=PA(2:29,:);%数据的标准化N=size(PA);forj=1:N(1,2)PAHminmax=minmax(PA(:,j)');fori=1:N(1,1)PA(i,j)=(PA(i,j)-PAHminmax(1,1))/(PAHminmax(1,2)-PAHminmax(1,1));endend%以每天的从0时计数起，每隔十五分钟作为输入P=PA(1:21,:);%以5月31的间隔十五分钟的发电量作为目标向量T=PA(22:N(1,1),:);%创建一个BP神经网络，每一个输入向量的取值范围为[0,1]，隐含层有22个神经元，%输出层有一个神经元，隐含层的激活函数为tansig，输出层的激活函数为%logsig，%训练函数为梯度下降函数，即标准学习算法fori=1:21a(i,1)=0;a(i,2)=1;endnet=newff(a,[21,7],{'tansig','logsig'},'traingd');net.trainParam.epochs=25000;net.trainParam.goal=0.01;%设置学习速率为0.1LP.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;fori=1:n(1,1)forj=1:n(1,2)a(k,1)=T1(i,j);b(k,1)=T0(i,j);k=k+1;endendN=size(a);%绘制误差图plot(1:N(1,1),b(:,1),1:N(1,1),a(:,1),'r')gridonlegend('实际值','预测值')xlabel('时点x'),ylabel('发电功率y');title('PA-5发电功率实时函数图像')图13-5PA风电机功率5.31~6.6日预测值根据建立的输入层和输出层数学模型，对PA利用MATLAB编程求解（以5月31日为例）：clc,clearload('data.mat')PA=PA(2:29,:);%数据的标准化N=size(PA);forj=1:N(1,2)PAHminmax=minmax(PA(:,j)');fori=1:N(1,1)PA(i,j)=(PA(i,j)-PAHminmax(1,1))/(PAHminmax(1,2)-PAHminmax(1,1));endend%以每天的从0时计数起，每隔十五分钟作为输入P=PA(1:21,:);%以5月31的间隔十五分钟的发电量作为目标向量T=PA(22,:);%创建一个BP神经网络，每一个输入向量的取值范围为[0,1]，隐含层有22个神经元，输出层有一个神经%元，隐含层的激活函数为tansig，输出层的激活函数为logsig，训练函数为梯度下降函数，即标准学习算法fori=1:21a(i,1)=0;a(i,2)=1;endnet=newff(a,[21,1],{'tansig','logsig'},'traingd');net.trainParam.epochs=30000;net.trainParam.goal=0.01;%设置学习速率为0.1LP.lr=0.1;%训练网络net=train(net,P,T);%预测5月31的发电量数据T1=sim(net,P);%预测值%PA.5月31日发电量真实值T0=PA(22,:);%预测值与实际值的误差fori=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)),'-*')gridonxlabel('时点x'),ylabel('发电功率误差y');title('PA-5发电功率随时点误差变化图像')%函数绘图figure(2)plot(1:N(1,2),T0,1:N(1,2),T1,'r')legend('实际值','预测值')gridonxlabel('时点x'),ylabel('发电功率y');title('PA-5发电功率随时点变化图像')图13-7风电功率预测图图13-8风电功率误差图通过MATLAB编程，可得5月31日内任一时点任一状态的概率，MATLAB代码如下：clc,clear;loadPAPA=PA(2:29,:);yt=PA(22:28,:);n=length(yt);%增长率波动情况forj=1:7fori=1:n-1y(j,i)=(yt(j,i+1)-yt(j,i))/yt(j,i);endendt=1:n-1;y1=-0.5;y2=-0.3;y3=0;y2=0.3;y3=0.5;plot(t,y(1,:),t,y1,'-r',t,y2,'-r',t,y3,'-r')gridonxlabel('时点t'),ylabel('发电功率增长率y');title('PA-5发电功率随时点变化图像')%构造马尔科夫链模型k=1;forj=21:22forl=1:nzsj(1,k)=PA(j,l);k=k+1;endendn1=size(zsj);b=zeros(1,n1(1,2)-1);%存放状态值%210-1-2%'2表示快速下降、1表示缓慢下降、0表示基本保持不变、1表示缓慢上升、2表示快速上升'c=zeros(1,25);%表示相邻状态发生的概率j=1;fori=1:n1(1,2)-1a(j,i)=(zsj(j,i+1)-zsj(j,i))/zsj(j,i);ifa(j,i)>=0.3b(j,i)=2;elseif(a(j,i)<0.3)&&(a(j,i)>0)b(j,i)=1;elseifa(j,i)==0b(j,i)=0;elseifa(j,i)>-0.3&&a(j,i)<0b(j,i)=-1;elseifa(j,i)<-0.3b(j,i)=-2;endend%c（1,25）表示相邻状态发生的概率统计fori=1:n1(1,2)-2if(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;endendd=zeros(1,5);i=1;%统计某一个状态到另一个状态的次数和forj=1:25if(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);elsed(i,5)=d(i,5)+c(i,j);endende=zeros(5,5);%状态概率矩阵i=1;forj=1:25if(j<6)if(d(i,1)==0)e(5,j)=0;elsee(1,j)=c(i,j)/d(i,1);%该2状态转移到另一个状态的概率endelseif(j>5&&j<11)if(d(i,2)==0)e(5,j-5)=0;elsee(2,j-5)=c(i,j)/d(i,2);%该1状态转移到另一个状态的概率endelseif(j>10&&j<16)if(d(i,3)==0)e(5,j-10)=0;elsee(3,j-10)=c(i,j)/d(i,3);%该0状态转移到另一个状态的概率endelseif(j>15&&j<21)if(d(i,4)==0)e(5,j-15)=0;elsee(4,j-15)=c(i,j)/d(i,4);%该-1状态转移到另一个状态的概率endelseif(d(i,5)==0)e(5,j-20)=0;elsee(5,j-20)=c(i,j)/d(i,5);%该-2状态转移到另一个状态的概率endendendf=b(:,n1(1,2)-1);%最后的增长模式g=zeros(1,5);m=1;if(f(1,1)==2)%预测h=[10000]*e;fork=1:96h=h*e;yuce(m,:)=h;m=m+1;endelseif(f(1,1)==1)g(1,:)=[01000];h=g(1,:)*e;fork=1:96h=h*e;yuce(m,:)=h;m=m+1;endelseif(f(1,1)==0)g(1,:)=[00100];h=g(1,:)*e;fork=1:96h=h*e;yuce(m,:)=h;m=m+1;endelseif(f(1,1)==-1)g(1,:)=[00010];h=g(1,:)*e;fork=1:96h=h*e;yuce(m,:)=h;m=m+1;endelseif(f(1,1)==-2)g(1,:)=[00001];h=g(1,:)*e;fork=1:96h=h*e;yuce(m,:)=h;m=m+1;endendw=size(yuce);fori=1:w(1,1)max1=max(yuce(i,:));forj=1:w(1,2)ifyuce(i,j)==max1;liehao(i,1)=j;%1表示快速下降、2表示缓慢下降、3表示基本不变、4表示缓慢上升、5表示快速上升endendend图13-9A机组增长率变化趋势图分别绘制PA、PB、PC、PD与P4、P58的相对误差变化图，Maltab程序代码如下：clc,clear,closeallfigure(1)%计算PA的预测误差与P4的预测误差的规律load('data.mat')%PA的神经网络模拟后的预测相对误差值%P4的神经网络模拟后的预测相对误差值n1=size(WA);k=1;fori=1:n1(1,2)forj=1:n1(1,1)A(k,1)=WA(j,i);%5.31-6.5合并于一列k=k+1;endendl=1;fori=1:n1(1,2)forj=1:n1(1,1)w4(l,1)=W4(j,i);%5.31-6.5合并于一列l=l+1;endendt=n1(1,1)*n1(1,2);plot(1:t,A(:,1),'r',1:t,w4(:,1),'b')gridontitle('PA的预测误差与P4的预测误差随时间变化图');legend('PA的预测误差','P4的预测误差')xlabel('时点x'),ylabel('发电功率预测误差y');figure(2)%计算PA的预测误差与P58的预测误差的规律load('data.mat')%PA的神经网络模拟后的预测相对误差值%P4的神经网络模拟后的预测相对误差值n1=size(WA);k=1;fori=1:n1(1,2)forj=1:n1(1,1)A(k,1)=WA(j,i);%5.31-6.5合并于一列k=k+1;endendl=1;fori=1:n1(1,2)forj=1:n1(1,1)w58(l,1)=W58(j,i);%5.31-6.5合并于一列l=l+1;endendt=n1(1,1)*n1(1,2);plot(1:t,A(:,1),'r',1:t,w58(:,1),'b')gridontitle('PA的预测误差与P58的预测误差随时间变化图');legend('PA的预测误差','P58的预测误差')xlabel('时点x'),ylabel('发电功率预测误差y');figure(3)%计算PA的预测误差与P58的预测误差的规律subplot(2,1,1)load('data.mat')%PA的神经网络模拟后的预测相对误差值%P4的神经网络模拟后的预测相对误差值n1=size(WA);k=1;fori=1:n1(1,2)forj=1:n1(1,1)A(k,1)=WA(j,i);%5.31-6.5合并于一列k=k+1;endendl=1;fori=1:n1(1,2)forj=1:n1(1,1)w58(l,1)=W58(j,i);%5.31-6.5合并于一列l=l+1;endendt=n1(1,1)*n1(1,2);plot(1:t,A(:,1),'r',1:t,w58(:,1),'b')gridontitle('PA的预测误差与P58的预测误差随时间变化图');legend('PA的预测误差','P58的预测误差')xlabel('时点x'),ylabel('发电功率预测误差y');subplot(2,1,2)holdon%计算PA的预测误差与P4的预测误差的规律load('data.mat')%PA的神经网络模拟后的预测相对误差值%P4的神经网络模拟后的预测相对误差值n1=size(WA);k=1;fori=1:n1(1,2)forj=1:n1(1,1)A(k,1)=WA(j,i);%5.31-6.5合并于一列k=k+1;endendl=1;fori=1:n1(1,2)forj=1:n1(1,1)w4(l,1)=W4(j,i);%5.31-6.5合并于一列l=l+1;endendt=n1(1,1)*n1(1,2);plot(1:t,A(:,1),'r',1:t,w4(:,1),'b')gridontitle('PA的预测误差与P4的预测误差随时间变化图');legend('PA的预测误差','P4的预测误差')xlabel('时点x'),ylabel('发电功率预测误差y');图13-10PA和P4、P58相对误差变化图图13-11PB和P4、P58相对误差变化图图13-12PC和P4、P58相对误差变化图图13-13PD和P4、P58相对误差变化图设置隐层神经元个数为15，利用MATLAB求解（以PA在5月31日为例）：clc,clearload('data.mat')PA=PA(2:29,:);%数据的标准化N=size(PA);forj=1:N(1,2)PAHminmax=minmax(PA(:,j)');fori=1:N(1,1)PA(i,j)=(PA(i,j)-PAHminmax(1,1))/(PAHminmax(1,2)-PAHminmax(1,1));endend%以每天的从0时计数起，每隔十五分钟作为输入P=PA(1:21,:);%以5月31的间隔十五分钟的发电量作为目标向量T=PA(22,:);%实际值，应预测的值T=T';n=size(P);k=1;fori=1:n(1,1)forj=1:n(1,2)P1(k,1)=P(i,j);k=k+1;endend%SolveanAutoregressionTime-SeriesProblemwithaNARNeuralNetwork%P1-feedbacktimeseries.targetSeries=tonndata(T,false,false);%CreateaNonlinearAutoregressiveNetworkfeedbackDelays=1:6;hiddenLayerSize=15;net=narnet(feedbackDelays,hiddenLayerSize);%openlooporclosedloopfeedbackmodes.[inputs,inputStates,layerStates,targets]=preparets(net,{},{},targetSeries);%SetupDivisionofDataforTraining,Validation,Testingnet.divideParam.trainRatio=70/100;net.divideParam.valRatio=15/100;net.divideParam.testRatio=15/100;%TraintheNetwork[net,tr]=train(net,inputs,targets,inputStates,layerStates);%TesttheNetworkoutputs=net(inputs,inputStates,layerStates);errors=gsubtract(targets,outputs);performance=perform(net,targets,outputs)%ViewtheNetworkview(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)图13-16PA的实时预测图图13-17均方差图长期预测MATLAB程序如下（以P4为例）：clc,clearloadP4P4=P4(2:29,:);%数据的标准化N=size(P4);forj=1:N(1,2)P4Hminmax=minmax(P4(:,j)');fori=1:N(1,1)P4(i,j)=(P4(i,j)-P4Hminmax(1,1))/(P4Hminmax(1,2)-P4Hminmax(1,1));endend%以每天的从0时计数起，每隔十五分钟作为输入P=P4(1:21,:);%以5月31的间隔十五分钟的发电量作为目标向量T=P4(22:28,:);%实际值，应预测的值n=size(P);k=1;fori=1:n(1,1)forj=1:n(1,2)P1(k,1)=P(i,j);k=k+1;endendm=1;n1=size(T);fori=1:n1(1,1)forj=1:n1(1,2)T1(m,1)=T(i,j);m=m+1;endend%SolveanAutoregressionTime-SeriesProblemwithaNARNeuralNetwork%ScriptgeneratedbyNTSTOOL%CreatedSunNov2717:30:25CST2011%%Thisscriptassumesthisvariableisdefined:%%P1-feedbacktimeseries.targetSeries=tonndata(T1,false,false);%CreateaNonlinearAutoregressiveNetworkfeedbackDelays=1:6;hiddenLayerSize=15;net=narnet(feedbackDelays,hiddenLayerSize);%PreparetheDataforTrainingandSimulation%ThefunctionPREPARETSpreparestimeseriesdataforapartic