一种面向水电站水库的随机优化调度方法与流程

本发明涉及水库调度技术领域，特别是涉及一种面向水电站水库的随机优化调度方法。

背景技术

水库优化调度研究兴起于20世纪60年代，迄今为止，已经在优化调度理论、以及调度模型构建和求解上取得了一系列丰硕的研究成果。然而，这些成果很少应用到水电站的实际调度中。究其原因，主要是由于水库调度面临着径流不确定性的难题。优化调度模型常采用这样两种方式考虑径流不确定性：①隐随机方式，即假定未来径流轨迹是历史径流的重演，通过历史径流场景反映未来径流的不确定性特征。②显随机方式，即基于非确定性的水文时间序列及其概率分布描述随机入流。其中，随机动态规划(sdp)是应用最为广泛的显随机优化方法。

sdp通过径流状态转移概率来最终量化径流不确定性。伴随着贝叶斯理论的发展和应用，其状态转移概率可根据径流预报信息不断地滚动更新，出现了大量的考虑径流预报不确定性的sdp模型—贝叶斯sdp(bayesiansdp，bsdp)模型。然而，目前状态转移概率主要通过以频率代替概率的经验估计方法获得，其计算精度完全依赖于实测/人工生成径流样本。当实测样本有限或者人工生成径流样本不能保留原始样本的特征时，计算结果将失真，会降低随机优化调度的效益。而且，随着径流不确定性被考虑到sdp中，bsdp不仅要考虑径流自身的转移规律，还要计算预报径流和实测径流的转移概率。此时，当实测径流样本不够时，采用随机模拟模型丰富径流样本的方法将变得不可行，因为同时产生人工径流系列和其预报径流系列，而不改变原始径流特性和其与径流预报系列的复杂相关关系，是十分困难的。

技术实现要素：

发明目的：为解决现有技术的不足，本发明提供一种面向水电站水库的随机优化调度方法。

技术方案：为实现上述发明目的，采用以下技术方案：

一种面向水电站水库的随机优化调度方法，包括以下步骤：

(1)建立径流及其预报不确定性的量化模型；

(2)提出径流及其预报不确定性的理论估计方法；

(3)构建不同径流描述方式下的随机动态规划模型，提取调度规则，以指导水电站水库实际运行。

进一步的，所述步骤(1)包括以下步骤：

(11)为了描述径流自身随机性，建立径流先验状态转移概率

将径流视为一个简单的一阶马尔科夫过程，并以径流状态转移概率p^t-1ij来描述径流自身的随机性，即p^t-1ij＝p[qt∈j|qt-1∈i]，它表示已知时段t-1的径流qt-1为等级i时，时段t的径流qt为等级j的概率；

由于径流状态转移概率p^t-1ij是在未获得任何径流预报信息的基础上得到的，将其称之为径流先验状态转移概率；

(12)根据径流预报不确定性，建立径流状态转移的似然概率，并进一步建立径流的可预测性概率

用似然概率p^tjk表示径流预报的不确定性，即p^tjk＝p[q^ft∈k|qt∈j]，它表示时段t实测径流值qt处于等级j时，预报径流值q^ft处于等级k的概率；

当获得了径流先验状态转移概率和似然概率之后，利用全概率公式得到径流的可预测性概率p^tjl，即p^tjl＝p[q^ft+1∈l|qt∈j]，它表示时段t实际径流处于等级j时，t+1时段预测径流处于等级l的概率；径流的可预测性概率计算公式为：

(13)获取径流状态转移的后验概率

根据获取到的新的预报信息，利用贝叶斯定理，推导径流状态转移的后验概率；这样，随着径流预报信息不断滚动更新，不断地对径流先验状态转移概率进行修正，径流状态转移的后验概率的计算公式为：

进一步的，所述步骤(2)包括以下步骤：

(21)以实测径流样本和预报径流样本为数据基础，分别建立各时段实测径流和预报径流的边缘分布；

(22)以copula函数为工具，建立相邻时段实测径流的联合分布和同一时段预报径流和实测径流的联合分布；

假设随机变量x,y的边缘分布分别为：

fx(x)＝p(x≤x)(3)；

fy(y)＝p(y≤y)(4)；

x,y的联合分布为：

fx,y(x,y)＝p(x≤x,y≤y)(5)；

根据sklar定理：令f为一个二维分布函数，其边缘分布为u＝fx(x)，v＝fy(y)；则存在一个二维copula函数c，使得对任意x,有：

fx,y(x,y)＝c(u,v)＝c(fx(x),fy(y))(6)；

采用archimedeancopula家族的两种常见的copula函数，即gumbel-hougaardcopula或者frankcopula，拟合随机变量x,y的联合分布；

gumbel-hougaardcopula函数公式为：

frankcopula函数公式为：

其中，θ为copula函数的结构相关参数，可根据其与kendall秩相关系数的关系直接计算；

当随机变量x，y为相邻时段的实测径流，或者为同一时段的实测径流和预报径流时，利用上述copula函数，分别建立相邻时段实测径流的联合分布，以及同一时段预报径流和实测径流的联合分布；

(23)采用条件概率公式，计算描述径流自身随机性的先验状态转移概率和描述径流预报不确定性的似然概率

已知随机变量x,y的联合分布fx,y(x,y)，给定x1<x≤x2，y≤y2的条件概率公式为：

式中，c(*)为一种给定的copula函数；u＝fx(x)、v＝fy(y)为x、y的边缘分布；ui＝fx(xi)，vi＝fy(yi)，i＝1，2；

此时，给定x1<x≤x2，y1<y≤y2的条件概率计算公式为：

p(y∈(y1,y2]|x∈(x1,x2])＝p(y≤y2|x1<x≤x2)-p(y≤y1|x1<x≤x2)(10)；

若给定随机变量x为时段t的实测径流qt，且qt∈(x1,x2]，随机变量y为时段t+1的实测径流qt+1，且qt+1∈(y1,y2]，则t+1时段径流的先验状态转移概率可通过式(10)计算；

若给定随机变量x为时段t的实测径流qt，且qt∈(x1,x2]，随机变量y为时段t的预报径流q^ft，且q^ft∈(y1,y2]，则t时段径流的似然概率也可通过公式(10)计算。

(24)基于公式(2)的贝叶斯定理和公式(1)的全概率公式，计算径流的后验状态转移概率和可预测性概率。

进一步的，所述步骤(3)中水电站水库随机优化调度以考虑径流不确定性下的期望发电量最大为目标，其目标函数表示为：

其中，f为调度期t内的最大期望发电量；bt(st-1,qt,st)为时段t水电站的发电出力；st-1和st分别为时段t水库初、末库容；qt为时段t径流；eqt为期望算子；△t为调度时段；

模型需满足水电站水库调度的水量平衡约束、库容约束、出库流量约束、发电出力约束、水库特征曲线约束。

进一步的，所述步骤(3)中，当径流预报精度很低时，考虑径流预报信息反倒会误导调度决策，采用无预报sdp模型提取调度规则；当径流预报精度较高时，采用耦合径流自身随机性和预报不确定性的贝叶斯sdp模型提取调度规则；特殊情况下，当前时段的径流预报信息完全准确时，采用完美预报sdp模型提取调度规则；三种模型分别为：

(a)无预报sdp模型

首先利用径流不确定性的理论估计方法计算先验状态转移概率，然后利用公式(12)所示的递推方程逆向递推计算，获得能指导水电站水库实际运行的调度规则；

无预报sdp模型，不考虑当前时段的径流预报信息，只考虑径流自身的随机转移规律；t时段径流由上一时段的径流状态qt-1确定，水库调度决策由实测径流qt-1和初始库容st-1共同决定；无预报sdp模型的递推方程为：

其中，j为实测径流等级指标；p(qt∈j|qt-1)为t时段的径流先验状态转移概率；bt(st-1,qt∈j,st)为t时段初、末库容分别为st-1和st，径流qt∈j时的即时效益；ft(st-1,qt-1)为给定初始状态为st-1,qt-1情况下，时段t至t的最大期望效益；

(b)bsdp模型

首先利用径流不确定性的理论估计方法计算先验状态转移概率、似然概率，并通过贝叶斯定理和全概率公式得到径流的后验状态转移概率和可预测性概率，然后利用公式(13)所示的递推方程逆向递推计算，获得能指导水电站水库实际运行的调度规则；

bsdp模型除了考虑径流自身的随机转移规律之外，还利用贝叶斯定理将径流预报不确定性以似然概率的形式考虑到递推方程中；t时段径流由前一时段径流状态qt-1和本时段的径流预报q^ft共同确定，水库调度决策由实测径流qt-1、预报径流q^ft和初始库容st-1决定；此时bsdp模型的递推方程为：

其中，j，k分别为实测径流和预报径流的等级指标；q^ft和q^ft+1分别为t时段和t+1时段的径流预报值；p(qt∈j|qt-1,q^ft)为t时段径流的后验状态转移概率；p(q^ft+1∈k|qt∈j)为t+1时段的径流的可预测性概率；

(c)完美预报sdp模型

首先利用径流不确定性的理论估计方法计算先验状态转移概率，然后利用公式(14)所示的递推方程逆向递推计算，获得能指导水电站水库实际运行的调度规则；

完美预报sdp模型假定当前时段有准确的径流预报信息，t时段的径流预报值等于实测值，即q^ft＝qt；t时段的水库调度决策由实测径流qt和初始库容st-1共同决定；完美预报sdp模型的递推方程为：

其中，j为实测径流等级指标；p(qt+1∈j|qt)为t+1时段的径流先验状态转移概率；bt(st-1,qt,st)为t时段初、末库容分别为st-1和st，径流为qt时的即时效益；ft(st-1,qt)为给定初始状态为st-1,qt情况下，时段t至t的最大期望效益。

有益效果：与现有技术相比，本发明具有以下优点：

(1)突破了传统经验估计方法在量化径流不确定性时，完全受制于实测径流样本的缺陷。传统经验估计方法的状态转移概率计算完全依赖于历史径流样本，采用以频率代替概率的经验估计方法获得。令p(qt+1∈j|qt∈i)、p(q^ft∈j|qt∈i)分别表示t时段实测径流为等级i时，t+1时段先验状态转移概率和t时段的似然概率。由于历史径流样本有限，传统方法往往存在两个缺点：1)当t时段无历史实测径流数据处于等级i时先验状态转移概率和似然概率将无法计算；2)当t+1时段无历史实测径流处于等级j，或者t时段无历史预报径流处于等级j时，先验状态转移概率和似然概率矩阵中会出现大量的0值。由于水文现象的随机性和模糊性，依赖历史径流样本计算的状态转移概率的0值，有时并不合理。理论估计方法通过建立实测径流和预报径流的边缘分布，实现了径流的插补延长，克服了传统方法的缺点。

(2)计算更加简便。理论估计方法只需要计算相邻时段实测径流，以及同一时段实测和预报径流的联合分布，即可得到先验状态转移概率和似然概率的显示表达式。然后，基于贝叶斯定理和全概率公式，可计算径流的后验状态转移概率和可预测性概率。本发明推荐采用archimedeancopula函数拟合径流的联合分布，其只含有一个结构相关参数θ，且θ可根据其与kendall秩相关系数的关系直接计算。而经验估计方法虽然计算思路简单，但是计算量较大。特别是当历史径流样本量不够时，还需要先生成人工径流样本，增加了计算量。

(3)避免了生成人工径流样本。理论估计方法通过建立径流的边缘分布可实现对径流的插补和延长，避免了借助其它径流模拟模型生成人工径流样本，降低了引入径流随机模拟模型带来的不确定性。而且，随着预报不确定性被考虑在sdp中，不仅要计算径流自身前后时段的状态转移概率，还要计算预报径流和实测径流的状态转移概率。此时，要同时产生人工径流系列和其预报径流系列，而不改变原始径流特性和其与预报径流系列的相关关系特征，基本上是不可能的。

(4)数据存储空间大大减小。经验估计方法需要存储每个时段的状态转移概率，所需数据存储空间随着径流等级和状态变量的增加呈指数型增长。而理论估计方法只需要存储各时段的径流边缘分布参数和联合分布参数，即可利用公式实时计算状态转移概率。以边缘分布采用皮尔逊ⅲ型分布、联合分布采用archimedeancopula函数为例，理论估计方法每一个时段只需存储边缘分布的三个参数(均值、cv、cs)和copula函数的一个结构相关参数(θ)，大大的节省了数据存储空间。

(5)提高随机优化调度的效益。理论估计方法克服了经验估计方法量化径流不确定性的缺陷，降低了由于引入径流随机模拟模型带来的不确定性，可提高径流不确定性的量化精度，进而提高随机优化调度效益。

附图说明

图1为本发明方法流程图；

图2为径流不确定性分析流程图；

图3为nf-sdp模型递推过程图；

图4为bsdp模型递推过程图；

图5为pf-sdp模型递推过程图。

具体实施方式

下面结合附图和具体实施例对本发明的技术方案进行详细说明。

本发明的一种面向水电站水库的随机优化调度方法，如图1所示，首先建立了径流及其预报不确定性量化模型，然后以copula函数为理论基础，提出了径流及其预报不确定性的理论估计方法，最后考虑实际径流预报水平，建立了不同径流描述方式下的sdp模型，以提取调度规则，指导水电站水库的实际运行。具体为：

(1)建立径流及其预报不确定性描述模型

(11)为了描述径流自身随机性，建立径流先验状态转移概率

本发明将径流视为一个简单的一阶马尔科夫过程，并以径流状态转移概率p^t-1ij来描述径流自身的随机性，即p^t-1ij＝p[qt∈j|qt-1∈i]，它表示已知时段t-1的径流qt-1为等级i时，时段t的径流qt为等级j的概率。

由于径流状态转移概率p^t-1ij是在未获得任何径流预报信息的基础上得到的，将其称之为径流先验状态转移概率。

(12)根据径流预报不确定性，建立径流状态转移的似然概率，并进一步建立径流的可预测性概率；

本发明用似然概率p^tjk表示径流预报的不确定性，即p^tjk＝p[q^ft∈k|qt∈j]，它表示时段t实测径流值qt处于等级j时，预报径流值q^ft处于等级k的概率。似然概率的大小与预报水平息息相关，当预报精度较高时，预报径流值和实测径流值取相同等级的概率越高。因此，预报径流和实测径流的似然概率矩阵可以反映预报误差信息。

当获得了径流先验状态转移概率和似然概率之后，可利用全概率公式得到径流的可预测性概率p^tjl，即p^tjl＝p[q^ft+1∈l|qt∈j]，它表示时段t实际径流处于等级j时，t+1时段预测径流处于等级l的概率。径流的可预测性概率计算公式为：

(13)获取径流状态转移的后验概率

由于径流状态转移概率p^t-1ij是在未获得任何径流预报信息的基础上得到的，将其称之为径流先验状态转移概率。一旦获取到新的预报信息，可利用贝叶斯定理，推导径流状态转移的后验概率。这样，随着径流预报信息不断滚动更新，可不断地对先验状态转移概率进行修正，以提高径流描述精度。径流状态转移的后验概率的计算公式为：

(2)提出径流及其预报不确定性的理论估计方法

本发明首先建立各时段实测径流和预报径流的边缘分布；然后基于copula函数，建立相邻时段实测径流，以及同一时段实测和预报径流的联合分布；然后，根据公式推导直接得到状态转移概率的计算公式。以此为基础，提出了一套完整的径流及其预报不确定性的理论估计方法流程。如图2所示，主要包括以下子步骤：

(21)以实测径流样本和预报径流样本为数据基础，分别建立各时段实测径流和预报径流的边缘分布；

根据《水利水电工程设计洪水计算规范》(sl44-2006)规定，水文频率曲线的线型一般采用皮尔逊ⅲ型(p-iii型)曲线，特殊情况，经分析论证后可采用其它线型。因此，本发明推荐采用p-iii型分布建立径流的边缘分布，实际中经论证也可采用其它分布。

p-iii型分布的概率密度函数为：

式中：γ(α)为α的伽玛函数；α、β、a0分别为p-iii型分布的形状、尺度和位置参数，α>0，β>0。这三个参数与总体样本的三个统计参数均值变差系数cv、偏态系数cs可以相互推求，具有下列关系：

(22)以copula函数为工具，建立相邻时段实测径流的联合分布和同一时段预报径流和实测径流的联合分布

假设随机变量x,y的边缘分布分别为：

fx(x)＝p(x≤x)(5)；

fy(y)＝p(y≤y)(6)；

x,y的联合分布为：

fx,y(x,y)＝p(x≤x,y≤y)(7)；

根据sklar定理：令f为一个二维分布函数，其边缘分布为u＝fx(x)，v＝fy(y)；则存在一个二维copula函数c，使得对任意x,有：

fx,y(x,y)＝c(u,v)＝c(fx(x),fy(y))(8)；

archimedeancopula是一类重要的copula函数，由于其构造方便、容易使用等优点，在水文领域应用广泛。本发明推荐采用archimedeancopula家族的两种常见的copula函数，即gumbel-hougaardcopula或者frankcopula，拟合随机变量x,y的联合分布；

gumbel-hougaardcopula函数公式为：

frankcopula函数公式为：

其中，θ为copula函数的结构相关参数，可根据其与kendall秩相关系数的关系直接计算；

当随机变量x，y为相邻时段的实测径流，或者为同一时段的实测径流和预报径流时，可利用上述推导过程，分别建立相邻时段实测径流的联合分布，以及同一时段预报径流和实测径流的联合分布。

(23)采用条件概率公式，计算描述径流自身随机性的先验状态转移概率和描述径流预报不确定性的似然概率

已知随机变量x,y的联合分布fx,y(x,y)，给定x1<x≤x2，y≤y2的条件概率公式为：

式中，c(*)为一种给定的copula函数；u＝fx(x)、v＝fy(y)为x、y的边缘分布；ui＝fx(xi)，vi＝fy(yi)，i＝1，2；

此时，给定x1<x≤x2，y1<y≤y2的条件概率计算公式为：

p(y∈(y1,y2]|x∈(x1,x2])＝p(y≤y2|x1<x≤x2)-p(y≤y1|x1<x≤x2)(12)；

若给定随机变量x为时段t的实测径流qt，且qt∈(x1,x2]，随机变量y为时段t+1的实测径流qt+1，且qt+1∈(y1,y2]，则t+1时段径流的先验状态转移概率可通过式(12)计算；

若给定随机变量x为时段t的实测径流qt，且qt∈(x1,x2]，随机变量y为时段t的预报径流q^ft，且q^ft∈(y1,y2]，则t时段径流的似然概率也可通过公式(12)计算。

(24)基于公式(2)的贝叶斯定理和公式(1)的全概率公式，计算径流的后验状态转移概率和可预测性概率。

(3)构建不同径流描述方式下的sdp模型，提取调度规则，指导水电站水库实际运行

水电站水库随机优化调度以考虑径流不确定性下的期望发电量最大为目标，其目标函数表示为：

其中，f为调度期t内的最大期望发电量；bt(st-1,qt,st)为时段t水电站的发电出力；st-1和st分别为时段t水库初、末库容；qt为时段t径流；eqt为期望算子；△t为调度时段；

模型应当满足以下约束条件：

水量平衡约束：

st＝st-1+(qt-rt)△t(14)；

库容约束：

st,min≤st≤st,max(15)；

流量约束：

rt,min≤rt≤rt,max，qfd,t≤qfd,max(16)；

rt＝qfd,t+qqs,t(17)；

发电出力约束：

nt＝a*qfd,t*△ht(18)；

nt≤nmax(19)；

水位～库容曲线约束：

zt＝f(st)，st＝f^-1(zt)(20)；

尾水位～流量曲线约束：

zdr,t＝g(rt),rt＝g^-1(zdr,t)(21)；

其中，rt表示时段t水库的放水流量，m³/s；rt,max,rt,min分别表示时段t水库的放水流量上、下限t,m³/s；st,max,st,min分别表示时段t水库的库容上、下限，m³；qfd,t,qqs,t分别表示时段t水电站的发电流量和弃水流量，m³/s；qfd，max表示时段t水电站的发电流量上限，m³/s；nt,nmax分别表示时段t水电站的发电出力和其装机容量，kw.h；a表示出力有效系数；△ht表示时段t发电水头，m；zt,zdr,t分别表示时段t水库的水位和尾水位，m；f(*)表示水位-库容关系约束；g(*)表示尾水位-流量关系约束。

根据是否考虑径流预报信息及其不确定性，本发明构建了3种sdp模型。当径流预报精度很低时，考虑径流预报信息反倒会误导调度决策，本发明推荐采用无预报sdp模型提取调度规则；当径流预报精度较高时，本发明推荐采用耦合径流自身随机性和预报不确定性的bsdp模型提取调度规则。特殊情况下，当前时段的径流预报信息完全准确时，推荐采用完美预报sdp模型提取调度规则。根据实际径流预报水平，可分别建立如下三种实施例：

1)实施例1-无预报sdp模型

无预报sdp(noforecastingsdp,nf-sdp)模型，不考虑当前时段的径流预报信息，只考虑径流自身的随机转移规律。t时段径流由上一时段的径流状态qt-1确定，水库调度决策由实测径流qt-1和初始库容st-1共同决定。nf-sdp模型的递推方程为：

其中，j为实测径流等级指标；p(qt∈j|qt-1)为t时段的径流先验状态转移概率；bt(st-1,qt∈j,st)为t时段初、末库容分别为st-1、st，径流qt∈j时的即时效益；ft(st-1,qt-1)为给定初始状态为st-1,qt-1情况下，时段t至t的最大期望效益。nf-sdp模型的递推过程见附图3。

实施例1首先利用本发明提出的径流不确定性的理论估计方法计算先验状态转移概率，然后利用公式(22)所示的递推方程逆向递推计算，获得能指导水电站水库实际运行的调度规则。实际调度中，利用前一时段的实测径流等级和初始库容作为调度判别指标，实时确定本时段的末库容状态。

2)实施例2-bsdp模型

bsdp除了考虑径流自身的随机转移规律之外，还利用贝叶斯定理将径流预报不确定性以似然概率的形式考虑到递推方程中。t时段径流由前一时段径流状态qt-1和本时段的径流预报q^ft共同确定，水库调度决策由实测径流qt-1、预报径流q^ft和初始库容st-1决定。此时bsdp模型的递推方程为：

其中，j，k分别为实测径流和预报径流的等级指标；q^ft和q^ft+1分别为t时段和t+1时段的径流预报值；p(qt∈j|qt-1,q^ft)为t时段径流的后验状态转移概率；p(q^ft+1∈k|qt∈j)为t+1时段的径流的可预测性概率。bsdp模型的递推过程见附图4。

实施例2首先利用本发明提出的径流不确定性的理论估计方法计算先验状态转移概率、似然概率，并通过贝叶斯定理和全概率公式得到径流的后验状态转移概率和可预测性概率，然后利用公式(23)所示的递推方程逆向递推计算，获得能指导水电站水库实际运行的调度规则。实际调度中，利用前一时段的实测径流等级、本时段的预报径流等级，以及初始库容状态作为调度判别指标，实时确定本时段的末库容状态。

3)实施例3-完美预报sdp模型

完美预报sdp(perfectforecastingsdp，pf-sdp)模型假定当前时段有准确的径流预报信息。t时段的径流预报值等于实测值，即q^ft＝qt。t时段的水库调度决策由实测径流qt和初始库容st-1共同决定。pf-sdp模型的递推方程为：

其中，j为实测径流等级指标；p(qt+1∈j|qt)为t+1时段的径流先验状态转移概率；bt(st-1,qt,st)为t时段初、末库容分别为st-1和st，径流为qt时的即时效益；ft(st-1,qt)为给定初始状态为st-1,qt情况下，时段t至t的最大期望效益。pf-sdp模型的递推过程见附图5。

实施例3首先利用本发明提出的径流不确定性的理论估计方法计算先验状态转移概率，然后利用公式(24)所示的递推方程逆向递推计算，获得能指导水电站水库实际运行的调度规则。实际调度中，利用本时段的实测径流等级和初始库容状态作为调度判别指标，实时确定本时段的末库容状态。