电力系统潮流计算C语言程序及说明.docx
《电力系统潮流计算C语言程序及说明.docx》由会员分享，可在线阅读，更多相关《电力系统潮流计算C语言程序及说明.docx（32页珍藏版）》请在冰豆网上搜索。
电力系统潮流计算C语言程序及说明
实验目的
根据所给的电力系统，编制潮流计算程序，通过计算机进行调试，最后完成一个切实可行的电力系统计算应用程序。
通过自己设计电力系统计算程序使同学对电力系统分析有进一步理解，同时加强计算机实际应用能力的训练。
程序计算原理
1、概述
应用计算机进行电力系统计算，首先要掌握电力系统相应计算的数学模型；其次是运用合理的计算方法；第三则是选择合适的计算机语言编制计算程序。
建立电力系统计算的相关数学模型，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式。
该数学模型的建立往往要突出问题的主要方面，即考虑影响问题的主要因素，而忽略一些次要因素，使数学模型既能正确地反映实际问题，又使计算不过于复杂。
运用合理的计算方法，就是要求所选用的计算方法能快速准确地得出正确结果，同时还应要求在解算过程中占用内存少，以利提高计算机的解题规模。
选择合适的语言编写程序，就是首先确定用什么计算机语言来编制程序；其次是作出计算的流程图；第三根据流程图用选择的语言编写计算程序。
然后上机调试，直到语法上无错误。
本程序采用C语言进行编程。
所编制的程序难免存在逻辑错误，因此先用一个已知结果的系统作为例题进行计算。
用程序计算的结果和已知结果相比较，如果结果相差甚远就要逐步分析程序的计算步骤，查出问题的出处；如果结果比较接近，则逐步分析误差来源；直到结果正确为止。
2、电力系统潮流计算的程序算法
潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压（幅值及相角）、网络中的功率分布及功率损耗等。
目前计算机潮流计算的方法主要有牛顿-拉夫逊算法和PQ分解法。
牛顿-拉夫逊算法是数学上求解非线形方程组的有效方法，具有较好的收敛性，曾经是潮流计算中应用比较普遍的方法。
PQ快速分解法是从牛顿-拉夫逊算法演变而来的，是将纯数学的牛顿-拉夫逊算法与电力系统具体特点相结合并进行简化与改进而得出的。
PQ快速分解法比牛顿-拉夫逊算法大大提高了计算速度和节省了内存，故而本程序以PQ快速分解法进行潮流计算。
1）形成节点导纳矩阵
（1）自导纳的形成
对节点i其自导纳Yii是节点i以外的所有节点都接地时节点i对地的总导纳。
显然，Yii应等于与节点i相接的各支路导纳之和，即
式中，yi0为节点i与零电位节点之间的支路导纳；yij为节点i与节点j之间的支路导纳。
（2）互导纳的形成
对节点i与节点k之间的互导纳是节点i、k之间的支路导纳的负值，即
不难理解
。
若节点i和k没有支路直接相连时，便有Yik=0
（3）含变压器支路的处理
若节点p、q间接有变压器，如下图所示，则可作出其∏型等值电路为：
图1变压器∏型等值电路
则p、q的自导纳和节点间的互导纳分别为
2）计算不平衡功率△P、△Q并形成修正方程式
对每一个PQ节点或每一个PV节点都可以根据下列公式计算出有功功率增量△P
而对于每一个PQ节点还可以根据下面的公式计算出无功功率增量△Q
在有功功率增量和无功功率增量不满足如下约束条件时
利用PQ分解法则可以形成如下修正方程
3）利用因子表法求解修正方程
在电网计算中经常遇到这样的问题，对方程组需要反复多次求解，而每次求解仅改变常数项F，系数矩阵保持不变。
按照一般的高斯消去法，对每一改变的常数项，形成包括常数项及系数矩阵在内的增广矩阵，然后消去回代求出其解。
可以看出，每次对增广矩阵中A矩阵元素的消元都是重复的，为了避免这种重复，我们把对相同的系数矩阵重复进行的消去与对不同的常数项进行的消去分开进行，因此对系数矩阵的消去只需进行一次，并在消去的过程中将对常数项进行消去运算的运算因子保存下来，形成所谓因子表，这就是因子表法。
因为因子表记录了高斯消去法对常数项进行消去的全部信息，利用它便可对不同常数项进行消去，形成上三角矩阵，最后求出全部未知数。
在使用PQ分解法时，其系数矩阵是在迭代过程中保持不变的，所以为了节省内存和缩短运算时间我们采取了因子表法。
同时由于电网的节点导纳矩阵是稀疏阵和对称阵，于是我们可以采取只保存系数矩阵的上三角阵来使运算更为简化。
若线性方程组一般形式如下：
其中
称为系数矩阵，
称为未知数向量，
称为常数项向量。
将矩阵A的元素进行如下处理：
得到因子表
其中
；
再利用因子表进行前代过程，求出每次迭代后的常数项。
其前代公式是：
求得向量
；
再由因子表与前代得到的向量F，得到方程组
求解出此方程即可得到线性方程组的解向量
。
4）多次迭代最终求得V和
以及全线路功率
利用上面所介绍的方法求解修正方程组
可以求得
和
。
再利用
求得每次迭代后的结果。
多次迭代当其满足约束条件
和
时，迭代结束。
迭代结束后即可得到各节点的V和
，再根据V、
来计算PV节点的无功功率Q和平衡节点的功率以及网络中的功率分布。
PV节点及平衡节点无功功率计算公式为：
平衡节点有功功率计算公式为：
以下图所标示的正方向，输电线路功率的计算公式如下：
图2支路功率计算
对其进行实部虚部进行分解可得P、Q计算公式为：
程序及说明
1、主要变量说明
1）结构体类型说明
（1）节点功率结构体
structNodetype
{
floatP,Q;
};
其中，P为节点的有功功率，Q为无功功率。
节点功率不区分负荷功率和发电机功率，其值为本节点连接的各支路输入功率及节点所接负荷、发电机功率之和，且规定功率流入节点为正，流出为负。
详细说明参见下一章“算例及结果”的第二节“源数据格式说明”。
（2）线路参数结构体
structLinetype
{
floatG,B,B0,k;
};
其中，G、B为线路的导纳和容纳；B0为线路的考虑变压器Π型等值电路后的对地充电容纳的一半Bc/2；k为折算到标准变压器支路后的变压器变比。
详细说明参见下一章“算例及结果”的第二节“源数据格式说明”。
2）变量说明
表2程序主要变量说明
主要变量
类型
含义
Node
int
系统总节点数
NP
int
PV+PQ节点数，即非平衡节点数
NQ
int
PQ节点数
Num
int*
原始节点编号与程序表示编号映射数组
No
structNodetype*
节点功率数组
V
float*
节点电压有效值数组
Dlta
float*
节点电压相角数组
Y
structLinetype**
线路参数矩阵
BP、BQ
float**
有功、无功简化雅克比矩阵B’、B”
count
unsignedint
PQ迭代次数
eP、eQ
const
有功、无功迭代精度控制
kp、kq
int
有功、无功迭代结束标志
dP、dQ
float*
有功、无功不平衡量数组
3、程序流程图
图3程序主流程图
图4迭代部分流程图
4、程序代码
/*FUNCTION:
POWERFLOW*/
/*WRITTENBY:
HUANG&YANG&TONG*/
/*LASTEDITED:
2008-11-24*/
#include
#include
/***宏定义***/
#defineeP0.00001
#defineeQ0.00001
#defineY_（i,j）（*（*（Y+i）-i+j））
#defineYij（*（Yi+j））
#defineYji（*（*（Y+j）-j+i））
#definePjiYji.G*cos（tmp）+Yji.B*sin（tmp）
#definePijYij.G*cos（tmp）+Yij.B*sin（tmp）
#defineQjiYji.G*sin（tmp）-Yji.B*cos（tmp）
#defineQijYij.G*sin（tmp）-Yij.B*cos（tmp）
/***结构体变量定义***/
structNodetype/*节点功率*/
{
floatP,Q;
};
structLinetype/*线路类型*/
{
floatG,B,B0,k;
};
/***子函数声明***/
voidin_node（）;/*读节点功率*/
voidin_line（）;/*读线路参数*/
voidB_Form（）;/*生成BP、BQ矩阵*/
voidfactor（）;/*求因子表*/
voidsolve（float**B,float*X,intN）;/*解方程组*/
voidPrtNode（）;/*打印节点参数*/
voidErrorMsg（intFlag）;/*错误提示信息*/
/***全局变量声明***/
intNode;/*节点数*/
int*Num;/*保存原始节点序号*/
intNP,NQ=0;/*PV+PQ、PQ节点数*/
structNodetype*No;/*节点数据*/
structLinetype**Y;/*线路参数*/
float**BP,**BQ;/*有功、无功简化雅克比矩阵B*/
float*V;/*节点电压有效值*/
float*Dlta;/*节点电压相角值*/
unsignedintcount=0;/*迭代计数*/
inti,j,k,m;/*通用下标值*/
floattmp;/*临时数据暂存*/
char*Type;/*节点类型*/
FILE*in,*out;/*输入、输出文件指针*/
/******************主函数******************/
/**↓****↓****↓**主函数**↓****↓****↓**/
intmain（void）
{
intkp=1,kq=1;/*P、Q精度判断：
1-不满足，0-满足精度要求*/
float*dP,*dPi,*dQ,*dQi;/*ΔP、ΔQ*/
floatDltai;
structLinetype*Yi;
structNodetype*Noi;
floattP,tQ;
if（（in=fopen（"Data.txt","r"））==NULL）ErrorMsg
（1）;
if（（out=fopen（"out.txt","w"））==NULL）ErrorMsg
（2）;
in_node（）;/*读取节点参数并统计PQ、PV节点数*/
in_line（）;/*读取线路参数并形成Y矩阵*/
B_Form（）;/*形成B（BP&BQ）矩阵*/
factor（）;/*求B因子式（仍保存于BP&BQ）*/
for（i=0;idP=（float*）malloc（sizeof（float）*NP）;/*ΔP*/
dQ=dP;/*ΔQ*//*ΔP、ΔQ不同时存在，故而可共用空间*/
loop:
/****开始迭代****/
if（kp==0&&kq==0）gotoloopEnd;
count++;/*迭代次数加一*/
if（count==65535）ErrorMsg（99）;/*不收敛，退出*/
kp=0;
for（i=0;i{
dPi=dP+i;
Yi=*（Y+i）-i;
Dltai=*（Dlta+i）;
*dPi=0;
for（j=0;j{
tmp=Dltai-*（Dlta+j）;/*tmp即δij*/
if（i>j）*dPi+=*（V+j）*（Pji）;
else*dPi+=*（V+j）*（Pij）;
}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
*dPi*=*（V+i）;
*dPi=（*（No+i））.P-*dPi;/*求得ΔPi*/
if（fabs（*dPi）>0x8fffffff）ErrorMsg（99）;/*不收敛，退出*/
if（fabs（*dPi）>eP）kp=1;/*有不满足精度的ΔP即令kp=1*/
*dPi/=*（V+i）;/*求得常数项ΔPi/Vi*/
}
if（kp==0）gotoloopQ;
solve（BP,dP,NP）;
for（i=0;iloopQ:
if（kp==0&&kq==0）gotoloopEnd;
kq=0;
for（i=0;i{
dQi=dQ+i;
Yi=*（Y+i）-i;
Dltai=*（Dlta+i）;
*dQi=0;
for（j=0;j{
tmp=Dltai-*（Dlta+j）;/*tmp即δij*/
if（i>j）*dQi+=*（V+j）*（Qji）;
else*dQi+=*（V+j）*（Qij）;
}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
*dQi*=*（V+i）;
*dQi=（*（No+i））.Q-*dQi;/*求得ΔQi*/
if（fabs（*dQi）>0x8fffffff）ErrorMsg（99）;/*不收敛，退出*/
if（fabs（*dQi）>eQ）kq=1;/*有不满足精度的ΔQ即令kq=1*/
*dQi/=*（V+i）;/*求得常数项ΔQi/Vi*/
}
if（kq==0）gotoloop;
solve（BQ,dQ,NQ）;
for（i=0;igotoloop;/*无功迭代，则必定需要下一轮回迭代判断*/
loopEnd:
/****迭代结束****/
free（dP）;/*释放内存空间*/
/****计算PV节点和平衡节点的无功功率Q****/
for（i=NQ;i{
Noi=No+i;
Yi=*（Y+i）-i;
Dltai=*（Dlta+i）;
for（j=0;j{
tmp=Dltai-*（Dlta+j）;/*tmp即δij*/
if（i>j）（*Noi）.Q+=*（V+j）*（Qji）;
else（*Noi）.Q+=*（V+j）*（Qij）;
}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
（*Noi）.Q*=*（V+i）;
}
/****计算平衡节点的有功功率P****/
i=NP;
Noi=No+i;
Dltai=*（Dlta+i）;
for（j=0;j{
tmp=Dltai-*（Dlta+j）;/*tmp即δij*/
（*Noi）.P+=*（V+j）*（Pji）;
}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
（*Noi）.P*=*（V+i）;
/****输出最终结果****/
fprintf（out,"\n\n【潮流计算结果（节点）】（迭代次数k=%3d）\n",count-1）;
PrtNode（）;
/****计算全部线路功率****/
fprintf（out,"\n\n【潮流计算结果（线路）】\n"）;
fprintf（out,"线路PQ\n"）;
for（k=0;k{
i=*（Num+k）;
Yi=*（Y+i）-i;
Dltai=*（Dlta+i）;
Noi=No+i;
for（m=0;m{
j=*（Num+m）;
if（j==i）continue;
tmp=Dltai-*（Dlta+j）;/*tmp即δij*/
if（j{
if（Yji.B==0）continue;/*若Bij=0，则节点i、j无直接联系*/
tP=*（V+j）*（Pji）;
tP=*（V+i）*Yji.G-tP;
tP*=*（V+i）;
tQ=-*（V+j）*（Qji）;
tQ-=*（V+i）*（Yji.B-Yji.B0/Yji.k）;
tQ*=*（V+i）;
}
else
{
if（Yij.B==0）continue;/*若Bij=0，则节点i、j无直接联系*/
tP=*（V+j）*（Pij）;
tP=*（V+i）*Yij.G-tP;
tP*=*（V+i）;
tQ=-*（V+j）*（Qij）;
tQ-=*（V+i）*（Yij.B-Yij.B0）;
tQ*=*（V+i）;
}
fprintf（out,"S[%d,%d]=（%10.6f,%10.6f）\n",k+1,m+1,-tP,-tQ）;
}
}
fclose（out）;
system（"cmd/cstartout.txt"）;
return（0）;
}
/**↑****↑****↑**主函数**↑****↑****↑**/
/******************主函数******************/
/****************子函数：
读节点数据****************/
voidin_node（）
{
structNodetype*Noi;/*临时变量*/
fscanf（in,"%d%d",&Node,&k）;/*读取节点数Node*/
NP=Node-1;/*PV+PQ节点数，即非平衡节点数目*/
Num=（int*）malloc（sizeof（int）*Node）;/*开Node个空间，每节点一个*/
V=（float*）malloc（sizeof（float）*Node）;/*电压*/
Dlta=（float*）malloc（sizeof（float）*Node）;/*电压相角*/
No=（structNodetype*）malloc（sizeof（structNodetype）*Node）;/*节点功率*/
j=1;
while（k!
=0）/*若k=0，表明节点数据读取完毕*/
{
switch（k）
{
case1:
k=NQ;NQ++;break;/*NQ统计PQ节点个数*/
case2:
k=NP-j;j++;break;/*从NP-1个空间倒着保存PV节点*/
case3:
k=NP;break;/*平衡节点*/
default:
ErrorMsg（3）;
}
Noi=No+k;
fscanf（in,"%d%f%f%f%f",&i,&（*Noi）.P,&（*Noi）.Q,V+k,Dlta+k）;
i--;/*节点编号减一，以和程序表达方式兼容*/
*（Num+i）=k;/*第i个Num元素中存放i节点在No中的下标*/
fscanf（in,"%d",&k）;/*读取节点类型*/
}
if（NQ+j!
=Node）ErrorMsg（4）;/*检验节点数据个数是否够Node个*/
fprintf（out,"【节点参数表】\n"）;
PrtNode（）;
fprintf（out,"总节点：
%d\nPQ节点：
%d\nPV节点：
%d\n",Node,NQ,NP-NQ）;
}
/************子函数：
读线路数据，并形成节点导纳矩阵Y************/
voidin_line（）
{
structLinetype*Yi;
floatR,X,k,B;
m=sizeof（structLinetype）;
Y=（structLinetype**）malloc（m*Node）;/*先开Node行,每一个节点一行*/
for（i=0;i{/*即以上三角存储Y矩阵*/
*（Y+i）=（structLinetype*）malloc（m*（Node-i））;
Yi=*（Y+i）-i;
for（j=i;j}
while（feof（in）==0）/*文件指针到文件末*/
{
fscanf（in,"%d%d%f%f%f%f",&i,&j,&R,&X,&k,&B）;
i--;j--;
i=*（Num+i）;/*转换节点号为该节点在程序中的储存编号*/
j=*（Num+j）;
（*（*（Y+i）））.B+=B;/*将对地充电导纳累加到自导纳*/
（*（*（Y+j）））.B+=B;
if（k!
=1.0）
{
X*=k;R=0;
tmp=（1-k）/X;
（*（*（Y+i）））.B+=tmp;/*将变压器的对地充电容纳累加到自导纳*/
（*（*（Y+j）））.B+=-（tmp/k）;
B=tmp;
k=-k;
}
if（i>j）{tmp=i;i=j;j=tmp;k=1/k;B*=k;}
Yi=*（Y+i）-i;/*以Yi代替*（Y+i）-i，简化表达式并避免重复计算*/
Yij.B0=B;/*保存ij0、ji0对地充电电容到Bij0*/
Yij.k=k;/*且有B0ji=B0ij/k*/
tmp=R*R+X*X;
R/=tmp;
X/=tmp;
Yij.G=-R;/*生成互导纳*/
Yij.B=X;
（*（*（Y+i）））.G+=R;/*将线路互导纳累加到自导纳*/
（*（*（Y+i）））.B+=-X;
（*（*（Y+j）））.G+=R;
（*（*（Y+j）））.B+=-X;
}
fclose（in）;
fprintf（out,"\n【节点导纳矩阵Y】\n"）;
for（k=0;k{
i=k;
i=*（Num+i）;/*查取第i节点在程序中存储序号*/
for（j=0;jfor（m=k;m{
j=*（Num+m）;/*查取第j节点在程序中存储序号*/
if（ielsefprintf（out,"（%10.6f,%10.6f）",Y_（j,i）.G,Y_（j,i）.B）;
}
fprintf（out,"\n"）;
}
}
/****************子函数：
生成BP、BQ矩阵****************/
voidB_Form（）
{
f
升级会员 