一种含光伏发电的发电系统可靠性评估方法

1.本发明属于电力系统分析技术领域，具体涉及一种含光伏发电的发电系统可靠性评估方法，尤其涉及一种基于改进拉丁超立方抽样(latin hypercube sampling，lhs)的含光伏发电的发电系统可靠性评估方法。


背景技术：

2.随着传统化石能源的逐渐枯竭和大气污染、温室效应等环境问题的日益突出，光伏发电在世界范围内得到了快速发展。然而，光伏出力具有强烈的随机性和不确定性，未来高比例光伏发电的并网增加了电力系统运行分析的复杂性，因此有必要对包含光伏发电的电力系统可靠性评估问题开展研究。
3.发电系统可靠性评估方法主要分为解析法和蒙特卡洛方法两大类，其中解析法的计算量随系统规模的增大呈指数增长，同时存在着新能源发电系统解析模型建模过程复杂的问题，因此基于蒙特卡洛方法的可靠性评估方法在含新能源的电力系统可靠性评估中得到了更广泛的应用。针对传统蒙特卡洛方法存在的抽样样本大、收敛速度慢的问题，研究人员提出了拉丁超立方抽样(latin hypercube sampling，lhs)、交叉熵抽样等方差减小技术。lhs方法通过分层抽样在不改变输出随机变量样本分布的情况下，使得产生的样本分布更加均匀，因而能够在抽取少量样本的情况下提高抽样效率；而交叉熵抽样(cross entropy，ce)方法则在保持原样本数学期望不变的条件下，以交叉熵最小为目标寻找最优重要抽样概率密度分布，即通过改变样本空间的概率分布使得抽取的样本以更大的概率落在关注的参数取值区域，从而缩减了样本数量。由此可见，lhs方法和交叉熵抽样方法从不同的角度提高了系统状态的抽样效率，将交叉熵理论与lhs抽样相结合，能够进一步提高电力系统可靠性评估的抽样效率。


技术实现要素：

4.针对上述问题，本发明提出一种含光伏发电的发电系统可靠性评估方法，在实现光伏出力随机生产模拟和建立系统等值负荷模型的基础上，采用交叉熵方法求解发电机组的近似最优概率密度分布，然后基于所求得的近似最优概率分布函数，采用lhs方法对发电机组状态进行分层抽样，从而加速了含光伏发电的电力系统可靠性评估的计算速度。
5.为了实现上述技术目的，达到上述技术效果，本发明通过以下技术方案实现：
6.一种含光伏发电的发电系统可靠性评估方法，包括：
7.步骤1：参数初始化；
8.步骤2：基于光伏发电统计数据，模拟产生光伏发电出力序列；
9.步骤3：将原始负荷序列与所述光伏发电出力序列进行点对点修正，获取修正后的负荷序列；
10.步骤4：对修正后的负荷序列进行线性化分或者聚类，建立等值负荷的多状态模型；
11.步骤5：使用交叉熵重要抽样法计算发电系统中发电机的最优不可用度；
12.步骤6：根据所述发电机的最优不可用度，利用lhs方法生成发电系统发电机状态序列和对应的似然比序列；
13.步骤7：基于所述发电系统发电机状态序列、对应的似然比序列和等值负荷的多状态模型，计算出含光伏发电的发电系统可靠性指标。
14.可选地，所述光伏发电出力序列的产生方法包括以下步骤：
15.步骤2
‑
1，读取光伏发电统计数据中的光伏出力历史小时序列数据，记为样本数组a，其样本数为n
pv
，舍去样本数据中的0，剩余的数据构成样本数组x
pv
，其样本数为n；
16.步骤2
‑
2，假设数组x
pv
中的n个样本变量为{x
pv,1
,x
pv,2
,...,x
pv,n
}，样本变量x
pv
最小值为a，最大值为b，概率密度函数为f(x
pv
)，则光伏出力的非参数核密度估计为其中，j1是样本数n的索引，带宽σ为标准差，k(
·
)表示核函数，且
17.步骤2
‑
3，将变量x
pv
的定义域[a,b]等分为任意的m个点，用步骤2
‑
2中的非参数核密度估计法对这m个点分别求出其对应的概率密度估计值并计算得各点的概率估计值为：
[0018][0019]
其中i
pv
＝1,2,...,m；
[0020]
步骤2
‑
4，以步骤2
‑
3中各点的概率估计值计算处的累积概率估计值：
[0021][0022]
其中，j
pv
是i
pv
的索引，的最小值记为c，的最大值记为d；
[0023]
步骤2
‑
5，以步骤2
‑
4中各个估计值为间隔点，将累积概率区间分为m
‑
1份，在每一个小区间内，利用三次样条插值法得到该区间上样本值x
pv
关于累积概率估计值y
pv
的三次多项式：其中i
m
＝1,2,
…
,m
‑
1；
[0024]
步骤2
‑
6，从小到大依次取其中i1＝1,2,...,n，若则将代入步骤2
‑
5中对应区间的多项式求得所需样本值否则
[0025]
步骤2
‑
7，将步骤2
‑
1中舍去的(n
pv
‑
n)个0数据与步骤2
‑
6中求得的n个样本值排列成新的数组z，z的样本数是n
pv
，其数据为从小到大排列的数据，数组z即为模拟的光伏出力数据；
[0026]
步骤2
‑
8，对步骤2
‑
7中模拟的光伏出力序列z按原始数据的位置特性进行排序，得
到光伏电站时序出力序列z’。
[0027]
可选地，所述步骤2
‑
8中排序的具体过程为：
[0028]
步骤2
‑8‑
1，对步骤2
‑
1中的样本数组a标记位置，对应的位置序列矩阵记为b＝[1,2,
…
,n
pv
]
t
；
[0029]
步骤2
‑8‑
2，将样本数组a中的数据从小到大排列得到新的数组a’，对应的位置矩阵变为b’＝[i2,
…
,n
pv
,
…
,j2]
t
，其中，i2∈{1,2,...,n
pv
}，j2∈{1,2,...,n
pv
}；
[0030]
步骤2
‑8‑
3，新建一个n
pv
×
2的矩阵c，将步骤2
‑
7中模拟的光伏出力序列z作为矩阵c的第一列，位置矩阵b’作为矩阵c的第二列，得到新的矩阵c；
[0031]
步骤2
‑8‑
4，对步骤2
‑8‑
3中的矩阵序列b’按从小到大的顺序排列，则矩阵c中对应的z变为z’，z’即为符合原始数据特性的光伏电站时序出力数据。
[0032]
可选地，所述修正后的负荷序列的具体获取方法包括：
[0033]
将原始负荷序列逐小时减去光伏发电对应时刻出力序列，得到修正后的等值负荷序列。
[0034]
可选地，所述多状态等值负荷模型包括：负荷水平序列和负荷水平概率序列其中，wp
iw
为负荷水平wl
iw
对应的负荷概率，iw＝1,2,...,n
w
；n
w
为多级负荷水平总数，
[0035]
可选地，所述发电系统中发电机的最优不可用度的计算方法包括以下步骤：
[0036]
步骤5
‑
1，交叉熵参数初始化：交叉熵重要抽样样本规模n；分位数ρ；平滑因子p；最大迭代次数n
kmax
；给定m台发电机组出力和该机组对应的不可用度向量u、初次迭代不可用度向量v0＝u，其中v0为1
×
m阶向量；迭代计算次数k＝1；
[0037]
步骤5
‑
2，生成n个系统状态序列，具体方法为：产生n个[0,1]区间的随机数向量，并与k
‑
1次迭代计算使用的系统发电机最优不可用度序列v
k
‑1进行比较，若随机数则令否则令其中表示第n
s
个随机数向量中第个元素，表示第n
s
个系统发电机状态向量中的第个元素；
[0038]
步骤5
‑
3，根据步骤5
‑
2得到的系统发电机状态序列，计算n个系统发电机状态下的系统发电机出力总量，将n个系统发电机出力总量按从小到大顺序进行排序，得到升序排列的系统发电机供电总量序列；
[0039]
步骤5
‑
4，根据步骤5
‑
3得到的升序排列的系统发电机出力总量序列，更新可安全供电总负荷l
k
的值，确定系统在n个状态下的切负荷标志序列；
[0040]
步骤5
‑
5，根据步骤5
‑
2得到的系统发电机状态序列和此次迭代计算使用的发电机最优不可用度序列v
k
‑1，计算n个系统状态似然比序列；
[0041]
步骤5
‑
6，根据步骤5
‑
5得到的系统发电机状态似然比序列，计算此次(即第k次)系统发电机最优不可用度序列v
k
；
[0042]
步骤5
‑
7，判断可安全供电总负荷l
k
是否达到指定的平均负荷p
l
，亦或迭代次数k是否达到最大值n
kmax
，若满足条件，即l
k
＝p
l
或k＝n
kmax
时，转入步骤6；否则令k＝k+1，并返回步骤5
‑
2。
[0043]
可选地，所述步骤5
‑
4中计算系统切负荷标志序列的具体步骤包括：
[0044]
步骤5
‑4‑
1，记步骤5
‑
3得到的升序排列的系统发电机供电总量序列为p；若p
t
>p
l
，则令l
k
＝p
t
；否则，令l
k
＝p
l
；从而计算得到第k次迭代的安全供电总负荷l
k
，其中ρ为分位数，p
t
为系统发电机供电总量序列p中的第t个元素，其中t＝ρ
×
n；p
l
为给定平均负荷，计算公式为
[0045]
步骤5
‑4‑
2，记步骤5
‑
3得到的未排序前的系统发电机供电总量序列为s
l
；记第n
s
个系统发电机状态序列为系统在n个系统发电机状态下的切负荷标志序列为h；若则否则其中表示系统在发电机状态下的供电总量，表示状态下系统对应的切负荷标志值。
[0046]
可选地，所述步骤5
‑
5中所述似然比的计算公式为：
[0047][0048]
其中，w表示系统发电机状态似然比序列，为发电系统状态的似然比，表示第个系统发电机状态序列，表示第个系统发电机状态序列中第台发电机状态值，m为系统发电机数，表示系统第台发电机的初始不可用度值，为系统第个发电机在第k
‑
1次迭代计算使用的不可用度值；其中
[0049]
可选地，步骤5
‑
6中所述v
k
的计算公式为：
[0050][0051]
引入平滑因子p，按式v
′
kj
′
＝pv
kj
′
+(1
‑
p)v
(k
‑
1)j
′
计算，得到v
k
，其中v
k,j
′
表示第k次迭代计算得到的发电机不可用度序列v
k
中第j
′
台发电机的不可用度值，为系统状态x
i
′
的似然比，x
i
′
表示第i
′
个系统发电机状态序列，x
i
′
j
′
表示第i
′
个系统发电机状态序列第j
′
台发电机状态值，h(x
i
′
)表示系统状态x
i
′
的切负荷标志值，其中i
′
＝1,2,...,n，j
′
＝1,2,...,m。
[0052]
可选地，所述发电系统发电机状态序列和对应的似然比序列的生产方法包括以下步骤：
[0053]
步骤6
‑
1，lhs参数初始化：拉丁超立方法采样规模n2；系统发电机台数m；
[0054]
步骤6
‑
2，根据步骤5得到的系统发电机最优不可用度序列v
*
，计算系统发电机状态原始采样矩阵r；
[0055]
步骤6
‑
3，逐行随机生成m
×
n2阶系统发电机状态排列矩阵，它的每一行由[1,n2]范围内的n2个不重复的正整数构成，整数的值代表该系统发电机状态采样矩阵对应行对应元素的排列位置；
[0056]
步骤6
‑
4，根据步骤6
‑
2得到的原始采样矩阵r，计算系统发电机状态行相关系数矩阵ρ；
[0057]
步骤6
‑
5，将步骤6
‑
4得到的相关系数矩阵ρ分解，构造行相关性更小的系统发电机状态排列矩阵q
s*
；
[0058]
步骤6
‑
6，将步骤6
‑
2得到的系统发电机状态原始采样矩阵r按照步骤6
‑
5得到的排列矩阵q
s*
重新排列，生成m
×
n2阶系统发电机状态采样矩阵；
[0059]
步骤6
‑
7，根据系统发电机最优不可用度序列v
*
，结合步骤6
‑
6得到的采样矩阵，产生新的m
×
n2阶系统发电机状态矩阵wx；
[0060]
步骤6
‑
8，根据步骤6
‑
7产生的系统发电机状态矩阵wx，计算1
×
n2阶系统发电机状态似然比序列ww。
[0061]
可选地，所述步骤6
‑
2中所述系统发电机状态原始采样矩阵r的计算公式为：
[0062][0063]
其中，n2为拉丁超立方抽样方法的采样规模数，e和λ都为定常数，r
j*z
表示矩阵r的第j
*
行第z列元素，v
*j*
为发电机不可用度序列v
*
中第j
*
台发电机的不可用度值，z＝1,2,...,n2，j
*
＝1,2,...,m；
[0064]
可选地，所述步骤6
‑
4中所述系统发电机状态相关系数矩阵ρ的计算公式为：
[0065][0066]
其中，m为系统发电机数；ρ
i
″
j
″
表示相关系数矩阵ρ中第i
″
行与第j
″
行的相关系数，r
i
″
q
表示矩阵r的第i
″
行第q列元素，r
j
″
q
表示矩阵r的第j
″
行第q列元素，表示矩阵r的第i
″
行平均值，表示矩阵r的第j
″
行平均值。
[0067]
可选地，所述步骤6
‑
5中相关系数矩阵ρ的分解方法为：
[0068]
采用cholesky分解方法，将相关系数矩阵ρ分解为dd
t
的矩阵乘积形式，其中d为非奇异下三角矩阵；再根据步骤6
‑
3生成的排列矩阵q
s
，生成一个行相关性更小的矩阵g，计算公式为再根据矩阵g每行元素值的大小，将元素用对应整数编号替换，生成系统发电机状态排列矩阵q
s*
。
[0069]
可选地，所述含光伏发电的发电系统可靠性指标的计算方法包括：
[0070]
步骤7
‑
1，令系统发电机状态序列索引wc＝0；
[0071]
步骤7
‑
2，令负荷水平索引wk＝0；
[0072]
步骤7
‑
3，wk＝wk+1；wc＝wc+1；
[0073]
步骤7
‑
4，从步骤6中生成的系统发电机状态矩阵wx中逐次选取系统发电机状态序列wx
wc
，逐次选取负荷水平wl
wk
；wx
wc
为发电机状态矩阵wx中第wc个序列，wl
wk
为负荷水平序列wl中第wk个元素，其对应的负荷概率为wp
wk
；
[0074]
步骤7
‑
5，更新第wk个负荷水平下系统可靠性指标的指示函数；
[0075]
步骤7
‑
6，判断wk<n
w
是否成立，若满足条件，则转入步骤7
‑
7；否则，转入步骤7
‑
3；
[0076]
步骤7
‑
7，更新多级负荷水平下的系统可靠性指标及方差系数；
[0077]
步骤7
‑
8，若步骤7
‑
7所计算的方差系数满足收敛条件，则输出发电系统可靠性指标；否则，转入步骤7
‑
2。
[0078]
可选地，所述步骤7
‑
5中，更新第wk个负荷水平wl
wk
下的系统可靠性指标指示函数的具体方法为：对步骤7
‑
4中的系统状态序列wx
wc
进行系统状态分析，判断系统状态wx
wc
是否产生切负荷；若产生切负荷，则指标指示函数为其中sd为切负荷量，其大小为负荷水平wl
wk
与系统m台发电机在状态wx
wc
下发电总量的差值；若没有产生切负荷，即系统m台发电机在状态wx
wc
下发电总量大于负荷水平wl
wk
，则指标指示函数为f
lolp
(wx
wc
)＝0，f
eens
(wx
wc
)＝0。
[0079]
可选地，所述步骤7
‑
7更新多级负荷水平下的系统可靠性指标及方差系数的计算公式为：
[0080][0081][0082][0083]
与现有技术相比，本发明的有益效果：
[0084]
本发明采用分段三次样条插值函数克服了光伏出力的累积概率分布函数的反函数无法写出具体表达式时不能应用lhs方法的问题，在保证lhs方法抽样精度和效率都足够高的同时，能准确地模拟出光伏发电出力序列。
[0085]
本发明以减小重要抽样概率密度函数与理论上最优重要抽样概率密度函数的交叉熵最小为目标函数，求取近似最优的重要抽样概率分布，在保证可靠性指标期望值近似不变的基础上，提高了稀有事件的出现概率，大大减小了电力系统可靠性评估抽样次数，提高了抽样效率；在此基础上进一步应用lhs方法，以交叉熵寻优得到的最优重要抽样概率分布进行分层抽样，在提高稀有事件概率的基础上保证了样本点区间全覆盖，避免了数据出现截尾现象，从而进一步提高了系统状态的抽样效率。本发明提供的方法充分利用了交叉熵抽样和lhs抽样方法的优点，有效提高了电力系统可靠性评估的计算速度。
附图说明
[0086]
为了使本发明的内容更容易被清楚地理解，下面根据具体实施例并结合附图，对本发明作进一步详细的说明，其中：
[0087]
图1是本发明实施例的发电系统可靠性评估方法的流程示意图；
[0088]
图2是使用交叉熵重要抽样法，寻优求取系统发电机最优不可用度的流程示意图；
[0089]
图3是使用拉丁超立方抽样，求取系统发电机状态序列和似然比序列流程示意图；
[0090]
图4是计算含光伏发电的发电系统可靠性指标流程示意图。
具体实施方式
[0091]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明的保护范围。
[0092]
下面结合附图对本发明的应用原理作详细的描述。
[0093]
本发明中提供了一种基于改进lhs的含光伏发电的发电系统可靠性评估方法，具体包含的步骤有：参数初始化；模拟产生光伏发电出力序列；修正原始负荷序列；建立等值负荷多状态模型；使用交叉熵重要抽样法计算发电系统中发电机的最优不可用度；利用lhs方法生成发电系统发电机状态序列和对应的似然比序列；计算含光伏发电的发电系统可靠性指标。本发明在实现光伏出力随机生产模拟和建立系统等值负荷模型的基础上，通过将交叉熵重要抽样和lhs抽样法相结合，有效地提高了系统随机抽样的效率，加速了含光伏发电的电力系统可靠性评估的计算速度。本发明中所提及的系统均指电力系统。
[0094]
如图1所示，所述基于改进lhs的含光伏发电的发电系统可靠性评估方法，包含步骤如下：
[0095]
步骤1：参数初始化；
[0096]
步骤2：基于光伏发电统计数据，模拟产生光伏发电出力序列；
[0097]
在本发明实施例的一种具体实施方式中，所述步骤2具体包含如下步骤：
[0098]
步骤2
‑
1，读取光伏发电统计数据中的光伏出力历史小时序列数据，记为样本数组a，数组a的长度记为n
pv
；将样本数据中的0舍去，剩余的数据构成数组x
pv
，数组x
pv
的长度为n；
[0099]
步骤2
‑
2，假设数组x
pv
中的n个样本变量为{x
pv,1
,x
pv,2
,...,x
pv,n
}，样本变量x
pv
最小值为a，最大值为b，概率密度函数为f(x
pv
)，则光伏出力的非参数核密度估计为
[0100]
其中，j1是样本数n的索引，带宽σ为标准差，k(
·
)表示核函数，且
[0101]
步骤2
‑
3，将变量x
pv
的定义域[a,b]等分为任意的m个点，用步骤2
‑
2中的非参数核密度估计法对这m个点分别求出其对应的概率密度估计值并计算得各点的概率估计值为；
[0102]
[0103]
其中，i
pv
＝1,2,...,m；
[0104]
步骤2
‑
4，以步骤2
‑
3中各点的概率估计值计算处的累积概率估计值：
[0105][0106]
其中，j
pv
是i
pv
的索引，的最小值记为c，的最大值记为d；
[0107]
步骤2
‑
5，以步骤2
‑
4中各个估计值为间隔点，将累积概率区间分为m
‑
1份，在每一个小区间内，利用三次样条插值法得到该区间上样本值x
pv
关于累积概率估计值y
pv
的三次多项式：其中i
m
＝1,2,
…
,m
‑
1；
[0108]
步骤2
‑
6，从小到大依次取其中i1＝1,2,...,n，若则将代入步骤2
‑
5中对应区间的多项式求得所需样本值否则
[0109]
步骤2
‑
7，将步骤2
‑
1中舍去的(n
pv
‑
n)个0数据与步骤2
‑
6中求得的n个样本值排列成新的数组z，z的样本数是n
pv
，其数据为从小到大排列的数据，数组z即为模拟的光伏出力数据；
[0110]
步骤2
‑
8，对步骤2
‑
7中模拟的光伏出力序列z按原始数据的位置特性进行排序，得到光伏电站时序出力序列z’；
[0111]
所述步骤2
‑
8中排序的具体过程为：
[0112]
步骤2
‑8‑
1，对步骤2
‑
1中的样本数组a标记位置，对应的位置序列矩阵记为b＝[1,2,
…
,n
pv
]
t
；
[0113]
步骤2
‑8‑
2，将样本数组a中的数据从小到大排列得到新的数组a’，对应的位置矩阵变为b’＝[i2,
…
,n
pv
,
…
,j2]
t
，其中，i2∈{1,2,...,n
pv
}，j2∈{1,2,...,n
pv
}；
[0114]
步骤2
‑8‑
3，新建一个n
pv
×
2的矩阵c，将步骤2
‑
7中模拟的光伏出力序列z作为矩阵c的第一列，位置矩阵b’作为矩阵c的第二列，得到新的矩阵c；
[0115]
步骤2
‑8‑
4，对步骤2
‑8‑
3中的矩阵序列b’按从小到大的顺序排列，则矩阵c中对应的z变为z’，z’即为符合原始数据特性的光伏电站时序出力数据。
[0116]
步骤3：将原始负荷序列与步骤2生成的光伏出力序列进行点对点修正，获取修正后的负荷序列；
[0117]
其中，修正方法是：将原始负荷序列逐小时减去光伏发电对应时刻出力序列。
[0118]
步骤4：对修正后的负荷序列进行线性化分或者聚类，建立等值负荷的多状态模型；
[0119]
其中，多状态等值负荷模型中包括负荷水平序列负荷水平概率序列其中，wp
iw
为负荷水平wl
iw
对应的负荷概率，iw＝1,2,...,n
w
；n
w
为多级负荷水平总数，
[0120]
步骤5：使用交叉熵重要抽样法计算系统中发电机的最优不可用度；
[0121]
如图2所示，在本发明实施例的一种具体实施方式中，所述步骤5具体包含如下步骤：
[0122]
步骤5
‑
1，交叉熵参数初始化：交叉熵重要抽样样本规模n；分位数ρ；平滑因子p；最大迭代次数n
kmax
；给定m台发电机组出力和该机组对应的不可用度向量u、初次迭代不可用度向量v0＝u，其中v0为1
×
m阶向量；迭代计算次数k＝1；
[0123]
步骤5
‑
2，生成n个系统状态序列，具体方法为：产生n个[0,1]区间的随机数向量，并与k
‑
1次迭代计算使用的系统发电机最优不可用度序列v
k
‑1进行比较，若随机数则令否则令其中表示第n
s
个随机数向量中第个元素，表示第n
s
个系统发电机状态向量中的第个元素。
[0124]
步骤5
‑
3，根据步骤5
‑
2得到的系统发电机状态序列，计算n个系统发电机状态下的系统发电机出力总量，将n个系统发电机出力总量按从小到大顺序进行排序，得到升序排列的系统发电机供电总量序列；
[0125]
步骤5
‑
4，根据步骤5
‑
3得到的升序排列的系统发电机出力总量序列，更新可安全供电总负荷l
k
的值，确定系统在n个状态下的切负荷标志序列；
[0126]
其中，计算系统切负荷标志序列的具体做法为：
[0127]
步骤5
‑4‑
1，记步骤5
‑
3得到的升序排列的系统发电机供电总量序列为p；若p
t
>p
l
，则令l
k
＝p
t
；否则，令l
k
＝p
l
；从而计算得到第k次迭代的安全供电总负荷l
k
，其中ρ为分位数，p
t
为系统发电机供电总量序列p中的第t个元素，其中t＝ρ
×
n；p
l
为给定平均负荷，计算公式为
[0128]
步骤5
‑4‑
2，记步骤5
‑
3得到的未排序前的系统发电机供电总量序列为s
l
；记第n
s
个系统发电机状态序列为系统在n个系统发电机状态下的切负荷标志序列为h；若则否则其中表示系统在发电机状态下的供电总量，表示状态下系统对应的切负荷标志值。
[0129]
步骤5
‑
5，根据步骤5
‑
2得到的系统发电机状态序列和此次迭代计算使用的发电机最优不可用度序列v
k
‑1，计算n个系统状态似然比序列；
[0130]
其中，似然比的计算公式：
[0131][0132]
其中，w表示系统发电机状态似然比序列，为系统状态的似然比，表示第个系统发电机状态序列，表示第个系统发电机状态序列中第台发电机状态值，m为系统发电机数，表示系统第台发电机的初始不可用度值，为系统第个发电机在第k
‑
1次迭代计算使用的不可用度值；其中
[0133]
步骤5
‑
6，根据步骤5
‑
5得到的系统发电机状态似然比序列，计算此次(即第k次)系
统发电机最优不可用度序列v
k
；
[0134]
其中，v
k
的计算公式：并引入平滑因子p，按式v
′
kj
′
＝pv
kj
′
+(1
‑
p)v
(k
‑
1)j
′
计算，得到v
k
，其中v
k,j
′
表示第k次迭代计算得到的发电机不可用度序列v
k
中第j
′
台发电机的不可用度值，w(x
i
′
；u,v
k
‑1)为系统状态x
i
′
的似然比，x
i
′
表示第i
′
个系统发电机状态序列，x
i
′
j
′
表示第i
′
个系统发电机状态序列第j
′
台发电机状态值，h(x
i
′
)表示系统状态x
i
′
的切负荷标志值，其中i
′
＝1,2,...,n，j
′
＝1,2,...,m。
[0135]
步骤5
‑
7，判断可安全供电总负荷l
k
是否达到指定的平均负荷p
l
，亦或迭代次数k是否达到最大值n
kmax
，若满足条件，即l
k
＝p
l
或k＝n
kmax
时，转入步骤6；否则令k＝k+1，并返回步骤5
‑
2。
[0136]
步骤6：根据所述发电机的最优不可用度，利用lhs方法生成发电系统发电机状态序列和对应的似然比序列；
[0137]
如图3所示，在本发明实施例的一种具体实施方式中，所述步骤6具体包含如下步骤：
[0138]
步骤6
‑
1，参数初始化；拉丁超立方法采样规模n2；系统发电机台数m；
[0139]
步骤6
‑
2，根据步骤5得到的系统发电机最优不可用度序列v
*
，计算系统发电机状态原始采样矩阵r；
[0140][0141]
其中，n2为拉丁超立方抽样方法的采样规模数，e和λ都为定常数，r
j*z
表示矩阵r的第j
*
行第z列元素，v
*j*
为发电机不可用度序列v
*
中第j
*
台发电机的不可用度值，z＝1,2,...,n2，j
*
＝1,2,...,m；
[0142]
步骤6
‑
3，逐行随机生成m
×
n2阶系统发电机状态排列矩阵，它的每一行由[1,n2]范围内的n2个不重复的正整数构成，整数的值代表该系统发电机状态采样矩阵对应行对应元素的排列位置；
[0143]
步骤6
‑
4，根据步骤6
‑
2得到的原始采样矩阵r，计算系统发电机状态行相关系数矩阵ρ；
[0144]
其中，计算相关系数矩阵ρ的公式为：
[0145][0146]
其中m为系统发电机数；ρ
i
″
j
″
表示相关系数矩阵ρ中第i
″
行与第j
″
行的相关系数，r
i
″
q
表示矩阵r的第i
″
行第q列元素，r
j
″
q
表示矩阵r的第j
″
行第q列元素，表示矩阵r的第
i
″
行平均值，表示矩阵r的第j
″
行平均值。
[0147]
步骤6
‑
5，将步骤6
‑
4得到的相关系数矩阵ρ分解，构造行相关性更小的系统发电机状态排列矩阵q
s*
；
[0148]
其中，相关系数矩阵分解方法为cholesky分解方法。这种方法将相关系数矩阵ρ分解为dd
t
的矩阵乘积形式，其中d为非奇异下三角矩阵；再根据步骤6
‑
3生成的排列矩阵q
s
，生成一个行相关性更小的矩阵g，计算公式为再根据矩阵g每行元素值的大小，将元素用对应整数编号替换，生成系统发电机状态排列矩阵q
s*
。
[0149]
步骤6
‑
6，将步骤6
‑
2得到的系统发电机状态原始采样矩阵r按照步骤6
‑
5得到的排列矩阵q
s*
重新排列，生成m
×
n2阶系统发电机状态采样矩阵；
[0150]
步骤6
‑
7，根据系统发电机最优不可用度序列v
*
，结合步骤6
‑
6得到的采样矩阵，产生新的m
×
n2阶系统发电机状态矩阵wx；
[0151]
其中，生成系统发电机状态序列的方法与步骤5
‑
2中所提到的方法相同。
[0152]
步骤6
‑
8，根据步骤6
‑
7产生的系统发电机状态矩阵wx，计算1
×
n2阶系统发电机状态似然比序列ww；
[0153]
其中，计算似然比序列的公式与步骤5
‑
5中的相同，应注意带入的是系统发电机最优不可用度序列v
*
，发电机状态为步骤6
‑
7计算得到的系统发电机状态。
[0154]
步骤7：基于所述发电系统发电机状态序列、对应的似然比序列和等值负荷的多状态模型，计算出含光伏发电的发电系统可靠性指标。
[0155]
如图4所示，在本发明实施例的一种具体实施方式中，所述步骤7具体包括如下步骤：
[0156]
步骤7
‑
1，令系统发电机状态序列索引wc＝0；
[0157]
步骤7
‑
2，令负荷水平索引wk＝0；
[0158]
步骤7
‑
3，wk＝wk+1；wc＝wc+1；
[0159]
步骤7
‑
4，从步骤6中生成的系统发电机状态矩阵wx中逐次选取系统发电机状态序列wx
wc
，逐次选取负荷水平wl
wk
；wx
wc
为发电机状态矩阵wx中第wc个序列，wl
wk
为负荷水平序列wl中第wk个元素，其对应的负荷概率为wp
wk
；
[0160]
步骤7
‑
5，更新第wk个负荷水平下系统可靠性指标的指示函数；
[0161]
其中，更新第wk个负荷水平wl
wk
下的系统可靠性指标指示函数的具体方法：对步骤7
‑
4中的系统状态序列wx
wc
进行系统状态分析，判断系统状态wx
wc
是否产生切负荷；若产生切负荷，则指标指示函数：
[0162]
f
lolp
(wx
wc
)＝ww(wx
wc
)
×
wp
wk
[0163]
f
eens
(wx
wc
)＝8760
×
sd
×
ww(wx
wc
)
×
wp
wk
[0164]
其中，sd为切负荷量，其大小为负荷水平wl
wk
与系统m台发电机在状态wx
wc
下发电总量的差值；若没有产生切负荷，即系统m台发电机在状态wx
wc
下发电总量大于负荷水平wl
wk
，则指标指示函数为f
lolp
(wx
wc
)＝0，f
eens
(wx
wc
)＝0。
[0165]
步骤7
‑
6，判断wk<n
w
是否成立，若满足条件，则转入步骤7
‑
7；否则，转入步骤7
‑
3；
[0166]
步骤7
‑
7，更新多级负荷水平下的系统可靠性指标及方差系数；
[0167]
其中，更新多级负荷水平下的系统可靠性指标及方差系数的计算公式为：
[0168][0169][0170][0171]
步骤7
‑
8，若步骤7
‑
7所计算的方差系数满足收敛条件，则输出发电系统可靠性指标；否则，转入步骤7
‑
2。
[0172]
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。