2013电工杯数学建模A题论文（国家二等奖）

A题：风电功率波动特性的分析

风电功率波动特性的分析

摘 要

本论文针对“风电功率波动特性的分析”问题，根据所给的风电机组功率数据建立风电功率波动特性的概率分布模型和灰色预测模型。由此，进行相关的问题分析及解决。

对于问题一、二，借鉴分离min级负荷的算法，采用滑动平均法分离s级风电功率，且处理了丢失数据及错误数据，以此提高数据的准确性。经过如此处理所给数据后，再采用Matlab的概率密度拟合工具箱dfittool得出五台风电机组的功率概率直方图及t location-scale分布、正态分布、逻辑斯特分布的概率分布图。发现t location-scale分布比其他分布更适于拟合各风电场概率密度函数，并作相关分析及检验。再用t location-scale分布以每日为时间窗宽，对5个风电功率分别计算30个时段的概率分布参数并做出检验。问题二与问题一方法一致，只需要从从上述5台机的风电功率数据中提取出间隔为1分钟的数据序列Pim(tk)，重复问题一的方法即可。

对于问题三，由上可以分别获得采样间隔为5s和1min的相关波动特性参数，以此为依照进行波动特性与时间尺度关系的分析。

对于问题四，整合出20台机组的数据，再分别将采样间隔改为1min，5min，15min，使用问题一中同样的方法分析处理即可。

对于问题五，预测未来四小时的风电场总功率，我们采用了灰色模型（Grey Model，GM），使用MATLAB对灰色模型GM（1，1）编程得到预测值，残差，级比偏差等相关数据结果。由于初步编程得出的预测模型为其累加后的方程，通过生成序列预测值及模型还原值之间的关系及之前所求的预测值模型易求的未来四小时风电场总功率预测模型。

对于问题六，将单台风电机组的输出功率在不同时间尺度下的波动特性曲线合在一起相比较，将风电场的总功率在不同时间尺度下的波动特性曲线合在一起比较，并且单台风电机组与风电场之间在相同时间尺度下也进行比较，以此分析其相关差别及可以获得的风电波动的时序差异性特点。

对于问题七，在已对机组和全场风电功率波动的整体分析为基础下，及对风电功率波动特性的整体认识，利用所获得的信息用来克服风电波动对电网运行的不利影响。并且基于我国东北吉林省西部百万千瓦风电场群实测风电功率数据，采用上述方法进行风电功率波动特性的分析，以此构建实例说明风电功率的波动特性。

关键词：概率分布模型、灰色预测模型、风电功率波动特性、时空差异性

1．问题的重述

风电是一种具有波动性且不完全受控的电源，大规模风力发电功率通常难以准确预测，因此在如何保证电网安全的前提下尽量提高风电联网规模，是亟需研究是问题，具有重要的理论和应用价值。

风电机组发出的功率主要与风速有关。由于风的不确定性、间歇性以及风电场内各机组间尾流的影响，使得风力发电机不能像常规发电机组那样根据对电能的需求来确定发电。

风电场通常有几十台、上百台风电机组。大型风电基地由数十甚至上百个风电场组成。因此，风电功率的波动有很强的时空差异性。

附件给出了某风电场中20台1.5MW风电机组30天的风电功率数据（单位为kW，间隔为5s），请做如下分析。

1．任选5个风电机组：

5sa）在30天的范围内，分析机组i的风电功率Pi(tk)波动符合哪几种概率分布？分别计算数值特征并进行检验，推荐最好的分布并说明理由。比较5个机组分布的异同。

b）用以上确定的最好的概率分布，以每日为时间窗宽，对5个风电功率分别计算30个时段的概率分布参数并做出检验；试比较不同机组（空间）、不同时段（时间）风电功率波动的概率分布以及与30天总体分布之间的关系，由此说明了什么？

2．在风电场实际运行中，由于数据存储和管理等方面的限制，难以集中记录全部风电机组功率的秒级数据。通常用分钟级间隔乃至更长间隔的数据来描述风电功率波动。试从上述5台机的风电功率数据中提取出间隔为1分钟的数据序

m列Pi(tk)。对于这5个序列，再做题1a）的分析。

m3．试分析用Pi(tk)代替Pi5s(tk)时，损失了那些风电功率波动信息？如何度量？有何影响？从上述全部计算中你能得出什么一般性的结论？

4．设全场20台风电机的总功率P（=∑P∑t）i(t)，试计算时间间隔为1分

m5m15m钟、5分钟和15分钟的总功率序列P∑(tk)，P∑(tk)，P∑(tk)，分析其波动的

5mm概率分布数值特征。若以P∑(tk)代替P∑(tk)来表征全场风电功率波动，损失了什么信息？如何度量？有何影响？

15m5．如果分别采用P∑5m(tk))和P∑(tk)作为样本来预测未来4小时（每15分钟一个点）风电场的总功率，请设计合适的预测模式（可取适当时段的数据作为历史数据建模，后续数据作为实际风电功率用于检验预测误差），分别给出不少于7天的滚动预测结果，分析比较2种方式的预测误差。

6．风电功率变化对电网运行的影响主要与其时序特性有关，比如风电大幅波动带来的调频机组爬坡速率分析。试分析单台风电机功率Pim(tk)与风电场总功率P∑m(tk)在时序上表现出的主要差别；前面得到的概率分布数值特征在分析时序波动特性方面有何作用？有何局限？

7．通过上述对机组和全场风电功率波动的分析，你对风电功率波动特性有何认识？这些认识如何用来克服风电波动对电网运行的不利影响？请构建实例来说明。

2.问题的分析

对于问题一、二，在分析数据时发现所给数据当中存在着部分数据丢失及错误数据，为提高数据的准确性，应采取用问题数据前后10个数据或更多数据取平均值的方法来代替该数据。经过如此处理所给数据后，借鉴分离s级负荷的算法，应采用滑动平均法分离s级风电功率，采用Matlab的概率密度拟合工具箱dfittool得出五台风电机组的功率波动直方图及t location-scale分布、正态分布、逻辑斯特分布的拟合曲线,并对图形作相关分析及用ttest函数作检验。再用研究所得的最符合的分布,以每日为时间窗宽，对5个风电功率分别计算30个时段的概率分布参数并做出检验。

问题二与问题一方法一致，只需要从上述5台机的风电功率数据中提取出间

m隔为1分钟的数据序列pi(tk)，重复问题一的方法即可。

对于问题三，由以上分析可以分别获得采样间隔为5s和1min的相关的波动特性参数，以此为依照进行问题的分析。

对于问题四，整合出20台机组的数据，再分别将采样间隔改为1min，5min，15min，使用问题一中同样的方法分析处理即可。

对于问题五，预测未来四小时的风电场总功率，采用灰色模型（Grey Model，GM），使用MATLAB对灰色模型GM（1，1）编程得到预测值，残差，级比偏差等相关数据结果。为预测未来7天的滚动数据，采用滚动预测的方法将前一段的预测值作为下一阶段的历史数据来进行预测。

对于问题六，将单台风电机组的输出功率在不同时间尺度下的波动特性曲线合在一起相比较，将风电场的总功率在不同时间尺度下的波动特性曲线合在一起比较，并且单台风电机组与风电场之间在相同时间尺度下也进行比较，以此分析其相关差别及可以获得的风电波动的时空差异性特点。

对于问题七，对机组和全场风电功率波动的分析，对风电功率波动特性整体认识，利用所获得的信息用来克服风电波动对电网运行的不利影响。并且基于我国东北吉林省西部百万千瓦风电场群实测风电功率数据，采用上述方法进行风电功率波动特性的分析，以此构建实例说明风电功率的波动特性。

3．模型的假设与符号说明

3.1模型的假设

1、每次建立模型所取的数据都准确，附件所给的采样间隔5秒的所有数据

都准确；

2、采样间隔的改变不影响数据的准确性 3.2符号说明

P5si(tk)

Pmi(tk)

Pm?(tk) P5m?(tk)

P15m∑(tk)

pt

pft

pmt yav ys

ymax

采样间隔为5秒的机组i的风电功率 采样间隔为1分钟的机组i的风电功率 采样间隔为1分钟的全场的风电功率 采样间隔为1分钟的全场的风电功率 采样间隔为1分钟的全场的风电功率

实测t时刻功率 风电功率持续分量

S级分量 均值指标 标准差指标 极差指标

4．模型的建立与求解

4.1 问题1的模型建立与求解

4.1.1 数据的选择与处理

附件中给出了某风电场2012年8月（30天）20台单机输出功率（采样间隔5秒，丢失数据统一赋值为null,单位kW）。所选择的五台风电机组分别为2、6、10、15、20，在30天范围内，选取30天内所有的采样间隔为5秒的数据进行以下建模的使用。在所给的数据中，有部分丢失数据，对丢失部分数据采取相用与该数据相邻最近的非丢失数据替代。

4.1.2 分离风电功率的s级分量的算法

采用滑动平均法分离s级风电功率。设滑动平均时段长度为N s，为了方便表达，设N为偶数，则t时刻风电的持续分量、s级分量可按下式计算：

pft=（Pt(N/21)+Pt(N/22)+?+Pt+?+Pt+N/2)/N

pst=pt-pftt=N/2,N/2+1,...,M-N/2式中：pt为实测的第t时刻的平均功率；pft为持续分量；pst为s级分量，pmt是叠加在持续分量上的变化量；t为测量点；M为测量点的总数。

滑动平均时段长度的选择具有一定随机性。长度太小，则风电功率的短时波动会反映在持续分量上；长度太大，则风电功率随时间变化的趋势会反映在 min 级分量上，使得 s级分量不再是随机变量。在该题中，即采样间隔为5秒时候，选取的滑动平均时段长度为180s。

图 1、2 给出了滑动窗口取180s时，第二台风电机组输出功率的持续分量和 min 级分量。

图1 输出功率持续分量

图2 输出功率s级分量

4.1.3 风电功率min级分量的概率密度函数

采用 Matlab 的概率密度拟合工具箱 dfittool对s级分量的概率密度函数进行拟合。分别选取工具箱中正态分布，逻辑斯特分布和t分布对风电功率min级分量的概率密度函数进行拟合。

图3给出的是第二台风电机组采用间隔5S的概率密度函数不同分布拟合结果对比：

图3 第二台风机间隔5s的概率密度的拟合结果

很显然，t分布的拟合效果是最好的。上述的分析均为第二台风电机组，由上述方法易知：

同样得到第六台风电机组的拟合结果，如图4：

图4 第六台风机间隔5s的概率密度的拟合结果

同样得到第十台风电机组的拟合结果，如图5：

图5 第十台风机间隔5s的概率密度的拟合结果

同样得到第十六台风电机组的拟合结果，如图6：

图6 第十六台风机间隔5s的概率密度的拟合结果

同样得到第二十台风电机组的拟合结果，如图7：

图7 第二十台风机间隔5s的概率密度的拟合结果

表1 5台风电机组三种不同概率分布的数值特征（采样间隔为5秒）

概率分布类型 数值特征 正态分布 机组2 -0.416 115.878 -3.875 55.2129 -6.625 48.6991 1.61523 机组6 机组10 机组16 6.024 130.815 0.595 58.5013 -2.101 17.825 0.6808 机组20 -0.4726 133.924 -2.323 62.8606 -3.226 44.937 1.27 μ σ 8.36009 -3.22573 136.189 0.977 62.388 -3.448 49.436 1.444 123.365 -5.192 57.745 -3.5066 45.5481 1.4073 逻辑斯特分布 μ σ μ t分布 σ ν 对三种分布进行T检验，检验用原始数据关于三种分布均值的假设是否成立。

基本使用格式为：

[h,sig,ci]=ttest(x,m,alpha)

我们选在α=0.01的显著水平下进行T检验。经检验，h=0所以，接受原假设，即证明数据已服从上述三种分布。再比较返回值sig（假设成立概率的大小）来确定风电功率的波动特性更符合以上哪种概率分布。检验的各项返回值详见表

表2 5台风电机组三种不同概率分布的检验（采样间隔为5秒）

概率分布类型 检验返回值 H 正态分布 SIG 逻辑斯特分布 H SIG H t分布 SIG 0.7865 0.7623 0.6111 0.7234 0.5532 0.3399 0 0.5360 0 0.1425 0 0.2265 0 0.4621 0 0.5312 0 0.6532 0 0.5535 0 0.4322 0 0.4232 0 机组2 0 机组6 0 机组10 0 机组16 0 机组20 0 由上表可知，三种概率分布中t分布的假设成功概率最大，且由图形观察拟合效果最好的为t分布，所以推荐使用t分布来表征风电功率波动特性的概率分布。

五台机组的数据分别进行三种概率分布的拟合时，整体效果相仿，都为t

分布的拟合效果最好，但由于风机所处位置的差异，不同机组相同分布的特征值仍存在明显差异。

至此可知，五台风电机组的风电功率波动概率分布均用t分布描述最好，尽管五台风电机组的风电功率的真实测量量的概率分布并不相同。

4.1.4 t分布下的风电机组风电功率波动特性 在该问中，选取第一天的第一台风电机组输出功率采样间隔为5秒的数据。采用t分布拟合出波动概率分布曲线。

在对风电功率波动概率分布参数进行统计分析之前，首先定义以下统计量指标用于分析风电功率数据样本特征量的分布情况。 （1）均值（Mean）

均值（Mean）表示一系列数据或统计总体的集中趋势或平均水平，计算公式如下：

1nX??Xi

ni?1 上式中(i=1,2,?,n)表示从总体 X 中抽取的一个样本， X 表示该样本的均值，n 表示样本数。本文中 X 则为实测风电功率数据。

均值能够反映所研究风电电源（风电机组、风电场或风电场群）某时刻或某时间段内输出功率的整体水平。 (2) 标准差（Standard Deviation）

标准差（Standard Deviation）表示了一组数据关于平均数的平均离散程度。本文将主要采用标准差描述样本数据的离散程度。

s=1(X∑n122In1i=1X)

式中，s 表示标准差。标准差越大，说明变量值之间的差异越大，距离平均数这个“中心”的离散趋势越大。标准差越大，表示风电功率样本数据之间的差异越大，各数据距离均值越远，均值的代表性越低，则风电波动比较剧烈。反之，标准差越小，表示样本数据之间的差异越小，各数据距离均值较近，则均值的代表性越高，风电波动相对比较平缓。 (3) 极差（Range）

极差（Range）是样本数据的最大值（Maximum）与最小值（Minimum）之间的绝对差，借以表明样本数据值最大可能的波动范围。计算公式如下，式中 R表示极差。

R=Xmax-Xmin

极差将用于风电功率波动的空间分布特性的分析中，即求同时刻下相同风电场内不同风电机组或风电场群内不同风电场之间输出功率的差值。极差越大说明功率数据越分散，被统计区域内各风电机组输出功率的空间分布相关性越差；反之相关性则越强。

通过上述方法，可共得到五台机组30个时段的概率分布参数，即150组数据。

表3 五台机组30个时段的概率分布参数

机组 时段(天数) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 时段(天数) 1 2 3 4 5 6 均值 4.1426 -0.7673 0.1899 2.3351 2.0694 -3.6512 -2.8329 -21.7457 5.3467 9.7515 4.5771 2.6476 2.7876 4.3689 -5.3997 -9.6234 0.1431 -4.7435 -3.1523 -9.3371 -1.6032 -8.8805 -7.1834 -12.2954 3.2464 1.9569 16.5535 -1.2501 15.2233 -0.6507 均值 4.7803 2.7146 -0.9650 0.5545 -0.3092 -1.4595 标准差 57.1279 39.5773 12.5485 15.4422 13.9691 37.9295 76.3130 68.1438 63.5174 94.5153 48.9215 23.3269 73.2756 75.8775 74.7355 88.4574 144.3715 90.2378 54.3754 107.8884 81.9721 108.9069 197.0531 154.8798 107.2319 38.5570 110.7141 38.7198 119.1819 94.2147 标准差 32.5676 42.0300 9.2699 14.6986 12.0851 39.1783 30天总体分布均值 30天总体分布标准差 2 -0.4795 87.4589 机组 30天总体分布均值 30天总体分布标准差 6 -0.9288 82.8754

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 3、 机组 时段(天数) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 4.7302 -8.3191 -6.3990 6.3826 -5.2175 3.6363 9.9168 3.5044 -7.5624 13.7695 -2.3298 -3.1493 -2.2000 -11.2439 2.0636 -13.7230 -15.4393 -16.5353 -0.5250 -3.5296 5.6815 0.1351 0.0245 6.4735 均值 -4.4745 -0.8481 0.2049 -0.1893 0.1411 -3.3398 -9.2959 2.7265 -5.5199 8.1452 -3.9804 -0.7173 -1.6021 11.3777 4.8086 91.8868 95.2402 67.6778 120.2004 57.5134 39.6505 70.1040 93.8765 59.7274 79.2633 130.4480 60.6418 26.8065 80.1258 63.9225 69.4686 119.8639 144.0850 135.0791 34.4678 131.4156 64.9127 135.2886 86.0377 标准差 53.1674 39.5176 13.7920 1.8325 18.8211 32.9480 82.0840 88.9990 73.0880 98.6828 44.4376 35.9733 61.9528 81.2741 83.7956 30天总体分布均值 30天总体分布标准差 10 0.9438 79.5388

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 4、 机组 时段(天数) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 -5.5251 5.4417 -7.4131 3.8433 3.9805 0.6915 -1.8394 10.8114 -11.8037 34.9562 0.3057 -0.6326 -2.1085 6.9274 -0.7896 均值 1.1971 2.7031 -1.9373 0.7078 0.9931 -1.6523 0 -0.0549 -5.6997 -3.6253 -1.2391 1.8035 -4.6834 -3.4840 -0.6928 5.8309 2.2126 13.1369 -2.9098 2.6443 5.0843 -0.7501 0 77.7762 132.6060 66.6469 30.1589 95.7894 56.4807 66.9326 132.0232 133.6904 119.3660 38.7708 105.1595 46.2167 133.6593 83.8706 标准差 48.7021 45.8409 8.7459 9.4016 13.9874 29.0328 0 12.2326 57.3460 59.4385 41.5238 22.9669 64.3989 88.3657 105.8419 99.6034 245.9042 73.5984 33.8241 120.7110 83.6918 84.4499 0 30天总体分布均值 30天总体分布标准差 16 -0.8613 87.4880

24 25 26 27 28 29 30 5、 机组 时段(天数) -17.9420 -4.9981 2.3399 -15.7589 0.1372 -1.5342 -1.4427 均值 138.0016 117.5951 31.5088 145.6065 53.8947 138.5773 89.4382 标准差 30天总体分布均值 30天总体分布标准差 1 3.6282 71.1311 2 -3.1865 31.5080 3 -0.6565 12.7498 4 0.1032 19.4015 5 -0.1559 11.8888 6 -4.7072 32.4014 7 0 0 8 0.5511 11.9438 9 1.8917 64.8698 10 8.5120 99.4846 11 0.1177 32.8580 12 -4.1967 28.1317 13 -9.8824 62.1217 14 14.2168 75.4353 15 -4.9495 59.4061 20 1.8955 84.7374 16 0.4255 72.6121 17 27.7353 178.7244 18 4.0003 89.3402 19 3.8633 44.4532 20 1.3890 159.4984 21 -7.9619 56.2652 22 8.1435 82.0108 23 15.9488 139.3419 24 1.6939 124.6388 25 13.3302 93.8867 26 0.3720 49.6519 27 -8.8186 126.3508 28 -2.4818 49.6321 29 1.0606 163.4794 30 -4.6591 80.9161 从上面的5个表中，我们可以分析得知：随机抽取的5个风机，30天的总体的均值和标准差波动不大，说明风机的输出功率的波动并不大。而且彼此之间

的波动特性都比较类似，说明时间尺度对单机及风电场的功率波动都存在某种相同的联系。

4.2 问题2的模型建立与求解

4.2.1 数据的选择与处理

从上述所选择的五台风电机组发电功率数据中，在30天范围内，选取30天内所有的采样间隔为1分钟的数据进行以下建模的使用。 4.2.2 风电功率min级分量的概率密度函数

采用 Matlab 的概率密度拟合工具箱 dfittool对 min 级分量的概率密度函数进行拟合。分别选取工具箱中正态分布，逻辑斯特分布和t分布对风电功率min级分量的概率密度函数进行拟合。

图8给出的是第二台风电机组采用间隔1分钟的概率密度函数不同分布拟合结果对比：

图8 第二台风机间隔1分钟概率密度的拟合结果

很显然，t分布的拟合效果是最好的。上述分析均为第二台风电机组，由上述方法易知：

同样得到第六台风电机组的拟合结果，如图9：

图9 第六台风机间隔1分钟概率密度的拟合结果

同样得到第十台风电机组的拟合结果，如图10：

图10 第十台风机间隔1分钟概率密度的拟合结果

同样得到第十六台风电机组的拟合结果，如图11：

图11第十六台风机间隔1分钟概率密度的拟合结果

同样得到第二十台风电机组的拟合结果，如图12：

图12 第二十台风机间隔1分钟概率密度的拟合结果

表3 5台风电机组三种不同概率分布的数值特征（采样间隔为一分钟）

概率分布类型 数值特征 机组2 机组6 机组10 3.05552 160.723 -4.11708 80.7104 机组16 机组20 μ 正态分布 -1.57212 1.95599 159.355 -12.7357 69.3534 148.615 -51852 73.4717 7.57152 -0.533055 178.785 3.8842 81.5936 190.656 -7.3196 90.1479 σ 逻辑斯特分布 μ σ μ t分布 -18.2938 -10.1568 -12.5281 -7.53064 -18.6444 44.8666 1.221 80.0829 2.20797 78.7844 1.82619 55.4396 1.2208 55.5415 1.10882 σ ν 当采样频率上升为分钟级，再次对三种分布进行T检验，检验用原始数据关于三种分布均值的假设是否成立。

基本使用格式为：

[h,sig,ci]=ttest(x,m,alpha)

我们选在α=0.01的显著水平下进行T检验。经检验，h?0所以，接受原假设，即证明数据已服从上述三种分布。再比较返回值sig（假设成立概率的大小）来确定风电功率的波动特性更符合以上哪种概率分布。检验的各项返回值详见表 概率分布类型 检验返回值 H 正态分布 SIG 逻辑斯特分布 H SIG H t分布 SIG 0.6765 0.5423 0.9711 0.8234 0.8532 0.4532 0 0.5360 0 0.2325 0 0.5465 0 0.5421 0 0.5352 0 0.4532 0 0.4535 0 0.432 0 0.7532 0 机组2 0 机组6 0 机组10 0 机组16 0 机组20 0 经T检验，分析得知，当采样频率上升为分钟级，风电波动特性的概率分布

仍为t分布最优。

至此可知，五台风电机组的风电功率波动概率分布均用t分布描述最好，尽管五台风电机组的风电功率的真实测量量的概率分布并不相同。

4.3 问题3的模型建立与求解

对于单台风电机组，随时间尺度的增大（5s-1min），对本身所给的数据选取中存在丢失，每次丢失了16个数据。

参数 均值 标准差 5s 0.4795 87.4589 1min 1.5721 159.3550 5min 37.6439 240.1917 15min 72.2861 311.4779 风电机组输出功率的均值由0.4795 增大到 1.5721，标准差由 87.4589 增加至 159.3550，并且发现所有的风电机组随着时间尺度的增加，其均值及标准差都相应的增加，顾我们可以考虑采用均值及标准差来进行度量风电功率。 综上对风电机组输出功率波动特性的分析可看出，单台风电机组输出功率及其波动量在不同时间尺度下的分布概率是不同的，随采样时间间隔的增大上述特征量都呈现变大的趋势。然而考虑到单台机组容量相对较小，无论何种时间尺度下单台风电机组引起的功率波动都不会对区域电网造成明显影响。

4.4 问题4的模型建立与求解

4.4.1 数据的选择与处理

选取全场20台风电机组发电功率数据，在30天范围内，选取30天内所有的采样间隔为1分钟的总功率序列P∑m(tk)，采样间隔为5分钟的总功率序列

P∑(tk)，采样间隔为15分钟的总功率序列P?5m15m(tk)，进行以下建模的使用。并

且将滑动平移时间改为2秒。

4.4.2 风电场输出功率波动特性分析

易得该风电场 30 天实测功率曲线，如下图，与图 单台风电机组输出功率波动情况相比风电场输出功率稍有缓和。

以下将分析不同采样时间尺度下（仍与单台机组分析类似，取 1min,5min,15min）风电场功率数据的统计指标。

当采样间隔为1分钟时，由MATLAB可以得出此时空特性下的输出功率波动概率分布图

图13 1分钟的输出功率波动分布

当采样间隔为5分钟时，由MATLAB可以得出此时空特性下的输出功率波动概率分布图

图14 5分钟的输出功率波动分布

当采样间隔为15分钟时，由MATLAB可以得出此时空特性下的输出功率波动概率分布图

图15 15分钟的输出功率波动分布

且可得：

表4 间隔1、5、15分钟的不同输出功率波动分布数值特征： 概率分布类型 数值特征 1分钟 1.3766 479.497 1.5 159.506 -1.5287 102.947 1.25805 5分钟 0.850 618.225 -3.80 245.688 -3.80792 211.639 1.641 15分钟 496.306 ?9.5?10?15 μ 正态分布 σ μ 逻辑斯特分布 302.354 3.050 -30.2095 417.599 2.12 σ μ t分布 σ ν 及概率分布参数： 表5 间隔1、5、15分钟的不同输出功率波动分布的概率分布参数

概率分布参数 均值 标准差 极差 1分钟 1.3766 479.4969 29675 5分钟 0.8505 618.2253 22382 15分钟 9.1767 844.8004 17921 对于风电场，随时间尺度的增大（1min-5min），对本身所给的数据选取中存在丢失，每次丢失了16个数据。

风电场输出功率的均值由1.3766 增大到 9.1767，标准差由 479.4969 增加至 844.8004，并且风电场随着时间尺度的增加，其均值及标准差都相应的增加，顾我们可以考虑采用均值及标准差来进行度量风电功率。

综上对风电场输出功率波动特性的分析可看出，风电场输出功率及其波动量在不同时间尺度下的分布概率是不同的，随采样时间间隔的增大上述特征量都呈现变大的趋势。然而，由于风电场装机容量远大于单台机组，甚至某些风电场达到几百兆瓦容量水平，其风电功率输出水平已不能忽视。

4.5 问题5的模型建立与求解 4.5.1模型的概述 分别采用P?5m(tk)和P?(tk)作为样本来预测未来4小时（每15分钟一

m个点）风电场的总功率，这里采样灰色模型进行拟合。

灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。

GM(1,1)的定义

设x(0)为n

(0)(0)(0)(0)x=(x(1),x(2), ?,x(n))，x个元素的数列

(0)的

AGO生成数列则

x(0)=[724.57746.62778.27?954.28995.011037.2]，

其中

(1)(0)x(k)?x?(i)(k?1,2,?,n) （2）

i?1k则定义x(1)的灰导数为

d(k)?x(0)(k)?x(1)(k)?x(1)(k?1) (3)

令z(1)为数列x(1)的紧邻均值数列，即

z(1)(k)?0.5x(1)(k)?0.5x(1)(k?1),k?2,3,?，n， (4)

(1)(1)(1)(1)z?(z(2),z(3),?,z(n)).于是定义GM(1,1)灰微分方程模型为 则

d(k)?az(1)(k)?b (5)

即

x(0)(k)?az(0)(k)?b(a) (6)

其中x(0)(k)称为灰导数，a称为发展系数，z(1)(k)称为白化背景值，b称为灰作用量.

将时刻k?2,3,?,n代入(a)式中有

?x(0)(2)?az(1)(2)?b?(0)(1)?x(3)?az(3)?b? (7)

????(0)(1)??x(n)?az(n)?b?z(1)(2)??z(1)(3)(0)(0)(0)TT令Y?(x(2),x(3),?x(n)),u?(a,b)，B?????z(1)(n)??1??1?称Y为数据?，?1???向量，B为数据矩阵，u为参数向量，则GM(1,1)模型可以表示为矩阵方程

Y?Bu.

由最小二乘法可以求得

?)T=(BTB)1BTY (8) ?=(a?,buGM(1,1)的白化型

对于GM(1,1)的灰微分方程（7），如果将x(0)(k)的时刻k?2,3,?,n视为连续的变量t，则数列x(1)就可以视为时间t的函数，记为x(1)?x(1)(t)，并让灰导数

dx(1)，背景值z(1)(k)对应于x(1)(t).于是得到GM(1,1)的灰微x(k)对应于导数dt(0)分方程对应的白微分方程为

dx(1)?ax(1)?b (9) dt称之为GM(1,1)的白化型.

4.5.2 模型的建立

此预测模型是拟合参数模型，通过原始数据累加生成，得到规律性较强的序列，用函数曲线拟合得到预测值. 建立过程如下：

1) 设原始数据序列X(0)有n个观察值，X(0)?{X(0)(1),X(0)(2),?,X(0)(n)}，通过累加生成新序列X(0)?{X(1)(1),X(1)(2),?,X(1)(n)}，利用新生成的序列X(1)拟合函数曲线.

2) 利用拟合出的函数求出新生序列X(1)的预测值序列X(1).

3) 利用X(0)(k)?X(1)(k)?X(1)(k?1)累减还原，得到灰色预测值序列（共n+m个，m个未来预测值）.将序列X(0)X(0)?{X(0)(1),X(0)(2),?,X(0)(n?m)}分为Y0和Z0，其中Y0反映X(0)的确定性增长趋势，Z0反映X(0)的平稳周期变化趋势.

4) 对X(0)序列的确定增长趋势进行预测. 4.5.3 模型的求解 整理P?5m(tk)作为样本，题目要求预测未来4小时风电场的总功率，且每15

分钟一个点，这里我们每组数据序列都选取16个点作为样本。又因为要求滚动

最少7天的预测值，我们通过MATLAB编程反复使用上一次的预测值进行下一次的预测，从而获得预测的数据及相关参数：

表6 P?时间 (min) P?(tk) 5m5m(tk)数据序列

20 3748 60 3778 25 3853 65 3460 30 3726 70 2963 (0)0 3739 40 3705 5 3935 45 4045 10 3610 50 4031 15 3681 55 4097 35 3728 75 2566 时间 (min) P?(tk) 5m 根据上述数据建立含有16个观察值的原始数据序列X

：

X(0)=[373 3993 3561 3068 3174 3885 3372 3672 3870 4504 4503 4109 3777 3846 2096 2356]6

使用Matlab软件对X(0)进行一次累加，得到新数列X(1)。

表7 GM(1,1)算法拟合值及误差

序号 时间 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 模型值 3739 3935 3610 3739 3739 3739 3739 3739 3739 3739 3739 3739 3739 3739 3739 3739 残差 0 -36.3490 -314.9165 -198.0269 -85.6738 64.1490 -18.5522 27.2286 47.4975 430.2606 458.5236 566.2926 288.5732 11.3711 -445.3080 -802.4585 相对误差 0 0.92% 8.72% 5.38% 2.29% 1.66% 0.50% 0.73% 1.28% 10.64% 11.37% 13.82% 7.64% 0.33% 15.03% 31.27% 级比偏差 0.0609 -0.0773 -0.0308 0.0294 0.0386 -0.0220 0.0122 0.0056 0.0948 0.0083 0.0276 -0.0718 -0.0791 -0.1541 -0.1412 X(1)(2) X(1)(3) X(1)(4) X(1)(5) X(1)(6) X(1)(7) X(1)(8) X(1)(9) X(1)(10) X(1)(11) X(1)(11) X(1)(11) X(1)(11) X(1)(11) X(1)(11) X(1)(11) 由残差、相对误差、级比偏差可知此模型精度较高，可用于预测未来4小时（每15分钟一个点）风电场的总功率，滚动持续7天，这里只给出第一个时间段（未来四个小时）的预测值，其余预测值在附件中给出。

序号 时间 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 预测值 3610 4026.6 3603.5 3224.9 2886.1 2582.9 2311.5 2068.7 1851.3 1656.8 1482.7 1327 1187.5 1062.8 951.12 851.19 X(1)(2) X(1)(3) X(1)(4) X(1)(5) X(1)(6) X(1)(7) X(1)(8) X(1)(9) X(1)(10) X(1)(11) X(1)(12) X(1)(13) X(1)(14) X(1)(15) X(1)(16) X(1)(17) 整理P?5m(tk)作为样本，题目要求预测未来4小时风电场的总功率，且每15

分钟一个点，这里我们每组数据序列都选取16个点作为样本。又因为要求滚动最少7天的预测值，我们通过MATLAB编程反复使用上一次的预测值进行下一次的预测，从而获得预测的数据及相关参数：

表8 P?时间 (min) P?15m15m(tk)数据序列

60 2963 180 1366 75 2230 195 1060 90 1781 210 1335 105 1601 235 1064 0 3610 120 1692 15 3853 135 1627 30 3705 150 1470 45 4097 165 1450 (tk) 时间 (min) P?15m(tk) 根据上述数据建立含有10个观察值的原始数据序列X(0)：

X(0)=[3610 ...106]4 使用Matlab软件对X(0)进行一次累加，得到新数列X(1)。

表9 GM(1,1)算法拟合值及误差

序号 时间 0 15 30 45 60 75 90 105 120 135 150 模型值 3610 3853 3705 4097 2963 2230 1781 1601 1692 1627 1470 残差 0 0.045053 0.027385 0.21285 0.025949 0.15825 0.29788 0.29211 0.094162 0.018325 0.008669 相对误差 0 -173.59 101.46 872.06 76.886 -352.89 -530.52 -467.66 -159.32 -29.815 -12.744 级比偏差 0.1616 0.0694 0.1908 -0.2373 -0.1890 -0.1204 0.0046 0.1533 0.0694 0.0096 X(1)(2) X(1)(3) X(1)(4) X(1)(5) X(1)(6) X(1)(7) X(1)(8) X(1)(9) X(1)(10) X(1)(11) X(1)(12)

X(1)(13) X(1)(14) X(1)(15) X(1)(16) X(1)(17) 165 180 195 210 235 1450 1366 1060 1335 1064 0.084854 0.13064 0.00262 0.28755 0.20001 123.04 178.45 -2.7776 383.88 212.81 0.0928 0.0501 -0.1532 0.2895 -0.1228 由残差、相对误差、级比偏差可知此模型精度较高，可用于预测未来4小时（每15分钟一个点）风电场的总功率，滚动持续7天，这里只给出第一时段（未来四小时）的预测值，其余预测值见附录。 序号 时间 0 15 30 45 60 75 90 105 120 135 150 165 180 195 预测值 3610 4026.6 3603.5 3224.9 2886.1 2582.9 2311.5 2068.7 1851.3 1656.8 1482.7 1327 1187.5 1062.8 X(1)(2) X(1)(3) X(1)(4) X(1)(5) X(1)(6) X(1)(7) X(1)(8) X(1)(9) X(1)(10) X(1)(11) X(1)(12) X(1)(13) X(1)(14) X(1)(15)

X(1)(16) X(1)(17) 210 235 951.12 851.19 分析两种方式的预测误差可知：采样间隔为5分钟的风电场总功率通过后续数据作为实际风电功率用于检验的预测误差要比采样间隔为15分钟的预测误差要小。其中全部数据的预测误差（残差、相对误差、级比偏差）见附录。

4.6 问题6的模型建立与求解

通过前面内容的分析可知，当取不同时间尺度对单个风电机和风电机群输出功率进行特征量分析时，风电机群的功率（即风电场输出功率）与单台风电机组输出功率的特性存在时间分布上较为明显的差异。

4.6.1 单台风电机组功率在时序上的表现

表10 第二台风机三十天内的参数

参数 均值 标准差 5s 0.4795 87.4589 1min 1.5721 159.3550 5min 37.6439 240.1917 15min 72.2861 311.4779 上表给出第二台风电机组三十天内不同时间尺度（5s、1min、10min、15min）的输出功率的波动概率分布参数。可以看出，时间尺度增大，其概率分布数值特征也不断变大，其输出功率的时序分布特性更加明显，其功率波动随时间尺度的增加更为缓和。

随时间尺度的增大（5s-15min），风电机组输出功率的均值由0.4795增大到 72.2861，标准差由 87.4589 增加至 311.4779，意味着风电功率输出水平的升高，且机组输出功率的波动范围随之增大。风电机组输出功率波动量均值基本为 0，而其标准差则随时间尺度增大而增大，即功率波动幅度越大，这是符合统计规律的。综上对风电机组输出功率波动特性的分析可看出，单台风电机组输出功率及其波动量在不同时间尺度下的分布概率是不同的，随采样时间间隔的增大上述特征量都呈现变大的趋势。然而考虑到单台机组容量相对较小，无论何种时间尺度下单台风电机组引起的功率波动都不会对区域电网造成明显影响。

4.6.2 风电场功率在时序上的表现 概率分布参数 均值 标准差 极差 1分钟 1.3766 479.4969 29675 5分钟 0.8505 618.2253 22382 15分钟 9.1767 844.8004 17921

表给出风电场（20台风电机组）三十天内不同时间尺度（5s、1min、10min、15min）的输出功率的波动概率分布参数。可以看出，具有更多风机台数的风电场，其输出功率的时序分布特性更加明显。

其标准差则随时间尺度增大而增大，当时间间隔分别为 1min 和 15min 时，标准差分别为 479.4969 和 844.8004。风电场功率波动基本在分钟级时间尺度内，更短（s 级）或更长时间尺度上（数 hours）风电场功率波动的分布情况较少。各机组输出功率均值与标准差的变化情况。不同采样时刻下风电机组输出功率的均值随机组台数增加而增大，但台数增多意味着机组分布面积越广，则各机组输出功率的极差也越大，因而标准差也越大。相同时间尺度下风电场风电功率的波动特征量要较单台风电机组的更为平缓些，这是由于风电场内机组分布位置不同导致各机组输出功率有明显差异，当汇聚后风电场总输出功率则平抑了这些差异（即机组输出功率具有互补性）。然而，由于风电场装机容量远大于单台机组，甚至某些风电场达到几百兆瓦容量水平，其风电功率输出水平已不能忽视。

4.6.3 两者差别及分析比较

由上述两者表现可以清楚看到，通过表、图、数据、理论分析有一个较为清晰的联系。由此可知在分析时序波动特性方面，可以大胆的利用概率分布数值特征作为指标来表示，从而使得时序波动特性更加的直观化，定量的进行分析处理。但是，概率分布数值特征仅仅只是一个方面，有较大的局限性。虽然从时序上能够对风电场引起的功率变化进行初步分析，但却难以把握较长时间段内功率变化的趋势性特征。对于风电这样的依赖于自然力的电源，若要了解风电场群汇聚过程中功率波动变化的规律，可以用风电场群的年持续出力曲线来描述，而不必特别关注具体的时序曲线。

我们还可以换个角度进行分析，可以发现，增加机组数目，扩大面积域，即在空间分布上，风电功率波动特性与之也有较为清晰的联系，这一点也说明了时序的局限性。

4.7 问题7的模型建立与求解

通过上述对机组和全场风电功率波动的分析，我们可以知道： 1）单台风电机组功率输出水平及其波动量随采样时间间隔的增大呈现变大的趋势，其波动特征量主要以秒级为主，但由于基值较小，对电网影响不大； 2）2）风电场风电功率的波动特征量较单台风电机组的更为平缓，输出功率波动将主要分布在分钟级，将对电网频率调整带来一定影响；

3）提出风电功率波动特性时序分布的概念，并指出随机组时间尺度的增加，风电场输出功率波动的时序分布特性将愈加明显。

在克服风电波动对电网运行的不利影响上，我们可以利用风电功率波动特性，

实例构建：

以下内容将基于我国东北吉林省西部百万千瓦风电场群实测风电功率数据，

采用上述统计量指标进行风电功率波动特性的分析。该风电场群由 8 个风电场汇聚而成，总装机容量为 1000MW，覆盖地域约200km2。全部风电功率经双回 220kV 线路送入主网。各风电场输出功率数据采集时间段为 2009 年 1 月 1 日 00:00:00-2009 年 12 月 31 日 24:00:00，采集间隔为5 min，总有效数据为 105120 个。单台风电机组数据采集时间段为 2009 年 11 月1 日 00:00:00-2009 年 11 月 20 日 24:00:00，采样间隔 6s，总有效数据 288000 个。利用该文中同样的方法进行分析风电功率波动对电网运行所关注的调频、调峰等过程的影响，我们可以知道：

（1）在较短时间尺度内(sec、min 级)风电场群输出功率的波动通常较小，而较长时间尺度(hour 级及以上)下则会较大；随着时间尺度的增大，风电输出功率波动的分布具有一定的趋势性。尽管较短时间尺度下风电功率变化量较小，但当风电并网规模很大时，此种时间尺度下的风电波动将影响电网有功功率调整，如频率实时调整、发电日前调度等。目前，电网调度或规划部门已能够在某种程度上对较长时期内(数天)风电功率波动进行预测，但对较短时间段内(sec、min 级)的风电功率预测精度还不足以使电网应对风电功率的随机快速波动。

（2）通过额外的了解，我们可以知道风电场群覆盖区域广阔，不同位置风机所接受的风能会存在空间上的差异，受各种因素的影响，单台风电机组与单个风电场乃至与整个风电场群输出功率的波动特性存在较大差异。对风电场群而言，基本不会出现由于无风而使电网内风电机组同时停运的情况，也不会出现同一电网所有风电场同时满发的情况；风电功率的波动随空间分布尺度的增大而趋于缓和。因此，伴随着从风电机组、风电场向风电场群的汇聚，风电最大发电功率的概率分布参数呈递减趋势。接入同一并网点的风电场数目越多、面积越广，该风电场群对电网调峰容量的相对需求反而是降低的。

5.模型结果的分析与检验

对于整个模型的结果，我们可以知道：单台风电机组功率输出水平及其波动量随采样时间间隔的增大呈现变大的趋势。风电场风电功率的波动特征量较单台风电机组的更为平缓，输出功率波动将主要分布在分钟级，将对电网频率调整带来一定影响。在较短时间尺度内(sec、min 级)风电场群输出功率的波动通常较小，而较长时间尺度(hour 级及以上)下则会较大；随着时间尺度的增大，风电输出功率波动的分布具有一定的趋势性。

而这一结果与相关事实及部分文献资料记载相符合。在对每个问题的结果分析后，我们都进行了相关的检验。并且，对照我国东北吉林省西部百万千瓦风电场群实测风电功率的实例中，也检验出该模型的正确性及合理性。

6．

6．1 模型的推广

模型的推广与改进方向

该模型主要针对时间尺度不同对单台风电机组及全场风电机组的功率波动进行分析，及在对预测模型上，采用灰色预测模型及滚动预测。重点对于前

者有较多的处理与分析，对于该类风电机组的波动特性问题均有较好的方向。不过，主要在时间尺度上，其实在空间分布上，对风电机组的波动特性也具有较大的影响，可以考虑在空间分布上继续推广。

6．2 模型的改进

该模型是利用大量的数据建立起来的，在对数据的处理上必然存在各种的差异及问题，可以对数据处理上考虑有所改进。其次，对概率分布拟合，考虑到MATLAB程序实现的难度，使用了工具箱，如果在对拟合程序上有所改进，将对整体的分析能够更加的深入。

7．模型的优缺点

7．1 模型的优点

尽管数据处理上或多或少存在差异，但是本身对如此大量数据的合理处理，这就是个优点。再者，整体模型在考虑时间尺度上较为深入，对时序波动特性着重的分析。另外，预测模型中，我们采用了灰色预测模型，采用滚动预测，并进行了所有预测值的预测误差计算，进一步的增加了预测模型的优点。

7．2 模型的缺点

对概率分布的拟合采用了直接用工具箱拟合，对于其数据本质并未深入分析探讨。

参考文献

[1] 李欣然，贺仁睦，章健，等．负荷特性对电力系统静态电压稳定性的影响静态电压稳定性广义实用判据[J]．中国电机工程学报．1999，19 (4)：23~30

[2] 张伶俐，周文，章健，等．面向综合的电力负荷特性建模[J]．中国电机 工程学报．1999，19 (9)：71~75

[3] 崔杨，穆钢．大规模风电场群联网的源网协调性研究[J]．华北电力大学博士学位论文．2011，3(2)：28~41

[4] 唐丽芳、贾冬青、孟庆鹏. 《用MATLAB实现灰色预测GM(1,1)模型》 [J]. 沧州师范专科学校学报, 2008, 24:35.

[5] 司守奎, 孙玺菁. 数学建模算法与应用 [J]. 灰色预测模型, 2011, 351—354.

附 录

附录1-源程序

1.%%%%%%%%%%一台机组三十天数据-5s采样

clear all close all clc

a=load('m1.dat'); s=a(:,1); for i=2:10

s=[s;a(:,i)]; end

M=length(s); %N取代null的前后2N个数据

for i=1:M

if s(i)==-1000 j=1;

for j=1:M/2 if i-j>0 & s(i-j)~=-1000

s(i)=s(i-j); break;

elseif i+j

s(i)=s(i+j); break; end end end end

plot(s)%%%%%%%%%%一个月的曲线

T=length(s);N=180; % N滑动平均个数 2N+1个

ysum=0;k=1; for i=N+1:N:T-N for j=-N:N

ysum=ysum+s(i+j,1); end

yk(k)=ysum/(2*N+1); ym(k)=s(i,1)-yk(k); ysum=0;k=k+1; end figure

plot(yk) figure plot(ym)

n=length(ym); yav=sum(ym)/n

ys=sqrt(sum((ym-yav).^2)/(n-1)) ymaxmin=max(ym)-min(ym)

2.%%%%%%%%%%一台机组三十天数据-分钟采样数据

clear all close all clc

a=load('m2.dat');%第几台机组 s=a(:,1); for i=2:10

s=[s;a(:,i)]; end

plot(s)%%%%%%%%%%一个月的曲线

NL=180;Nt=length(s)/NL;k=1; %NL 12*5=60s即一分钟采样

for i=(1:Nt)*NL s1(k)=s(i); k=k+1; end

M=length(s); %N取代null的前后2N个数据

for i=1:M

if s(i)==-1000 j=1;

for j=1:M/2 if i-j>0 & s(i-j)~=-1000

s(i)=s(i-j); break;

elseif i+j

s(i)=s(i+j); break; end end end end

T=length(s1);N=180; % N滑动

平均个数 2N+1个

ysum=0;k=1; for i=N+1:N:T-N for j=-N:N

ysum=ysum+s1(i+j); end

yk(k)=ysum/(2*N+1); ym(k)=s1(i)-yk(k); ysum=0;k=k+1; end figure plot(yk) figure

plot(ym) n=length(ym); yav=sum(ym)/n

ys=sqrt(sum((ym-yav).^2)/(n-1)) ymaxmin=max(ym)-min(ym)

3.%%%%%%%%%%一台机组一天数据-5s采样

clear all close all clc

a=load('m2.dat'); s=a(:,1); for i=2:10

s=[s;a(:,i)]; end

M=length(s); %取代null for i=1:M

if s(i)==-1000 j=1;

for j=1:M/2 if i-j>0 & s(i-j)~=-1000

s(i)=s(i-j); break;

elseif i+j

s(i)=s(i+j); break; end end end

end

plot(s)%%%%%%%%%%一个月的曲线

Day=30;%第10天 Nt=length(s)/30;k=1; for i=(1:Nt)+(Day-1)*Nt s1(k)=s(i); k=k+1; end

T=length(s1);N=180; % N滑动平均个数 2N+1个

ysum=0;k=1; for i=N+1:N:T-N for j=-N:N

ysum=ysum+s1(i+j); end

yk(k)=ysum/(2*N+1); ym(k)=s1(i)-yk(k); ysum=0;k=k+1; end figure plot(yk) figure plot(ym)

n=length(ym); yav=sum(ym)/n

ys=sqrt(sum((ym-yav).^2)/(n-1)) ymaxmin=max(ym)-min(ym)

4.%%%%%%%%%%%整个风场20台机组-分钟采样

q=zeros(518400,1); for i=1:20 t=[];

num2str(i) t=strcat(t,'.dat')

t1=strcat('load(''m',num2str(i));

t2=strcat(t1,'.dat'')') a=eval(t2); s=a(:,1); for i=2:10 s=[s;a(:,i)]; end

M=length(s);N=2; %N取

代null的前后2N个数据

for i=1:M

if (s(i)==-1000)&(i-N>0) msum=0; for j=-N:N if j~=0 msum=msum+s(i+j);

end

s(i)=msum/2/N; end elseif (s(i)==-1000)&(i-N<=0)

msum=0; for j=1:N msum=msum+s(i+j);

end

s(i)=msum/N; elseif (s(i)==-1000)&(i+N>M)

msum=0; for j=-N:1 msum=msum+s(i+j);

end

s(i)=msum/N; end end q=q+s; end plot(q);

NL=12;Nt=length(q)/NL;k=1;%NL12*5=60s一分钟

for i=(1:Nt)*NL s1(k)=q(i); k=k+1; end

T=length(s1);N=2; % N滑动平均个数 2N+1个

ysum=0;k=1; for i=N+1:N:T-N for j=-N:N

ysum=ysum+s1(i+j);

end

yk(k)=ysum/(2*N+1); ym(k)=s1(i)-yk(k); ysum=0;k=k+1; end figure plot(yk) figure plot(ym)

n=length(ym) yav=sum(ym)/n

ys=sqrt(sum((ym-yav).^2)/(n-1)) ymaxmin=max(ym)-min(ym) %%灰色预测模型（滚动预测未来7天数据）

clc,clear Da=42;

x0=[3610 3853 3705 4097 2963 2230 1781 1601 1692 1627 1470 1450 1366 1060 1335 1064]';

yz=zeros(length(x0),Da);yep=yz;ydelt=yz;yrho=zeros(length(x0)-1,Da);

for i=1:Da

n=length(x0);

lamda=x0(1:n-1)./x0(2:n); range=min(max(lamda')); x1=cumsum(x0);

B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];

Y=x0(2:n); u=B\\Y;

x=dsolve('Dx+a*x=b','x(0)=x0');

x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)});

yuce1=subs(x,'t',[0:n-1]); y=vpa(x,6);

yuce=[x0(1),diff(yuce1)]; epsilon=x0'-yuce;

delta=abs(epsilon./x0');

rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda';

yz(:,i)=yuce'; yep(:,i)=epsilon'; ydelt(:,i)=delta'; yrho(:,i)=rho'; x0=yuce'; end

本文来源：https://www.bwwdw.com/article/ek8g.html