王渤权,王丽萍,唐 勇,刘明浩,李传刚
(1.华北电力大学 可再生能源学院,北京 102206;2.国家电网公司国家电力调度控制中心,北京 100031)
我国具有十分丰富的水能资源,位居世界的第一位,多利用水能资源从而使水电厂多发电,以节约煤炭资源的消耗,这一举动在当今资源短缺的时代有着重要的意义,这也是电力系统对水电部门提出的要求[1]。而水电站发电量的多少主要取决于径流过程,目前水电站工作人员多采用径流预报值来做发电计划的制定,然而,由于河川径流的随机性、地区性,常常出现预报值不够精准,预报误差较大的情况,由此便产生一定的出力误差,影响着发电计划的制定,从而导致被迫修改发电计划。由于目前国内外尚未有对由于预报误差而引起的出力误差进行相关研究,为方便水电站工作人员充分掌握由于预报级别以及预报误差产生的出力误差信息,以便能够及时地对发电计划进行制定与修改,本文对预报值、预报误差以及出力误差之间的关系进行探讨,建立三者的三维联合分布模型。对于联合分布的求解,国内外学者进行了大量的工作,费永法[2]基于事件积原理,采用一种解析算法计算出不同河流洪水遭遇的概率问题,一定程度上解决了多元随机变量的联合分布问题。戴昌军、梁忠民[3]通过将正态化变换的Moran方法、传统经验频率方法以及多维转化为一维的方法应用在水文领域中进行比较分析,结果表明Moran方法具有较好的拟合度。王文圣、丁晶[4]采用核估计理论来构造一种多变量的非参数模型,并且将其应用于水文学中,为随机水文学的发展提供了新的思路。郑红星、刘昌明[5]采用了经验频率的方法对南水北调中不同的水文区域降水遭遇问题进行了计算与分析。以上专家学者们对多维随机变量的联合分布的研究均做出了一定的贡献,但这些方法对于求解预报值、预报误差以及出力误差的联合分布均存在一些不足之处:以上方法对于边缘分布函数要求多为正态分布并且有些方法对于两者的相关性为线性要求,而对于预报误差与出力误差而言,其分布函数不一定为正态分布,尤其对于出力误差分布的研究甚少,其分布未知,并且二者存在着复杂的数学关系,呈现出非线性的关系,故不适合采用上述方法进行求解计算。最大熵法是依据样本信息来推求未知分布的一种推断方法,能够客观地反映事物的本质,可以拟合任意形式的分布,文献[6]证明对于预报值以及预报误差的分布特征均可应用最大熵法进行描述,并且文献[7]指出运用最大熵方法拟合出来的最大熵分布函数较正态分布函数能够更好地反映预报误差的特点,而Copula函数[8-10]是通过边缘分布函数与相关性结构两方面构造的函数,并且适用于任何边缘分布函数,对其变量间相关性的要求也比较广泛,线性非线性均可,具有更强的灵活性与适应性。
基于此,本文将最大熵理论与Copula函数理论结合,建立三维联合分布模型来求解预报值、预报误差以及出力误差的三维联合分布,在此基础上结合贝叶斯理论进而计算出在一定预报径流值与预报误差的情况下其出力误差的条件概率的大小,为水电站工作人员调度管理提供理论依据。
针对出力误差分布未知的情况,利用最大熵理论中可以求解出任意形式的分布,以及Copula函数中适用于任何形式边缘函数的特点,将最大熵理论与Copula函数理论进行结合,从而建立三维联合分布模型对预报值、预报误差以及出力误差的联合分布进行求解。
本文建立的三维联合分布模型形式如下所示:
P(X1≤x1,X2≤x2,X3≤x3)=C(u1,u2,u3)=
φ-1[φ(u1;θ)+φ(u2;θ)+φ(u3;θ)]
(1)
对于上述模型的求解,主要是针对两个未知参数λij和θ的估计,其值的好坏直接影响着模型的准确程度,此处采用两阶段估计法对其进行求解,即先计算λij,之后由其估计值再计算出θ值。其主要计算步骤如下:
(1)对于λij的估计采用非线性规划法求解,计算公式如下:
(3)
其中:
f(xi)=exp(λi0+∑λijxji)
式中:f(xi)为最大熵概率密度分布函数;mij为第i个变量第j阶原点矩。
(2)利用得到的λij的值采用极大似然法[11]进行对θ值的估计,计算公式如下:
(4)
由此得到的λij和θ值代入式(1)即得到预报值、预报误差以及出力误差的三维联合分布函数。
模型能否很好的反映样本数据的分布特征,描述各变量之间的相关关系主要取决于对模型的检验评价,为此,本文分别采用P-P图和均方根误差法(RMSE)[12-14 ]对其进行检验,其中P-P图是用来检验数据是否服从某一指定分布的一种方法,将变量的累积比例与相应的指定分布的累积比例点绘在一张图中,并且计算回归系数,其值越接近于1,说明指定分布的拟合效果越好。而RMSE是用来衡量理论分布频率与经验频率之间的差距,其计算公式为:
(5)
式中:MSE为均方差;N为样本容量;pc为由模型得出的理论频率;p0为经验频率;i为样本序号。
由式(5)可知,其值越小,模型的拟合效果越好。
模型的整体计算流程图如图1所示。
图1 三维联合分布模型计算流程图Fig.1 The calculation flow chart of three-dimensional joint distribution model
本文选取某水电站A为例进行实例分析,分别利用该流域10年(2004-2013年)的每日12个时段的预报流量值以及实际值计算出预报出力值和实际出力值,进而得到与预报误差相应的出力误差,以此为数据样本进行计算。
由于预报值、预报误差与水电站出力误差存在着复杂的数学联系,出力误差不仅和预报值大小以及预报误差有关,还和水位因素有关。由于同一时期预报值的流量级别大致为一个范围,在整个调度期水位范围也大致相同,而不同时期水位范围却不尽相同,并且对于不同时期预报误差与出力误差也具有不同的相关关系,故采用时期的划分来反映水位的特点。因此,本文分为枯水期、过渡期和汛期3个时期进行探讨,按照当地流域资料划分的月份如表1所示。
表1 时期划分表Tab.1 Period division table
将划分好的时期对预报值、预报误差以及出力误差进行相关性分析,进而确定三维联合分布模型中生成元的类型,其中,弱相关采用AMH Copula[15]生成元:
φ(u)=ln{[1-θ(1-u)]/u}
正相关采用Gumbel-Hougaard Copula生成元:
φ(u)=-(lnu)θ
当呈现出负相关时,采用Frank Copula生成元:
φ(u)=-ln[(e-θu-1)/(e-θ-1)]
其不同时期三者之间的秩相关系数如表2所示。
表2 各变量间秩相关性系数Tab.2 The rank correlation coefficient among variables
由表2计算结果可知,在三个时期预报值与预报误差、预报值与出力误差的秩相关系数在0.1左右,均呈现出弱相关性。而预报误差与出力误差三个时期的秩相关系数在枯水期最大,为0.623;过渡期次之,为0.420;而在汛期二者的秩相关系数较弱,为0.251。通过分析上述变量间相关性可知,预报值与预报误差,预报值与出力误差呈现弱相关性,而预报误差与出力误差呈现出正相关性,故本文分别选取AMH Copula生成元以及Gumbel-Hougaard Copula生成元代入式(1)建立两种形式的模型:模型1与模型2,并且对两种模型进行检验,从而选出最优模型。
对模型1与模型2中的参数λi和θ利用式(2)~(4)进行估计,得到结果如表3~表6所示。
(1)清万树《词律》中载《卜算子》调9体[注]柳永与张先的名为《卜算子》的词作,实则是《卜算子慢》,且柳永之作为《卜算子慢》正体。[1]118-121
表3 各时期预报值拉格朗日因子Tab.3 Lagrangian multipliers of forecast value in each period
表4 各时期预报误差拉格朗日因子Tab.4 Lagrangian multipliers of forecast error in each period
表5 各时期出力误差拉格朗日因子Tab.5 Lagrangian multipliers of output error in each period
由计算出的参数λi和θ即可得到三维联合分布模型1和模型2,为了验证模型的可行性,采用χ2检验法检验模型的边缘分布的拟合程度,判断其是否通过检验。由于模型1与模型2所求边缘分布相同,故只针对一种模型进行检验,显著性水平取0.05,其检验结果如表7所示。
表6 两种模型中的θ值Tab.6 Theta values of two models
表7 χ2值检验结果Tab.7 The test results of χ2
由表7可以看出,预报值、预报误差以及出力误差分布的χ2值分别为3.21、5.51和6.17,均小于检验临界值11.07、9.48、9.48。通过检验要求,说明采用最大熵法拟合出的预报值分布、预报误差分布以及出力误差分布满足检验要求,是可行的、合理的。
针对两种形式的模型,为进一步选择出拟合度较好的一种,采用P-P图以及均方根误差法来评价两种模型优劣,从而确定最终的三维联合分布模型形式。
本文利用样本数据计算经验频率,再分别用模型1与模型2计算相应的理论频率,其中经验频率计算公式如式(15)所示:
P(x1,x2,x3)=P(X1≤x1,X2≤x2,X3≤x3)=
(6)
式中:∑Ni为样本系列中满足条件的总的数据对数;N为样本系列的数据对数。
得到的经验频率与模型1、模型2计算出相应的理论频率,共同点绘到图中,如图2~图7所示。
图2 模型1枯水期P-P图Fig.2 The P-P figure of model 1 in dry period
图3 模型2枯水期P-P图Fig.3 The P-P figure of model 2 in dry period
图4 模型1过渡期P-P图Fig.4 The P-P figure of model 1 in transitional period
图5 模型2过渡期P-P图Fig.5 The P-P figure of model 2 in transitional period
图6 模型1汛期P-P图Fig.6 The P-P figure of model 1 in flood period
图7 模型2汛期P-P图Fig.7 The P-P figure of model 2 in flood period
为了进一步评价两个模型的优劣性,分别按照式(5)计算二者的均方根误差值进行比较,所得结果如表8所示。
表8 两种模型在不同时期的均方根误差值Tab.8 The RSME values of two models in different periods
由表8可以看出,采用模型1函数计算得到的三个时期的均方根误差分别为0.024 1、0.030 7、0.022 5,而模型2的均方根误差分别为0.019 8、0.024 6、0.025 1,在枯水期以及过渡期模型2均方根误差较小,模型2更优;在汛期模型1的均方根误差较小,模型1更优。故在枯水期以及过渡期采用模型2来描述其三维联合分布,在汛期采用模型1来描述。
由此便得到在不同时期预报值、预报误差以及出力误差三者的联合分布,对于不同时期而言,其预报误差以及出力误差的二维联合分布如图8~图10。
图8 枯水期二维联合分布图Fig.8 Two-dimensional joint distribution figure in dry period
图9 过渡期二维联合分布图Fig.9 Two-dimensional joint distribution figure in transitional period
图10 汛期二维联合分布图Fig.10 Two-dimensional joint distribution figure in flood period
由图8~图10可以看出,在不同时期,预报误差与出力误差的二维联合分布呈现出不同的形状,说明预报误差与出力误差随着时期的不同而呈现出不同的关系,分时期进行研究是有必要并且合理的。在此基础上,运用所求出的三维联合分布模型结合贝叶斯理论便可以求得在预报值与预报误差一定时,水电站出力误差的条件概率大小,使得水电站工作人员方便制定与修改发电计划,更好地进行规划与管理工作。
由贝叶斯理论可知:
P(A|B)=P(AB)/P(B)
(7)
式中:P(A|B)为在事件B发生的情况下发生A的条件概率;P(AB)为事件A和B同时发生的概率;P(B)为事件B发生的概率;此处事件A表示出力误差;事件B表示预报值与预报误差。
对于连续型的概率密度函数来说,一个点的概率值为0,一个范围的概率才有意义,故首先对各个时期的预报值、预报误差进行级别划分,以此来进行概率计算,其分类结果如表9、表10所示。
表9 预报值分类结果 m3/sTab.9 Classification results of forecast value
表10 预报误差与出力误差级别划分 %Tab.10 Level classification of forecast error and output error
由上述分类结果结合贝叶斯统计理论即可得到出力误差的条件概率分布,其计算公式如式(17)所示:
(8)
式中:X、Y、Z分别表示为预报值,预报误差以及出力误差。
将三维联合分布模型代入式(8)便可计算出各个时期的出力误差条件概率大小,如:当预报值为6级,预报误差为3级时,计算得到各个时期的出力误差条件概率分布,其结果如表11所示。
表11 出力误差条件概率分布Tab.11 Conditional probability distribution of output error
由表11可以看出,当预报误差一定时,在不同时期呈现出不同的出力误差概率值。其中枯水期出力误差为3级时的概率为90%,在过渡期出力误差为4级时的概率值最大,为93.7%,而对于汛期,出力误差主要为6级,其概率为98.3%,造成这种不同主要是因为导致出力误差的原因除了预报误差外还与来流大小有着重要的关系,来流的大小主要从两方面影响着出力误差的结果,一是发电流量的多少,二是水库水位的高低;在汛期,由于此级别下的来水已足够大,导致水电站接近满发,预报误差对出力误差的影响较弱,故此时对出力误差基本没有影响,出力误差在6级处较为集中。此外,由于出力误差除了和当前时段的预报误差有关,还有当前水位与之前的预报误差有着重要的关系,故导致同一时期相同预报值与预报误差下出力误差出现不同级别的概率情况。此结果可为水电站工作人员在预报值与实际值产生一定偏差时修改发电计划提供科学的参考依据。
针对水电站运行中由于来水偏差产生出力误差而被迫修改发电计划的情况,本文对入库流量预报值、预报误差及其相应的出力误差三者之间的相关性进行了分析,并且基于最大熵理论与Copula函数理论建立三维联合分布模型,求解出三者之间的联合分布。在此基础上,利用该模型结合贝叶斯理论求解出不同时期不同预报级别与预报误差情况下水电站出力误差的条件概率分布。结果表明,在枯水期和过渡期,预报误差与出力误差相关性较强,并且成正相关,汛期各变量之间的相关性较弱;所建立的三维联合分布模型与样本数据拟合较好,能够充分反映样本数据的分布特征,各项指标均满足检验要求,用该模型表述预报值、预报误差以及出力误差的联合分布是可行的、合理的;此外,由该模型求解出的水电站出力误差条件概率不同时期呈现出不同的分布状态,符合实际情况,能够有效反映不同时期不同预报值与预报误差下其出力误差的特点。该结果可为水电站工作人员由于预报来水偏差而及时修改发电计划提供重要的指导信息,为水电站管理工作提供科学的参考依据。
由于影响水电站出力误差的因素有很多,除了预报值大小,预报误差以外,还有输电线效率,意外事故,人为因素等相关内容,本文仅考虑了一般情况,即预报值、预报误差对出力误差的影响,在下一阶段,将考虑对上述其他因素进行探讨,进一步完善对出力误差分布的研究。
□
[1] 宋萌勃,岳延兵,陈吉琴,等.水库调度与管理[M]. 郑州:黄河水利出版社,2013.
[2] 费永法. 一种计算洪水条件概率的方法[J]. 水文,1989,(1):18-23.
[3] 戴昌军,梁忠民. 多维联合分布计算方法及其在水文中的应用[J]. 水利学报,2006,(2):160-165.
[4] 王文圣,丁 晶.水文学确定性和不确定性方法的耦合初探[J].水文,2014,(2):3-7.
[5] 郑红星,刘昌明. 南水北调东中两线不同水文区降水丰枯遭遇性分析[J]. 地理学报,2000,(5):523-532.
[6] 李宪东. 基于最大熵原理的确定概率分布的方法研究[D]. 北京:华北电力大学,2008.
[7] 刁艳芳,王本德,刘 冀.基于最大熵原理方法的洪水预报误差分布研究[J].水利学报,2007,(5):591-595.
[8] 方 彬,郭生练,肖 义,等.年最大洪水两变量联合分布研究[J]. 水科学进展,2008,(4):505-511.
[9] 张 翔,冉啟香,夏 军,等. 基于Copula函数的水量水质联合分布函数[J].水利学报,2011,(4):483-489.
[10] CHOW DHARY H, ESCOBAR L A, SINGH V P. Identification of suitable copulas for bivariate frequency analysis of flood peak and volume data[J]. Hydrology Reasearch, 2011,42(2/3):193-215.
[11] 谢 华,罗 强,黄介生. 基于三维copula函数的多水文区丰枯遭遇分析[J].水科学进展,2012,(2):186-193.
[12] 康 玲,何小聪. 南水北调中线降水丰枯遭遇风险分析[J]. 水科学进展,2011,(1):44-50.
[13] Vandenberghe S,Verhoest N E C,Onof C, et al. A comparative copula-based bivariate frequency analysis of observed and simulated storm events: A case study on Barlett-Lewis modeled rainfall[J]. Water Resource Research, 2011,47:W07529.
[14] SONG S B, Vijay P Singh. Frequency analysis of droughts using the Plackett copula and parameter estimation by genetic algotirhm[J]. Stochastic Enviromental Research and Risk Assessment,2010,doi:10.1007/s00477-0101-0364-5.
[15] 李 计,李 毅,宋松柏,等.基于Copulas函数的多维干旱变量联合分布[J].自然资源学报,2013,(2):312-320.