基于冒泡排序的风储联合发电系统储能元件优化运行方法与流程

本发明涉及一种基于冒泡排序的风储联合发电系统储能元件优化运行方法，属于风电控制技术领域。

背景技术：

由于风速的波动性和不确定性，使得大规模风电并网在节能减排的同时，也给电力系统带来一定的风险。目前，随着储能元件的成本降低及储能技术的进一步提升，近年来越来越多的研究人员致力于将储能技术应用于改善风电的波动性。现已有相关研究将储能技术接入包含风电场的发电系统，并进行可靠性评估。

一般包含风电场和储能元件的发电系统可靠性分析主要分为系统内各元件状态抽取、风电场各项参数、储能元件的参数及运行策略。目前关于元件状态抽取主要分为解析法和模拟法，模拟法又称蒙特卡洛方法，蒙特卡洛方法主要分为序贯蒙特卡洛方法和非序贯蒙特卡洛方法。关于风电场各项参数的研究，主要是指风速的产生、风机的位置、风电场间的影响因子、尾流效应、允许风电最大接入比例等。关于储能元件部分，目前研究多集中于储能元件的充放电速度、最大放电深度、容量、接入位置以及储能策略。

目前已有相关文献对储能元件的运行方式即储能策略进行研究，但是大多数文件的储能策略是传统储能策略或者满充满放策略，并且往往将多个储能元件等效为一个或两个大型储能元件，忽略元件的故障停运，并且考虑的情况较少。少量文献考虑了风电接入比例，制定相应的储能策略，但是一般情况下忽略在风电场提供的电能超过风电接入比例，并且不满足系统负荷情况下的溢出风电，导致这部分电量的损失。极少数文献考虑到了这部分风电，但是仍是按照元件原始编号对这部分风电进行存储，并未按照储能元件内部剩余电量的优先顺序对储能策略进行优化。

技术实现要素：

本发明所要解决的技术问题是克服现有技术的缺陷，提供一种基于冒泡排序的风储联合发电系统储能元件优化运行方法，将所有处于故障状态的储能元件容量清零并使之不参与充放电过程，按照剩余的储能元件在这一时刻内部存储电量的由高到低的顺利进行充电，在风电场提供的电能超过风电接入比例，并且不满足系统负荷情况下，使风电利用率进一步提升。

为解决上述技术问题，本发明提供一种基于冒泡排序的风储联合发电系统储能元件优化运行方法，包括以下步骤：

1)参数初始化，所述参数包括：时间参数t、仿真次数n、系统年均缺电量eens、系统缺电概率lolp和年均缺电时间lolt，初始化结果为：t为一年8760小时，n为1000次，eens、lolp、lole初始值均为0；

2)根据ieee-rbts可靠性测试系统，运用序贯蒙特卡洛方法模拟得到常规机组的“故障—运行”状态序列；将常规机组的“故障—运行”状态序列乘以常规机组的额定功率，得到一年内每台常规机组每小时的出力；将每小时各个常规机组的出力求和，得到整个rbts系统一年的时序出力；

3)根据风电场所在区域的历史风速数据，基于最小信息准则确定风电场的自回归滑动平均模型即arma模型阶数；

4)根据所述步骤3)得到的风电场的arma模型阶数及风电场所在区域的历史风速数据，在matlab中调用armax函数模拟产生风电场的风速数据；

5)根据风机的故障率、修复率，运用序贯蒙特卡洛方法模拟得到风机的“故障—运行”状态序列；

6)根据步骤4)的风速数据和步骤5)的状态序列，结合风机的风速—功率特性曲线，得到风电场一年时间内的功率序列；

7)根据步骤2)得到的常规机组的时序出力gc(t)和步骤6)得到的风电场的功率序列gw(t)，系统负荷为pl(t)，考虑风电允许接入比例得到每个时刻下的系统过剩电量δe(t)和过剩风电电量δgw(t)；

8)根据第m台储能元件的故障率和修复率，运用序贯蒙特卡洛方法模拟得到储能元件的“故障—运行”状态序列；

9)利用步骤8)的状态序列中状态值是否为0判断当前时刻是否为故障状态，若为0表示该储能元件处于故障状态，则当前时刻的储能元件的电量清零；

10)当t时刻风电和常规机组出力不满足条件a时，根据系统过剩电量δe(t)，确定t时刻储能系统实际的充放电量；转入步骤12)；

11)当t时刻风电和常规机组出力满足条件a时，根据系统过剩电量δe(t)和过剩的风电电量δgw(t)，确定t时刻储能系统实际的充放电量；转入步骤12)；

12)将确定的储能元件充电电量或放电电量δeb(i,t)和t时刻剩余电量eb(i,t)相加即可得到下一时刻各储能元件的剩余电量δeb(i,t+1)；

13)若执行步骤10)或者步骤11)后，系统过剩电量δe(t)<0，则此时刻系统处于缺电状态，缺电量记为se(t)，且se(t)＝-δe(t)，同时系统缺电时间加1，lole＝lole+1；

14)令时间递增1，即t＝t+1，重复步骤10)到步骤13)，直到t>8760；

15)计算一年的缺电量eens，其值为一年内所有se(t)的总和；

16)计算仿真叠加缺电量eens＝eens+eens，重复步骤2)到步骤16)，每循环一次，仿真次数加1，直到仿真次数达到n，其中，eens的初始值为0；

17)计算系统年均缺电量eens＝eens/n，系统年均缺电小时数lolt＝lole/n，系统缺电概率lolp＝lolt/8760。

前述的步骤2)，步骤5)和步骤8)中，运用序贯蒙特卡洛方法模拟得到常规机组/风机/储能元件的“故障—运行”状态序列，具体步骤如下：

2-1)，在[0,1]之间产生随机数u1、u2；

2-2)，根据第i台常规机组/风机/储能元件的故障率λi和修复率μi，得到一次平均无故障工作时间一次平均修复时间mttf＝-ln(u2)/μi；

2-3)，重复步骤2-2)，直到该常规机组/风机/储能元件所有的mttr和mttf之和超过一年8760个小时；

2-4)，将8760个小时中所有mttr时间内常规机组状态记为1，表示常规机组正常运行；mttf时间内常规机组状态记为0，表示常规机组故障，得到一年内继续的常规机组/风机/储能元件的“故障—运行”状态序列；

2-5)，重复步骤2-2)到2-4)，直到获取rbts系统内所有常规机组/风机/储能元件的状态序列。

前述的步骤6)风电场一年时间内的功率序列，具体包含以下步骤：

6-1)，根据切入风速vci和额定风速vr，计算特性参数a、b、c：

6-2)根据步骤4)中模拟出的风速数据，得到正常运行下，每台风机每小时的功率pt，将pt按照时间顺序排列即可得到风机一年的功率序列p，

其中，vt表示风速，vco表示切出风速，pw为风机的额定功率；

6-3)将步骤5)得到的每台风机的“故障—运行”状态序列点乘功率序列p，即可得到考虑风机故障状态的功率序列p'；

6-4)将每小时内各风机的功率序列p'求和即可得到整个风电场每小时的实际功率。

前述的步骤7)中，求解系统过剩电量δe(t)和过剩风电电量δgw(t)，具体包含以下步骤：

7-1)考虑“风电超出允许接入比例而系统吸收允许接入的风电后仍处于缺电状态”这一条件，将此条件记为条件a：

7-2)风电和常规机组出力不满足条件a时，系统过剩电量δe(t)如下，δe(t)正值表示过剩的电量，δe(t)为负表示系统实际缺电的电量：

7-3)当风电和常规机组出力满足条件a时，系统过剩电量δe(t)和过剩风电电量δgw(t)如下：

前述的步骤10)中，储能系统实际的充放电量，具体求解过程如下：

10-1)参数初始化，储能系统内共有m台储能元件，每台储能元件的额定容量为ee，t时刻第i台储能元件剩余电量为eb(i,t)，初始时刻第i台储能元件剩余电量eb(i,0)为0，即所有储能元件初始时刻的剩余电量均为0，每小时的最大充放电功率为pm；

10-2)根据步骤7-2)得到系统t时刻的过剩的电量或者缺电的电量δe(t)，

若δe(t)为正，则此时储能系统吸收此过剩的电量，求解第i台储能元件的充电电量δeb(i,t)，具体过程如下：

10-2-1)根据步骤9)，判断第i台储能元件在t时刻的状态，若为故障状态，则i＝i+1，重新执行步骤10-2-1)，若为正常状态，则执行步骤10-2-2)；

10-2-2)对于第i台储能元件，该设备最大吸收的电量δeb0(i,t)如下：

10-2-3)考虑每小时的最大充放电功率pm，那么每台储能元件每小时的实际充电电量如下：

10-2-4)对于第i+1个储能元件，因为前一个储能元件吸收了电量，因此修正δe(t)，

δe(t)＝δe(t)-δeb(i,t)

10-2-5)令i递增1，即i＝i+1，重复步骤10-2-2)到10-2-4)，直至δe(t)为零或i>m，即所有储能元件充电结束；

若δe(t)为负，则此时储能系统补充系统此时的缺电电量，求解第i台储能元件的放电电量δeb(i,t)，具体过程如下：

10-2-a)，根据步骤9)，判断第i台储能元件在t时刻的状态，若为故障状态，则i＝i+1，重新执行步骤10-2-a)，若为正常状态，则执行步骤10-2-b)；

10-2-b)对于第i台储能元件，该设备最大放电的电量δeb0(i,t)如下：

10-2-c)考虑每小时的最大充放电功率pm，那么每台储能元件每小时的实际放电电量如下：

10-2-d)对于第i+1个储能元件，因为前一个设备放出了电量，因此修正δe(t)，

δe(t)＝δe(t)-δeb(i,t)

10-2-e)令i递增1，即i＝i+1，重复步骤10-2-b)到10-2-d)，直至δe(t)为零或i>m，即所有储能元件放电结束。

前述的步骤11)中，储能系统实际的充放电量，具体求解过程如下：

11-1)参数初始化，储能系统内共有m台储能元件，每台储能元件的额定容量为ee，t时刻第i台储能元件剩余电量为eb(i,t)，初始时刻第i台储能元件剩余电量eb(i,0)为0，即所有储能元件初始时刻的剩余电量均为0，每小时的最大充放电功率为pm；

11-2)根据步骤7-3)得到系统t时刻的系统过剩电量δe(t)和过剩的风电电量δgw(t)；根据冒泡排序的算法，将t时刻的储能元件按照此刻储能元件剩余电量eb(i,t)的多少从高到低排序，并进行编号，得到储能元件的优先放电顺序；

优先放电排序方法如下：

11-2-1)初始化储能元件的放电顺序数组s，数组内有m个元素，m个元素的初始值s(1)——s(m)依次为1至m；

11-2-2)建立一个新的储能元件剩余电量数组e，将原始剩余电量数组eb中所有值赋值给此e数组；

11-2-3)设置两层循环，第一层循环1≤i≤m-1，第二层循环1≤j≤m-i，i，j初始值均设置为1；

11-2-4)第一层循环，初始值i＝1；

11-2-5)第二层循环，对任意一个i值，其初始值j＝1；

11-2-6)对于第j个储能元件，如果其存储的电量e(j)小于第j+1个储能元件存储的电量e(j+1)，则交换剩余电量数组e中两者的电量，同时交换放电顺序数组s中两者的编号，然后转入步骤11-2-7)；如果e(j)不小于e(j+1)，则直接转入步骤11-2-7)；

11-2-7)令j递增1，即j＝j+1，若满足j≤m-i，则转入步骤11-2-6)，否则转入步骤11-2-8)；

11-2-8)令i递增1，即i＝i+1，若满足i≤m-1，则转入步骤11-2-5)，否则，转入步骤11-2-9)；

11-2-9)循环结束，得到优先放电顺序数组s，s(1)存放了最大剩余电量的设备序号，s(2)存放了剩余电量排第二的设备序号，依此类推，直到第m个元素即s(m)存放了最少剩余电量的设备编号；

11-3)根据确立的优先放电顺序数组s，确立每台储能元件的放电电量δeb(k,t)，具体过程如下：

11-3-1)令i＝1；

11-3-2)令k＝s(i)；

11-3-3)根据步骤9)，判断第k台储能元件在t时刻的状态，若为故障状态，则令i＝i+1，重新执行步骤11-3-2)；若为正常状态，则执行步骤11-3-4)；

11-3-4)对于第k台储能元件，该设备最大放电电量δeb0(k,t)如下：

11-3-5)考虑每小时的最大充放电功率pm，那么每台储能元件每小时的实际放电电量如下：

11-3-6)对于第s(i+1)台储能元件，因为前一个设备放出了电量，因此修正δe(t)，

δe(t)＝δe(t)-δeb(k,t)

11-3-7)令i递增1，即i＝i+1，重复步骤11-3-2)到11-3-6)，直至δe(t)为零或i>m，即所有设备放电结束；

11-4)若i>m，直接转步骤12)；否则剩余未参与放电的储能元件可吸收过剩的风电电量δgw(t)，确立每台储能元件的充电电量δeb(k,t)，具体过程如下：

11-4-1)令k＝s(i)；

11-4-2)根据步骤9)判断第k台储能元件在t时刻的状态，若为故障状态，则令i＝i+1，重新执行步骤11-4-1)，若为正常状态，则执行步骤11-4-3)；

11-4-3)对于第k台储能元件，该设备最大吸收的电量δeb0(k,t)如下：

11-4-4)考虑每小时的最大充放电功率pm，那么每台设备每小时的实际充电电量如下：

11-4-5)对于第s(i+1)台储能元件，因为前一个设备吸收了电量，因此修正δgw(t)，

δgw(t)＝δgw(t)-δeb(i,t)

11-4-6)令i递增1，即i＝i+1，重复步骤11-4-1)到11-4-5)直至δgw(t)为零或i>m，即所有设备充电结束。

本发明的有益效果是：

本发明考虑风电机组的故障模型、储能元件的故障模型、风电允许接入比例等因素。在风电场提供的电能超过风电接入比例，并且不满足系统负荷情况下，本发明方法可以避免传统策略按照储能元件原始顺序进行充放电所带来的电量损失，相比传统方法能够向系统提供相同或者更多的充电储能元件，因此能够使风电利用率进一步提升，提高了发电系统的可靠性。

附图说明

图1是本发明方法的流程图；

图2是本发明中计算储能元件优先放电顺序的流程图。

具体实施方式

下面对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

如图1所示，本发明的基于冒泡排序的包含风电场的发电系统储能元件运行方法，包括以下步骤：

步骤1，参数初始化，参数包括：时间参数t、仿真次数n、系统年均缺电量eens、系统缺电概率lolp和年均缺电时间lolt，初始化结果为：t为一年8760小时，n为1000次，eens、lolp、lole初始值为0。

步骤2，根据ieee-rbts可靠性测试系统，运用序贯蒙特卡洛方法模拟得到常规机组的“故障—运行”状态序列，结合每台常规机组的额定功率，得到rbts系统一年的时序出力。

具体包含以下步骤：

步骤2-1，在[0,1]之间产生随机数u1、u2。

步骤2-2，根据第i台常规机组的故障率λi和修复率μi，得到一次平均无故障工作时间一次平均修复时间mttf＝-ln(u2)/μi。

步骤2-3，重复步骤2-2，直到该常规机组所有的mttr和mttf之和超过一年8760个小时。

步骤2-4，将8760个小时中所有mttr时间内常规机组状态记为1，表示常规机组正常运行；mttf时间内常规机组状态记为0，表示常规机组故障，得到一年内继续的“故障—运行”状态序列。

步骤2-5，重复步骤2-2到2-4，直到获取rbts系统内所有常规机组的状态序列。

步骤2-6，将常规机组的状态序列乘以常规机组的额定功率，得到一年内每台常规机组每小时的出力。

步骤2-7，将每小时各个常规机组的出力求和，得到整个rbts系统一年的时序出力。

步骤3，根据风电场所在区域的历史风速数据，基于最小信息准则(aic原则)确定风电场的自回归滑动平均模型即arma模型阶数。

步骤4，根据步骤3所得到的风电场的arma模型阶数及风电场所在区域的历史风速数据，在matlab中调用armax函数模拟产生风电场的风速数据。

步骤5，根据风机的运行参数(故障率、修复率)，运用序贯蒙特卡洛方法模拟得到风机的“故障—运行”状态序列，具体过程同步骤2-1到步骤2-5。

步骤6，根据步骤4和步骤5，结合风机的风速—功率特性曲线，得到风电场一年时间内的功率序列。

具体包含以下步骤：

步骤6-1，根据切入风速vci和额定风速vr，计算特性参数a、b、c：

步骤6-2，根据步骤4中模拟出的风速数据，得到正常运行下，每台风机每小时的功率pt，将pt按照时间顺序排列即可得到一年的功率序列p，

其中，vt表示风速，vco表示切出风速，pw为风机的额定功率。

步骤6-3，将步骤5得到的每台风机的“故障—运行”状态序列点乘功率序列p，即可得到考虑风机故障状态的功率序列p'。

步骤6-4，将每小时内各风机的功率序列p'求和即可得到整个风电场每小时的实际功率。

步骤7，根据步骤2得到的常规机组的时序出力gc(t)和步骤6得到的风电场的时序出力gw(t)，系统负荷为pl(t)，考虑风电允许接入比例得到每个时刻下的系统过剩电量δe(t)和过剩风电电量δgw(t)，具体包含以下步骤：

步骤7-1：考虑一种特殊情况，“风电超出允许接入比例而系统吸收允许接入的风电后仍处于缺电状态”这一条件的数学表达式如式下，将此条件记为条件a：

步骤7-2：风电和常规机组出力不满足条件a时，系统过剩电量δe(t)如下，δe(t)正值表示过剩的电量，δe(t)为负表示系统实际缺电的电量：

步骤7-3：当风电和常规机组出力满足条件a时，系统过剩电量δe(t)和过剩风电电量δgw(t)如下：

步骤8，根据第m台储能元件的故障率和修复率，运用序贯蒙特卡洛方法模拟得到储能元件的“故障—运行”状态序列，具体过程同步骤2-2到步骤2-5。

步骤9，利用步骤8的状态序列中状态值是否为0判断当前时刻是否为故障状态。若为0表示该储能元件处于故障状态，则当前时刻的储能元件的电量清零。

步骤10，当t时刻风电和常规机组出力不满足条件a时，那么储能系统根据步骤7-2中系统过剩电量δe(t)，确定t时刻储能系统实际的充放电量。

具体过程如下：

步骤10-1，参数初始化。储能系统内共有m台储能元件，每台储能元件的额定容量为ee，t时刻第i台储能元件剩余电量为eb(i,t)，初始时刻第i台储能元件剩余电量eb(i,0)为0，即所有储能元件初始时刻的剩余电量均为0，每小时的最大充放电功率为pm。

步骤10-2，根据步骤7-2得到系统t时刻的过剩的电量或者缺电的电量δe(t)，若δe(t)为正，则此时储能系统吸收此过剩的电量，以此确立每台机组的充电电量，δeb(i,t)为正，表示第i台储能元件的充电电量。

具体过程如下：

步骤10-2-1，根据步骤9，判断第i台储能元件在t时刻的状态，若为故障状态，则i＝i+1，重新执行步骤10-2-1，若为正常状态，则执行步骤10-2-2。

步骤10-2-2，对于第i台储能元件，该设备最大吸收的电量δeb0(i,t)如下：

步骤10-2-3，考虑每小时的最大充放电功率pm，那么每台储能元件每小时的实际充电电量如下：

步骤10-2-4，对于第i+1个储能元件，因为前一个储能元件吸收了部分电量，因此修正δe(t)，

δe(t)＝δe(t)-δeb(i,t)

步骤10-2-5，令i递增1，即i＝i+1，重复步骤10-2-2到10-2-4，直至δe(t)为零或i>m，即所有储能元件充电结束。

若δe(t)为负，则此时储能系统补充系统此时的缺电电量，因此确立每台机组的放电电量，δeb(i,t)为负，表示该储能元件的放电电量，具体过程如下：

步骤10-2-a，根据步骤9，判断第i台储能元件在t时刻的状态，若为故障状态，则i＝i+1，重新执行步骤10-2-a，若为正常状态，则执行步骤10-2-b。

步骤10-2-b，对于第i台储能元件，该设备最大放电的电量δeb0(i,t)如下：

步骤10-2-c，考虑每小时的最大充放电功率pm，那么每台储能元件每小时的实际放电电量如下：

步骤10-2-d，对于第i+1个储能元件，因为前一个设备放出了部分电量，因此修正δe(t)，

δe(t)＝δe(t)-δeb(i,t)

步骤10-2-e，令i递增1，即i＝i+1，重复步骤10-2-b到10-2-d，直至δe(t)为零或i>m，即所有储能元件放电结束。

步骤11：当t时刻风电和常规机组出力满足条件a时，那么储能系统根据步骤7-3中系统缺电量δe(t)和过剩的风电电量δgw(t)，确定t时刻储能系统实际的充放电量。

具体过程如下：

步骤11-1，参数初始化。储能系统内共有m台储能元件，每台储能元件的额定容量为ee，t时刻第i台储能元件剩余电量为eb(i,t)，初始时刻第i台储能元件剩余电量eb(i,0)为0，即所有储能元件初始时刻的剩余电量均为0，每小时的最大充放电功率为pm。

步骤11-2，根据步骤7-3得到系统t时刻的缺电的电量δe(t)和过剩的风电电量δgw(t)。

步骤11-3，根据冒泡排序的算法，将t时刻的储能元件按照此刻储能元件剩余电量eb(i,t)的多少从高到低排序，并进行编号，得到储能元件的优先放电顺序。

排序方法参见图2，具体包含以下步骤：

步骤11-3-1，初始化储能元件的放电顺序数组s，数组内有m个元素，m个元素的初始值s(1)——s(m)依次为1至m。

步骤11-3-2，建立一个新的储能元件剩余电量数组e，将原始剩余电量数组eb中所有值赋值给此e数组，此步骤目的在于保证以下各步操作对原始剩余电量数组eb不造成改变。

步骤11-3-3，冒泡排序算法中有两层循环，第一层循环1≤i≤m-1，第二层循环1≤j≤m-i。i，j初始值均设置为1。第一层循环的意义在于进行m-1次第二层循环，可选出m-1个最大值。第二层循环的意义是在m-i个设备中选择出最大值，保证前一次得出的最大值不包含在当前的比较值当中。m-1次选择最大值后，还剩余1个元素值可不必比较。

步骤11-3-4，第一层循环，初始值i＝1。

步骤11-3-5，第二层循环，对任意一个i值，其初始值j＝1。

步骤11-3-6，对于第j个储能元件，如果其存储的电量e(j)小于第j+1个储能元件存储的电量e(j+1)，则交换剩余电量数组e中两者的电量，同时交换放电顺序数组s中两者的编号，然后转入步骤11-3-7；如果e(j)不小于e(j+1)，则直接转入步骤11-3-7。

步骤11-3-7，令j递增1，即j＝j+1。若满足j≤m-i，则转入步骤11-3-6，否则转入步骤11-3-8。

步骤11-3-8，令i递增1，即i＝i+1。若满足i≤m-1，则转入步骤11-3-5，否则，转入步骤11-3-9。

步骤11-3-9，循环结束，得到优先放电顺序数组s。s(1)存放了最大剩余电量的设备序号，s(2)存放了剩余电量排第二的设备序号，依此类推，直到第m个元素即s(m)存放了最少剩余电量的设备编号。

步骤11-4，根据确立的优先放电顺序数组s，确立每台储能元件的放电电量δeb(k,t)，以补充系统缺电量，δeb(k,t)为负值表示放电电量。

具体过程如下：

步骤11-4-1，i＝1。

步骤11-4-2，k＝s(i)。

步骤11-4-3，根据步骤9，判断第k台储能元件在t时刻的状态，若为故障状态，则令i＝i+1，重新执行步骤11-4-2。若为正常状态，则执行步骤11-4-4。

步骤11-4-4，对于第k台储能元件，该设备最大放电电量δeb0(k,t)如下：

步骤11-4-5，考虑每小时的最大充放电功率pm，那么每台储能元件每小时的实际放电电量如下：

步骤11-4-6，对于第s(i+1)台储能元件，因为前一个设备放出了部分电量，因此修正δe(t)，

δe(t)＝δe(t)-δeb(k,t)

步骤11-4-7，令i递增1，即i＝i+1，重复步骤11-4-2到11-4-6，直至δe(t)为零或i>m，即所有设备放电结束。

步骤11-5，若i>m，直接转入步骤12。否则剩余未参与放电的储能元件可吸收过剩的风电电量δgw(t)，确立每台储能元件的充电电量δeb(k,t)，δeb(k,t)为正值，表示充电电量。

具体过程如下：

步骤11-5-1，k＝s(i)。

步骤11-5-2，根据步骤9，判断第k台储能元件在t时刻的状态，若为故障状态，则令i＝i+1，重新执行步骤11-5-1，若为正常状态，则执行步骤11-5-3。

步骤11-5-3，对于第k台储能元件，该设备最大吸收的电量δeb0(k,t)如下：

步骤11-5-4，考虑每小时的最大充放电功率pm，那么每台设备每小时的实际充电电量如下：

步骤11-5-5，对于第s(i+1)台储能元件，因为前一个设备吸收了部分电量，因此修正δgw(t)，

δgw(t)＝δgw(t)-δeb(i,t)

步骤11-5-6，令i递增1，即i＝i+1，重复步骤11-5-1到11-5-5直至δgw(t)为零或i>m，即所有设备充电结束。

步骤12，将步骤10或步骤11确立的储能元件充电电量或放电电量δeb(i,t)和t时刻剩余电量eb(i,t)相加即可得到下一时刻各储能元件的剩余电量δeb(i,t+1)。

步骤13，根据步骤10或步骤11得到系统t时刻的缺电量。若执行步骤10或者步骤11后，δe(t)<0，则说明此时刻系统处于缺电状态，缺电量记为se(t)，且se(t)＝-δe(t)，同时系统缺电时间lole＝lole+1。

步骤14，令t＝t+1，重复步骤10到步骤13，直到t>8760。

步骤15，计算一年的缺电量eens，其值为一年内所有se(t)的总和。

步骤16，计算eens＝eens+eens(初始eens为0)，重复步骤2到步骤16，每循环一次，仿真次数加1，直到仿真次数超过n。

步骤17，计算系统年均缺电量eens＝eens/n，系统年均缺电小时数lolt＝lole/n，系统缺电概率lolp＝lolt/8760。

本发明首先对常规机组进行建模，得到了常规机组的发电出力数据。其次结合arma模型和风速—功率特性曲线，考虑风机的故障率，模拟了风电场的时序出力数据。然后根据这两种发电数据结合时序负荷数据，获取了系统的过剩电量以及过剩的风电量。然后基于冒泡排序的原理，根据储能元件内部剩余电量，确立优先放电的顺序，对过剩的风电电量进行最大化的存储，进一步提高的风电的利用率，改善了系统的可靠性。

实施例

为了进一步的说明本发明的优势，建立一种实际系统的可能情况：设系统此时负荷为200mw，常规机组提供170mw，风电为30mw，风电允许接入比例为10％，则过剩风电量为10mw，负荷此时缺电量为10mw。系统内含有5台电池，每台电池的容量为5mwh。考虑两种储能系统运行过程：1)按照5台电池的物理序号，确定每台电池的充放电量。2)按照本发明方法，即5台电池剩余电量由高到低排序，确定每台电池的充放电量。电池的原有电量以及在两种策略作用下，各电池对缺电系统的放电量和吸收过剩风电的电量如表1、表2所示。

表1物理序号下各储能元件电量

表2冒泡排序序号下各储能元件电量

由表1、表2可看出，按照冒泡排序所确立的各储能元件充放电量比按照物理序号确立的充放电量多吸收过剩的风电4mw，能进一步提高风电的利用率。同时，排序后的储能系统剩余电量也高于物理序号下系统的剩余电量，相当于提高了储能系统下一次放电的可用容量，能有效提高系统可靠性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。