一种基于CU三角分解求取电力系统节点阻抗矩阵的方法.pdf

一种基于CU三角分解求取电力系统节点阻抗矩阵Z的方法，属于电力系统分析计算领域，主要步骤如下：形成节点导纳矩阵Y；对Y阵进行快速CU三角分解；用CWk＝Ek求取Wk阵；用UZk＝Wk求取Zk阵对角元Zkk及以上的非对角元素；根据对称性求取对角元Zkk以左的非对角元素；写Z阵数据到数据文件。本发明利用了比LDU三角分解法计算效率更高的CU三角分解法来求取Z阵元素；利用Y阵以及C、U阵元素的对称关系快速形成C、U因子阵的合成阵；利用Z阵元素的对称性及Zk阵元素的计算顺序，省去对Z阵下三角元素的回代求取而直接获得，大大提高了形成Z阵速度。用本发明方法对IEEE-30、-57、-118节点系统验算，与传统的LDU三角分解法相比，计算速度可提高约33％。

权利要求书
1.  一种基于的CU三角分解求取电力系统节点阻抗矩阵的方法，其特征包括以下步骤：
步骤1：读取各支路数据文件；
步骤2：形成节点导纳矩阵Y；
步骤3：利用Y阵以及C、U阵元素的对称关系快速形成C、U因子阵的合成阵；
步骤4：根据方程CUZk＝Ek对方程CWk＝Ek求取Wk阵；通过方程UZk＝Wk求取Zk阵对角元Zkk及以上的非对角元素；
步骤5：根据对称性求取对角元Zkk以左的非对角元素；
这种根据Z阵元素的对称性得到对角元Zkk以左元素的计算方式可减少50％非对角元素的计算。也可不求Z阵下三角的非对角元素，以进一步加快计算速度。
步骤6：将Z阵写入数据文件。

说明书一种基于CU三角分解求取电力系统节点阻抗矩阵的方法
技术领域
本发明属于电力系统分析计算领域，涉及一种基于CU三角分解求取电力系统节点阻抗矩阵的方法。
背景技术
节点阻抗矩阵Z在电力系统中应用十分广泛且具有重要的作用。传统的求解Z阵的方法有支路追加法、导纳矩阵Y消元求逆法、LDU三角分解法等。在常用的各种传统方法中，由于LDU三角分解法的计算速度最快，因而使用最多，其特点是利用了适于求解常系数线性方程的三角分解法，在对Y阵进行LDU三角分解后，可将一个对n×n阶Z阵元素的求解分成对n个列矩阵Zk元素的求解。
实际上，在各种三角分解法中，A＝LR或A＝CU三角分解法形成因子阵所需计算变量的总数最少，约为1.5n2或最多2n2，各需求解2个方程组，各有1个中间矩阵变量Wk。而LDU三角分解法形成因子阵所需计算元素的总数为3n2，远远多于LR或CU三角分解法，且需求解3个方程组，有2个中间矩阵变量Wk、Hk。因此LDU三角分解法的计算效率远比LR或CU三角分解法的计算效率要低。一般情况下在求取同类方程时LDU三角分解法的计算速度比LR或CU三角分解法的计算速度约慢30％以上。因此用LDU三角分解法求取Z阵元素并非是对三角分解法应用的最佳选择。而另一方面，CU三角分解法在求取同类方程时其计算速度比LR三角分解法的计算速度又快约3％。因此用三角分解法求取Z阵元素时，CU三角分解法应该是最佳选择。
各种传统三角分解法是独自建立各自的因子阵，因此不容易发现各个因子阵元素之间的相互关系并加以利用。在计算各个因子阵的元素时是按图3“┘”(反L)的方式从左上角到右下角或者按图4“行”的方式从上到下一个一个地计算所有元素，这2种元素的计算方式均无法利用Y阵以及因子阵元素的对称关系。这2种元素的计算方式对所有元素均是利用计算公式一次计算完成，这里称其为公式法。通过分析可以发现，公式法也不利于对元素计算过程的改进。因此各种传统三角分解法形成因子阵的计算效率很低。
此外一般传统的LDU三角分解法在求解Z阵元素的过程中未考虑利用Z阵元素的对称性，因此对非对角元素存在大量的重复计算，因此其计算时间并不理想。
传统的LDU三角分解法Z阵元素的计算顺序为：Z1,…,Zk,…,Zn，过程如图5所示。
发明内容
本发明的目的是克服现有方法的不足之处，提供一种基于CU三角分解求取电力系统节点阻抗矩阵的方法。
本发明是通过以下技术方案实现的。
一种基于CU三角分解法快速求取电力系统节点阻抗矩阵的方法，其特征包括以 下步骤：
步骤1：读取各支路数据文件；
步骤2：形成节点导纳矩阵Y；
步骤3：利用Y阵以及C、U阵元素的对称关系快速形成C、U因子阵的合成阵；
步骤3具体实施过程如下：
根据YZk＝Ek，令Y＝CU，得CUZk＝Ek。再将CUZk＝Ek进一步分解为CWk＝Ek，UZk＝Wk二个方程。
形成C、U二个因子阵的合成阵时，其元素的计算方式对形成因子阵的速度有较大的影响，而本发明所采用CU三角分解法充分利用了各元素之间的关系可以快速形成C、U因子阵的合成阵，主要特点如下：
(1)以左下式4阶Y阵为例，根据C、U阵结构的特点，可建立右下式的4阶合成阵。合成阵中，C、U阵各元素与Y阵元素的关系如中下式所示。
Y11Y12Y13Y14Y21Y22Y23Y24Y31Y32Y33Y34Y41Y42Y43Y44c11u12u13u14c21c22u23u24c31c32c33u34c41c42c43c44]]>
c11=Y11u12=Y12/c11u13=Y13/c11u14=Y14/c11c22=u23=u24=c21=Y21Y22-c21u12(Y23-c21u13)/c22(Y24-c21u14)/c22c32=c33=u34c31=Y31Y32-c31u12Y33-(c31u13+c32u23)[Y34-(c31u14+c32u24)]/c33c42=c43=c44=c41=Y41Y42-c41u12Y43-(c41u13+c42u23)Y44-(c41u14+c42u24+c43u34)]]>
(2)形成C、U因子阵的合成阵不但节省存贮单元，而且有利于充分展示和利用Y阵以及C、U阵元素的相应关系和合成阵中各元素之间的对称关系，还有利于理解和应用本发明方法中提出的过程法C、U阵元素的计算方式。
(3)先按“Г”(倒L)方式确定其对角元以右非零元素的位置，再根据对称性确定对角元以下非零元素的位置。
(4)然后按类似高斯消元的方式分步计算“Г”所包含的上三角中的所有元素，而下三角的所有元素均根据对称性获得，这可减少约50％元素的计算，但多一点赋值语句。这种分步计算的方式与传统的三角分解法中公式法的计算完全不同，可将其称为 过程法。
(5)上述分析表明本发明方法中是按过程法形成C、U因子阵，其各个cij、uij元素不是一次、而是多次分步计算完成，特别是uij元素的计算。
(6)过程法计算过程示意图如图6所示。如对各行对角元以左的某个元素进行类似消元计算时仅计算其同行对角元及以右的元素、而同行对角元以左的所有元素均不被计算(阴影部分)，这些元素均是通过对称性得到。如分别对ck1和ck2进行类似消元计算时，均只要分步计算ckk～ukn的元素，而对角元ckk以左的元素均不计算，可直接由对称性获得。
步骤4：根据方程CUZk＝Ek对方程CWk＝Ek求取Wk阵；通过方程UZk＝Wk求取Zk阵对角元Zkk及以上的非对角元素；
步骤4具体实施过程如下：
(1)对方程CWk＝Ek求解Wk阵的过程
(2)对方程UZk＝Wk求解Zk阵的过程
根据(1)中求出的Wk阵的元素，利用UZk＝Wk方程回代过程求取相应的Zk阵元素。本发明方法中对Zk阵的计算顺序为：Zn,…,Zk,…,Z1，且在计算各个Zk阵中，仅计算对角元Zkk及其以上的元素。
步骤5：根据对称性求取对角元Zkk以左的非对角元素；
由于步骤4中在计算各个Zk阵时仅计算对角元Zkk及其以上的元素，如果需要可根据Z阵元素的对称性得到对角元Zkk以左的元素。这种计算方式可减少50％非对角元素的计算。如果不求取Z阵下三角的非对角元素，步骤5可省略，则计算速度还将大大提高。
本发明方法中Z阵元素的求取过程如图7所示，其中带上标“′”的元素表示根据对称性得到的非对角元素。
步骤6：将Z阵写入数据文件以备后续程序使用。
考虑到程序的结构化，形成Z阵程序到此结束，而所形成的Z阵数据文件可由下一个程序调用执行。
本发明方法与传统的LDU三角分解法相比主要具有以下四点优势：
(1)利用了比LDU三角分解法计算效率更高的CU三角分解法求取Z阵元素。
CU三角分解法需仅求解2个方程组和1个中间矩阵变量Wk，第一个方程组计算变量的个数为n×n，第二个方程组计算变量的个数为n(n+1)/2，因此其计算元素的总数约为1.5n2。即使加上根据对称性得到下三角元素的个数n(n-1)/2，其计算元素的总数也仅约为2n2。因此，本发明方法的计算变量总数约为传统LDU三角分解法的1/2或2/3，且少求取1个中间矩阵变量Hk，因此计算效率更高。
(2)利用合成阵不但节省存贮单元，且有利于理解和应用本发明的过程法。
形成C、U因子阵的合成阵不但节省存贮单元，而且有利于充分展示和利用Y阵以及C、U阵元素的相应关系和合成阵中各元素之间的对称关系，还有利于理解和应 用本发明方法中提出的过程法C、U阵元素的计算方式。
(3)利用了Y阵以及C、U阵元素的对称关系快速形成C、U因子阵。
本发明方法则按“Г”方式，不但能够根据对角元以右非零元素的位置确定对角元以下非零元素的位置，而且可按类似消元方式分步计算其同行对角元及以右的元素，即分步计算“Г”方式下包含上三角中所有元素的内容，而下三角中的所有元素根据对称性获得。这样可减少因子阵形成过程中约50％元素的计算。
(4)利用Z阵元素的对称性及Zk阵元素的计算顺序，仅计算Zk阵对角元Zkk及其以上的非对角元素，而利用对称性得到对角元Zkk以左的非对角元素，这样可减少回代过程50％非对角元素的计算。如果不求取Z阵下三角的非对角元素，则计算速度还将大大加快。
因此，鉴于上述4点，在用Y阵求取Z阵时，本发明方法与传统的LDU三角分解法相比计算速度将大大提高。
附图说明
图1为传统法的LDU三角分解法求取Z阵元素的计算流程图。
图2为本发明方法求取Z阵元素的计算流程图。
图3为传统三角分解法计算各个因子阵的元素时按“┘”(反L)的方式从左上角到右下角计算所有元素的示意图。
图4为传统三角分解法计算各个因子阵的元素时按“行”的方式从上到下计算所有元素的示意图。
图5为传统的LDU三角分解法Z阵元素的计算过程示意图。
图6为过程法计算过程示意图。
图7为本发明方法中Z阵元素的求取过程示意图。
具体实施方式
本发明将通过以下实施举例作进一步说明。
实施例1。
以n×n阶节点系统为例，分别比较传统的LDU三角分解法和本发明方法求取整个Z阵元素的计算过程。比较结果如表1所示。
表1传统的LDU三角分解法和本发明方法求解整个Z阵元素计算过程的比较

根据表1可以看出：
(1)传统的LDU三角分解法是将一列列Zk阵元素全部求解出来，其过程可分解为对LWk＝Ek，DHk＝Wk，UZk＝Hk三个方程组进行求解。这三个方程组各有n个方程，每个方程均要求解n个变量，其计算变量总数为3n2个，而中间矩阵变量有2个。
(2)本发明方法求解Z阵元素仅需求解Zk阵对角元Zkk及其以上的元素。其过程可分解为对CWk＝Ek，UZk＝Wk二个方程组进行求解。
求解方程CWk＝Ek需计算n×n个变量，求解方程UZk＝Wk需计算变量个数为n(n+1)/2，其计算变量的总数约为1.5n2，即使加上根据对称性得到下三角元素的个数n(n-1)/2，其计算变量的总数也仅约为2n2。
因此，本发明方法计算变量的总数约为传统LDU三角分解法的1/2或2/3，且少求取1个中间矩阵变量Hk，因此计算效率更高。
(3)此外根据形成因子阵时元素的计算公式，即使对于同样的因子阵形成过程，在形成LDU因子阵时，计算每个因子阵元素所需的元素总数也远远多于形成CU因子阵时所需的元素总数。
(4)传统方法形成因子阵时要计算所有的元素，并没有考虑利用元素的对称性，计算过程的重复率较高。而本发明方法不但利用对称性来确定对角元及以右元素与对角元以下元素的值和非零元素的位置，而且在计算过程中可利用对称性减少约50％元素的计算。
实施例2。
分别用传统的LDU三角分解法(图1)、以及本发明方法(图2)对IEEE-30、-57、-118节点系统的Y阵求取整个Z阵，比较其在“分解”、“回代”及“分解+回代”过程中的平均计算时间。计算结果如表2所示。
表2二种方法求取IEEE各个系统整个Z阵元素在“分解”、“回代”及“分解+回代”过程中平均计算时间的比较

T1：传统LDU三角分解法“分解”过程的平均迭代时间
T2：本发明方法“分解”过程的平均迭代时间
T′1：传统LDU三角分解法“回代”过程的平均迭代时间
T′2：本发明方法“回代”过程的平均迭代时间
T″1：传统LDU三角分解法“分解+回代”过程的平均迭代时间
T″2：本发明方法“分解+回代”过程的平均迭代时间
从表2可以看出，本发明方法与传统LDU三角分解法相比有以下几点不同：
(1)本发明方法的“分解”过程计算时间减少了约48％；
(2)本发明方法的“回代”过程计算时间减少了约29％；
(3)本发明方法的“分解+回代”过程计算时间减少了约33％；
以上计算结果表明本发明方法在求取Z阵元素时与传统LDU三角分解法相比，不仅提供了一种求取Z阵元素的新方法，且在“分解”、“回代”、“分解+回代”过程的计算速度上都具有较大优势。
本方法可以采用任何一种编程语言和编程环境实现，这里采用C++编程语言，开发环境是Visual C++。