潮流计算软件设计.docx
- 文档编号:8998113
- 上传时间:2023-05-16
- 格式:DOCX
- 页数:47
- 大小:274.70KB
潮流计算软件设计.docx
《潮流计算软件设计.docx》由会员分享,可在线阅读,更多相关《潮流计算软件设计.docx(47页珍藏版)》请在冰点文库上搜索。
潮流计算软件设计
电力系统潮流软件设计
1.原始数据的输入
目前计算机的速度和计算方法已经使我们能很快得到计算结果,但是上机以前的准备工作却非常耗费时间,而且也容易出错。
因此在程序设计时,必须尽可能减轻上机前的准备工作,尽可能利用计算机代替人工繁琐的工作。
在这里,原始数据的填写格式是很关键的一个环节,它与程序使用的方便性和灵活性有着直接的关系。
原始数据输入格式的设计,主要是从使用的角度出发,原则是简单明了,便于修改。
输入格式的简单明了就可以减轻数据填写的工作量,并减少或避免程序使用者在填写数据时发生错误。
电力系统潮流计算往往需要进行多种运行方式的调整和比较,因此在数据格式上考虑计算过程中修改数据的方便性就显得非常重要。
以下所介绍的本潮流程序中所用到的6个信息:
N:
系统节点的总数;
Nb:
系统中支路数,即输电线路条数、变压器数的总和;
Ng:
发电机节点总数;
Nl:
负荷节点总数;
V0:
系统平均电压,在迭代过程中,以它作为电压的初值;
eps:
迭代收敛所要求的精确度。
潮流计算程序所需要的原始数据,分别归纳为下几个结构体:
支路数据结构体Branch_Type,有5个数据成员,对应于每条支路的5个数据。
该结构体数组定义为:
structBranch_Type
{
inti,j;
doubleR,X,YK;
}Branch[400];
当支路为输电线路时,这5个数据成员分别表示:
i:
输电线路一端的节点号;
j:
输电线路另一端的节点号;
R:
输电线路的电阻;
X:
输电线路的电抗;
Y0:
输电线路充电电容的容纳。
(a)(b)
图1-1网络支路的等值电路
如图1-1(a)所示。
当支路为变压器支路时,这5个数据成员分别表示:
i:
变压器一端的节点号;
j:
变压器另一端的节点号,这两个节点号有一个带有负号,作为变压器支路的标志;
R:
变压器的电阻;
X:
变压器的电抗(RT和XT都是归算到变压器标准变比侧的数值);
Y0:
变压器的非标准变比(设在节点号为负的一侧);
变压器的模拟电路如图1-1(b)所示。
发电机节点和负荷节点数据分别定义为结构体数组,结构体有4个数据成员,其内容是相同的。
发电机节点数据定义为:
structGenerator_Type
{
doubleP,Q;
inti;
doubleV;
}Generator[50];
负荷节点数据定义为:
structLoad_Type
{
doubleP,Q;
inti;
doubleV;
}Load[300];
对于发电机节点,P、Q填正号;对于负荷节点,P、Q填负号。
对于发电机节点和负荷节点,若为PQ节点,这些数据成员分别表示:
P:
节点的有功功率;
Q:
节点的无功功率;
i:
节点的编号;
V:
该节点正常运行的电压。
若节点为PV节点时数据成员分别表示:
P:
节点的有功功率;
Q:
节点无功功率的上限;
i:
节点的编号;
V:
节点需要维持的电压,负号是PV节点的标志。
在节点数据输入计算机后,为了提高计算效率,应统计PV节点的总数Npv,并形成PV节点结构体数组。
PV节点结构体定义为:
structPVNode_Type
{
doubleV;
inti;
}PVNode[100];
每个PV节点有两个数据,第一个数据为PV节点的给定电压V,第二个数据为相应的节点号i。
在形成PV结构体数组的同时,把发电机节点或负荷节点数组中V前面负号去掉。
2.数据优化
在实用潮流计算程序中,对数据的输入次序应不加以限制,这样便于数据的填写和修改。
输入以后,在计算机内对数据再进行排队和整理。
2.1负荷节点的优化
为了简化程序,需要对负荷节点的顺序进行优化,优化后各节点按照节点编号的顺序进行排列,该部分程序的实现框图为图2-1。
图2-1负荷节点优化框图
其中定义了一个中间工作结构体数组LoadN[300],该结构体定义为:
structLoad_Type
{
doubleP,Q;
inti;
doubleV;
}LoadN[300];
用来存放优化后的节点参数。
2.2发电机节点的优化
由于发电机节点的优化与负荷节点优化相似,因此只给出程序框图(图2-2)。
图2-2发电机节点优化框图
2.3对支路数据的排队整理优化
为了使支路数据的排列方式适合形成导纳矩阵的上三角部分,整理过后支路数据按以下次序排列:
(1)两端节点号应把小号排在前边,大号排在后边。
(2)各支路按其小节点号的顺序排列。
实现这部分的框图为图2-3。
首先对框图中的符号做一介绍:
图2-3支路数据优化框图
BranchN[400]:
为一个中间工作结构体数组,该结构体定义为:
structBranchN_Type
{
inti,j;
doubleR,X,YK;
}BranchN[400];
用来存放优化后的支路数据。
m[k]:
用来存放小节点号为i支路的另一节点号。
s[k]:
用来存放该支路原先编号。
该函数首先将原始数据中节点按照
(1)的要求优化,然后按照
(2)的要求对支路顺序进行再次优化。
2.4对支路特殊情况的优化
当某些支路的电阻R大于电抗X的2倍时,影响到PQ分解法的收敛性,这时程序将自动把支路分成两条支路(增加一个节点),满足PQ分解法要求电阻小于电抗的条件。
由图2-4构成的程序来实现这部分功能。
对框图中变量的说明:
a:
用来存放原始数据中支路的个数。
i:
为计数变量。
该函数的实现思想是:
当检查出某个支路满足上述情况时,则在该支路中增加一个节点,节点号为N+1;增加的支路号为Nb+1,原支路仅保留电阻,新支路参数为原支路的电抗,而对地支路不变。
图2-4对支路特殊情况的优化框图
2.5平衡节点的优化
为了简化程序,去掉一些判断,要求把平衡节点排在最后,即作为第N个节点,同时还要求这个节点为负荷节点。
如果平衡节点没有负荷,则该节点的负荷功率填零。
这样保证了节点N既是发电机节点,又是PV节点,又是负荷节点。
图2-5平衡节点优化框图
图2-5用来实现这一功能。
对框图中变量的说明:
j:
用来存放原始数据中平衡节点的节点号。
t[i]:
为一个标识数组。
t[i]=1表示第i条支路的小节点号为j,t[i]=2表示第i条支路的大节点号为j。
s[i]:
也为一个标识数组。
s[i]=1表示第i条支路的小节点号为N,s[i]=2表示第i条支路的大节点号为N。
该函数的设计思想是:
先把平衡节点找出(如图2-5A部分),然后与最后一个节点进行交换(如图2-5B、C部分),同时保持原有支路数据形式不变(如图2-5D部分)。
3.稀疏导纳矩阵的形成
3.1基本公式
当电力系统中i、j两点间输电线路的阻抗为zij时,节点i、j之间互导纳为
Yij=-
=-yij(3-1)
式中:
yij是阻抗zij的倒数,即输电线串联支路的导纳;Yij是导纳矩阵中i行j列的非对角元素。
由于导纳矩阵的对称性,一般
Yij=Yji
支路i、j对导纳矩阵中i、j两行对角元素的影响可表示为如下的增量:
△Yii=△Yjj=
=yij(3-2)
这里导纳矩阵对角元素Yii和Yij也就是节点i、j的自导纳。
当节点i连有导纳为Y0i的接地支路时,它对导纳矩阵的影响仅仅使i行对角元素增加如下的分量:
△Yii=Y0i(3-3)
(a)(b)
图3-1变压器的等值电路
当i、j两节点间的支路是非标准变比的变压器时[如图3-1(a)所示],我们可用П型等值电路来模拟[见图3-1(b)所示],因此i、j之间互导纳可按下式计算:
Yij=-
=-
yij
i、j节点自导纳分别有如下的增量:
△Yii=
yij(3-5)
△Yjj=yjj(3-6)
3.2稀疏导纳矩阵的处理
电力系统的导纳矩阵不仅具有对称性,而且具有稀疏性。
当i、j节点之间没有直接联系时,导纳矩阵中非对角元素Yij及Yji应为零。
由于导纳矩阵的对称性,在计算机中可以只存储其上三角部分或下三角部分,在本程序中,储存导纳矩阵的上三角部分及对角元素,因此,其中每个非对角元素Yij的下标都满足i 对角元素按节点编号顺序存放在对角元素结构体数组中: structYii_Type { doubleG,B; }Yii[300]; 其中G存放对角元素的实部,B存放对角元素的虚部,每个数组的元素个数与系统节点数相等。 由于上三角矩阵中非对角元素和系统中不接地支路一一对应,非对角元素的个数等于网络中不接地支路数Nb。 为了节约内存和提高计算速度,在计算机内存中只储存非零元素,把非零非对角元素“挤实”在一起。 为了识别非对角元素的行号和列号,我们在每个元素后存放相应的列下标。 非对角元素结构体数组定义为: structYij_Type { doubleG,B; intj; }Yij[400]; 按照这样的排列,取一个互导纳,就可以同时把该元素的列号取出来,为了判断该元素的行号,需要借助于数组NYseq[300]。 数组NYseq[300]按导纳矩阵行号的顺序存放各行非对角元素的首地址(事实上,存放的是各行第一个非对角元素在导纳矩阵非零非对角元素中的顺序号)。 本程序中还定义导纳矩阵中各行非对角元素的个数NYsum[300]。 一般我们有 NYsum[i]=NYseq[i+1]-NYseq[i](3-7) 由于非对角元素是逐行向下排列的,所以就很容易判断出各非对角元素的行号。 按照对支路原始数据处理的要求,可以得出支路排列的顺序和导纳矩阵非对角元素的排列顺序完全一样。 因此只要顺序取出支路数据,按照式(3-1)求出倒数取负号之后,连同该支路的大节点号(即列下标)顺序送入Yij数组,就形成了导纳矩阵的上三角部分。 3.3导纳矩阵形成过程及框图 在本程序中,适应P-Q分解法的需要,导纳矩阵分为两步形成。 第一步只用不接地支路构成导纳矩阵,不考虑接地支路(包括变压器非标准变比)的影响。 这里同时形成两个导纳矩阵,即通常意义上的系数矩阵Y和不考虑输电线路电阻的系统导纳矩阵Y ,以适应BX法的要求,这两个导纳矩阵实际上只是半成品。 Y主要用来形成BX法所要求的第一个因子表,当该因子表形成后,就在半成品的基础上把接地支路及变压器非标准变比的影响加进去,形成完整的系统导纳矩阵。 其中Y 用来形成BX法第二个因子表,而Y将在整个迭带求解过程及线路潮流计算过程中发挥作用。 只考虑不接地支路构成导纳矩阵的程序框图如图3-2所示。 整个形成过程需要把不接地支路扫描一遍,对每条不接地支路作两方面的工作。 首先把阻抗求倒数并取负号后连同大节点号送到导纳矩阵非对角元素数组Yij(对应于Y)和Yij1(对应于Y )中形成非对角元素,然后把阻抗的倒数累加到该支路两端节点的自导纳上去[见式(3-2)]。 为了累加形成对角元素,在计算开始时应对数组Yii和Yii1清零[见图3-2中⑴、⑵框]。 ⑵框中NYsum为临时工作数组,定义为NYsum[300],在其中累计导纳矩阵各行非对角元素的个数,因此也需要预先清零。 由于两个导纳矩阵的结构是相同的,共用一个NYsum数组。 图3-2形成不接地支路的导纳矩阵 对不接地支路的扫描用一个循环语句来完成[图中⑶框]。 ⑷框把支路的有关数据送进中间工作单元,因为支路为变压器节点号可能为负,所以在这里对节点号取绝对值。 在⑸框中,把阻抗的倒数即支路导纳放到中间工作单元Gij、Bij中,而把支路电抗的倒数放到中间工作单元b_ij中。 ⑺框判断支路是否为变压器支路。 若为变压器支路,则导纳需除以变压器非标准变比后再取负号送到导纳矩阵非对角元素数组Yij和Yij1中;否则直接取负号送到导纳矩阵非对角元素数组中[⑻框]。 ⑼框向Yij和Yij1数组送列号。 这样就把支路阻抗数据变成了导纳矩阵的上三角部分。 在⑽~⑿框中根据式(3-2)累计有关节点的自导纳。 在⒀框中统计小节点号的不接地支路数目,这也就是导纳矩阵上三角部分每行非对角元素的个数。 至此,完成了一条不接地支路的处理;当循环由1做到Nb时形成了只考虑不接地支路的导纳矩阵。 ⒁~⒃框是由NYsum数组根据式(3-7)形成NYseq数组。 3.4追加接地支路的程序框图 图3-3追加对地支路框图 追加接地支路包括两部分内容,即追加对地电容支路和考虑变压器非标准变比的影响,其程序框图如图3-3所示。 整个计算过程需要对支路数据再进行依次扫描,扫描是由一个循环语句来控制[图中⑴框]。 在⑵框中把支路有关数据送入中间工作单元。 ⑵框根据节点号i、j的符号判断所取的支路是输电线路还是变压器支路。 当i、j中任一个为负时,为变压器支路,否则为输电线路。 当所取的支路为输电线路是,转入⑿~⒁框,向相应的节点累计自导纳部分。 当所取支路为变压器支路时,转入⑷~⑾框。 在⑷框中判断非标准变比设在支路的哪一侧。 如1章中所述,当i<0时,非标准变比就设在i侧,否则设在j侧。 由图3-1(b)及式(3-5)可知,在非标准变比侧自导纳应累计 yij,但在形成不接地支路的导纳矩阵时,该点自导纳累计了 ,因此需要再追加累计( -1) =(1- )Yij。 图中⑹~⑼框就是完成这些运算的。 在非标准变比侧自导纳应累计yij,在形成不接地支路的导纳矩阵时,该点自导纳同样累计了 ,因此需要再追加累计(K-1) =(1-K)Yij,在⑽框和⑾框中完成这些运算。 这样,顺次把支路数据扫描、处理一遍,就形成了描述网络的完整的导纳矩阵。 4.稀疏系数矩阵线性方程式的求解 4.1修正方程式的解法及计算公式 在P-Q分解法潮流计算的迭代过程中,需要反复求解修正方程式 △P/V=B △θV(4-1) △Q/V=B △V(4-2) 如前所述,这两个方程的系数矩阵(B 和B )在迭代过程中保持不变,只要求对不断变化的常数项(即误差项△P/V、△Q/V)求解出相应的修正量△θV及△V。 在这种情况下,可以先将系数矩阵进行三角分解,然后只要用分解出的三角矩阵(或因子表)对不同常数项进行前代及回代的运算,即可得到要求的修正量。 系数矩阵的三角分解可以利用递推公式求得,也可以利用高斯消去法求得,这两种方法在运算量及内存量上都是等效的。 本程序利用高斯消去法对系数矩阵进行三角分解并形成因子表的计算方法。 在P-Q分解法潮流程序中,为了在迭代过程中轮流求解式(4-1)及式(4-2)需要形成两个因子表,为此,可以把式(4-1)及式(4-2)统一为如下的形式: B△X=△I(4-3) 在P-θ迭代时,式中B即B ,△X为△θV,△I为△P/V;在Q-V迭代时,式中B为B ,△X为△V,△I为△Q/V。 以下简单归纳以下形成因子表及常数项进行前代及回代运算的有关公式,对于式(4-3)的系数矩阵进行B三角分解以后,可以得到以下形式的因子表: (4-4) 因子表中上三角部分元素组成了上三角矩阵U: U= 其中元素 Uij=Bij(i)(i 因子表中对角元素组成了对角矩阵D: D= 其中元素 Dii= (4-6) 式(4-5)中Bij(i)为因子表中上三角部分i行j列元素,其上标(i)表示此元素由原来系数矩阵元素Bij经过i次运算得来。 这i次运算中包括i-1次消去运算及一次规格化运算: Bij(k)=Bij(k-1)-Bik(k-1)Bkj(k)=Bij(k-1)-Bik(k-1)Ukj (k=1,2,…,i-1;j=k+1,k+2,…,n-1)(4-7) Bij(i)= (j=i+1,i+2,…,n-1)(4-8) 由于系数矩阵B为对称矩阵,在因子表中不需要保留其下三角部分。 在形成因子表的过程中需要用到下三角部分的元素Bik(k-1)(k Bik(k-1)=Bki(k)Bkk(k-1)=Uki (4-9) 利用式(4-6)~(4-9)就可以由对称系数矩阵B有规律的计算出因子表的所有元素。 对式(4-3)来说,当系数矩阵为对称矩阵时,其具体计算公式为: △Ij(i)=△Ij(j-1)-Uij△Ii(i-1)(j=2,3,…,n-1;i=1,2,…,j-1)(4-10) △Ij(j)=Djj△Ij(j-1)(4-11) 显然,顺次取因子表各元素参加一次运算就完成了前代过程的计算。 回代过程计算公式为: △xi=△Ii(i)- (j=n-1,n-2,…,1)(4-12) 4.2因子表形成程序框图 由于式(4-3)中的系数矩阵和导纳矩阵B具有相同的结构和稀疏性,因而在一般情况下因子表的上三角矩阵U也是稀疏的。 但由于在消去过程中可能会增加一些新的非零元素,所以U和B的结构不完全相同。 为了节约内存和提高计算速度,在因子表中将不保留零元素,而把上三角矩阵中非零元素紧凑的排列在一起,对每个上三角矩阵U的非零元素都附一个列下标,并对其中每一行都给出非零元素的数目。 以下框图(图4-1)来说明形成因子表的具体过程。 首先对框图中符号作一说明。 flag: 是一个标识变量,当flag=1时,该框图所表示的程序形成第一个因子表(与B 相对应)。 当flag=2时形成第二个因子表(与B 相对应)。 flag在该子程序调用前根据主程序 图4-1形成因子表框图 要求置1或置2,作为参数传个子程序。 n_pv: 为PV节点数组的计数变量,在形成第二个因子表时,用它可以判断系数矩阵中应该去掉哪些列和哪些行。 i: 为因子表正在形成的行号变量。 j: 为列下标变量。 n_u: 是因子表上三角矩阵元素计数变量。 i_above: 为消去行号计数变量(是被消行号计数变量,因子表按行消去过程中,依次取行上面的到行)。 i_pv: 为PV节点的节点号变量。 Btemp: 为临时变量。 count: 为临时计数变量。 形成因子表所需要的原始数据可由以下数组取得,这些数组的定义及内容详见1章和3章: Load: 负荷功率数组。 PVNode: PV节点数组。 Yij: 导纳矩阵非对角元素数组。 NYseq: 导纳矩阵各行非对角元素的首地址数组。 Yii: 导纳矩阵对角元素数组。 形成的因子表将放在以下数组中: NUsum: 存放因子表上三角矩阵各行非对角元素数。 D: 存放因子表的对角元素。 U: 存放因子表上三角矩阵元素,定义为结构体数组: structU_Type { doublevalue; intj; }U1[300],U2[300]; 其中存放该元素的数值,存放该元素的列下标。 最后,在框图中还有一个很重要的工作数组B[N],形成因子表的运算主要在这个数组中进行。 现在分别介绍框图中各个部分的作用。 图4-1中A、B、C三部分的作用是“传递”,通过这三部分的工作把系数矩阵中待消行(i行)的元素按其下标稀疏的排列在工作数组B中。 D框的作用是把工作数组B中的待消行元素按照式(4-7)进行消去运算。 最后通过E框的工作把数组B中的元素按式(4-8)进行格式化,并把数组B中非零元素搜集起来,紧密的排列到U数组中去。 这样,从第一行(i=1)做到第N-1行(i=N-1),就形成了完整的因子表。 整个计算过程是一行号为循环变量的。 以下将详细讨论图中各细框的工作情况。 首先介绍当flag=1时即形成第一因子表时的工作情况。 A部分把工作数组B从i+1到N-1单元全部充零,并把导纳矩阵对角元素的虚部送到B数组的第i个单元。 B部分把导纳矩阵非对角元素的虚部按其列下标送到B数组的相应单元中去。 在flag=1时,程序不执行C部分中的运算,直接转入到D部分。 至此,系数矩阵的第i行元素已稀疏的按其列下标排列在B数组中,以下将在D部分中按行对工作数组(即待消行i)中的元素进行消去运算。 一般的说,在形成因子表第i行各元素时,工作数组应该与i-1行以前已形成的各行因子表元素进行消去运算。 因此,在D部分中安排了消去行号i_above=1到i_above=i-1的循环过程。 该框开始时,对n_u赋值1,为顺序取用因子表上三角矩阵已形成各元素作好准备。 由于因子表及系数矩阵B都具有稀疏的特性,消去过程比较复杂,以下用图所示的例子来说明。 图4-2中因子表的第一行及第二行已经形成。 为了清楚起见,在图4-1中已将这两行因子表元素展开排列,实际上在数组U中它们是密集排列的(见图4-3)。 因子表第一行 因子表第二行 待消的第三行 图4-2形成因子表时的消去过程 在图4-2所示的例子中,第一行有三个非对角元素,第二行有三个非对角元素。 由式(4-5)可知,图中各元素的具体意义为 U13=B13 (1),U15=B15 (1),…,U28=B28 (2) 图4-3因子表上三角矩阵的存放形式 工作数组中B33为系数矩阵第三行的对角元素,B36及B38为非对角元素。 首先讨论因子表第一行与工作数组的消去过程。 根据式(4-7),可以写出 B3j (1)=B3j-B31 B1j (1) 式中: 下标应该由待消行号3开始,即j=3,4…。 因为系数矩阵是对称矩阵,如上所述,在因子表中不必保留下三角部分,因此上式中B31还必须利用式(4-9)求出: B31= B13 (1)= U13 得到后,就可以进行消去运算。 先从,即对角元素起 B33 (1)=B33-B31 B13 (1)=B33-B31 U13(4-13) 当j=4时: B34 (1)=B34-B31 B14 (1)=B34-B31 U14(4-14) 但在图4-2所示的例子中是零元素,因此上式变为 B34 (1)=B34 这样,式(4-14)所表示的运算对工作数组(即待消行)中的元素没有任何影响。 因此,在第一行与第三行进行消去运算时,只需要从因子表第一行中列下标为3的元素开始,顺次取以后的元素(见图4-3),按照其列下标与工作数组中相应的元素进行消去运算。 在本例中,除了按式(4-13)进行计算以外,还应进行以下两次运算: B35 (1)=B35-B31 U15 (4-15) B36 (1)=B36-B31 U16 式(4-15)中,B35为零(见图4-2),即系数矩阵中B35为零元素,但与第一行进行消去运算后变成了非零元素,因而在工作数组(即待消行)中出现了一个注入元素。 现在讨论第二行对第三行的消去运算。 根据式(4-7),消去过程的计算公式为 B3j (2)=B3k (1)-B32 (1) B2j (2)(4-16) 式中: B32 (1)= B23 (2)= U23 由图4-2可知,在该例中U23=0,因此B32 (1)也等于零,这样式(4-16)变为 B3j (2)=B3j (1) 也就是说,由于U23是零元素,因此第二行不必对第三行进行消去运算。 在一般情况下,当工作数组中待消行号为i,而第i (i 行对i行的消去运算过程。 在图4-1的D部分中,为了判断i 行是否需要对工作数组中的待消行进行消去运算,必须对因子表i 行中元素逐个检查其列下标是否等于待消行的行号i[D部分中⑴]。 当查到列下标等于i的元素时,首先按照式(4-9)求出Bi,i (i -1)[D部分中⑵],然后用该行剩余元素(包括列下标为i的
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 潮流 计算 软件设计
![提示](https://static.bingdoc.com/images/bang_tan.gif)