数学建模论文题目：A题风电功率预测参赛队员信息：队员1 队员2 队员3 姓名张旭汪宝施灵卫专业物理学物理学电气工程及其自动化邮箱1786428369@q 1026253546@q404842219@qq.com风电功率预测问题摘要：风电接入电网时，大幅度的风电功率变化对电网的功率平衡调节会产生不利影响.对风电功率进行准确预测，电力调度部门就可以根据预测的功率安排调度计划，保证电网功率平衡和运行安全.通过建立数学模型解决风电功率预测问题.首先对P A ,P B ,P C ,P D ,P 58进行统计分析，分别求出A 、B 、C 、D 单机每天的功率和，用MATLAB 分别绘制出单机每天的总功率和随天数变化的四组曲线，四组曲线变化趋势相同，具有相同的周期T=7天.比较P A P B ,P C ,P D ,P 58随天数变化曲线，以T=7天为一个周期，建立前一个周期内功率变化函数模型，预测下一个周期内的功率.建立三种函数模型.模型一：利用MATLAB 对一周内七天的曲线进行拟合，得到一个三次曲线 = ，以7天为周期，预测下个周期功率.模型二:建立灰色系统模型.选一周内7天的28个数据进行累加得到新的数据序列 ，j=28利用灰色系统理论进行求解 ，得到七天的功率满足的函数关系，根据函数变化趋势预测功率. 模型三：建立自回归移动平均模型ARIMA(p,d,q)，通过 作图发现时间序列是非平稳的，采取适当差分则会显示出平稳序列的性质，对序列分析求解得到参数，最后得到线性函数 根据函数变化趋势，预测下一周内的功率变化.关键词：时间序列 灰色系统 非线性拟合 ARIMA(p,d,q)d ct b a t t +++23()t f ()})28(,),2(),1({)0()0()0(0x x x x =()()()()j j jx x ∑=101()()t x 0()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ332211（一）问题重述日前预测是预测明日24小时96个时点的风电功率数值.实时预测是滚动地预测每个时点未来4小时内的16个时点的风电功率数值.某风电场由58台风电机组构成.(每十五分钟预测一个点)问题1：风电功率实时预测及误差分析。

具体要求：1）采用不少于三种预测方法（至少选择一种时间序列分析类的预测方法）；2）预测量：a．P A, P B, P C, P D；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：试分析风电机组的汇聚对于预测结果误差的影响。

问题3：进一步提高风电功率实时预测精度的探索。

（二）问题分析通过分别对P A,P B,P C,P D,P58每天功率曲线进行比较，得到P A,P B,P C,P D,P58具有相同的周期T=7；对前一个周期内的数据进行分析，预测下一个周期内的功率变化曲线.问题一：利用上一个周期内的数据建立回归曲线模型，灰色系统模型，ARIMA(p,d,q) 三种模型，对功率进行预测并计算精确度.问题二：由于预测每个单机时存在一定误差，所以当单机汇聚时会使误差增大.问题三：由于功率受到风速、风向、空气密度等因素的影响，所以功率变化比较复杂. 为提高预测精度，利用一个周期内的数据均值对ARIMA(p,d,q)求解.（三）模型假设1.各个单机功率互不影响2.每个单机运行的外部条件相同3.功率变化随天数存在周期4.各个单机的周期变化天数相同（四）符号说明T 周期 t时刻 a 、b 、c 、d 三次曲线参数为某天的数据序列 ()()t x 0为t 时刻功率为将某天数据累加得到新序列为自回归系数 为移动平均系数（五）模型建立与求解分别对A 、B 、C 、D 单机每天的功率和，画出功率和随天数变化的曲线，A 、B 、C 、D 的图像基本相同，下图为A 的功率和随天数变化的曲线θθθp⋅⋅⋅21,()x0()x1ϕϕϕp ⋅⋅⋅21,每天的总功率随天数变化具有周期性，对图形进行分析可以得到周期为7天.画出每天P A ,P B ,P C ,P D ,P 58随时间的变化的曲线，也以7天为一个周期，建立数学模型模拟前一个周期的功率变化曲线，根据曲线的变化趋势，预测下一个周期的功率变化曲线.其中模型一直接用的函数的周期性.对模型二,将时间t 代入函数表达式，得到功率.模型三与模型二相同.模型一：回归曲线模型对P A 进行拟合：拟合三次回归模型，该回归模型的回归系数为 为了方便引入矩阵写成矩阵利用最小二乘法对 进行估计，即使 的值最小.P A的回归系数的估计值 5月31号部分真实值与预测值真实249 107.0 49.8 56.5 155.2 41.9 385.6 11.6 153.1 预测294 301 307 312 317 322 325 329 332 6月1号到6月6号部分真实值与预测值真实25.2 152.5 66.6 343.9 0.37 187.3 123.9 18.7 17.9 预测337 340 342 343 343 342 340 338 335[]T d c b a =β()()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=672....21f f f Y ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=1672672672 (122211112)32323X []Td c b a =ββX Y =β()()ββX Y X Y T--[]6.29857.0002.0109.16-⨯=-β真实56.3 377.4 435.9 154.6 380.4 68.2 49.6 26.2 -1.3 预测332 329 326 323 320 318 317 317 319 真实68.6 325.2 219.6 690 836.4 131.5 预测319 322 323 332 340 350总偏差 =128562 回归平方和 =165248 相关系数R=1.35月31号到6月6号真实值（实线）与预测值（虚线）的比较对P B 进行拟合：P B的回归系数的估计值 5月31号部分真实值与预测值真实224.5 305.1 65.1 65.3 99.7 210.4 406.4 25.4 150 预测294 301 308 314 320 325 330 334 338 6月1号到6月6号部分真实值与预测值真实11.2 155.3 55.1 398.9 -0.38 148 112.6 39.1 47.7 预测344 349 353 352 358 359 359 360 360 真实46.7 387.2 488.7 90.1 47.34 114.2 102.2 37.7 -0.75[]4.29463.0002.0101.26-⨯=-β291∑⎪⎭⎫ ⎝⎛-=-y y i T ss 291^∑⎪⎭⎫ ⎝⎛-=-y y i Rss预测359 359 360 361 362 364 368 372 378 真实78.2 308.1 168.7 848.9 853.8 139.3 预测372 385 394 405 433 451R=1.355月31号到6月6号真实值（实线）与预测值（虚线）的比较对P C 进行拟合：PC的回归系数的估计值 5月31号部分真实值与预测值真实141.2 59.5 22.6 188.3 210.7 631.6 36 208.8 预测296 303 308 314 319 323 327 330 6月1号到6月6号部分真实值与预测值真实28.5 82.3 49.7 352.0 -0.56 144.1 139.9 40.2 79.3 预测333 337 341 342 343 342 340 338 335 真实126.6 403.5 502.8 358.3 152.2 94.9 97.8 15 -0.75 预测331 327 323 318 314 310 306 303 301[]6.29656.0002.0108.16-⨯=-β291∑⎪⎭⎫⎝⎛-=-y y i T ss 291^∑⎪⎭⎫⎝⎛-=-y y i Rss真实79.4 280.8 230.5 852 849.5 144.0 预测299 299 301 304 309 315=266889 =148807 R=0.56 5月31号到6月6号真实值（实线）与预测值（虚线）的比较对P D 进行拟合：P D的回归系数的估计值 5月31号部分真实值与预测值真实160.5 105.4 49.0 31.7 197.8 338.8 521.8 34.9 87.5 预测255 265 275 284 292 300 306 313 318 6月1号到6月6号部分真实值与预测值真实-2.4 202.7 112.2 243.9 -0.75 107.2 234.9 98.0 40.2 预测327 334 339 342 345 345 345 345 344 真实43.2 252.9 382.2 318 115.2 72.5 51.2 32.5 -1.68 预测343 342 341 341 342 344 348 353 360 真实22.6 242.8 200.1 688.9 844.8 97.5[]1.25592.0003.0101.36-⨯=-β291∑⎪⎭⎫ ⎝⎛-=-y y i T ss 291^∑⎪⎭⎫ ⎝⎛-=-y y i Rss预测369 381 396 431 434 458 = 216118 =133551 R=0.62 5月31号到6月6号真实值（实线）与预测值（虚线）的比较对 进行拟合：的回归系数的估计 5月31号部分真实值与预测值真实1254 565 536 786 366 266 1423 598 935 预测1040 1171 1199 1225 1249 1270 1289 1306 1321 6月1号到6月6号部分真实值与预测值真实486 742 545 869 523 985 1237 879 985 预测1345 1363 1375 1382 1384 1382 1378 1370 1361 真实445 652 658 546 672 698 745 836 925 预测1351 1341 1332 1324 1317 1314 1314 1319 1328 真实856 622 769 678 985 1423 预测1343 1364 1393 1330 1475 1530p 4p 4[]114068.2009.0109.86-⨯=-β291∑⎪⎭⎫⎝⎛-=-y y iTss 291^∑⎪⎭⎫⎝⎛-=-y y i Rss=198775 =156899 R=0.79 5月31号到6月6号真实值（实线）与预测值（虚线）的比较对P 58进行拟合:P 58的回归系数估计值 5月31号部分真实值与预测值真实11634.8 9585.7 2462.6 2779 9054 11771 22250 1136.6 预测14300 14699 15065 15400 15705 15980 16227 16447 6月1号到6月6号部分真实值与预测值真实372.3 8154.4 4295.1 16014.2 -36 4749 6218.2 2171.1 预测16946 17178 17331 17414 17438 17411 17344 17245 真实1836 3239 16986 21741.8 10245.94 5045.6 4452 3661 预测17126 16994 16861 16735 16626 16544 16498 16499 真实1626 76.3 3225 15812 11165 33088 36925 7549 预测16555 16676 16872 17153 17528 18006 18598 19313=300691005 =357103273 R=1.2 5月31号到6月6号真实值（实线）与预测值（虚线）的比较.[]441043.169.34118.01016.1⨯-⨯=-β291∑⎪⎭⎫ ⎝⎛-=-y y i T ss 291^∑⎪⎭⎫ ⎝⎛-=-y y i R ss 291∑⎪⎭⎫ ⎝⎛-=-y y i T ss 291^∑⎪⎭⎫ ⎝⎛-=-y y i Rss由图形可知，模型能够较好的反映一周内的功率变化趋势，在实际值比较平稳的阶段，能够较好的预测功率.但是由于影响功率变化因素比较多，所以预测值与实际值存在很大的误差.要准确预测实际功率，要充分考虑各种因素的影响，记录超大量的数据，对数据进行分析，计算量就大大增加模型二：灰色理论模型等间隔的选取7天内的28个数据，组成原始数据序列 .对原始数据序列 对数据进行累加得到新的序列 设 满足一阶常微分方程 （1） 其中a ，u 是常数，并且满足的解为 （2）对等间隔取样的离散值（ ）,则为（3）通过最小二乘法来估计常数a ，u 的值()1x (1)(1)dx ax u dt+=(1)(1)00()t t x x t ==当时()()()()()a u e a u t x t x t t a +⎥⎦⎤⎢⎣⎡-=--001110=t ()()()()a u e a u x k x ak +⎥⎦⎤⎢⎣⎡-=+-1111()()()()}283,2,1{101⋅⋅⋅⋅==∑i j x i xi})28(,),2(),1({)0()0()0()0(x x x x=()0x因 作初值用，故将 代入方程（2），用差分代替微分，又因等间隔取样， 故得到类似有 于是由（2）得到把 移到右边，并写成向量的数量级形式得到下面的方程式（4）由于涉及到累加列 的两个时刻的值，因此取前后时刻的平均代替更为合理，将 替换为 将（4）写成矩阵表达式（5）则（5）式的矩阵形式为 （6） 方程组（6）的最小二乘估计()()11x ()()()()()()28,3,2111x x x ⋅⋅⋅⋅,1)1(=-+=∆t t t (1)(1)(1)(1)(0)(2)(2)(2)(1)(2),x x x x x t∆=∆=-=∆(1)(1)(0)(0)(3)()(3),...,().xx N x x N t t∆∆==∆∆(0)(1)(0)(1)(0)(1)(2)(2),(3)(3),..............................()().x ax u x ax u x N ax N u ìï+=ïïïï+=ïïíïïïïï+=ïïî)()1(i ax (0)(1)(0)(1)(0)(1)(2)[(2),1](3)[(3),1]()[(),1]a xx u ax x u a xN x N u ⎧⎡⎤=-⎪⎢⎥⎣⎦⎪⎪⎡⎤⎪=-⎪⎢⎥⎨⎣⎦⎪⎪⎪⎡⎤=-⎪⎢⎥⎪⎣⎦⎩t x ∆∆)1()1(x )()1(i x )()(i x i ()()1[()(1)],(2,3,...,).2i i x i x i i N +-=(1)(1)(0)12(1)(1)(0)12(1)(1)(0)12[(2)(1)]1(2)[(3)(2)]1(3).1[()(1)]1()x x x a x x x u x N x N x N ⎡⎤⎡⎤-+⎢⎥⎢⎥-+⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-+-⎢⎥⎣⎦⎣⎦ BU y =yB B B ua U T T 1)(ˆˆˆ-=⎥⎦⎤⎢⎣⎡=用MATLAB 求解得到 将 代入方程则预测函数 ，将t=24,25....48代入预测函数，预测下一周的功率.5月31号到6月6号部分真实值（实线）与预测值（虚线）的曲线由图看出实际功率随天数呈现周期性变化.对P B 进行求解：用MATLAB 求解得到 预测下一周的功率.5月31号到6月6号部分真实值（实线）与预测值（虚线）的曲线⎥⎦⎤⎢⎣⎡⨯=⎥⎥⎦⎤⎢⎢⎣⎡=-5932.351106853.26^^^u a U ^^,u a ()()()()a u e a u x k x ak +⎥⎦⎤⎢⎣⎡-=+-1111()()()()()()k x k x k x 11011-+=+⎥⎦⎤⎢⎣⎡⨯=⎥⎥⎦⎤⎢⎢⎣⎡=-81.148109549.45^^^u aU用MATLAB 求解得到 预测下一周的功率.5月31号到6月6号部分真实值（实线）与预测值（虚线）的曲线对P D 进行求解：用MATLAB 求解得到 预测下一周的功率.得到5月31号到6月6号部分真实值与预测值的曲线预测下一周的功率.得到5月31号到6月6号部分真实值与预测值的曲线所有的点基本上是一个值，处于中间范围内，但是我们还是找得到这几天的功率变化趋势，即呈现周期性变化.⎥⎦⎤⎢⎣⎡⨯=⎥⎥⎦⎤⎢⎢⎣⎡=-76.237102457.24^^^u aU ⎥⎦⎤⎢⎣⎡⨯=⎥⎥⎦⎤⎢⎢⎣⎡=-36.24110997.95^^^u a U用MATLAB 求解得到 预测下一周的功率.5月31号到6月6号部分真实值（实线）与预测值（虚线）的曲线对P 58进行求解：用MATLAB 求解得到 预测下一周的功率.5月31号到6月6号部分真实值（实线）与预测值（虚线）的曲线而且由预测值与真实值比较，预测的曲线基本呈现一条直线，根本原因是点较多，离散型太强，导致无法得到可靠的曲线.⎥⎦⎤⎢⎣⎡⨯=⎥⎥⎦⎤⎢⎢⎣⎡=-7224.979106159.14^^^u aU ⎥⎦⎤⎢⎣⎡⨯⨯=⎥⎥⎦⎤⎢⎢⎣⎡=-44^^^103908.1101208.2u aU模型三：建立自回归移动平均模型ARIMA(p,d,q)通过作图容发现，时间序列是非平稳的，而当通过差分后就会显示出平稳序列的性质.为了方便,没有考虑随机项，令q=0.对序列进行一阶差分 得到新的序列t tX Y ∇=利用公式 计算序列的自相关 系数 ,当 距离0越远时，说明次序列具有很强的相关性.对 进行零均值化得到新的序列t Z ，然后对序列t Z 样本的协方差函数，利用 得到样本的协方差估计，利用计算样本的自相关系数.我们建立的是三阶线性函数.通过公式求解模型参数.用MATLAB 进行编程计算.随机选择5月30号的等间隔12个数据，建立三阶线性函数，用此函数预测下一周内的功率.其中 p=6 对A P 参数为依次-0.08227 -0.4727 0.12 0.31 0.12 0.026代入线性方程 得到预测函数. 5月31号部分真实值（实线）与预测值（虚线）的曲线K γK γ1t t t X X X -∇=-t Y 0ˆˆˆk k γργ=∑-+=kk t t k Z Z 121121γ⎪⎪⎪⎪⎪⎭⎫⎝⎛⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⋅⋅⋅111115111151621ρρρρρρρρϕϕφ∑∑⎪⎭⎫⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-=---+-1212301Y Y Y Y Y Y t kk t t k γ()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ332211B P 参数依次为-0.1588 0.1328 -0.186 0.2759 0.18 -0.4368代入线性方程 得到预测函数. 5月31号部分真实值（实线）与预测值（虚线）的曲线C P 的参数依次是-0.3553 -0.1772 -0.126 -0.24 -0.03 -0.103代入线性方程 得到预测函数.5月31号部分真实值（实线）与预测值（虚线）的曲线D P 参数依次是-0.50 -0.1088 0.0012 0.2374 0.1735 0.1643代入线性方程 ，得到预测函数. ()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ332211()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ332211()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ3322115月31号部分真实值（实线）与预测值（虚线）的曲线4P 的参数依次是-0.167 0.0737 -0.1078 0.3998 0.0326 -0.13代入线性方程 得到预测函数. 5月31号部分真实值（实线）与预测值（虚线）的曲线的参数依次是-0.1581 0.068 -0.0692 0.3666 0.0358 -0.2807代入线性方程 得到预测函数. 5月31号部分真实值（实线）与预测值（虚线）的曲线()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ33221158P ()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ332211利用ARIMA(p,d,q)模型对5月31号的风电功率进行预测，模型的预测结果曲线和实际功率曲线很接近.但在变化规律发生改变时，预测曲线通常不能提前变化，具有延时性.同时为了提高预测的精度，所取得12组数据可以是一个周期内对应数据的平均值，依次来提高预测精度.对问题（一）：由模型一预测A P 准确率为 0.9707 合格率为1 准确率为0.8525 合格率0.93 由模型二预测A P 准确率为 0.9750 合格率为1 准确率为0.8621 合格率为0.90 由模型三预测A P 准确率为0.9936 合格率为1 准确率为0.8941 合格率为0.93并且我们能够通过比较预测值与真实值的拟合好坏，得到模型三更加精确.模型一与模型二通常只是能够反映数据的整体变化趋势，并不能很精确的预测数据.模型三能够更精确的预测数据，同时也会表现出滞后性.总体来说ARIMA(p,d,q)更加优越.对问题（二）：把预测汇聚后的准确率与预测单机的准确率进行比较，汇聚后的准确率较小.由于对每台单机的功率进行预测时，都存在一定的误差.所以当单机汇聚在一起，通过单机功率估计总功率时，会造成误差的积累，所以误差会增大.同时由于风电功率不稳定，所以电力调度部门会对风电总功率进行一定调整，也会减小汇聚后准确率.4P 4P 4P对问题（三）：风电功率受到多种因素的影响如风向、风速、空气密度等.当其中一个因素发生变化时,风电功率都会发生较大的变化.实际情况是风向、风速、空气密度具有随机性，所以风电功率也具有随机性.在建立模型三时，选取的是5月30号的12组数据，选取的数据具有随机性.为了使数据具有普遍性，选取从5月24号到5月30号每天相同时刻数据.相同时刻的数据取平均值，然后按照原来的顺序组成一组新数据.利用这组数据建立ARIMA(p,d,q)模型，预测单机A 风电功率. 用MATLAB 对ARIMA(p,d,q)模型参数求解得到参数写成矩阵为代入线性方程 得到预测函数.5月31号部分真实值（实线）与预测值（虚线）的曲线从上图中可以看出,取平均值后，得到的ARIMA(p,d,q)预测的A 的单机风电功率非常精确，除去极个别值预测值与真实值有较小的误差外，其余值基本相同.所以取7天内对应数据平均值，可以更好地反映数据的特征，建立的模型能更加精确预测风电功率.[]0621.01965.00132.02386.04761.04549.0-----()()()()P T p t t t t X X X X X ----+⋅⋅⋅+++=ϕϕϕϕ332211（六）模型评价模型优点：（1）回归曲线模型和灰色理论模型可以描述出风电功率的整体变化趋势，由回归曲线模型和灰色理论模型得到风电功率变化具有很大随机性.（2）ARIMA(p,d,q)可以在回归曲线模型和灰色理论模型基础上，更加精确地预测风功率.模型缺点：（1）这三种模型不适用于对长远数据的预测（2）模型的建立没有考虑到外部因素的影响（七）模型推广这三种模型都可以应用到很多的预测领域，其中ARIMA(p,d,q)可以对没有较强规律的数据进行较为准确预测.（八）参考文献[1] 韩中庚，数学建模方法及其应用[M]，北京：高等教育出版社，2005[2] 韩中庚，数学建模竞赛——获奖论文精选与点评，北京：科学出版社，2007（九）附录function [a,c]=xulie(X) filename=[X '.dat'];b=load(filename);n=96;m=28;a=zeros(m,n);c=zeros(m,5);for i=1:mfor k=1:na(i,k)=sum(b(i,1:k));endt=1:96;figure('visible','off');p=polyfit(t,a(i,:),4);c(i,:)=p;f=polyval(p,t);plot(t,a(i,:),'k',t,f,'--k');fnn=[X num2str(i) '.jpg'];saveas(gcf,fnn)close(gcf)endfunction [a,c]=liejie(X,m,n) filename=[X '.dat'];b=load(filename);a=[];for k=m:na=cat(2,a,b(m,:)); endt=1:672;figure('visible','off');p=polyfit(t,a,3);c=p;f=polyval(p,t);plot(t,a,'k',t,f,'--k');fnn=[X 'lianj' '.jpg'];saveas(gcf,fnn)close(gcf)function xkp=jian(lj,X,Y) [d,~]=lj(X,22,28);e=zeros(1,28);for t=1:28e(t)=d(t*24);endfn1=[X '.dat'];fn2=[Y '.mat'];b=load(fn1);s=load(fn2,Y);c=s.(Y);x1=b(15,1);u=c(2);a=c(1);xkp=zeros(1,28);for k=29:56xkp(k-28)=(x1-u/a)*exp((-a)*k)-(x1-u/a)*exp((-a)*(k-1)); endfigure('visible','off');figure(1);t=1:28;plot(t,e,'k',t,xkp,'--k');grid on;fnn=[X 'sz' '.jpg'];saveas(gcf,fnn);close(gcf);。