一种确定电力系统状态转移矩阵稳态特征值的方法及系统与流程

1.本发明涉及电力系统技术领域，并且更具体地，涉及一种确定电力系统状态转移矩阵稳态特征值的方法及系统。


背景技术：

2.电力系统针对次同步振荡和低频振荡等不同的电力系统场景分析，需要计算电力系统状态空间系统矩阵的特征值，而现代电力系统由于其接入规模变大、接入的分布式电源种类繁多和各种电力电子器件，使得电力系统的状态矩阵高达上万阶，其矩阵大小超出qr法计算全部特征值的求解能力。
3.实际分析中，电力系统的小干扰稳定只需要求解不稳定状态特征值、阻尼比最小特征值或者实部最大特征值即可，这些特征值亦可以成为关键特征值。因此，高阶电力系统的关键特征值的快速求取方法是进行小干扰稳定分析的可行和有效方法。对于特征值的振荡模态的分析方法，现有的研究工作主要集中在如何高效以及可靠地计算高阶电力系统的关键特征值。在数值计算中高阶电力系统关键特征值的计算方法是部分特征值算法，而这种针对大规模矩阵特征值问题的部分特征值算法是数值计算领域的研究重点之一。


技术实现要素：

4.针对上述问题，本发明提出了一种确定电力系统状态转移矩阵稳态特征值的方法，包括：
5.获取目标电力系统状态空间方程的状态转移矩阵；
6.针对目标电力系统稳态，确定cayley变换参数，根据cayley变换参数对状态转移矩阵进行cayley矩阵变换，获取矩阵t；
7.根据矩阵t，获取初始krylov分解；
8.将初始krylov分解变换为krylov
‑
schur分解，并根据krylov
‑
schur分解计算ritz值残差；
9.根据ritz值残差，对ritz值分类，获取收敛的ritz值，根据收敛的ritz值获取重构矩阵；
10.根据重构矩阵确定状态转移矩阵的稳态特征值。
11.可选的，获取初始krylov分解，具体为：计算矩阵t的任意阶的arnoldi分解，确定arnoldi分解过程，根据arnoldi分解过程建立初始krylov子空间，根据初始krylov子空间获取初始krylov分解。
12.可选的，方法还包括：计算矩阵t上任意阶的hessenberg矩阵的schur分解，根据hessenberg矩阵的schur分解，将初始krylov分解变换为krylov
‑
schur分解。
13.可选的，对ritz值分类，包括：将ritz值分为已经收敛的ritz值集合、截取阶段保留的未收敛的ritz值集合和截取阶段丢失的ritz值集合。
14.可选的，方法还包括：若ritz值的模值高于1且为收敛，对ritz值迭代截取直至
ritz值收敛。
15.本发明还提出了一种确定电力系统状态转移矩阵稳态特征值的系统，包括：
16.初始单元，获取目标电力系统状态空间方程的状态转移矩阵；
17.矩阵变换单元，针对目标电力系统稳态，确定cayley变换参数，根据cayley变换参数对状态转移矩阵进行cayley矩阵变换，获取矩阵t；
18.分解单元，根据矩阵t，获取初始krylov分解；
19.第一计算单元，将初始krylov分解变换为krylov
‑
schur分解，并根据krylov
‑
schur分解计算ritz值残差；
20.第二计算单元，根据ritz值残差，对ritz值分类，获取收敛的ritz值，根据收敛的ritz值获取重构矩阵；
21.输出单元，根据重构矩阵确定状态转移矩阵的稳态特征值。
22.可选的，分解单元获取初始krylov分解，具体为：计算矩阵t的任意阶的arnoldi分解，确定arnoldi分解过程，根据arnoldi分解过程建立初始krylov子空间，根据初始krylov子空间获取初始krylov分解。
23.可选的，第一计算单元还用于：计算矩阵t上任意阶的hessenberg矩阵的schur分解，根据hessenberg矩阵的schur分解，将初始krylov分解变换为krylov
‑
schur分解。
24.可选的，第二计算单元对ritz值分类，包括：将ritz值分为已经收敛的ritz值集合、截取阶段保留的未收敛的ritz值集合和截取阶段丢失的ritz值集合。
25.可选的，第二计算单元还用于：确定当ritz值的模值高于1且为收敛时，对ritz值迭代截取直至ritz值收敛。
26.本发明结合krylov子空间算法，只需一次cayley谱变换，即可以求出高阶电力系统稳定分析的关键特征值；
27.本发明提升了稳态特征值的提升计算速度以及计算精度，能更快速观测电力系统的状态变化；且节省了计算量，提升了计算速度、降低了计算空间需求。
附图说明
28.图1为本发明方法的流程图；
29.图2为本发明方法实施例的流程图；
30.图3为本发明方法实施例arnoldi分解过程；
31.图4为本发明方法实施例的全维特征值对比图；
32.图5为本发明系统的结构图。
具体实施方式
33.现在参考附图介绍本发明的示例性实施方式，然而，本发明可以用许多不同的形式来实施，并且不局限于此处描述的实施例，提供这些实施例是为了详尽地且完全地公开本发明，并且向所属技术领域的技术人员充分传达本发明的范围。对于表示在附图中的示例性实施方式中的术语并不是对本发明的限定。在附图中，相同的单元/元件使用相同的附图标记。
34.除非另有说明，此处使用的术语(包括科技术语)对所属技术领域的技术人员具有
通常的理解含义。另外，可以理解的是，以通常使用的词典限定的术语，应当被理解为与其相关领域的语境具有一致的含义，而不应该被理解为理想化的或过于正式的意义。
35.本发明提出了一种确定电力系统状态转移矩阵稳态特征值的方法，如图1所示，包括：
36.针对目标电力系统，获取目标电力系统状态空间方程的状态转移矩阵；
37.针对目标电力系统稳态，确定cayley变换参数，根据cayley变换参数对状态转移矩阵进行cayley矩阵变换，获取矩阵t；
38.根据矩阵t，获取初始krylov分解；
39.将初始krylov分解变换为krylov
‑
schur分解，并根据krylov
‑
schur分解计算ritz值残差；
40.根据ritz值残差，对ritz值分类，获取收敛的ritz值，根据收敛的ritz值获取重构矩阵；
41.根据重构矩阵确定状态转移矩阵的稳态特征值。
42.其中，获取初始krylov分解，具体为：计算矩阵t的任意阶的arnoldi分解，确定arnoldi分解过程，根据arnoldi分解过程建立初始krylov子空间，根据初始krylov子空间获取初始krylov分解。
43.其中，计算矩阵t上任意阶的hessenberg矩阵的schur分解，根据hessenberg矩阵的schur分解，将初始krylov分解变换为krylov
‑
schur分解。
44.其中，对ritz值分类，包括：将ritz值分为已经收敛的ritz值集合、截取阶段保留的未收敛的ritz值集合和截取阶段丢失的ritz值集合。
45.其中，若ritz值的模值高于1且为收敛，对ritz值迭代截取直至ritz值收敛。
46.下面结合具体实施例对本发明进行进一步的说明：
47.实施例步骤如图2所示和图3所示，包括：
48.步骤1：获取电力系统状态空间方程的状态转移矩阵。
49.步骤2：针对电力系统所需分析稳定状态确定cayley变换所需参数。位移变换点需在所求关键特征值附近，以增大其主导性。
50.对于复常数α1和α2(α1≠α2)，cayley变换为：
51.t＝(a
‑
a1i)(a
‑
a2i)
‑1ꢀꢀꢀ
(1)
52.式(1)中a为高阶电力系统的状态矩阵，t为经过cayley变换的系统矩阵，i为单位矩阵，经过上式的谱变换，不难推到出其特征值计算为：
[0053][0054]
式(2)中λ
a
和λ
t
分别为高阶系统状态矩阵的特征值以及谱变换后的特征值。
[0055]
cayley变换的唯一变换点α1和α2连接线段的垂直平分线将复平面划分为左右两侧，将α1侧的特征值映射到单位圆以外，将α2侧的特征值映射到单位圆以内，以便求取其主导特征值，即原系统状态矩阵的关键特征值。针对某一特定范围则可以选取某个阻尼比ζ1以及α2＝a b
j
，α2阻尼比小于ζ1，则有：
[0056][0057]
根据式(3)的变换则可以将特定范围内的特征值映射为矩阵t的模值大于1的所有特征值。
[0058]
步骤3：计算矩阵t的m阶arnoldi分解，由arnoldi过程建立初始k维krylov子空间(k＝2
×
m)，得到初始krylov分解。若定义矩阵t关于向量u的m阶krylov子空间为：
[0059]
k
m
(t,u)＝span{u,tu,...,t
m
‑1u}t∈c
n
×
n
ꢀꢀꢀ
(4)
[0060]
则t关于向量u的m阶arnoldi分解形式如下：
[0061][0062]
式(5)中：u∈c
n
×
m
，h为k阶上hessenberg矩阵；ek为第k个元素为1的单位向量。
[0063]
arnoldi分解过程，如图3所示，具体如下：
[0064]
初始化随机向量u1，使得||u1||2＝1；k＝1；
[0065]
计算向量w＝tu
k
；
[0066]
对w向量以及各个特征向量v1,v2....v
k
正交化处理；计算w与v
i
的点积，得到h
k
的第i行第j列元素为再将v
i
的分量从w中移除，更新w为：w＝w
‑
h
ij
v
i
；
[0067]
计算h
k
的第j 1行第j列元素为：h
k
(j 1,j)＝||w||2；
[0068]
将w标准化得到u
k
的第j 1列向量：u
j 1
＝w/h
k
(j 1,j)；
[0069]
令k＝k 1，循环进行2)到6)直至k＝k0，则跳出arnoldi分解过程；
[0070]
步骤4：计算上hessenberg矩阵h
k
的schur分解，将krylov分解变换为krylov
‑
schur分解；计算酉矩阵q
k
将h
k
转化为schur形，得到：
[0071]
tu
′
k
＝u
′
k
s
k
u
k 1
b
′
k
ꢀꢀꢀ
(6)
[0072]
式中：u
′
k
＝u
k
q
k
；
[0073]
计算所有ritz值的残差，将公式(6)中的schur矩阵排序后，则其sk矩阵为上三角矩阵，其对角元素既是其特征值。对于s矩阵的第i个对角元素其对应的特征值与特征向量为(λ
i
,v
i
)，其ritz值残差定义为：
[0074]
res
j
＝||b(j)u||2ꢀꢀꢀ
(7)
[0075]
式(7)中b(j)为式(6)中b
′
k
的第j个元素。若残差res均小于设定收敛误差允许值，则判定ritz收敛。可以通过将b(j)置为零，实现krylov
‑
schur分解的收缩，当所有代求的特征值收敛后，则krylov
‑
schur迭代结束。
[0076]
步骤5：将得到的ritz值分为三个部分。已经收敛的ritz值集合(ω
r1
)、截取阶段保留的未收敛的ritz值集合(ω
r2
)、截取阶段丢失的ritz值集合(ω
r3
)。若所有模值大于1的ritz值都收敛时进入步骤7；否则继续进行步骤6；
[0077]
步骤6：对s
k
矩阵按对角元素模值大小进行排序，由子空间特征值算法的收敛判据可知，每个特征值的收敛速度由其主导性的强弱而定，而主导性则由各特征值的模值相对其他特征值的模值大小而定，所以特征值的收敛次序不完全按照模值由大到小收敛，而是与其相对模值大小有关。对步骤5中的集合进行如下排序：
[0078][0079]
随着迭代次数的增多，收敛的ritz值逐渐增多，在重启动过程中应适当丢弃尚未收敛的最小模值的ritz值，其值的集合为ω
r3
。
[0080]
进行截取操作，保留krylov
‑
schur分解得到的ritz集合ω
r1
∪ω
r2
得到：
[0081]
tu
′
m
＝u
′
ms
′
m
u
k 1
b
′
m 1
ꢀꢀꢀ
(9)
[0082]
若k<2
×
m，则重新设置子空间大小；再利用arnoldi过程得到k维的krylov分解；锁定收敛的ritz值，继续计算式(9)，返回步骤3。
[0083]
步骤7：所有模值大于1的ritz值都收敛后，按照公式(2)求取重构矩阵t的特征值，所获取特征值即对应关注模态。
[0084]
下面以具体算例对本发明进行验证，如下：
[0085]
步骤1：获取复杂电力系统状态空间方程的状态转移矩阵。
[0086]
系统状态空间方程的状态转移矩阵为2371阶，获得方式为：在powerfactory仿真软件中根据我国某省220kv及以上电网实际拓扑参数，搭建仿真模型，基于仿真模型导出状态转移矩阵。
[0087]
步骤2：针对电力系统所需分析稳定状态确定cayley变换所需参数。
[0088]
cayley变换的唯一变换点α1和α2连接线段的垂直平分线将复平面划分为左右两侧，将α1侧的特征值映射到单位圆以外，将α2侧的特征值映射到单位圆以内，以便求取其主导特征值，即原系统状态矩阵的关键特征值。针对某一特定范围则可以选取某个阻尼比ζ1以及α2＝a b
j
，α2阻尼比小于ζ1。
[0089]
设定关注模态为低频振荡模态及潜在低频振荡模态，即频率为0.1hz～2.5hz、阻尼比小于0.1的模态。
[0090]
因此，关注特征值为复平面斜率为
±
10所包围的区域外、虚部在0.628j～15.7j范围之间的特征值。位移变换点需在所求关键特征值附近，以增大其主导性，因此本算例确定位移变换点为1 9j，所需输入的参数为α2、设置阻尼比大小、矩阵的输入文件以及子空间初始大小的设置。
[0091]
根据上述参数，将2371阶高维状态转移矩阵变换为cayley变换阵t。
[0092]
步骤3：计算矩阵t的m阶arnoldi分解，由arnoldi过程建立初始k维krylov子空间(k＝2
×
m)，得到初始krylov分解。
[0093]
分配子空间大小k为100，因此m＝50。
[0094]
步骤4：将krylov分解变换为krylov
‑
schur分解，并计算ritz值残差。
[0095]
步骤5：将得到的ritz值分为三个部分。已经收敛的ritz值集合(ω
r1
)、截取阶段保留的未收敛的ritz值集合(ω
r2
)、截取阶段丢失的ritz值集合(ω
r3
)，继续进行步骤6直至所有模值大于1的ritz值都收敛。
[0096]
步骤6：对步骤5中的集合进行排序，适当丢弃尚未收敛的最小模值的ritz值，保留krylov
‑
schur分解得到的ritz集合ω
r1
∪ω
r2
得到锁定收敛的ritz值，获取重构矩阵。
[0097]
步骤7：求取重构矩阵特征值。
[0098]
通过上述步骤，本发明的关键特征值，以及与对原2371阶矩阵通过qr算法求取的
全维特征值对比如图4所示，显然，本发明方法能够对高维复杂电力系统状态转移矩阵求取所有关注模态对应的特征值。
[0099]
从计算效率而言，计算该2371阶系统的特征根，通过传统qr法平均需要12.1176秒；dgeev算法平均需要2.3292秒；而本发明方法仅需0.9286秒，本发明方法相比于传统qr算法计算时间缩短92.33％，计算效率提升13.05倍；相比于dgeev算法，计算时间缩短68.09％，计算效率提升1.51倍。
[0100]
本发明还提出了一种确定电力系统状态转移矩阵稳态特征值的系统200，如图5所示，包括：
[0101]
初始单元201，针对目标电力系统，获取目标电力系统状态空间方程的状态转移矩阵；
[0102]
矩阵变换单元202，针对目标电力系统稳态，确定cayley变换参数，根据cayley变换参数对状态转移矩阵进行cayley矩阵变换，获取矩阵t；
[0103]
分解单元203，根据矩阵t，获取初始krylov分解；
[0104]
第一计算单元204，将初始krylov分解变换为krylov
‑
schur分解，并根据krylov
‑
schur分解计算ritz值残差；
[0105]
第二计算单元205，根据ritz值残差，对ritz值分类，获取收敛的ritz值，根据收敛的ritz值获取重构矩阵；
[0106]
输出单元206，根据重构矩阵确定状态转移矩阵的稳态特征值。
[0107]
其中，获取初始krylov分解，具体为：计算矩阵t的任意阶的arnoldi分解，确定arnoldi分解过程，根据arnoldi分解过程建立初始krylov子空间，根据初始krylov子空间获取初始krylov分解。
[0108]
其中，第一计算单元204还用于：计算矩阵t上任意阶的hessenberg矩阵的schur分解，根据hessenberg矩阵的schur分解，将初始krylov分解变换为krylov
‑
schur分解。
[0109]
其中，对ritz值分类，包括：将ritz值分为已经收敛的ritz值集合、截取阶段保留的未收敛的ritz值集合和截取阶段丢失的ritz值集合。
[0110]
其中，第二计算单元205还用于：确定当ritz值的模值高于1且为收敛时，对ritz值迭代截取直至ritz值收敛。
[0111]
本发明结合krylov子空间算法，只需一次cayley谱变换，即可以求出高阶电力系统稳定分析的关键特征值；
[0112]
本发明提升了稳态特征值的提升计算速度以及计算精度，能更快速观测电力系统的状态变化；且节省了计算量，提升了计算速度、降低了计算空间需求。
[0113]
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd
‑
rom、光学存储器等)上实施的计算机程序产品的形式。本发明实施例中的方案可以采用各种计算机语言实现，例如，面向对象的程序设计语言java和直译式脚本语言javascript等。
[0114]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流
程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0115]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0116]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0117]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0118]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。