毕业论文-稀疏矩阵技术在短路电流计算中的应用

第页绪论稀疏矩阵的背景稀疏矩阵（sparsematrix）在本科教材中的定义是指非零元素占全部元素的百分比很小的矩阵。早在20世纪，C.F.高斯、C.G.J.雅可比等历史名人就已经研究过利用矩阵稀疏性的办法。后来到了20世纪50年代，有一些处理稀疏问题开始慢慢的出现在了例如技术线性规划和边值问题数值解之中。60年代初期，W.F.廷尼关于直接法的研究标志着现代稀疏矩阵技术的开端。如今，伴随着大规模电网的高速发展，电网的短路容量以及短路电流的相继迅速增大直接导致了短路故障成为电力系统最常见最突出的故障类型之一。不过我们实现短路电流早期预测的方法多种多样，例如采用小波算法，形态学相结合的方法，结合BP神经网络的方法等均取得一定的效果，但其计算所使用的数据量依然巨大。但是根据短路电流预测要求的实时性、快速性，采用稀疏矩阵来提高系统运行速度，增大存储空间具有明显优势。短路电流的预测改革开放以来，伴随着大规模电网的高速发展，现代电网的最大短路电流水平也在“水涨船高”，并且在我国局部地区的短路水平如今已接近所规定的额定值，因此，只有通过对短路电流进行实时且快速的预报威胁方能够保障我国电力系统安全稳定以及持续的发展。在短路电流实时监测过程中，出现的矩阵是大规模的非零稀疏矩阵，占大量存储空间，运行速度慢，不符合我们电网实时监测的快速要求，根据稀疏矩阵的特点，我们在短路电流计算中运用稀疏矩阵的重要目的是：用它来提高存储空间以及算法运行速度。与此同时，这对提高我国大规模电网运行的安全可靠性也是具有十分重要作用的。1.3短路电流早期预测研究现状，现代电网的最大短路电流水平也在“水涨船高”，并且在我国局部地区的短路水平如今已接近所规定的额定值，因此，只有通过对短路电流进行实时且快速的预报威胁方能够保障我国电力系统安全稳定以及持续的发展。在短路电流实时监测过程中，出现的矩阵是大规模的非零稀疏矩阵，占大量存储空间，运行速度慢，不符合我们电网实时监测的快速要求，根据稀疏矩阵的特点，我们在短路电流计算中运用稀疏矩阵的重要目的是：用它来提高存储空间以及算法运行速度。与此同时，这对提高我国大规模电网运行的安全可靠性也是具有十分重要作用的。1.4本课题研究内容及章节安排本文主要研究的是稀疏矩阵在短路电流计算中的应用，采用稀疏矩阵来提高系统运行速度，增大存储空间具有明显优势。本文的主要框架如下所述:第一章是绪论。简要介绍了一下稀疏矩阵的背景，实时短路电流预测的概况和基本理论，以及稀疏矩阵在短路电流预测中的意义和发展现状。第二章简要介绍短路电流计算的基本理论。从介绍短路电流计算的基本概念出发,紧接着分析了短路电流计算的原理以及步骤,最后总结了实时短路电流预测的方法。第三章主要介绍十字链表法如何在短路电流计算中起作用以及基本的原理分析。第四章主要是相关的算例分析与程序介绍。短路电流计算的基本理论2.1短路电流计算的介绍短路电流计算属于电网计算中比较基础的。正是因为配电网辐射形且支路参数r和x相差不大的结构特点，导致有时会出现r>x的情况。而输电网的短路电流计算在220KV以上的情况下可以考虑运用节点导纳矩阵,并且如果选用了稀疏技术，则可以很大程度上达到改良短路电流计算的效率的目的。不过因为从节点编号优化到形成复数因子表耗时且麻烦，因此该方法并不是最优良的首选。许多电力科学笔者在进行探寻适合配电网短路电流计算方面提出了自己方法[2～4]。如何寻找一种对于配电网特点的短路电流最为合适和简单的计算方法，对此进行了大量的国内外文献的阅读和研究，基于此作者运用叠加原理，在此基础上又结合了配电网潮流计算，两项计算与原理的结合形成新的短路电流计算的方法,将abc三相模型运用其中，而在正常运行下所出现在故障点处突然叠加一附加的注入电流被称之为短路故障，首先是进行一次回推前推计算，然后便可以算出各支路的电流以及节点电压。这样情况下的故障点短路电流便是附加的注入电流,而故障点短路电流的计算可借助两方面条件，首先是故障前的三相潮流计算结果，然后是故障边界条件。2.1故障点短路电流计算运用潮流计算的方法可以获得短路故障前的电压、电流,在这种方法的支撑下，首先可以知道的是短路点故障前的电压。如果对短路做一个假设，比作是突然叠加的一个负的电压源,参考下图1,据图1显示，我们可以将短路故障做一个分解，分解成为正常运行方式以及具有一个电压源的故障分量,故障点短路电流即是该电压源所产生的电流。必须考虑到中国10kV配电网比较特殊，是一个不接地系统,基于此，笔者旨在论述短路故障中三相短路故障与两相短路故障的计算。2.1.1短路分析根据下图1(a)所示,若计算故障点短路电流可根据下列算式:式中的短路电流相量分别有I·fa,I·fb,I·fc且等同于a相、b相、c相;从故障点至根节点的所有支路自阻抗的和分别是Zaa,Zbb,Zcc;Zab=Zba,Zbc=Zcb,Zac=Zca,V·a|0|,V·b|0|,V·c|0|为从故障点至根节点的所有支路互阻抗之和;即是故障前故障点处电压。2.2两相短路分析如图1(b),故障分量的边界条件是:式中V·fb,V·fc为故障后b相、c相对地电压;I·fb,I·fc为b相、c相故障短路电流。同时,根据戴维南定理，并结合电压降与电流的关系，便可以得出:由式(2)、式(3)联立求解可得:因此,两相短路时故障点短路电流为:2.3两相短路分析参照国内关于短路电流的计算方法有多种，其中运算曲线法是其中一种比较广泛应用的方法，由于我国的电机参数自身具有较强的实用性等特征。在使用运算曲线方法去计算电力系统短路电流情况下，等值电抗采用的是各台发电机的dx′′，而网络中符合则不计算在内，进行网络化简的前提是必须作出等值网络，由此可以得出的结果分别是含有发电机电动势节点，然后是短路点的简化网络，得出的结果是各电源对于短路点间所出现的转移阻抗。通过统一的功率基准值进行归算，便可以得出等值网络所有阻抗，把各电源以及所出现的短路点间的转移阻抗，将之分别进行归算到有额定容量的电源，由此各电源的计算电抗便可以容易得出。2.4IEC标准参照IEC标准在电力系统中的电源，短路时出现可瞬时逆变运行的静止换流器装置等，当然也可作为网络元件进行处理，但对于出现的稳态短路电流，此时的同步、异步电动机还有静止换流器都不会产生任何影响，也就是不会被看做网络元件。当然在网络元件的出现不同情况的前提下，IEC有着不同的规定，比如网络元件的阻抗方面，IEC标准对以下元件有明确并且严格的规定，包括电力系统电源、异步电动机、发电机等。电阻R和电抗X可以是所有网络元件的参数。在网络中唯一可以加到短路点的是一个等效电压源CUn/3，也就是在短路点处标称电压Un，并在电压Un前加了一个系数C，主要是考虑到电压的变化，有些因素可以不包括在内，比如设备与线路电容还有旋转电机的次暂态特性变化等，参照IEC的标准，计算短路电流可在两种情况下分别进行，即为发电机近端短路还有发电机远端短路。2.5ANSI标准在进行简化短路电流的计算，最恰当的方法便是ANSI标准通过一个等效电压源E与一个等值阻抗串联（均取标幺值）的等值网络。在XEI=K⋅（1）一式中，我们所熟悉的K的求出也是由系统中的短路点以及从短路点看进去的系统的X/R值所共同决定，而不是单看某一方面。2.62三种短路电流计算标准或方法间的比较在近、远端短路，涵盖故障前电压、元件阻抗、网路结构区分、变压器模型等概念，并结合短路电流计算方法，综合上述两种因素进行对比三种短路电流计算标准或方法间的差异，详见表1。后述将结合图1具体阐述如何运用简单电力系统的短路电流计算，并依此分析三种方法计算结果差异所体现的规律性。系统元件参数见表2。上述5个图表是关于具有某种假设前提下的短路电流的计算，而这种假设的前提是通过改变线路等值阻抗计算发电机与短路点间的不同的计算。得出的结果可以参看表3，结果之间的差异图可以参看图2~图5。面对这样一个系统，结合差异图进行对比，要想使三种标准计算的短路电流值之间的差异很小，前提是在发电机与短路点间的计算阻抗较小和较大时的情况。对图2~图5进行对比查看，如果计算电抗有一个额定的范围即0.2～0.8之间时，运算曲线法所得出的短路电流的初始值的计算结果，以及运用IEC标准对短路电流初始值的计算结果，两种结果相对比，之间的差异便一目了然，假设将计算的额定范围进行改变，如果在0.4～1.6之间时，运算曲线法对短路0.05后短路电流周期分量的计算结果与运用IEC标准至今存在着一定的差异。同时也不难看出，不论是水轮机式发电机还是汽轮机式发电机，运用IEC标准与使用ANSI标准同时所计算的短路电流值存在的很大的差异，前者趋于大，后者趋于小。在发电机有水轮机和汽轮机两种情况下，对于不同计电抗下三种方法，依次进行短路电流结果的计算，并且对计算结果进行差异分析，比较之后的分析的优点相对于其他结果有较强的普遍性和指导性。比如说在计算电抗很小时，与之一样的是三种结果的对比所出现的差异也是较小的，我们可以用比较简单方便的运算曲线法去对发电机的近端短路（xjs≤0.2）进行计算。在0.2＜xjs＜1.6前提下，运算曲线法和IEC标准两种普遍且适用的计算方法所得出的计算结果之间的差异还是比较大，当xjs=0.4时，运算曲线法和IEC标准两种计算法所算出的水轮机短路电流周期分量初始值，两者的差异非常明显，差异率达到了20.7%。对于出现此类状况IEC标准是最恰当，且符合实际的。当然，在选择计算标准的同时，电机的参数和电网类型也是非常重要的因素。如果xjs＞1.6时，并且会出现发电机远端短路的特殊情况，比较合适计算电流的方法便是运算曲线法。而且我们也可以从运算曲线法的制定过得到一些启发：结合统计的方法，在不同转移阻抗的前提下，对三种方法计算结果的差异进行统计，由此得出的短路电流差异曲线相比之下更精确且更具代表性。并且更为准确的划分了不同的标准适用于不同的电抗区间的计算。稀疏矩阵在短路电流计算中的应用3.1三元组表存储算法三元组表的类型定义如下:#ｄｅｆｉｎｅＭＡＸＳＩＺＥＮ/＊非零元素的个数最多为Ｎ,由用户定义＊/ｔｙｐｅｄｅｆｓｔｒｕｃｔ{ｉｎｔｒｏｗ,ｃｏｌ;/＊该非零元素的行下标和列下标＊/ＥｌｅｍｅｎｔＴｙｐｅｅ;/＊该非零元素的值＊/}Ｔｒｉｐｌｅ;ｔｙｐｅｄｅｆｓｔｒｕｃｔ{Ｔｒｉｐｌｅｄａｔａ[ＭＡＸＳＩＺＥ+1];/＊非零元素的三元组表。ｄａｔａ[0]未用＊/ｉｎｔｍ,ｎ,ｌｅｎ;/＊矩阵的行数、列数和非零元素的个数＊/}ＴＳＭａｔｒｉｘ;3.2十字链表简介3.2.1十字链表基本原理十字链表用于存储稀疏矩阵的数据结构，其本身是借助单链表为基础扩展形成的。通过图1所示，这种结构是属于十字链表储存结构中5维稀疏方阵式的。其中行链表的定义是单链表中每1个含有附加头结点保存着矩阵的每行数据；而每列数据与之一样的保存，其被称之为列链表。下图2所展现的正是十字链表的存储结点，value、row、col、right、down所代表的是链表中存储节点所包含的5个域，其适用的对象分别是存储矩阵元素的数值、所在行的行号和列号、行链表及列链表指针，将对象结点所在行/列中后继结点的地址用来保存行/列表指针。将表头结点设置为1个，将行链表的附加头结点通过列链表指针链接成单链表，然后进行单链表的链接即借助链表指针辅助行链表的附加头结点，这便是附加列链表；使用同样的方法将列链表进行链接，便形成了附加行链表。保存矩阵的行数和列数需要借助表头结点的row和col，row所代表的行链表附加头结点其作用是用来保存其所在行的行号，而col所代表的列链表附加头结点，与row的作用相似，但是它保存的对象是所在列的列号。所有链表结点都可以借助十字链表的表头结点进行访问，比如附加点接头，如何访问链接点所在行/列的矩阵元素可以通过行/列链表的附加头结点进行。另，从与存储结点所对应的矩阵中某个元素开始，其后继结点也能够进行访问，当然其所在列的后继结点也是可以借助访问的。总结一下，十字链表本身具备下列2个特点：1）结构灵活，扩展性强。插入和删除操作在十字链表的操作中十分便利且高效，并且矩阵的维数可以得到的扩展。2）检索方式灵活。能够在行、列2个方向上进行矩阵元素的检索，这也是十字链表非常有代表的特征，而且另外非常重要的是增加了新的辅助数组，并且对矩阵的主对角元进行直接便利地访问，而且任意矩阵元素所在的行和列的后继结点可以借助矩阵元素所在的结点进行检索。3.3三角分解及前代–回代计算中如何应用十字链表的作用3.3.1三角分解法在进行电力系统仿真计算的过程中，雅克比矩阵和电力网络的节点导纳矩阵都是必要的存储系统，而且系统的微分–代数方程和电力网络方程所求出的值也是十分必要的[11,17]。解决线性方程我们可以采用牛顿法求解系统微分–代数方程时的修正方程，当然线性方程也涵盖了电力网络方程，基于此，我们需要对电力系统仿真中须频繁求解线性方程组Ax=b(1)式中：A为n维矩阵；x、b为n维列向量。求方程式(1)，首先要做一个形成因子，而这个因子是第一要做的是对矩阵A做三角进行分解形成的，再运用前代−回代运算方式，便可以得出(1)的解。关于矩阵A的三角分解流程可以参看图3。参看图3所示，矩阵A中第i行第j列元素是aij，运用高斯消去法对矩阵A进行按列消去的这样一个过程便是矩阵A的三角分解过程，但是计算的下三角元素和主对角元在消去的过程中保留看下来。其中，分2个步骤去对循环1的第p步计算分：1）运用规格化运算来解决第p行主对角元之后的元素；2）消去运算也是为了消去第p列主对角元之下的元素而进行的。参看循环2，若以第p行的主对角线元素app作为开始，并且对该行主对角元之后的元素apj进行依次访问，首先进行规格化计算的是apj，接着转入循环3之中。在3中，对第p列和第j列位于第p行之下的元素进行检索，并运用消去运算对第j列位于第p行之下的元素进行运算。综合上述分析所知，直接访问主对角元是三家分解过程的首要进行的过程，并且顺序检索主对角元开始进行行方向和列方向的检索，并且关于从上三角元素开始进行的列方向顺序检索式非常必要的。由于注入元会产生在三角的分解过程中，所以插入因素的操作时非常必要的。更加快捷的直接访问主对角元的另一个途径是增加辅助数组的十字链表，进行行方向和列方向的顺序检索可以借助矩阵的任意元素进行开始工作，而且十字链表完成插入元素的操作是非常便利并且高效的，参看图4。图4的与图3两个的计算原理是相同，由于计算过程中所出现的十字链表的操作方式，并且完成从属于稀疏矩阵的非零元素的访问，更重要的是增加结点的时间点即出现注入元时。3.4前代–回代计算矩阵A的因子表中第i行第j列元素我们可以假设为wij，而第i行元素则是bi为右端向量b的第i行元素，未知向量x是xi，前代−回代过程中利用因子表求解Ax=b方程组可以参看图5。若借用因子表的主对角元以及下三角元素对右端向量的消去计算过程便是前代过程，而求取未知向量则采用因子表的上三角元素和消去后的右端向量的过程便是回代过程。在前代过程的循环1所列的第i步之中，首先进行一系列的规格化，即第i行主对角元wii对该行的右端向量元素bi所进行的，接着便转入循环2，在循环2中进行顺序检索，则检索所得出的因子表是第i列位于主对角元之下的元素，同时进行的消法运算时相应行上的右端向量元素所做的。代表回代过程的循环1的第i步中，第一步进行的是对bi赋给xi的消去运算，接着下接循环2，而因子表中第i行主对角元之后的元素wij需要循环2进行依次访问，而接下来的回代运算需要同时利用已完成求解的xj去进行。综上以上分析可得出，顺序检索是需要在主对角元开始进行行方向和列方向的双向进行的，我们必须要做的是采用直接访问主对角元的方法，而这是进行前代−回代计算过程。十字链表稀疏矩阵技术而变得更加的高效是在于能够从主对角元开始，对行方向和列方向的顺序检索，完成前代−回代计算因借助了。图6介绍了矩阵A所运用的十字链表存储格式的前代−回代流程，其中向量元素bi、xi是b[i]、x[i]数组元素的代表。3.5十字链表改进方法把所有储结点通过指针将其链接起来，这是十字链表本身具备的特点，当然也会使十字链表的检索方式更具灵活性，而且也会让插入和删除功能更具备便捷性，当然其会任意排放十字链表的存储结点在内存当中的物理位置，而这个位置明显表现出来不是连续分布。总管建构十字链表的整个过程，矩阵元素每增加1个新的，存放该元素的相关信息的关键就是再分配1个存储结点，且插入十字链表的结点将会是那个新分配的，基于上述，在内存空间通常会出现属于十字链表的各个结点散列分布其中，但是其弊端便是在稀疏矩阵的运算效率将会降低。笔者拟借助两种内存分配方式进行稀疏矩阵的计算，并将对比两者的计算效率：1）存储十字链表需要分配1块连续的内存空间；2）在内存空间中散列的这些十字链表的结点，其前提是根据需求增加一个矩阵元素才能分配1个存储点。试验结果表明采用第1种方式时稀疏矩阵的运算效率高于第2种方式的效率，通过运用稀疏矩阵的运算方法对第1和第2两种实验结果的对比，便可以看出第1种的效率要高于第二种，造成这种现象的原因是由于计算机的高速缓冲存储器的工作原理所造成的。在介于中央处理器(centralprocessingunit，CPU)和主存储器之间的高速、小容量存储器便是高速缓冲存储器。纵观整个计算机发展的历史，主存储器存取速度相对于CPU的操作速度，前者的速度相对来说的要慢得多，因而CPU出现的问题是自身具备的高速处理能力不能得到充分发挥，由此也会对计算机的工作效率产生非常大的影响。由于上述原因，普遍采用的是高速缓冲存储器进行缓和CPU与主存储器之间速度不匹配的矛盾。相对于主存储器的庞大的存储量来说，高速缓冲存储器存储容量相比来说就显得比较小，只有住储存器的百分之几，虽然容量是其局限，但是其存取速度能却能和CPU相埒。硬件本身也可以自动户必须办法告诉缓冲储存器与主存储器之间信息的调度和传送。依据程序局部性这一原理，有一个很大的可能是正在使用的主存储器中某一单元临近的单元将会被=被使用。高速缓冲存储器可以直接进行。参照内存分配的方式如2时，突出特征是在内存中十字链表的结点其分布并不连续，而其他的结点也并不一定存在与某个链表结点邻近的存储单元中，于是当访问某个元素的时候，、会出现在高速缓冲存储器中不存在该元素的结点，那么相对的是十字链表的其他结点也并不一定存在于即将被调入高速缓冲存储器中。于是可以得知，梳理整个稀疏矩阵运算过程，较之于在内存空间中把十字链表的结点进行散列分布，并将十字链表存储于连续的内存空间中，对矩阵元素的访问CPU可以借助高速缓冲存储器进行，基于此运算效率也可以得到较高的提升。此外，因为注入元会出现在稀疏矩阵三角分解过程中，倘若采取内存分配方式2，那么在注入元出现在三角分解过程中时，要把1个存储结点分配给该元素，于是，内存分配操作会频繁地出现在三角分解的过程，随之出现的是计算效率的降低。如果把足够数量的存储结点进行预先的分配的话，则已经分配好的结点可以在出现注入元时直接使用，有反复分配内存所产生的负担也可以有效避免。如前所述，十字链表分配内存可以采用的方式为1，提高效率的最佳途径也是在稀疏矩阵的运算过程中，将高速缓冲存储器本身具备高速的特点进行最好的发挥。那么由于三角分解过程中因为注入元的出现须频繁进行内存分配对计算机所造成的额外负担可以得到很好的避免。综合上述原因，对于十字链表的改进笔者给出了自己的方法：1）以十字链表进行存稀疏矩阵的存储，在十字链表的构建之前，第一需要进行的是非零元素在稀疏矩阵的个数，还要估算注入元在三角分解过程中可能会出现的个数，而且需要对十字链表所需的存储节点的个数惊醒判断，然后十字链表需要使用根据所需结点个数分配1块连续的存储空间。2）倘若出现注入元个数超过三角分解过程所预计的个数，则需要在注入元内保存分配存储结点，而这个时候需要这样的一个存储空间即分配1块连续的包含一定数量结点的存储空间。假如注入元出现在之后的计算过程中，则可以对已经分配的存储结点进行直接使用，当然内存分配的操作也可以避免频繁进行。3.6十字链表潮流方法3.6.1导纳矩阵的形成在常规的潮流计算中,要预先在程序中开出足够大的数组以达到形成导纳矩阵的目的。静态数组的缺点依旧不能通过采用了压缩存储技术来克服。在程序无法准确洞悉所需内存空间的大小的前提下,我们只得开出比较大的数组来，但是在这样节点数比较小的情况下内存空间浪费了许多许多;同理，在节点数大时,程序又无法进行处理。因此，十字链表的优点就体现在能够让所有结点所需内存都是动态申请的,其中包括了表头结点(其多少由节点数决定)以及非零元素结点(其多少由支路数决定)。在这个矩阵中,对角元的Data为节点自导纳,非对角元的Data为该非零元素对应节点与其表头结点对应节点之间的互导。在读入节点数之后,表头结点链表将被程序初始化。表头结点数目比节点数大一,最后一个表头结点链接的链表存储节点对地导纳。规则是按照每读入一个节点便在对应的表头结点后插入一个对角元,并将对角元指针指向它;同理，每读入一条支路,如果是一条新的支路,就在此支路的每一个节点对应的表头结点后都插入一个对应另一个节点的非零元素结点,然后对这两个节点的互导进行修改,否则直接修改互导。在所有互导形成之后,就可直接计算每个节点的自导了。总的来说，方法就是累加本表头结点后的链表中的所有非零元素的Data值再乘以-1。节点的对地导纳也可选择将这一部分计算在里面。root指针指向总表头结点也就是我们表头结点链表中的第一结点。这种链表的结构见图3所示。其中包含了网架的拓扑结构和一些有联系的数值。程序中这个矩阵是基本不变的。具有相对独立性便是原始数据存储在这个矩阵中的优点。3.6.2本文潮流计算采用PQ分解法来形成矩阵B1,B2通过把原导纳矩阵去掉松弛节点形成的矩阵B1,B1矩阵所不同的是节点的行号。导纳矩阵中节点的行号和原有节点的保持一致,而B1矩阵中节点的行号为之前导纳矩阵中去掉松弛节点后重新排成的节点号。这种节点号的对应关系必须有程序所记录下来。类比可知，原导纳矩阵去掉松弛节点和PQ节点后形成了矩阵B2。同理,程序也必须记录B2矩阵的节点号与原导纳矩阵节点号的对应。另外可以通过只访问表头结点链表中对地导纳所对应的链表并且按行号进行搜寻则可得到相应的对地导纳。在PQ分解法中，虽然B1,B2矩阵每次只需形成一次但是实际上如果使用牛顿法的话,那么每次都需要形成雅可比矩阵,在这个时候，我们的导纳矩阵的相对独立性就显得较有优势。十字链表法其存储结构提供了将数学上的解方程方法与网络分析相分离的基础，采取的是通过把网架拓扑结构与相关数值一起存储的办法。接着的下面的几个步骤就是纯数学理论研究方面就可以解决的问题了。3.6.3最优排序及其它一些问题在方法原理方面，最优排序与其它存储结构的是一致的,但是在注入元的处理方法上，两者便有出入了。一般而言，静态数组的压缩存储对注入元的处理方法通常是采用在每行末尾预留空间,相应的，这也带来了一些不可避免的缺陷，这些权限的来源也有一部分是因为在于预留空间大小的设定问题。然而在十字链表方法中,注入元的处理并不像以上的那么复杂，相反就十分方便了。对注入元的处理只需在排序时考虑,而没有必要在形成矩阵或在系统分析时就进行考虑。即使出现排序时产生了注入元的情况也只需在对应行各插入一个新的结点就可以轻松解决了。在B1,B2矩阵形成并排过序的前提下,只需要解方程即可。主要是因子表的形成与消元和回代。可以采用LU分解法求因子表。原理相同,只需注意十字链表的存储结构。形成了因子表并消元回代之后,方程解出了。然后按照流程,进行反复迭代。当出现功率偏差量<收敛精度的时候便可结束计算。3.7算例由于原始数据中没有计入第9号节点的对地导纳，则下面是根据IEEE14节点模型计算的导纳矩阵(由于篇幅关系，本文只列出了部分,把导纳矩阵的实部略过),。导纳矩阵虚部为:[1,1]-19.49807[1,5]4.234983[1,2]15.263086[2,2]-30.354715[2,5]5.193927[2,4]5.115839[2,3]4.781863[2,1]15.263086[14,14]-5.344014[14,13]2.314964[14,9]3.029051由上可知,导纳矩阵输出的数据是按照十字链表中的顺序输出的，并且要先输出对角元(访问对角元指针即可得)B2矩阵:[1,1]-38.63517[1,4]1.8555[1,3]4.889513[1,2]21.578554[2,2]-35.527538[2,1]21.578554[3,3]-19.549006[3,1]4.889513[3,4]9.090083[9,9]-5.344014[9,8]2.314964[9,4]3.029051可以看出，原导纳矩阵和B2中结点的对应关系如下:12345678945791011121314B2最优排序后:[1,1]-35.527538[1,7]21.578554[2,2]-8.497018[2,3]4.402944[3,2]4.402944[3,3]-14.768337[3,9]10.365394[9,3]10.365394[9,6]3.029051[9,7]1.8555[9,8]9.090083[9,9]-24.282507B2矩阵排序后结点号与原B2矩阵中结点号的对应关系:123456789718932456限于篇幅B1矩阵略。3.8结论一般的数组存储方法归类于静态数据结构并且在程序的说明部分须给出其类型定义或变量说明。通常数组依然是空间预留的一种实现,这个意思也就是说也就还是存在静态数组的些许不足之处。有时候malloc函数或new操作符动态常用在某些程序中以来申请数组,不过因为其十字链表归根到底动态数据结构的这一性质,往往使得它的scale大小在程序执行过程中是可以变动的。在程序执行时，十字链表中结点的数目呈现动态增长,与此同时，导纳矩阵的数目随着数据的输入逐步形成。存储的分配较为灵活,更充分的利用内存是使用十字链表存储的一大特色，也是一大优点。结点的插入操作导致形成了十字链表法中的导纳矩阵。对应于电力系统中节点、利用十字链表法来实现线路的增加或删除等的结点操作比静态数组要简洁且迅速。例如在添加线路的方法与读入数据时添加线路时，一般情况下是采取分开存储拓扑结构和网架数据在两个数组里的方法。与此同时，为达到访问拓扑结构和数据的目的再通过另外的一个数组作索引来实现。为将网架的拓扑结构与数据存储在一起我们大家通常采用十字链表的存储方法。因此可在潮流计算中将系统分析方法与数学处理方法这两种方法分离开来以求达到各种数据的处理模块化的目的。封装的十字链表作为一种自定义的数据类型,可使程序更便于编制和维护，具有极强的可读性。算例分析4.1算例分析十字链表数据：print：print函数可以把函数图形保存成图片：[plain]

\o"viewplain"viewplain\o"copy"copyminbnd

=

-4*pi;

maxbnd

=

4*pi;

t

=

minbnd:0.1*pi:maxbnd;

plot(t,

sin(t),

'g',

'Linewidth',

2);

line([minbnd,

maxbnd],

[0,

0]);

%绘制x轴

axis([-10,

10,

-2,

2])

%定义显示的坐标区间:x在(-10,10)之间，y在(-2,2)之间

grid

on;

title('sin(x)');

xlabel('x');

ylabel('sin(x)');

print('-dpng','sin.png');

%保存为png图片，在Matlab当前的工作目录下

如下：打开Matlab当前的工作目录下可以看到有sin.png图片了print('-dpng','sin.png')表示保存为png图片，文件名为sin.png，其中第一个参数可以是：-dbmp：保存为bmp格式-djpeg：保存为jpeg格式-dpng：保存为png格式-dpcx：保存为pcx格式-dpdf：保存为pdf格式-dtiff：保存为tiff格式fprintf：fprintf函数可以将数据按指定格式写入到文本文件中：[plain]

\o"viewplain"viewplain\o"copy"copydata

=

[5,

1,

2;

3,

7,

4];

[row,

col]

=

size(data);

for

i=1:row

for

j=1:col

fprintf('data(%d,

%d)

=

%d\n',

i,

j,

data(i,

j));

%直接输出到屏幕；类似于C语言的输出格式

end

end

fprintf(fid,format,data)中的fid表示由fopen函数打开的文件句柄，如果fid省略，则直接输出在屏幕上，format是字符串形式的输出格式，data是要输出的数据。其中format可以为：[plain]

\o"viewplain"viewplain\o"copy"copy%c

单个字符

%d

有符号十进制数（%i也可以）

%u

无符号十进制数

%f

浮点数（%8.4f表示对浮点数取8位宽度，同时4位小数）

%o

无符号八进制数

%s

字符串

%x

小写a-f的十六进制数

%X

大小a-f的十六进制数

输出到文件：[plain]

\o"viewplain"viewplain\o"copy"copydata

=

[5,

1,

2;

3,

7,

4];

[row,

col]

=

size(data);

%求出矩阵data的行数和列数

%加t表示按Windows格式输出换行，即0xOD

0x0A，没有t表示按Linux格式输出换行，即0x0A

fid=fopen('test.txt',

'wt');

%打开文件

for

i=1:row

for

j=1:col

fprintf(fid,

'%d

',

data(i,

j));

%类似于C语言的输出格式

end

fprintf(fid,

'\n');

end

fprintf(fid,

'This

is

a

string\n');

fprintf(fid,

'%X',

hex2dec('ABCD'));

fclose(fid);

%最后不要忘记关闭文件！

就会在Matlab当前的工作目录下生成test.txt文件从文件中读取：我们可以使用fscanf函数fscanf：[plain]

\o"viewplain"viewplain\o"copy"copy%加t的理由和上面一样

fid=fopen('d:\test.txt',

'rt');

%把数据读到data中。其中data是2*3的矩阵

data=fscanf(fid,

'%d',

[2,

3]);

s=fscanf(fid,

'%s');

d=fscanf(fid,

'%X');

%关闭文件

fclose(fid);

disp(data);

disp(s);

disp(d);

从输出的结果看，我们发现fscanf读取数据时会忽略空格，直到回车为止！disp：disp函数直接将内容输出在Matlab命令窗口中：[plain]

\o"viewplain"viewplain\o"copy"copy%单字符串输出：

disp('Hello

World!');

%不同类型数据输出：

num1

=

1;

num2

=

2;

disp([

num2str(num1),

'

+

',

num2str(num2),

'

=

',

num2str(num1+num2)]);

输出：HelloWorld!1+2=3致谢路漫漫其修远，四年的本科生活只是人生长河中一个短暂停顿的驻点，然而就在这短暂的时光中我却学到了很多宝贵的东西，不仅学会了如何独立思考进行科学研究，更重要的是明白了许多人生的哲理，让我受益终生。在本文的撰写、定稿过程中，何老师提出了很多宝贵的意见，在此表示诚挚的感谢。再次衷心地感谢何洪英教授，感谢您在我四年的本科生学习中给予的指导，提供的帮助。您严谨的治学风范和乐观的生活态度将是我一生学习的榜样！感谢曾经朝夕相处的朋友认识他们是我莫大的幸福；感谢班长陈峰泉在撰写论文过程中给予的帮助；感谢实验室的师兄师姐们给予我的支持和帮助。感谢我敬爱的父母，感谢他们一直以来对我的关爱和支持，没有他们无私的付出，便没有我今天的成长，他们是我在人生路上不断前进的动力。最后，向所有帮助和关心过我的老师和朋友，表示我最真挚的感谢！参考文献[1]KundurP．电力系统稳定和控制[M]．北京：中国电力出版社，2002：33-388．[2]王成山，王守相．分布式发电供能系统若干问题研究[J]．电力系统自动化，2008，32(20)：1-4,31．WangChengshan，WangShouxiang．Studyonsomekeyproblemsrelatedtodistributedgenerationsystems[J]．AutomationofElectricPowerSystems，2008，32(20)：1-4,31(inChinese)．[3]杨金刚，房大中，李传栋．中期电压稳定的并行仿真算法[J]．电网技术，2009，33(3)：8-14.YangJingang，FangDazhong，LiChuandong．Parallelsimulationformid-termvoltagestabilityanalysis[J]．PowerSystemTechnology，2009，33(3)：8-14(inChinese)．[4]吴红斌，丁明．用于电力系统暂态稳定仿真的可变步长牛顿法[J]．中国电机工程学报，2010，30(7)：36-41．WuHongbin，DingMing．Newtonmethodwithvariablestepsizeforpowersystemtransientstabilitysimulation[J]．ProceedingsoftheCSEE，2010，30(7)：36-41(inChinese)．[5]ChenJJ，CrowML．Avariablepartitioningstrategyforthemultiratemethodinpowersystems[J]．IEEETransonPowerSystems，2008，23(2)：259-266．[6]TinneyWF，WarlkerJW．Directsolutionsofsparsenetworkequationsbyoptimallyorderedtriangularfactorization[J]．ProceedingsoftheIEEE，1967，55(11)：1801-1809．[7]何银菊，宋玮，周庆捷，等．面向对象的电力系统潮流计算与静态安全分析[J]．电网技术，2001，25(8)：11-14．HeYinju，SongWei，ZhouQingjie，etal．Anobject-orientedpowersystemloadflowcalculationandstaticsecurityanalysis[J]．PowerSystemTechnology，2001，25(8)：11-14(inChinese)．[8]杜正春，牛振勇，方万良．基于分块QR分解的一种状态估计算法[J]．中国电机工程学报，2003，23(8)：50-55．DuZhengchun，NiuZhenyong，FangWanliang．AblockQRbasedpowersystemstateestimationalgorithm[J]．ProceedingsoftheCSEE，2003，23(8)：50-55(inChinese)．[9]王晓东，李乃湖，丁恰．基于稀疏技术的原对偶内点法电压无功功率优化[J]．电网技术，1999，23(3)：23-26．WangXiaodong，LiNaihu，DingQia．Aprimal-dualinteriorpointmethodforoptimalvoltage/reactivepowercontrolwithsparsitystructure[J]．PowerSystemTechnology，1999，23(3)：23-26(inChinese)．[10]王成山，王丹，郭金川，等．基于网式链表–双层结构的电力系统时域仿真算法[J]．电力系统自动化，2008，32(16)：6-10．WangChengshan，WangDan，GuoJinchuan，etal．Apowersystemtime-domainsimulationplatformbasedonnet-formchainanddoublelayeralgorithm[J]．AutomationofElectricPowerSystems，2008，32(16)：6-10(inChinese)．[11]张伯明，陈寿孙．高等电力网络分析[M]．北京：清华大学出版社，1996：20-30．[12]HakavikB，HolenAT．Powersystemmodellingandsparsematrixoperationsusingobject-orientedprogramming[J]．IEEETransonPowerSystems，1994，9(2)：1045-1051．第35卷第5期电网技术39[13]ZhuHJ，CaiZX．Object-orientedmodelingofsparsematrixoperationinpowersystemsoftware[C]//6thInternationalConferenceonAdvancesinPowerSystemControl,OperationandManagement．HongKong，China：TheInstitutionofElectricalEngineersHongKong(IEEHK)，2003：732-737．[14]毛安家，郭志忠．电力系统计算中的二维稀疏结构技术[J]．继电器，2001，29(1)：19-21．MaoAnjia，GuoZhizhong．The2-ordersparsestructureinpowersystemcalculations[J]．Relay，2001，29(1)：19-21(inChinese)．[15]朱凌志，安宁．基于二维链表的稀疏矩阵在潮流计算中的应用[J]．电网技术，2005，29(8)：51-54．ZhuLingzhi，AnNing．Applicationoftwo-dimensionalchaintablebasedsparsematrixinpowerflowcalculation[J]．PowerSystemTechnology，2005，29(8)：51-54(inChinese)．[16]殷人昆．数据结构(用面向对象方法与C++语言描述)[M]．北京：清华大学出版社，2007：43-181．[17]余贻鑫，陈礼义．电力系统的安全性和稳定性[M]．北京：科学出版社，1988：107-149．[18]高速缓冲存储器的作用和工作原理[EB/OL]．[2009-07-30]．/storage-systems-114478.htm．[19]SPARSE1.3：Asparselinearequationsolver[EB/OL]．[2008-05-05]．/IPRO/Software/Description/sparse1.3.html．[20]王成山，高毅，王丹，等．考虑直流系统开关特性控制的变步长仿真算法[J]．中国电机工程学报，2009，29(34)：16-21．WangChengshan，GaoYi，WangDan，etal．Variable-stepsimulationmethodconsideringHVDCcontrolswithswitchingcharacteristics[J]．ProceedingsoftheCSEE，2009，29(34)：16-21(inChinese)．[21]王丹，王成山．基于数值微分法求导的分布式发电系统仿真算法[J]．电力系统自动化，2009，33(17)：81-85．WangDan，WangChengshan．Astabilitysimulationmethodofdistributedgenerationsystembasedonnumericaldifferentiationmethod[J]．AutomationofElectricPowerSystems，2009，33(17)：81-85(inChinese)．[22]MurachM，VachranukunkietP，NagvajaraP，etal．OptimalreconfigurableHW/SWco-designofloadflowandoptimalpowerflowcomputation[C]//IEEEPowerEngineeringSocietyGeneralMeeting．Montreal，Canada：IEEE，2006：1-5．[23]黄永皓,夏清等著.电力市场理论研宄与应用[M].北京:中国电力出版社,2002.[24]于尔座,韩放,谢等.电力市场[M].北京:中国电力出版社,1998.附录要实现稀疏矩阵的三元组存储方式并不难，举个例子，每个元素有三个域分别是i,j,e。分别代表了该非零元的行号、列号以及值。十字链表的形式可以理解成每一行是一个链表，而每一列又是一个链表。如图所示：

通过上图可知，每个结点还需要分别存放它横向和纵向的下一个结点的地址。以形成一个类似十字形的链表的结构。那么每个结点的结构体定义也就油然而生了。

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?typedef

struct

OLNode

{

int

i,

j;

//行号与列号

ElemType

e;

//值

struct

OLNode

*right,

*down;

//指针域

}OLNode,

*OList;

这样我们对结点的插入与删除就要修改两个指针域。为了方便对结点的操作，我们要创建头指针或者头结点。.怎么创建头指针或者头结点？可以通过创建OLNode结构的结点形成一段连续的地址空间来指向某一行或者某一列中的结点。

或者创建指针数组，数组元素中存放的地址就是某一行或者某一列的第一个结点的地址。来下面分析下两种方法：

第一种方法会浪费大量的空间，而因为指针变量的空间都是4个字节，所以相对来说第二种节省空间。毫无疑问我们选择第二种，也就是创建头指针。那么第二种我们用什么来实现？是数组还是动态内存分配？如果用数组我们要预先定义行和列的最大值，显然这不是一个好主意，而动态内存分配的方法我们可以在用户输入了行数与列数之后分配相应的一段地址连续的空间。更为灵活，所以我们选择动态内存分配。[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?typedef

struct

{

OLink

*Rhead,

*Chead;

int

mu,

nu,

tu;

//

稀疏矩阵的行数、列数和非零元个数

}CrossList;

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?int

main(void)

{

CrossList

M;

CreateSMatrix(&M);

}

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?int

CreateSMatrix(CrossList

*M)

{

int

i,

j,

m,

n,

t;

int

k,

flag;

ElemType

e;

OLNode

*p,

*q;

if

(M->Rhead)

DestroySMatrix(M);

do

{

flag

=

1;

printf("输入需要创建的矩阵的行数、列数以及非零元的个数");

scanf("%d%d%d",

&m,

&n,

&t);

if

(m<0

||

n<0

||

t<0

||

t>m*n)

flag

=

0;

}while

(!flag);

M->mu

=

m;

M->nu

=

n;

M->tu

=

t;

return

1;

}

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?//创建行链表头数组

M->Rhead

=

(OLink

*)malloc((m+1)

*

sizeof(OLink));

if(!M->Rhead)

exit(-1);

//创建列链表头数组

M->Chead

=

(OLink

*)malloc((n+1)

*

sizeof(OLink));

if(!(M->Chead))

exit(-1);

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?for(k=1;k<=m;k++)

//

初始化行头指针向量;各行链表为空链表

M->Rhead[k]=NULL;

for(k=1;k<=n;k++)

//

初始化列头指针向量;各列链表为空链表

M->Chead[k]=NULL;

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?do

{

flag

=

1;

printf("输入第%d个结点行号、列号以及值",

k);

scanf("%d%d%d",

&i,

&j,

&e);

if

(i<=0

||

j<=0)

flag

=

0;

}while

(!flag);

p

=

(OLink)

malloc

(sizeof(OLNode));

if

(NULL

==

p)

exit(-1);

p->i

=

i;

p->j

=

j;

p->e

=

e;

[cpp]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?if(NULL==M->Rhead[i]

||

M->Rhead[i]->j>j)

{

//

p插在该行的第一个结点处

//

M->Rhead[i]始终指向该行的第一个结点

p->right

=

M->Rhead[i];

M->Rhead[i]

=

p;

}

[cpp]\o"viewplain"viewplain\o"copy"copyHYPERLINK"/