风电功率预测问题 全国三等奖论文.doc

风电功率预测问题摘要：现今风力发电主要利用的是近地风能。近地风具有波动性、间歇性、低能量密度等特点，因而风电功率也是波动的。为减少大幅度地风电功率波动对电网的功率平衡和频率调节带来的不利影响，本文讨论多种预测模型在所给数据的基础上对风电场的发电功率进行尽可能准确地预测。本文采用了“基于小波变换和神经网络的风电功率预测方法”、“灰色预测模型”、“基于ARMA的卡尔曼预测模型”三种方法对风电功率进行预测。根据从一风电场获知的58台风电机28天内的相关运行数据，针对风电功率的随机时间序列特性，初步建立时间序列模型预测风电功率，同时基于相同的数据建立BP模型和小波分析，运用MATLAB软件编程预测输出，进行误差分析。在满足预测精度要求的基础上对比三种模型的结果，选出最优的模型。为实现风电功率的进一步提高精度预测，提出一种基于小波分析法与滚动式时间序列法组合建模的优化算法。关键词：小波神经网络预测 灰色预测模型 卡尔曼预测模型 优化算法一、问题重述根据百度百科，“风”是“跟地面大致平行的空气流动，是由于冷热气压分布不均匀而产生的空气流动现象”。风能是一种可再生、清洁的能源，风力发电是最具大规模开发技术经济条件的非水电再生能源。现今风力发电主要利用的是近地风能。近地风具有波动性、间歇性、低能量密度等特点，因而风电功率也是波动的。大规模风电场接入电网运行时，大幅度地风电功率波动会对电网的功率平衡和频率调节带来不利影响。如果可以对风电场的发电功率进行预测，电力调度部门就能够根据风电功率变化预先安排调度计划，保证电网的功率平衡和运行安全。因此，如何对风电场的发电功率进行尽可能准确地预测，是急需解决的问题。根据电力调度部门安排运行方式的不同需求，风电功率预测分为日前预测和实时预测。日前预测是预测明日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）。问题1：风电功率实时预测及误差分析。请对给定数据进行风电功率实时预测并检验预测结果是否满足附件1中的关于预测精度的相关要求。具体要求：（1） 采用不少于三种预测方法（至少选择一种时间序列分析类的预测方法）；（2） 预测量：a PA, PB, PC, PD； bP4； cP58。（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.1模型的假设：1.所给数据真实可靠，排除人为因素干扰带来的误差；2.短期内不存在大的自然灾害，例如地震、海啸及台风等等；3.预测期间风电机组分布无重大调整。2.2符号说明：L为滞后延迟算子；Pt 为风电功率的时间序列；p为自回归的阶数；t为零均值的系统白噪声； q为移动平均的阶数；为测量附加噪声；为第i+1 个时刻的预测风速；为第i+1个时刻的实际风速。三、问题分析风电装机容量的逐年快速增长，使得风电在电网中所占的比例也逐年增大。由于风能的随机性，风力发电使电力系统中不确定因素增加，对电力系统的安全运行提出了新的挑战。因而限制风电的发展规模。预测风电场的出力成为解决此问题的有效途径。如果能对风电场风速和发电功率进行比较准确的预测，则有利于电力系统调度部门必要时提前调整调度计划，有效的降低风电对电网安全的影响。3.1问题1的分析： 本文采用了“基于小波变换和神经网络的风电功率预测方法”、“灰色预测模型”、“基于ARMA的卡尔曼预测模型”三种方法对风电功率进行预测。根据从一风电场获知的58台风电机28天内的相关运行数据，针对风电功率的随机时间序列特性，初步建立时间序列模型预测风电功率，同时基于相同的数据建立BP模型和小波分析，运用MATLAB软件编程预测输出，进行误差分析。在满足预测精度要求的基础上对比三种模型的结果，选出最优的模型。3.2问题2的分析: 根据第一问建立的预测模型进行分析对电机组功率（PA，PB，PC，PD）的相对预测误差与多机总功率（P4，P58）预测的相对误差之间的规律，我们发现多机总功率（P4，P58）数据相对很平稳，所以用上面三个模型预测可以判断出多机总功率（P4，P58）预测的相对误差要比电机组功率（PA，PB，PC，PD）相对预测误差小。因此，通过风电机组汇聚可以较好的降低对风电功率预测的相对误差。3.3问题3的分析：为实现风电功率的进一步提高精度预测，提出一种基于小波分析法与滚动式时间序列法组合建模的优化算法。思路流程图如图所示：图 3.1 思路流程图 四、模型的建立与求解4.1问题1的求解：4.1.1基于小波变换和神经网络的风电功率预测方法影响风电功率的主要因素是风速，因此根据风电功率的历史数据直接建立风电功率预测模型是不准确的。而在风速预测的基础上建立风速和风电功率的关系模型进而预测风电场的输出功率，事实上没有真正实现风电场的输出功率预测。本文采用基于小波变换和神经网络的风电功率预测方法。主题思路是将风电功率和风速分别进行小波分解，得到风电功率和风速的细节部分和概貌部分。由于分解后的风速序列能很好地反映与之相对应的风电功率序列1，所以将风电功率和风速的概貌序列输入神经网络，建立风电功率概貌序列的预测模型，将风电功率和风速的细节序列输入神经网络，建立风电功率细节序列的预测模型，再将两模型的预测结果重构，得到风电功率的预测结果。（1）、小波变换对数据的预处理对风电功率P进行单尺度小波分解，得到风电功率的细节部分Pd和概貌部分Pa。概貌部分反映风电功率波动的周期性规律，细节部分反映风电功率波动的细节规律。用db3小波对历史风电功率序列和历史风速序列进行单尺度分解并重构。（2）建立BP神经网络模型1）、神经网络概述在数据分析时常常要用到预测，但传统的预测方法如回归分析和时问序列分析在涉及研究变量繁且数据又非线性时，有非常大的局限性。神经网络的出现为处理非线性问题提供了一条新途径。因为神经网络具有大规模的并行、分布式存储和处理、自组织、自适应和自学习能力，特别适合处理同时考虑许多因素和条件不精确，模糊的信息处理问题。BP网络一般具有三层或三层以上，包括输入层、中间层(即隐含层，可为一层或多层)和输出层，其结构形式如下图所示：图4.1 BP网络结构形式图BP神经网络学习过程分正向传播和反向传播，即信息的正向传播和误差的反向传播两个过程组成。在正向传播过程中，输入层各神经元负责接收来自外界的输入信息，并传递给中间层各神经元：中间层是内部信息处理层，负责信息变换，根据信息变化能力的需求，中间层可以设计为单隐层或者多隐层结构；最后一个隐层传递到输出层各神经元的信息，经进一步处理后，完成一次学习的正向传播处理过程，由输出层向外界输出信息处理结果。当实际输出与期望输出不符时，进入误差的反向传播阶段。通过网络将误差信号沿原来的连接通路返回，同时修改各层神经元的权值直至达到期望目标。周而复始的信息正向传播和误差反向传播过程，是各层权值不断调整的过程，也是神经网络学习训练的过程，此过程一直进行到网络输出误差减小到可以接受的程度。对于超短期风电功率预测，前一时刻风电场的输出功率，表征了预测的起始值，作为神经网络的输入神经元。2）、构建BP神经网络步骤如下Step1：输入输出归一化将输入，输出数据进行归一化处理，将每组数据都变为-1,1之间的数，所涉及的函数,这个函数的算法为：Step2：构建前馈网络可以按所分的类别利用MATLAB语言的BP神经网络工具箱提供的现成函数和BP 神经网络类，参考BP神经网络基本结构图4.2 BP神经网络基本结构建立一个前馈的BP神经网络模型:其中,第一个输入变量为列向量，存储的样本输入数据归一化后的最小值和最大值，第二个输入变量是一个行向量，存储BP神经网络各层的节点数输入，第三个输入变量是单元数组，由若干字符串构成，每个字符串对应于该层的传输函数类型。由此定义，可以确定BP神经网络的输入输出层。模型的输入：5月10日至5月29日第i时刻风电功率的xi值。 模型的输出：5月11日至5月30日第i时刻风电功率的yi值Step3：参数训练利用上述输入输出数据调用train()函数对BP神经网络参数进行训练。该函数调用格式为：X,Y分别存储样本点的输入和输出数据；Y1 ,E分别返回由BP神经网络计算出的输出和误差矩阵。Step4：泛化仿真调用函数进行泛化仿真，再调用函数把泛化仿真后的数据还原为原始的数量级。其中，分别为已归一化的输出量的最小值和最大值。Step5:确定输入数据根据题意，风电场每隔15分钟采样一组数据，每天一共有96个采样时刻，将5月30日第i时刻的风电功率值输入BP神经网络中，得到5月31日同时刻的风电功率值，并令i=i+1,回到Step1重复此算法直到i=96。Step6:输出5月31日风电功率的日前预测值在本题中，BP神经网络的输入层，隐含层，输出层对应的传输函数类型分别为,。3）、建立模型BP 神经网络的训练过程为: 对于a、5月31日0时0分至5月31日23时45分时间段的预测， 先用2006年5月29日0时0分到5月30日23时30分的指标历史数据作为网络的输入,用5月29日4时0分到5月30日23时45分的指标历史数据作为网络的输出,组成训练集对网络进行训练,使之误差达到满意的程度,用这样训练好的网络进行预测.采用滚动预测方法进行预测:滚动预测方法是通过一组历史数据预测未来某一时刻的值,然后把这一预测数据再视为历史数据继续预测下去,依次循环进行,逐步预测未来一段时期的值。用5月30日20时0分到23时45分数据作为网络的输入,5月31日0时0分的预测值作为网络的输出. 接着用5月30日20时15分到31日0时0分的数据作为网络的输入,31日0时15分的预测值作为网络的输出.依次类推,这样就得到5月31日0时0分至5月31日23时45分的预测值。目前在BP 网络的应用中,多采用三层结构. 根据人工神经网络定理可知,只要用三层的BP 网络就可实现任意函数的逼近. 所以训练结果采用三层BP模型进行模拟预测. 模型训练误差为0.00001,隐层单元数选取8个,学习速率为0.05,动态参数0.6,Sigmoid参数0.9,最大迭代次数3000.运行3000次后,样本拟合误差等于0.00021。3）、预测值的生成将小波分解后得到的风电功率的细节部分Pd和Pa概貌部分分别作为神经网络的输入，得到预测值的细节部分和概貌部分，将两者重构得到下一时刻的预测风电功率。对于5月31日0时0分至5月31日23时45分，运用matlab编程得到预测值如下：4.1.2运用卡尔曼滤波法预测风电功率本文由 ARMA方程推导出卡尔曼滤波状态方程和测量方程，从而将预测问题转化到状态空间，并利用卡尔曼滤波法预测了风电功率，与单独的时间序列分析法相比，卡尔曼滤波法能够提高风电功率的预测精度，并在一定程度上解决了时间序列分析法的预测时延问题，对电力系统的安全 稳定 经济运行以及提高运行效益具有重要意义。（1）、 ARMA模型 ARMA （p ,q）为（L）Pt=（L）t式中： L为滞后延迟算子；Pt 为风电功率的时间序列； p为自回归的阶数；t为零均值的系统白噪声； q为移动平均的阶数。ARMA模型的平稳条件是滞后多项式（L）的根在单位圆外，可逆条件为（L）的根都在单位圆外 ARMA模型对数据平稳性有要求，要在平稳时间序列的大前提下建模，所以要用ARMA模型预测风电功率，首先要检验风电功率时间序列的平稳性 时间序列平稳性检验的常用方法为增广 Dickey-Fuller （ADF ） 检验。ADF 检验包括1个回归方程：左边为序列的一阶差分项，右边则为序列的一阶滞后项 、滞后差分项，实际操作中，滞后阶数p的选择可根据保证t 是白噪声过程的最小p值的标准进行选择，在每种情况下，单位根检验都对回归式中pt-1的系数进行检验，如果系数显著不为零，那么pt包含单位根的假设将被拒绝，pt序列即是平稳。将风电功率时间序列pt识别为 ARMA（p，q）结构， 鉴于模型 ARMA（p，q） 的识别具有很大的灵活性，为了得到最合理的模型，本文采取了定阶步骤，根据时间序列的自相关、偏相关函数分析图，对 p=1， 2， 3， 4；q=1， 2， 3， 4， 5 多组可行阶数进行了参数估计，对所有备选模型进行模型诊断，筛选出备选模型集 由于许瓦兹信息准则（SIC ）的强一致性，在理论层面上能够渐进地选择真实模型，所以计算备选模型集中所有模型的 SIC 考虑到模型满足可逆性和稳定性条件，最终确定p和q，依照经典时间序列分析的步骤，在完成模型阶数识别后，使用极大似然估计法获得模型的参数估计模型如下：对于风电机组A：对于风电机组B：对于风电机组C：对于风电机组D：对于四个风电机组的和：对于58个风电机组的和： 运用 ARMA 模型对2006年 5月31日的96个时点及5月31日0时0分到6月6日23点45分风电功率进行预测，得到原始风电功率和预测风电功率曲线如图 1 所示 由图 1 可以看出，模型的预测结果曲线和实际功率曲线在形状上比较接近，但在变化规律发生改变时，预测曲线常常不能提前变化，而是存在一个“时滞” ，所以，虽然 ARMA 模型基本能捕捉到风电功率数据的变化规律，但该预测存在明显的延时性，而且预测精度不高 （2）、卡尔曼滤波模型要实现卡尔曼滤波法预测风电功率，首先必须推导出正确的状态方程和测量方程 因已通过时间序列分析建立了风电功率时间序列的 ARMA 模型，故可将 ARMA模型转换到状态空间，建立卡尔曼滤波的状态方程和测量方程1。对于一般的 ARMA模型其方程为： pt=a1pt-1+a m pt-m+t+b1t-1 +b m-1t-m+1m=max（p，q+1）此时pt 是经过零均值化所得时间序列ARMA模型，可写成状态空间模型为Xt=AXt-1+BEt Yt=CXt式中：Xt=（pt ,pt-1 ,pt-m-1 ）´,Et=(t, t-1, t-m+1) ´, C=(1,0,0) 且 a*= （ a1，a2， ，am-1），b*= （1，b1 , bm-1）； Im-1为状态空间模型的自回归系数。将风电功率序列Pt进行零均值化。对于风电机组A：对于风电机组B:对于风电机组C：对于风电机组D：对于4个风电机组和：对于4个风电机组和：Zt是一个零均值的序列，且Zt满足：将该式转换到状态空间，可得到状态方程式：测量方程式：式中：为测量附加噪声，为建模方便，假设其为白噪声。考虑到收敛速度并参考工程习惯，初始值近似取为Z（0） =0，p（0） =10 I，对输入噪声的协方差矩阵和测量噪声的协方差矩阵取经验值，应用 Matlab 软件实现卡尔曼滤波法预测风电功率，可得状态变量Zt 的预测值，再令由Zt得到Pt，即可得到预测风电功率曲线如下图所示由图可以看出：该模型的预测结果曲线和实际功率曲线在形状上较真实值更为接近，在风电功率变化规律发生变化时，卡尔曼滤波法有时能捕捉到变化信息，模型给出的预测值 Pt滞后于实际值的概率相对降低，使“时滞”问题较为缓和，所以，应用卡尔曼滤波法预测风电功率不仅提高了预测精度，而且在一定程度上解决了时间序列分析法的预测时延问题。根据附件1中“风电场预测预报考核指标”，计算出预测结果的准确率为：82.3%（三）、基于灰色模型的风电功率预测本文用传统灰色模型 GM（1，1）2，预测超短期风电功率，具体方法：时序风电功率为对数列按照传统 GM（1,1）模型的算法进行一次累加生成数列 即 构造 GM（1,1）模型的一阶微分方程:式中模型系数 a 和 b 可以用最小二乘法求得，即其中 预测生成序列为 预测的结果序列还原得 本文以风电机组在2006年5月30日22时0分至23时45分总共8个风速数据为历史时序风速数列，即 n=8，对此风电机组采用 GM（1，1）灰色模型预测从31日0时0分时刻开始的后续的15min、30min、45min、60min进行一步预测、两步预测、三步预测和四步预测风速。接着， 跟踪风速的实时采集数据不断更新历史时序风速数列，即去掉最老数据，保留实时采集数据，继续完成后续时间的一步至四步风速预测。预测误差r为：是第i+1 个时刻的预测风速；是第i+1个时刻的实际风速。实际风电功率和一步至四步预测结果如图所示： 图1风功率数据PA图2风功率数据PB图 3风功率数据PC图 4风功率数据PD图 5风功率数据P58注：6月1日至6月6日的灰色预测分析图见附录一。从图和表的结果可以看出， 在不断更新下一个15 min风速数据时，进行一步超短时风速预测的误差最小，其准确度可以接受，这说明灰色模型在超短时风速预测中应用的可行性和有效性。4.2问题2求解 根据第一问建立的预测模型进行分析对电机组功率（PA，PB，PC，PD）的相对预测误差与多机总功率（P4，P58）预测的相对误差之间的规律，我们发现多机总功率（P4，P58）数据相对很平稳，所以用上面三个模型预测可以判断出多机总功率（P4，P58）预测的相对误差要比电机组功率（PA，PB，PC，PD）相对预测误差小，由于平均相对误差：数据PAPBPCPDP4P58平均相对误差%18.9722.5122.7321.9924.7413.51因此，通过上面的图像可以得出一个结论：通过风电机组汇聚可以较好的降低对风电功率预测的相对误差。三、问题3的求解对于同一个预测问题，可以有多种预测方法，由于考虑的角度、方式等不同，不同的预测方法包含有不同的信息，将这些方法组合，可增大信息量，提高预测精度。为进一步提高预测精度，结合风电场功率随机性大的特点，并考虑风电场功率高阶矩的特征，建立了基于最大信息熵的风电场功率组合预测模型，应用实例表明，该方法可以有效提高预测精度。2因此，将问题一的三种模型的预测结果组合出优于每个单项的预测结果，即组合预测能有效地提高预测精度。本文采用回归组合预测法。现有3个预测模型，关于Pt的3个预测值 （i=1,2,3） 构造组合预测模型为： 选择使其加权平均 最小。可得出的加权最小 二乘估计为： 其中： 这种组合权值的计算方法叫做加权最小二乘方法。最终解得Pt=0.32+1.2Pt（1）+0.41Pt（2）+0.71Pt(3) 利用公式，基于问题一中得到的三种模型的预测结果可以得到综合测量法的测量结果：对于时间范围 a. 5月31日0时0分至5月31日23时45分；注：6月1日至6月6日图像见附录二。结果分析：在实际的风电场风速及功率预测中，可以应用的预测方法也有多种，并且各有优劣和不同的适用场合，即没有任何一种方法能保证在任何情况下都能获得最好的预测结果。因此，采用组合预测方法进行风速及风电场功率预测，充分利用各种单一预测模型的有用信息，势必可以提高风电场功率预测的精度。结语：目前国内外预测风电功率的方法模型都在不断完善中。如果能够针对预测环境的多变性 、预测者的滞后性 、预测资料的准确性 、预测方法的适宜性等影响因素，考虑到风速、风向、气温、气压、太阳辐射、空气温度、空气密度、云层高度、降水量、风电机的高度、风电机扫掠面积、风电机组的分布等综合因素，采用精确的风电功率预测模型，结合更加科学的数据处理方法，建立相应的风电场管理制度将会对预测系统，相信风电功率的预测精度还能进一步的提高。但预测只是在掌握现有信息的基础上，依照一定的方法与规律对未来的事情进行测算，并不能真确的了解事物的未来结果。而在全球气候变暖的大背景下，极端天气事件发生的概率和频率呈现增多的趋势，这有事风电功率预测增加了一层不确定性。所以，风电功率预测精度会不断提高，但不能无限提高。参考文献1 王世谦，苏娟，杜松怀 基于小波变换和神经网络的短期风电功率预测方法 农业工程学报 第26 卷增刊2 125-126 2010 年12 月2 时庆华，高山，陈昊 2种风电功率预测模型的比较 智能电网 第23 卷第6 期 33-34 2011 年6 月3 李俊芳，张步涵，谢光龙，李 妍，毛承雄 基于灰色模型的风速-风电功率预测研究 电力系统保护与控制 第38 卷第 19 期 152-153 2010 年10 月1 日4 韩爽，杨勇平，刘永前. 风电场功率短期预测方法研究附录一：6月1日至6月6日的灰色预测分析图：6月1日：图 1风功率数据PA图 2风功率数据PB图 3风功率数据PC图 4风功率数据PD图 5风功率数据P586月2日图 2风功率数据PA图 3风功率数据PB图 4风功率数据PC图 5风功率数据PD图 6风功率数据P586月3日图 7风功率数据PA图 8风功率数据PB图 9风功率数据PC图 10风功率数据PD图 11风功率数据P586月4日图 12风功率数据PA图 13风功率数据PB图 14风功率数据PC图 15风功率数据PD图 16风功率数据P586月5日图 17风功率数据PA图 18风功率数据PB图 19风功率数据PC图 20风功率数据PD图 21风功率数据P586月6日图 22风功率数据PA图 23风功率数据PB图 24风功率数据PC图 25风功率数据PD图 26风功率数据P58附录二：附录三：%小波神经网络预测源程序clearclc%读入数据Datapa datapb datapc datapd datap58%构造输入输出及测试矩阵%aa(1,:)=datapa(1,:),datapa(2,:),datapa(3,:),datapa(4,:),datapa(5,:),datapa(6,:),datapa(7,:),datapa(8,:),datapa(9,:),datapa(10,:),datapa(11,:),datapa(12,:),datapa(13,:),datapa(14,:),datapa(15,:),datapa(16,:),datapa(17,:),datapa(18,:),datapa(19,:),datapa(20,:),datapa(21,:),datapa(22,:);%aa(1,:)=datapb(1,:),datapb(2,:),datapb(3,:),datapb(4,:),datapb(5,:),datapb(6,:),datapb(7,:),datapb(8,:),datapb(9,:),datapb(10,:),datapb(11,:),datapb(12,:),datapb(13,:),datapb(14,:),datapb(15,:),datapb(16,:),datapb(17,:),datapb(18,:),datapb(19,:),datapb(20,:),datapb(21,:),datapb(22,:);%aa(1,:)=datapc(1,:),datapc(2,:),datapc(3,:),datapc(4,:),datapc(5,:),datapc(6,:),datapc(7,:),datapc(8,:),datapc(9,:),datapc(10,:),datapc(11,:),datapc(12,:),datapc(13,:),datapc(14,:),datapc(15,:),datapc(16,:),datapc(17,:),datapc(18,:),datapc(19,:),datapc(20,:),datapc(21,:),datapc(22,:);%aa(1,:)=datapd(1,:),datapd(2,:),datapd(3,:),datapd(4,:),datapd(5,:),datapd(6,:),datapd(7,:),datapd(8,:),datapd(9,:),datapd(10,:),datapd(11,:),datapd(12,:),datapd(13,:),datapd(14,:),datapd(15,:),datapd(16,:),datapd(17,:),datapd(18,:),datapd(19,:),datapd(20,:),datapd(21,:),datapd(22,:);aa(1,:)=datap58(1,:),datap58(2,:),datap58(3,:),datap58(4,:),datap58(5,:),datap58(6,:),datap58(7,:),datap58(8,:),datap58(9,:),datap58(10,:),datap58(11,:),datap58(12,:),datap58(13,:),datap58(14,:),datap58(15,:),datap58(16,:),datap58(17,:),datap58(18,:),datap58(19,:),datap58(20,:),datap58(21,:),datap58(22,:); for i=1:1:2000 input(i,:)=aa(:,i:i+15);endfor i=1:1:96 input_test(i,:)=aa(:,i+2000:i+2015);end for i=1:1:2000 output(i,:)=aa(:,i+16);end for i=1:1:96 output_test(i,:)=aa(:,i+2016);end %初始化M=16;N=1;n=6; lr1=0.01;lr2=0.001;maxgen=100; Wjk=randn(n,M);Wij=randn(N,n);a=randn(1,n);b=randn(1,n); d_Wjk=zeros(n,M);d_Wij=zeros(N,n);d_a=zeros(1,n);d_b=zeros(1,n); y=zeros(1,N);net=zeros(1,n);net_ab=zeros(1,n); inputn,inputps=mapminmax(input');outputn,outputps=mapminmax(output');inputn=inputn'outputn=outputn'%培养for i=1:maxgen error(i)=0; for kk=1;size(input,1) x=inputn(kk,:);待添加的隐藏文字内容1 yqw=outputn(kk,:); for j=1:n for k=1:M net(j)=net(j)+Wjk(j,k)*x(k); net_ab(j)=(net(j)-b(j)/a(j); end temp=mymorlet(net_ab(j); for k=1:N y=y+Wij(k,j)*temp; end end error(i)=error(i)+sum(abs(yqw-y); for j=1:n temp=mymorlet(net_ab(j); for k=1:N d_Wij(k,j)=d_Wij(k,j)-(yqw(k)-y(k)*temp; end temp=d_mymorlet(net_ab(j); for k=1:M for l=1:N d_Wjk(j,k)=d_Wjk(j,k)+(yqw(l)-y(l)*Wij(l,j); end d_Wjk(j,k)=-d_Wjk(j,k)*temp*x(k)/a(j); end for k=1:N d_b(j)=d_b(j)+(yqw(k)-y(k)*Wij(k,j); end d_b(j)=d_b(j)*temp/a(j); for k=1:N d_a(j)=d_a(j)+(yqw(k)-y(k)*Wij(k,j); end d_a(j)=d_a(j)*temp*(net(j)-b(j)/b(j)/a(j); end Wij=Wij-lr1*d_Wij; Wjk=Wjk-lr1*d_Wjk; b=b-lr2*d_b; a=a-lr2*d_a; d_Wjk=zeros(n,M); d_Wij=zeros(N,n); d_a=zeros(1,n); d_b=zeros(1,n); y=zeros(1,N); net=zeros(1,n); net_ab=zeros(1,n); endend %预测x=mapminmax('apply',input_test',inputps);x=x'for i=1:96 x_test=x(i,:); for j=1:1:n for k=1:1:M net(j)=net(j)+Wjk(j,k)*x_test(k); net_ab(j)=(net(j)-b(j)/a(j); end temp=mymorlet(net_ab(j); for k=1:N y(k)=y(k)+Wij(k,j)*temp; end end yuce(i)=y(k); y=zeros(1,N); net=zeros(1,n); net_ab=zeros(1,n);endynn=mapminmax('reverse',yuce,outpu