CN105490310A - 一种基于通用分布的含风电电力系统实时动态经济调度方法 - Google Patents

Abstract

本发明公开了一种基于通用分布的含风电电力系统实时动态经济调度方法，通过统计历史数据，离线生成通用分布的形状参数表，表征在不同风电功率预测值下的风电功率概率分布；输入12个调度周期系统负荷预测数据和风电功率预测数据等；计算火电机组在实时调度12个调度周期内的火电出力可行域；使用通用分布表征风功率，使用等微增率准则进行12个时间周期的风电火电静态经济调度，得到风电和火电的静态经济调度曲线，将随机经济调度转化为确定性经济调度模型；使用等微增率准则对火电机组的出力进行修正，在保证火电机组爬坡限制并保证火电机组发电计划的稳定性的基础上，使火电出力成本最优。以某实际风电功率为基础的仿真算例验证了所提方法的可行性和合理性。本方法具有良好的推广价值和应用前景。

Description

一种基于通用分布的含风电电力系统实时动态经济调度方法

技术领域

本发明属于电力系统运行与控制领域，涉及一种基于通用分布的含风电电力系统实时动态经济调度方法。

背景技术

随着风力发电的大规模开发，风电对电网安全运行和调度控制等带来巨大挑战,一种将风电作为调度量的实时动态经济调度方法对电力安全生产至关重要。

经济调度问题是电力系统的基本问题，在风电大规模接入之前，由于火电机组良好的可控性及电力系统负荷成熟的预测方法，以往的经济调度通常是日前调度和实时手动调节的组合。风电的接入，风电的自然属性(随机性和波动性)成为电力系统的不确定电源，在风电渗透率不高的情况下和负荷一样作为电力系统扰动源。随着风电渗透率的不断增加，风电的随机针对电力系统的影响越来越大。由于日前风功率预测误差较大，而风电的超短期预测(1h)精度较高，故国内外的研究大都是通过日前机组组合、日前经济调度和实时经济调度配合，多尺度进行调度。在含风电调度方面，国内外学者进行了大量基础研究，可分为四类。

(1)基于动态规划法的动态经济调度。此类算法解析度通常较高，算法逻辑清晰，计算速度快，例如等微增率算法。不足之处为此类算法通常将风电作为净负荷处理，如前所述，这种处理方式在风电渗透率不高的风电发展初期较为适用。在风电渗透率达到一定规模后，此类算法的弊端凸显。

(2)基于启发式方法的动态经济调度。此类算法的数学模型相对复杂，传统的算法难以求解，故通常适用蚁群、粒子群等启发式算法求解。启发式算法虽然智能性高，可以保证求出解，但每次优化结果是不确定的，有事无法保证得到的解是全局最优解，部分启发式算法不适合用于大规模系统求解。

(3)基于鲁棒优化方法的动态经济调度。鲁棒优化的核心思想是考虑最极端的边界情况，达到优化方法的鲁棒性。即牺牲一定经济性的条件下，保证模型对极端情况的有效处理。该方法可以有效求解发电计划，但由于鲁棒优化固有的鲁棒特性，在经济性等方面会有一定的降低。

(4)基于随机优化方法的动态经济调度。随机优化算法为将风电功率作为随机变量处理，场景法也是随机优化方法的一种常用方法。随机优化是一种处理含不确定性优化问题的有效方法，目前已广泛应用于含不确定性的电力系统经济调度问题中。如何准确描述风电功率的不确定性、有效求解对应的优化模型是含风电电力系统随机经济调度的关键问题。此类算法物理概念较为明确，但在计算速度方面目前还有很大的提升空间，尤其是对于对计算速度有较高要求的实时调度层面。

动态规划法由于其计算速度快，在实时调度中具有很大的潜力，但动态规划法大都将风电视为净负荷处理。在随机优化中，如何合理描述风电功率的不确定性并将其运用于电力系统的经济调度与优化运行中具有重要意义。描述风电功率的分布常用正态分布。对风电功率分布描述的准确性与对应随机动态经济调度模型的解的准确性密切相关。但目前尚未出现相关具有实用价值的技术方案。

发明内容

本发明针对现有技术的缺陷，提供基于通用分布的电力系统实时动态经济调度技术方案。

本发明提供的技术方案是：

1、一种基于通用分布的含风电电力系统实时动态经济调度方法，其特征在于，包括以下步骤：

步骤1：输入日前系统负荷预测数据和日前风电功率预测数据、火电机组运行参数、历史统计数据，所述历史统计数据包括不同风电功率预测水平下实际风电功率的通用分布参数α,β,γ；

步骤2：定义p_i,t为第i台火电机组t时刻的出力，火电机组的总数为I，i＝1,2,...,I，w_j,t为第j个风电场t时刻的计划出力，风电场的总数为J，j＝1,2,...,J，T为时刻的总数，t＝1,2,...,T；

计算从当前时刻也就是当前计算周期到后面1小时的调度时间域每个调度周期火电机组出力可能达到的范围，称为火电机组在t调度周期的可行域，该1小时的调度时间域为12个调度周期：

${\underset{\‾}{\mathrm{\α}}}_{it}\≤p_{it}\≤{\stackrel{\‾}{\mathrm{\α}}}_{it}$ 式一

式中：p_it表示t时刻第i台火电机组的出力；α _it和是表示t时刻第i台火电机组的出力下限和上限；

使用通用分布表征风功率，使用等微增率准则进行12个调度周期的风电火电静态经济调度，得到风电和火电的经济调度曲线，基于以下的目标函数以及约束条件：

目标函数： $\begin{array}{c}\min :f=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{f}_t=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\left(f_{g,t}+f_{w,t}\right)\\ =\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\[\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{C}_{g,i,t}\left(p_{it}\right)+\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{C}_{w,j,t}\left(w_{jt}\right)\]\end{array}$ 式二

火电机组出力功率约束条件：式三

火电机组爬坡约束条件：-δ_i≤p_it-p_i(t+1)≤δ_i式四

风电场出力功率约束条件：式五

系统功率平衡约束条件：式六

式中：f_tf_g,tf_w,t分别表示t时刻系统总发电成本、火电发电成本和风电成本；w_jt表示t时刻第j个风电场的出力；I表示火电机组数量,J表示风电场数量；P_t ^L表示t时刻系统总的负荷预测值；和p _i表示第i台火电机组的出力上下限；δ_i表示第i台火电机组的爬坡能力；表示第j个风电场的风电场装机容量；

其中，所述火电燃料成本C_g,i,t采用以下方式计算，

$C_{g,i,t}\left(p_{it}\right)=\frac{1}{2}{a}_i{p}_{it}^2+b_i{p}_{it}+c_i$ 式七

其中，p_it为第i台火电机组t时刻的出力，a_i,b_i,c_i为第i台火电机组的燃料成本系数；

所述风电成本C_w,j,t采用以下方式计算，

C_w,j,t(w_j,t)＝d_jw_jt+C_un,j,t(w_j,t)+C_ov,j,t(w_j,t)式八

其中，w_jt为第j个风电场t时刻的计划出力，d_j为第j个风电场的运行成本系数；

所述平均低估成本C_un,j,t采用风电场的平均弃风成本，采用以下方式计算，

$C_{un,j,t}\left(w_{j,t}\right)=k_{un,j}{\∫}_{w_{j,t}}^{w_{\max ,j}}(w_{av,jt}-w_{jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式九

其中，k_un,j为第j个风电场对应的低估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定；

所述平均高估成本C_ov,j,t采用系统的平均备用成本，采用以下方式计算，

$C_{ov,j,t}\left(w_{j,t}\right)=k_{ov,j}{\∫}_0^{w_{j,t}}(w_{jt}-w_{av,jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式十

其中，k_ov,j为第j个风电场对应的高估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定；

步骤3：利用步骤2得到的12个调度周期各自的风电场出力，修正火电机组各个调度周期的可行域；在新的火电机组可行域内，使用等微增率准则对火电机组的出力进行修正，保证火电机组爬坡限制并保证火电机组发电计划的稳定性；

步骤4：利用步骤2得到的风电场出力和步骤4得到的火电机组出力，输出火电机组和风电场的计划出力；

而且，所述火电燃料成本C_g,i,t采用以下方式计算，

$C_{g,i,t}\left(p_{it}\right)=\frac{1}{2}{a}_i{p}_{it}^2+b_i{p}_{it}+c_i$ 式十一

其中，p_it为第i台火电机组t时刻的出力，a_i,b_i,c_i为第i台火电机组的燃料成本系数；

而且，所述风电成本C_w,j,t采用以下方式计算，

C_w,j,t(w_j,t)＝d_jw_jt+C_un,j,t(w_j,t)+C_ov,j,t(w_j,t)式十二

其中，w_jt为第j个风电场t时刻的计划出力，d_j为第j个风电场的运行成本系数；

而且，所述平均低估成本Cun,j,t采用风电场的平均弃风成本，采用以下方式计算，

$C_{un,j,t}\left(w_{j,t}\right)=k_{un,j}{\∫}_{w_{j,t}}^{w_{\max ,j}}(w_{av,jt}-w_{jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式十三

其中，k_un,j为第j个风电场对应的低估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定；

而且，所述平均高估成本C_ov,j,t采用系统的平均备用成本，采用以下方式计算，

$C_{ov,j,t}\left(w_{j,t}\right)=k_{ov,j}{\∫}_0^{w_{j,t}}(w_{jt}-w_{av,jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式十四

其中，k_ov,j为第j个风电场对应的高估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定。

本发明利用通用分布刻画了风电功率的不确定性，以此为基础，提出了考虑风电低、高估成本的实时动态经济调度技术方案，包括以风电场的历史风电功率数据为基础，利用通用分布模型拟合不同风电功率预测水平下实际风电功率的分布，计及风电预测误差带来的惩罚成本，建立了基于通用分布的考虑风电低、高估成本的实时动态经济调度的随机优化模型；通过计算火电机组在各个调度周期的可行域，将对应的动态经济调度模型转化为12个静态经济调度模型，得到12个调度周期的火电机组和风电场出力，并将随机优化模型转化为确定性优化模型；使用动态规划算法修正火电机组的出力，求解动态经济调度模型，得到火电机组和风电场的实时计划出力。经验证，可知本发明技术方案的有效性，具有良好的推广价值和应用前景。

附图说明

图1是本发明实施例的通用分布与正态分布对于拟合风功率分布的拟合效果比较。

图2是本发明实施例的基于通用分布模型应用的含风电电力系统调度结构示意图。

图3是本发明实施例的风火静态经济调度等微增率算法。

图4是本发明实施例的火电机组在调度时间域的各个调度周期的运行可行域。

图5是本发明实施例的火电机组基本发电计划。

图6是本发明实施例的火电机组发电计划的修正过程。

图7是本发明实施例的火电机组发电计划的修正等微增率算法。

图8是本发明实施例的IEEE10机39节点系统图。

图9是本发明实施例的日负荷曲线。

图10是本发明实施例的丹麦海上风电预测与实测值。

图11是本发明实施例的丹麦海上风电通用参数表。

图12是本发明实施例的风电调度值和实测值。

图13是本发明实施例的两种调度模式的系统成本比较。

图14是本发明实施例的不同发电机数目下的算法计算时间。

图15是本发明实施例的采用不同分布表征风功率的计算时间。

图16是本发明实施例的调度模型对于保证火电机组发电计划稳定性的结果比较分析。

图17是本发明实施例的采用不同方式表征风功率的随机性对系统频率的影响比较。

具体实施方式

为了使本发明实施例的目的、技术方案、优点更加清晰，下面将结合本发明实施例和附图来介绍本发明的技术方案。

本发明提供的技术方案是一种基于通用分布的考虑风电低、高估成本的实时动态经济调度方法，原理如下：

将并网风电场历史风电功率的预测和实测数据标幺化，根据风电功率的预测值的不同，对历史风电功率数据进行分箱，在不同风电功率预测水平下，利用通用分布函数拟合不同预测箱下实测风电功率的分布，得到对应的通用分布参数；

考虑系统有功功率平衡约束、火电机组和风电场的出力限制、火电机组爬坡速率约束，建立基于通用分布的计及风电低、高估成本的含风电电力系统的实时动态经济调度模型；

基于预测风电功率的随机经济调度模型，使用通用分布表征风功率，并计算各个调度周期火电机组的出力可行域，使用等微增率准则进行12个调度周期的风电火电静态经济调度，得到考虑风电低、高估成本的实时随机动态经济调度模型的风电和火电的经济调度曲线。

得到各个调度周期风电场出力计划后，基于通用分布的随机动态经济调度问题转化成确定性经济调度模型。利用步骤二得到的12个调度周期各自的风电场出力，修正火电机组各个调度周期的可行域。在新的火电机组可行域内，使用等微增率准则对火电机组的出力进行修正，保证火电机组爬坡限制并保证火电机组发电计划的稳定性。最终输出日前火电机组和风电场的计划出力曲线。

首先，为便于实施参考起见，介绍基于通用分布的随机优化模型：

一、通用分布模型和形状参数表的生成。

1.1、通用分布模型。

若连续型随机变量X服从一个形状参数为α、β和γ的通用分布，则记为：

X～V(α,β,γ)(11)

其中，形状参数α、β和γ满足：

α＞0,β＞0,-∞＜γ＜+∞(12)

通用分布的概率密度函数(ProbabilityDensityFunction，PDF)定义为：

$f\left(x\right)=\frac{{\mathrm{\α\βe}}^{-\mathrm{\α}(x-\mathrm{\γ})}}{{(1+e^{-\mathrm{\α}(x-\mathrm{\γ})})}^{\mathrm{\β}+1}}---\left(13\right)$

通用分布的累积分布函数(CumulativeDistributionFunction，CDF)定义为：

F(x)＝(1+e^-α(x-γ))^-β(14)

1.2、风电概率密度直方图。

1.2.1、绘制风电概率密度直方图。

统计风电场预测和实测历史数据，得到风电功率实测时间序列和预测时间序列以后，以风电场最大出力为基值，将数据转化为标幺值，供后文研究使用。

本文所涉及的几个关键名词的定义如下：

(1)风电预测功率等级：将风电预测功率的取值范围0～1p.u.等分为M个小区间，各区间的宽度为1/M，取值范围分别为[0,1/M],[1/M,2/M],…,[(M-1)/M,1]，每个区间称作一个“风电预测功率等级”，第m级(m＝1,2,…,M)对应的预测功率取值范围是[(m-1)/M,m/M]。

(2)第m个预测功率等级对应的风电实测功率原始数据集合：在原始数据集合中的每一个时段，风电场有一个预测功率值和一个实测功率值，二者组成数对。对于每个属于第m级预测功率等级的风电预测功率值，都有一个与之成对的实测功率值，所有这些实测功率值构成的集合称作“第m个预测功率等级对应的风电实测功率原始数据集合”。

值得一提的是，M为用户自定义的正整数，可用的历史数据越多，M取值越大，反之越小。另一方面，如果历史数据总量固定，M取值越大，每个预测功率等级对应的风电实测功率原始数据集合中样本数量越少，反之越多。对于每个预测功率等级对应的风电实测功率原始数据集合，如果该集合的样本数量过少，将降低统计结果的可信度，这就失去了统计的意义。因此，为保障后文中统计分析的效果，应使每个预测功率等级对应的风电实测功率原始数据集合中样本数量足够多。经试验，对于本文采用的风电场数据，取M＝25可以达到较好的统计效果。

使用直方图表征风电功率概率分布，概率密度直方图的绘制步骤如下：

Step1：收集数据。对于不同的风电预测功率等级，实际风电功率的概率分布不同。对于第m个预测功率等级，选取其对应的风电实测功率原始数据集合作为样本集合。

Step2：确定直方图的总宽度。传统做法是：由样本集合中的最大值与最小值之差求得极差，直方图的总宽度(即组距与组数的乘积)取略大于极差的数。由于风电出力在0～1p.u.范围内，为简化问题，本文假设直方图的跨度也在0～1p.u.范围内，于是总宽度就等于1p.u.。

Step3：确定直方图的组数、组距和各组界限。假设直方图的组数为N；每组的宽度为1/N，即组距；第n组(n＝1,2,…,N)的下界为(n-1)/N,上界为n/N。注意，组数的确定要适当，组数过少可能会使数据的分布规律失真，组数过多可能会影响数据分布规律的明显性。经试验，对于本文采用的风电场数据，取N＝50可以达到较好的统计效果。

Step4：确定各组频数(FrequencyCounts)、概率和概率密度。把样本集合中的数据按照大小划分到直方图各组，统计各组中数据的个数，即频数；各组频数除以总样本数，就得到各组概率，用prob_n(n＝1,2,…,N)表示；各组概率除以组距1/N，就得到各组的概率密度，即N·prob_n。

Step5：绘制概率密度直方图。横坐标为风电实际功率，纵坐标为概率密度，直方图每组对应一个矩形，矩形的宽度为组距，高度为各组的概率密度。

由于直方图是由真实数据统计而来，没有采用任何近似方法，故本文将使用直方图表征的风电功率概率分布称为“实际分布”。实际分布本质上是离散型分布，其PDF定义如式(5)所示：

$f_{act}\left(x\right)=\left\{\begin{array}{cc}{\mathrm{prob}}_n\&CenterDot;N& \left(\frac{n-1}{N}\&le;x<\frac{n}{N},n=1,2,\mathrm{...},N\right)\\ {\mathrm{prob}}_N\&CenterDot;N& \left(x=1\right)\end{array}\right.---\left(15\right)$

对直方图各组概率离散求和就可以得到实际分布的CDF，如式(6)所示：

$F_{act}\left(x\right)=\left\{\begin{array}{cc}\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{\mathrm{prob}}_i& \left(\frac{n-1}{N}\&le;x<\frac{n}{N},n=1,2,\mathrm{...},N\right)\\ 1& \left(x=1\right)\end{array}\right.---\left(16\right)$

1.2.2、曲线拟合得到通用分布形状参数表。

正态和贝塔分布近似模型的参数均是通过所统计的样本集合的平均值μ和标准差σ直接或间接推导求得。这使得正态和贝塔分布的PDF/CDF曲线形状与μ和σ强相关，而与实际分布的PDF/CDF曲线形状弱相关，导致正态和贝塔分布近似模型模拟实际分布的效果较差。

为了能够更精确地模拟实际分布，本文提出：当使用通用分布近似表征实际分布时，不通过μ和σ求取通用分布的形状参数α、β和γ，而是使用通用分布PDF/CDF曲线直接对实际分布PDF/CDF曲线进行最小二乘拟合得到α、β和γ，请见图1，通用分布拟合效果相比于其他分布能够更加准确的拟合风功率实际分布。由于实际分布是以风电预测功率为均值的条件概率分布，故对于不同的风电预测功率等级，通用分布的形状参数不同。对于某一预测时间尺度下的风电场，可以通过拟合各风电预测功率等级下的实际分布PDF或CDF，生成一个通用分布的形状参数表(Look-UpTable)。

需要强调的是，PDF拟合法和CDF拟合法得到的通用分布形状参数大小可能有一定差异。那么，到底应该使用哪种拟合方法获取通用分布的形状参数呢？在经济调度问题中，一些经典算法通常需要计算低/高估风电代价的偏导数并将约束条件线性化，低/高估风电代价的偏导数和机会约束是关于风电功率概率分布CDF或其逆函数的表达式。所以，CDF拟合法得到的通用分布形状参数更实用、更直接和更精确。当求解实际电力系统经济调度问题时，应尽量采用CDF拟合法求取通用分布的形状参数。

另外需要注意的是，最小二乘拟合解(通用分布的三个形状参数α、β和γ)并不唯一，主要取决于初值如何选取。如果随机选取每个风电预测功率等级条件下的通用分布形状参数的初值，可能使某个预测功率等级对应的拟合解与其它预测功率等级对应的拟合解在正负号或数量级方面产生较大差异。当然，虽然有些拟合解的正负号或数量级“特殊”，但只要相应的通用分布曲线能够较好地拟合实际分布曲线，这些看似“特殊”的拟合解就具有实用性，所选取的初值就是可以接受的。

二、基于通用分布的风电火电静态经济调度。

2.1、基于通用分布的风功率随机性表征。

基于通用分布模型应用的含风电电力系统调度结构示意图如图2所示，在调度计划制定之前，风电场提供实测功率历史数据，风电预测系统提供预测功率历史数据，采用上节所述的方法，通过使用通用分布CDF对实际分布CDF进行最小二乘拟合，得到通用分布形状参数表。在调度周期内，风电功率预测系统给定当前风电功率预测值，该预测值属于某一预测功率等级，从通用分布形状参数表中可查询当前预测功率等级条件下的形状参数，进而得到当前预测功率等级条件下的通用分布模型，然后将其应用到调度模型和算法中。最终得到调度系统向风电场和常规机组发送的调度计划指令。

如果风电场的地理布局、装机容量和风机类型等条件过去没有太大的变化，则为生成通用分布的形状参数表，图2中的统计分析环节应采用全部的风电预测和实测功率历史数据。为提高统计分析的精度，应定期(如每天或每星期)将最新的数据补充到历史数据库中。生成通用分布形状参数表的过程并不会影响调度系统计算和决策的效率，原因如下：第一，生成通用分布形状参数表的耗时很短，例如，在一台Core-i72.70GHz笔记本电脑中的Matlab2012a环境下耗时不到1s；第二，形状参数表无需实时更新，每天或每星期更新一次即可。

2.2、风电火电静态经济调度。

对于实时经济调度，其目标函数和约束条件如下：

目标函数：

$\begin{array}{c}\min :f=\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{f}_t=\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\left(f_{g,t}+f_{w,t}\right)\\ =\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\&lsqb;\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}{C}_{g,i,t}\left(p_{it}\right)+\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{C}_{w,j,t}\left(w_{jt}\right)\&rsqb;\end{array}---\left(2\right)$

火电机组出力功率约束：

${\underset{\&OverBar;}{p}}_i\&le;p_{it}\&le;{\stackrel{\&OverBar;}{p}}_i---\left(3\right)$

火电机组爬坡约束：

-δ_i≤p_it-p_i(t+1)≤δ_i(4)

风电场出力功率约束：

$0\&le;w_{jt}\&le;w_j^r---\left(5\right)$

系统功率平衡约束：

$\underset{i=1}{\overset{I}{\&Sigma;}}{p}_{it}+\underset{j=1}{\overset{J}{\&Sigma;}}{w}_{jt}=P_t^L---\left(6\right)$

式中：f_tf_g,tf_w,t分别表示t时刻系统总发电成本、火电发电成本和风电成本；p_it表示t时刻第i台火电机组的出力；w_jt表示t时刻第j个风电场的出力；I表示火电机组数量,J表示风电场数量；P_t ^L表示t时刻系统总的负荷预测值；和p _i表示第i台火电机组的出力上下限；δ_i表示第i台火电机组的爬坡能力；表示第j个风电场的风电场装机容量。

目标函数中风电成本展开为：

$\begin{array}{c}{f}_{w,t}=\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}\&lsqb;d_j{w}_{jt}+k_{un,j}{\&Integral;}_{w_{j,t}}^{w_{\max ,j}}\left(w_{av,jt}-w_{jt}\right)f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}\\ +k_{ov,j}{\&Integral;}_0^{w_{j,t}}\left(w_{jt}-w_{av,jt}\right)f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}\&rsqb;\end{array}---\left(17\right)$

式中：d_j是风电场j的运行费用；k_un,j是风电场j的低估代价费用；k_ov,j是风电场j的高估代价费用；w_av,jt是风电场j的在t时刻的实际可能出力。

风电成本等微增率：

${\mathrm{\&lambda;}}_{j,t}=\frac{\&part;f_{w,t}}{\&part;w_{jt}}=d_j-k_{un,j}+(k_{un,j}+k_{ov,j})F_{jt}\left(w_{jt}\right)---\left(18\right)$

根据(144)式的通用分布函数CDF，(18)式可化为：

${\mathrm{\&lambda;}}_{j,t}=c_{1t}+c_{2t}{(1+e^{-{\mathrm{\&alpha;}}_t\left(w_{jt}-{\mathrm{\&gamma;}}_t\right)})}^{-{\mathrm{\&beta;}}_t}---\left(19\right)$

根据火电成本的二次曲线，火电成本等微增率可按下式计算：

${\mathrm{\&lambda;}}_{i,t}=\frac{\&part;f_{g,t}}{\&part;p_{it}}=a_i{p}_{it}+b_i---\left(20\right)$

式中：a_i和b_i表示火电机组i的二次成本函数系数。

根据火电机组初始运行点，可通过式(3)(4)(6)计算出调度域所有调度周期火电机组的功率范围，即可行域，如图4。按照可行域的计算方法，在计算可行域上限的时候，风电功率取最小值，即0；在计算可行域下限的时候，风电功率取最大值，即风电场装机容量。暂不考虑可行域中的爬坡因素，火电机组可行域上下限可由下式表示：

${\underset{\&OverBar;}{\mathrm{\&alpha;}}}_{it}\&le;p_{it}\&le;{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_{it}---\left(1\right)$

在火电机组可行域(1)和风电场装机容量(5)内，可通过经典的等微增率计算调度域(1h)内12个调度周期(5min)各自的风电场和火电机组出力，如图3。

3、基于动态规划法的火电机组出力修正方法。

3.1、火电机组基本发电计划。

火电机组基本发电计划的引入是降低不同调度周期对于调度域同一调度周期发电计划的更改，保证发电计划的稳定性。故某一调度周期的基本发电计划由上一调度周期计算得到但未使用的11个调度周期的火电出力和此调度周期计算的最远调度周期的火电计划组成，如图5所示。

3.2、火电机组运行可行域和系统功率平衡修正。

风电火电静态经济调度得到调度域各调度时间风电发电计划后，火电机组运行可行域可以使用风电发电计划进行修正。同时，系统功率平衡也可修正为：

$\underset{i=1}{\overset{I}{\&Sigma;}}{p}_{it}=P_t^L-\underset{j=1}{\overset{J}{\&Sigma;}}{w}_{jt}=P_t^{NL}---\left(21\right)$

式中：P_t ^NL表示t调度周期下的净负荷(系统总负荷减去风电出力)。

3.3、火电机组出力修正。

火电机组发电计划可以通过在每调度周期基本发电计划加上增量表示：

$p_{it}=p_{it}^0+{\mathrm{\&Delta;p}}_{it},(t=1...T)---\left(22\right)$

注意到是基于上一计算周期的净负荷P_t ^NL，0求得的,基于本计算周期净负荷P_t ^NL得到p_it，故Δp_iT,…Δp_i1即表示两个计算周期对于同一调度周期的净负荷变化，体现在如下的功率平衡中：

$\underset{i}{\overset{I}{\&Sigma;}}{\mathrm{\&Delta;p}}_{it}-P_t^{NL}-P_t^{NL,0},(t=1...T)---\left(23\right)$

式中： $P_1^{NL}\&equiv;P_1^{NL,0}$

为计算Δp_it,定义爬坡率变化值：ΔD_it

ΔD_it＝Δp_it-Δp_i(t+1)(24)

ΔD_i1，ΔD_i2...ΔD_iT可以通过从较远的调度周期向较近的调度周期通过等微增率准则进行计算。在计算过程中，当计算到t＝s阶段，即计算ΔD_is，假设对于t>sΔD_it＝0。即假设：ΔD_i(s+1)≡0

这样，在修正过程中，只有ΔD_is(i＝1...I)为变量。在修正到每一个调度周期时，ΔD_is的值决定所有的Δp_it(t＝1...s)值：

${\mathrm{\&Delta;p}}_{it}=\underset{k=t}{\overset{s}{\&Sigma;}}{\mathrm{\&Delta;D}}_{ik}={\mathrm{\&Delta;D}}_{is}+\underset{k=t}{\overset{s-1}{\&Sigma;}}{\mathrm{\&Delta;D}}_{ik}---\left(25\right)$

当ΔD_is变化时,t＝s...1调度周期的火电机组发电计划根据式(25)平行移动，如图6。由于等微增率准则的计算速度非常快，故修正过程耗时很短。在每一个调度周期的修正过程中，等微增率准则的各项参数如图7。

其中：

${\mathrm{\&epsiv;}}_s=(P_s^{NL}-P_{s+1}^{NL})-(P_s^{NL,0}-P_{s+1}^{NL,0})---\left(26\right)$

$\begin{array}{c}{g}_s=\frac{a_i}{s}\left(\underset{m=1}{\overset{s-1}{\mathrm{\&Sigma;}}}\underset{t=m}{\overset{s-1}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;D}}_t+\underset{t=1}{\overset{s}{\mathrm{\&Sigma;}}}{p}_{it}^0\right)\\ =\frac{1}{s}\left(\left(s-1\right)g_{s-1}+a_i\left(\left(s-1\right){\mathrm{\&Delta;D}}_{s-1}+P_{is}^0\right)\right)\\ g_1=a_i{p}_{i1}^0\end{array}---\left(27\right)$

当s＝12,所有调度周期的ΔD_it计算完毕,修正后的火电机组出力即：

$p_{it}=p_{it}^0+\underset{k=t}{\overset{T}{\&Sigma;}}{\mathrm{\&Delta;D}}_{ik},(i=1...I,t=1...T)---\left(28\right)$

四、算例及仿真。

本发明以IEEE10机39节点系统标准算例为例进行了仿真计算。系统10台发电机中，G1为调频机组，不进行调度，G7、G9只相应日前调度出力曲线，不参与实时动态经济调度，如图8所示。风电使用丹麦海上风电实际数据，在8、15和39接入，设置风电运行成本24$/MWh，风电低估成本60$/MWh，风电高估成本120$/MWh。实时调度算法计算一天的机组调度出力，整个系统使用RTDS进行一天的仿真计算，7台BLO机组(即实时调度机组)的出力由实时调度计算结果决定，忽略系统负荷预测误差，系统一天负荷值见图9，系统一天风电预测和实测值见图10，通用分布参数见图11。比较使用通用分布将风电作为调度量和传统将风电作为净负荷两种调度模式的调度结果。

4.1、风电调度出力。

如图12，系统前三小时风电调度出力与预测出力的关系。根据分析，若风电低高估成本极大且相等，PDF曲线不偏轴时，则风电的调度量应在使低高估期望相等的点，即CDF＝0.5点，也即PDF最高点，故对于风电分布偏轴的PDF曲线，风电调度量应在尖峰偏0.5p.u.处。而对于将风电作为净负荷的算法，其PDF曲线相当于只在其调度值的箱内，故将风电作为调度量的结果应为风电预测值略靠近0.5p.u.。而实际系统运行过程中，低估成本要高于高估成本，从此因素考虑，风电调度值应略低于风电预测值。

综合如上两个因素，可以得到风电调度值与期望值的一个大致的规律：当风电低于0.5p.u.时，风电调度值在预测值上下；当风电高于0.5p.u.时，风电调度值低于风电预测值。图12结果证明了这个分析。5—105min时，风电调度值在预测值上下，105min—180min，风电调度值低于风电预测值。

4.2、系统运行成本。

如图13，系统前三小时风电作为调度量和风电作为净负荷两种调度模式的成本比较，可以看出，风电作为调度量相比于风电作为净负荷，可以大大减少由风电随机性带来的系统对风电低高估成本，减少总的系统发电成本。

4.3、调度模型计算时间。

如图14，火电机组数量与算法计算时间的关系，可以看出，调度算法在火电机组数量为7台时，计算时间只有0.097s，且计算时间和火电机组数量呈线性关系(区别于二次规划法，计算时间和火电机组数量的平方呈线性关系)。算法计算时间满足实时调度的要求。

如图15，通用分布和正态分布对于表征风功率，并将风功率作为调度量进行调度算法计算时间。正态分布也可以用来描述风电功率的随机性，但正态分布相比于通用分布的拟合效果较差。除此之外，虽然正态分布通过查表也可计算风电成本等微增率，但由于每次迭代计算中通用分布的解析表达式计算时间比正态分布的查表计算时间快几十倍，且通用分布由于有解析表达式，可以进行线性化减少迭代计算次数。故由于这两个方面，使用基于通用分布的含风电经济调度计算时间远远小于基于正态分布的含风电经济调度计算时间。

4.4、调度模型对于保证火电机组发电计划稳定性的分析。

在实时调度中，每个计算周期(5min)计算下后面一小时的发电计划，即后面12个调度周期(12×5min)的发电计划，但只使用下一个调度周期的发电计划。故每个调度周期的发电计划被重复计算了12次才被使用，假设现在是0min，需要计算5、10…60min这12个调度周期的发电计划。第60min的发电计划，在0、5…55min的计算周期均被计算，但只使用55min计算周期的发电计划结果。在这12个计算周期计算的第60min的计划，必须有一定的稳定性，尽量减少不必要的发电计划更改，否则会使发电机无所适从。

见图16，比较将风电作为净负荷进行调度的本文方法和经典的二次规划算法(将风电作为净负荷处理)，对于火电机组调度出力稳定性的比较。可以看出，9:35制定的9:40到10:35的发电计划，本文方法和二次规划法的对于相同的一台发电机其处理曲线相同。在9:40制定9:45到10:40的发电计划，本文方法对于9:45到10:35的发电计划和9:35制定的是相同的，即没有出现发电计划变动；而使用二次规划算法制定的发电计划在10:00到10:35发生了变动。在9:45，二次规划法制定的发电计划也出现了变动，而本文方法在三个计算时间对于相同调度周期的出力结果是相同的，即本文算法相比于传统算法能有效保证发电计划的稳定性。

4.5、采用不同方式表征风功率的随机性对系统频率的影响比较。

分别计算三种调度模式：1)将风电作为净负荷处理(即不进行调度)，2)使用正态分布表征风电功率并进行调度，3)使用通用分布表征风电功率并进行调度。三种调度模式均在相同的负荷值(见图9)、相同的风功率预测值和实测值(见图10)下进行一天的滚动计算，记录系统频率偏差超过某一设定值的次数，用以表征系统有功平衡情况。

见图17结果：1)将风电作为净负荷系统频率偏差最多，2)使用正态分布表征风电功率次之，3)使用通用分布表征风电功率最少。经过之前分析，由于风电预测的特性，风电作为调度量的结果应为风电预测值略靠近0.5p.u.，风电随机性对于系统有功平衡的影响最小。1)将风电作为净负荷，即调度量结果就是预测值，故最不利于系统有功平衡；2)使用正态分布表征风电功率有使风电预测值略靠近0.5p.u.的效果，但由于正态分布为对称分布，不具有通用分布的偏轴特性，所以其对系统有功平衡的保证能力强于将风电作为净负荷，但不如3)使用通用分布表征风电功率。

Claims (1)

1.一种基于通用分布的含风电电力系统实时动态经济调度方法，其特征在于，包括以下步骤：

步骤1：输入日前系统负荷预测数据和日前风电功率预测数据、火电机组运行参数、历史统计数据，所述历史统计数据包括不同风电功率预测水平下实际风电功率的通用分布参数α,β,γ；

步骤2：定义p_i,t为第i台火电机组t时刻的出力，火电机组的总数为I，i＝1,2,…,I，w_j,t为第j个风电场t时刻的计划出力，风电场的总数为J，j＝1,2,...,J，T为时刻的总数，t＝1,2,...,T；

计算从当前时刻也就是当前计算周期到后面1小时的调度时间域每个调度周期火电机组出力可能达到的范围，称为火电机组在t调度周期的可行域，该1小时的调度时间域为12个调度周期：

${\underset{\&OverBar;}{\mathrm{\&alpha;}}}_{it}\&le;p_{it}\&le;{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_{it}$ 式一

式中：p_it表示t时刻第i台火电机组的出力；α _it和是表示t时刻第i台火电机组的出力下限和上限；

使用通用分布表征风功率，使用等微增率准则进行12个调度周期的风电火电静态经济调度，得到风电和火电的经济调度曲线，基于以下的目标函数以及约束条件：

目标函数： $\begin{array}{c}\min :f=\underset{t=1}{\overset{T}{\&Sigma;}}{f}_t=\underset{t=1}{\overset{T}{\&Sigma;}}(f_{g,t}+f_{w,t})\\ =\underset{t=1}{\overset{T}{\&Sigma;}}\&lsqb;\underset{i=1}{\overset{I}{\&Sigma;}}{C}_{g,i,t}\left(p_{it}\right)+\underset{j=1}{\overset{J}{\&Sigma;}}{C}_{w,i,t}\left(w_{jt}\right)\&rsqb;\end{array}$ 式二

火电机组出力功率约束条件：式三

火电机组爬坡约束条件：-δ_i≤p_it-p_i(t+1)≤δ_i式四

风电场出力功率约束条件：式五

系统功率平衡约束条件：式六

式中：f_t，f_g,t，f_w,t分别表示t时刻系统总发电成本、火电发电成本和风电成本；w_jt表示t时刻第j个风电场的出力；I表示火电机组数量,J表示风电场数量；表示t时刻系统总的负荷预测值；和p _i表示第i台火电机组的出力上下限；δ_i表示第i台火电机组的爬坡能力；表示第j个风电场的风电场装机容量；

其中，所述火电燃料成本C_g,i,t采用以下方式计算，

$C_{g,i,t}\left(p_{it}\right)=\frac{1}{2}{a}_i{p}_{it}^2+b_i{p}_{it}+c_i$ 式七

其中，p_it为第i台火电机组t时刻的出力，a_i,b_i,c_i为第i台火电机组的燃料成本系数；

所述风电成本C_w,j,t采用以下方式计算，

C_w,j,t(w_j,t)＝d_jw_jt+C_un,j,t(w_j,t)+C_ov,j,t(w_j,t)式八

其中，w_jt为第j个风电场t时刻的计划出力，d_j为第j个风电场的运行成本系数；

所述平均低估成本C_un,j,t采用风电场的平均弃风成本，采用以下方式计算，

$C_{un,j,t}\left(w_{j,t}\right)=k_{un,j}{\&Integral;}_{w_{j,t}}^{w_{\max ,j}}(w_{av,jt}-w_{jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式九

其中，k_un,j为第j个风电场对应的低估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定；

所述平均高估成本C_ov,j,t采用系统的平均备用成本，采用以下方式计算，

$C_{ov,j,t}\left(w_{j,t}\right)=k_{ov,j}{\&Integral;}_0^{w_{j,t}}(w_{jt}-w_{av,jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式十

其中，k_ov,j为第j个风电场对应的高估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定；

步骤3：利用步骤2得到的12个调度周期各自的风电场出力，修正火电机组各个调度周期的可行域；在新的火电机组可行域内，使用等微增率准则对火电机组的出力进行修正，保证火电机组爬坡限制并保证火电机组发电计划的稳定性；

步骤4：利用步骤2得到的风电场出力和步骤4得到的火电机组出力，输出火电机组和风电场的计划出力；

而且，所述火电燃料成本C_g,i,t采用以下方式计算，

$C_{g,i,t}\left(p_{it}\right)=\frac{1}{2}{a}_i{p}_{it}^2+b_i{p}_{it}+c_i$ 式十一

其中，p_it为第i台火电机组t时刻的出力，a_i,b_i,c_i为第i台火电机组的燃料成本系数；

而且，所述风电成本C_w,j,t采用以下方式计算，

C_w,j,t(w_j,t)＝d_jw_jt+C_un,j,t(w_j,t)+C_ov,j,t(w_j,t)式十二

其中，w_jt为第j个风电场t时刻的计划出力，d_j为第j个风电场的运行成本系数；

而且，所述平均低估成本Cun,j,t采用风电场的平均弃风成本，采用以下方式计算，

$C_{un,j,t}\left(w_{j,t}\right)=k_{un,j}{\&Integral;}_{w_{j,t}}^{w_{\max ,j}}(w_{av,jt}-w_{jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式十三

其中，k_un,j为第j个风电场对应的低估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定；

而且，所述平均高估成本C_ov,j,t采用系统的平均备用成本，采用以下方式计算，

$C_{ov,j,t}\left(w_{j,t}\right)=k_{ov,j}{\&Integral;}_0^{w_{j,t}}(w_{jt}-w_{av,jt})f_j\left(w_{av,jt}\right){\mathrm{dw}}_{av,jt}$ 式十四

其中，k_ov,j为第j个风电场对应的高估成本系数，w_av,j,t为第j个风电场t时刻的实际可能出力，f_j(w_av,j,t)为第j个风电场在对应风电预测水平下实际可能出力的概率密度函数，表达形式为通用分布的概率密度函数，根据t时刻对应通用分布参数α,β,γ确定。