Matlab及C语言在潮流计算运用Word文件下载.docx
《Matlab及C语言在潮流计算运用Word文件下载.docx》由会员分享，可在线阅读，更多相关《Matlab及C语言在潮流计算运用Word文件下载.docx（29页珍藏版）》请在冰点文库上搜索。
）Δx°
+[f'
'
）（Δx°
）2]/2+…
若所取的|Δx°
|足够小，则含（Δx°
）²
的项及其余的一切高阶项均可略去，并使其等于零，即：
f（x¹
）≈f（x°
=0
故得
Δx°
=-f（x°
）/f'
从而
x¹
=x°
-f（x°
可见，只要f'
）≠0，即可根据上式求出第一次的修正估值x¹
，若恰巧有f（x¹
）=0，则方程的真实解即为x*=x¹
若f（x¹
）≠0，则用上述方法由x¹
再确定第二次的修正估值x²
如此反复叠代下去，直到求得真实解x*为止。
二，节点电压用极坐标牛顿—拉夫逊法潮流计算
节点的功率方程写成
其中
式中，
是两节点电压的相角差。
把节点功率表示为节点电压的幅值和相角的函数。
在有n个节点的系统中，假定第1～m号节点为节点，第m+1～n-1号节点为节点，第n号节点为平衡节点。
Un和θn是给定的，PV节点的电压幅值Um+1～Un-1也是给定的。
因此，只剩下n-1个节点的电压相角θ1，θ2，…，θn-1和m个节点的电压幅值U1，U2，…，Um是未知量。
实际上，对于每一个PQ节点或每一个PV节点都可以列写一个有功功率不平衡量方程式，一共包含了n-1+m个方程式，正好同未知数的数目相同，而比直角坐标形式的方程少了n-1-m个。
对于方程式可以写出修正方程式如下：
H是（n-1）×
（n-1）阶方阵，其元素为；
N是（n-1）×
m阶矩阵，其元素为；
K是m×
（n-1）阶矩阵，其元素为；
L是m×
m阶方阵，其元素为。
在这里把节点不平衡功率对节点电压幅值的偏导数都乘以该节点电压，相应地把节点电压的修正量都除以该节点的电压幅值，这样，雅可比矩阵的表达式就具有比较整齐的形式。
三，matlab及c语言运算过程的具体实例
例：
网络接线如图，各支路阻抗和各节点功率均以标幺值标于图中，其中节点2连接的实际是发额定功率的发电厂，设节点1的电压保持为1.06，用牛顿-拉夫逊法计算系统中的潮流分布。
0.02+j0.06
0.45+j0.15
0.08+j0.24
0.01+j0.03
0.4+j0.05
0.04+j0.12
0.06+j0.18
0.2+j0.2
G
5
2
1
3
4
-0.02+j0.16
四，程序运行结果：
Matlab程序运行结果如下：
Y=
3.7500-11.2500i-2.5000+7.5000i0-1.2500+3.7500i0
-2.5000+7.5000i10.8333-32.5000i-1.6667+5.0000i-1.6667+5.0000i-5.0000+15.0000i
0-1.6667+5.0000i12.9167-38.7500i-10.0000+30.0000i-1.2500+3.7500i
-1.2500+3.7500i-1.6667+5.0000i-10.0000+30.0000i12.9167-38.7500i0
0-5.0000+15.0000i-1.2500+3.7500i06.2500-18.7500i
count=
4
U=
1.05551.05721.03191.03311.0600
Angle=
0.04170.0055-0.0261-0.02360
Error_Angle_U=
1.0e-008*
-0.5980
-0.0366
0.1334
0.1379
-0.8207
-0.4688
-0.5584
-0.6216
S=
0.6000-0.1000i0.2000+0.2000i-0.4500-0.1500i-0.4000-0.0500i0.0677+0.1531i
Line_S=
00.2905-0.1333i00.2991-0.0163i0
-0.2956+0.1229i00.2185+0.0764i0.2028+0.0746i0.0764-0.0760i
0-0.2109-0.0812i0-0.0918-0.0132i-0.1389-0.0789i
-0.2931-0.0031i-0.1960-0.0786i0.0919+0.0130i00
0-0.0771+0.0758i0.1448+0.0773i00
Line_I=
00.2697+0.1377i00.2824+0.0273i0
-0.2790-0.1178i00.2071-0.0711i0.1922-0.0695i0.0719+0.0723i
0-0.2022+0.0840i0-0.0886+0.0151i-0.1326+0.0799i
-0.2836+0.0097i-0.1879+0.0806i0.0886-0.0146i00
0-0.0727-0.0715i0.1366-0.0729i00
C语言运行结果如下：
节点导纳矩阵：
3.750000+j-11.250000-2.500000+j7.500000-0.000000+j-0.000000-1.250000+j3.750000-0.000000+j-0.000000
-2.500000+j7.50000010.833333+j-32.500000-1.666667+j5.000000-1.666667+j5.000000-5.000000+j15.000000
-0.000000+j-0.000000-1.666667+j5.00000012.916667+j-38.750000-10.000000+j30.000000-1.250000+j3.750000
-1.250000+j3.750000-1.666667+j5.000000-10.000000+j30.00000012.916667+j-38.750000-0.000000+j-0.000000
-0.000000+j-0.000000-5.000000+j15.000000-1.250000+j3.750000-0.000000+j-0.0000006.250000+j-18.750000
迭代次数：
count=4
节点相角分别为:
θ=0.0417100.005480-0.026114-0.0235940.000000
节点电压分别为:
U=1.0555251.0571781.0319311.0330721.060000
节点功率为:
S=
0.600000+j-0.1000000.200000+j0.200000-0.450000+j-0.150000-0.400000+j-0.0500000.067700+j0.153099
请按任意键继续...
经过比较matlab和c语言运行的结果是一样的！
最终的潮流分布如下图：
0.0764-0.0760i
--0.1389-0.0789i
-0.0918-0.0132i
0.2905-0.1333i
0.2028+0.0746i
0.2956-0.1229i
-0.06+0.1i
0.2931+0.0031i
0.0677+0.1531i
0.2185+0.0764i
0.2991-0.0163i
0.2109+0.0812i
-0.092-0.013i
0.0771-0.0758i
-0.1448-0.0773i
0.1960+0.0786i
五，小结
在计算机潮流计算中，使用牛顿-拉夫逊法进行潮流计算，非常快捷。
在程序编写过程中，对公式的由来应该非常清楚，不然在编写过程中易出错，比如for循环的次数，一些细节问题都应该反复琢磨。
在程序调试过程中，应该分步调试，不要一开始就运行完整的程序。
对程序中出现的问题应该从运行结果中分析问题的所在，有什么问题解决什么问题，抓住问题的要害。
在潮流计算中运用matlab比运用c语言要方便好多，matlab里有大量的数学函数直接调用。
用c语言写的话每一个函数都要自己编写，写起来很痛苦,其中求逆矩阵就用了一百多行。
最后想说的是，想学没有什么不可以。
六，程序代码
Matlab程序：
%Writer:
ChenZhou
%Time:
2009年11月9日-----2009年11月15日
%编译环境：
matlabR2007b
clc
%说明：
为了使节点按照先PQ，再PV节点，最后平衡节点的次序编号，以便与公式对照，节点1与节点5对调。
%节点阻抗矩阵
Z=[0,0.04+0.12i,0,0.08+0.24i,0;
0.04+0.12i,0,0.06+0.18i,0.06+0.18i,0.02+0.06i;
0,0.06+0.18i,0,0.01+0.03i,0.08+0.24i;
0.08+0.24i,0.06+0.18i,0.01+0.03i,0,0;
0,0.02+0.06i,0.08+0.24i,0,0];
%求互导纳
fori=1:
forj=1:
ifZ（i,j）==0
y（i,j）=0;
else
y（i,j）=1/Z（i,j）;
end
end
%求导纳
ifi~=j
Y（i,j）=-y（i,j）;
fork=1:
Y（i,j）=sum（y（i,:
））;
end
%导纳矩阵
G=real（Y）;
B=imag（Y）;
%赋初值
U=[1,1,1,1,1.06];
Angle=[0,0,0,0,0];
%节点功率
S=[0.6-0.1i,0.2+0.2i,-0.45-0.15i,-0.4-0.05i,0];
P=real（S）;
Q=imag（S）;
Error_Power=ones（8,1）;
count=0;
%迭代次数
whilemax（Error_Power）>
1e-5
fori=1:
H（i,i）=0;
N（i,i）=0;
M（i,i）=0;
L（i,i）=0;
Error_P（i）=0;
Error_Q（i）=0;
Error_P（i）=Error_P（i）-U（i）*U（j）*（G（i,j）*cos（Angle（i）-Angle（j））+B（i,j）*sin（Angle（i）-Angle（j）））;
Error_Q（i）=Error_Q（i）-U（i）*U（j）*（G（i,j）*sin（Angle（i）-Angle（j））-B（i,j）*cos（Angle（i）-Angle（j）））;
Error_P（i）=Error_P（i）+P（i）;
Error_Q（i）=Error_Q（i）+Q（i）;
Error_Power=[Error_P,Error_Q]'
;
%求雅可比矩阵
%当i~=j时H,N,M,L如下
H（i,j）=-U（i）*U（j）*（G（i,j）*sin（Angle（i）-Angle（j））-B（i,j）*cos（Angle（i）-Angle（j）））;
N（i,j）=-U（i）*U（j）*（G（i,j）*cos（Angle（i）-Angle（j））+B（i,j）*sin（Angle（i）-Angle（j）））;
L（i,j）=-U（i）*U（j）*（G（i,j）*sin（Angle（i）-Angle（j））-B（i,j）*cos（Angle（i）-Angle（j）））;
M（i,j）=U（i）*U（j）*（G（i,j）*cos（Angle（i）-Angle（j））+B（i,j）*sin（Angle（i）-Angle（j）））;
%当i=j时H，N，M，L如下
H（i,i）=H（i,i）+U（i）*U（j）*（G（i,j）*sin（Angle（i）-Angle（j））-B（i,j）*cos（Angle（i）-Angle（j）））;
N（i,i）=N（i,i）-U（i）*U（j）*（G（i,j）*cos（Angle（i）-Angle（j））+B（i,j）*sin（Angle（i）-Angle（j）））;
M（i,i）=M（i,i）-U（i）*U（j）*（G（i,j）*cos（Angle（i）-Angle（j））+B（i,j）*sin（Angle（i）-Angle（j）））;
L（i,i）=L（i,i）-U（i）*U（j）*（G（i,j）*sin（Angle（i）-Angle（j））-B（i,j）*cos（Angle（i）-Angle（j）））;
N（i,i）=N（i,i）-2*（U（i））^2*G（i,i）;
L（i,i）=L（i,i）+2*（U（i））^2*B（i,i）;
J=[H,N;
M,L];
Error_Angle_U=-（inv（J）*Error_Power）;
Error_Angle=zeros（1,4）;
Error_U=zeros（1,4）;
Error_Angle（i）=Error_Angle_U（i）;
Error_U（i）=Error_Angle_U（i+4）*U（i）;
4
Angle（i）=Angle（i）+Error_Angle（i）;
U（i）=U（i）+Error_U（i）;
count=count+1;
%count
%平衡节点的注入功率
P（5）=U（5）*U（j）*（G（5,j）*cos（Angle（5）-Angle（j））+B（5,j）*sin（Angle（5）-Angle（j）））+P（5）;
Q（5）=U（5）*U（j）*（G（5,j）*sin（Angle（5）-Angle（j））-B（5,j）*cos（Angle（5）-Angle（j）））+Q（5）;
S（5）=P（5）+Q（5）*sqrt（-1）;
%节点电压用复数表示
Complex_U（i）=U（i）*exp（Angle（i）*sqrt（-1））;
%元件端功率
Line_S（i,j）=conj（Complex_U（i））*（conj（Complex_U（i））-conj（Complex_U（j）））*conj（y（i,j））;
%元件端电流
Line_I（i,j）=conj（Line_S（i,j））/conj（Complex_U（i））;
Y,count,U,Angle,Error_Angle_U,S,Line_S,Line_I
C语言程序：
此程序一共三个文件，分别为InverseMatrix.h，Functions.h，ChaoLiuJiSuan.c
主函数：
//Writer:
//Time:
//编译环境：
visualstdio2008
//说明：
为了使节点按照先PQ，再PV节点，最后平衡节点的次序编号，以便与公式对照，节点与节点对调。
#include<
stdio.h>
math.h>
#include"
InverseMatrix.h"
Functions.h"
voidmain（）
{
unsignedinti,j,count=0;
Z_to_Y（）;
printf（"
\n"
）;
for（i=0;
i<
5;
i++）
{
for（j=0;
j<
j++）
{
printf（"
%f+j%f"
Y_Re[i][j],Y_Im[i][j]）;
}
}
One_Times（）;
while（（max_array（Error_P）>
0.00001）||（max_array（Error_Q）>
0.00001））
One_Times（）;
count++;
count=%d"
count+1）;
θ="
for（i=0;
i++）
%f"
Angle[i]）;
}printf（"
U="
U[i]）;
Calculate_S（）;
\n节点功率为:
S=\n"
S_P[i],S_Q[i]）;
}
求逆矩阵的头文件（InverseMatrix.h文件）
intmax=8;
//==============雅可比矩阵=====================================
doubleJ[8][8]={0};
intswitchnum[2][8];
voidswitchrow（intr,intk）{/*将矩阵中的R,K行交换*/
inti;
doubletemp;
for（i=0;
max;
{
temp=J[k][i];
J[k][i]=J[r][i];
J[r][i]=temp;
voidswitchcol（intr,intk）{/*将矩阵中的R,K行交换*/
temp=J[i][k];
J[i][k]=J[i][r];
J[i][r]=temp;
intInverse_Matrix（）
inti,j,r,k;
for（k=0;
k<
max;
k++）
temp=0;
r=k;
for（j=k;
=max;
j++）/*找主元*/
if（temp<
fabs（J[r][j]））
temp=fabs（J[r][j]）;
r=j;
if（0==temp）/*非奇异矩阵,退出*/
return（0）;
if（r!
=k）
switchrow（r,k）;
switchnum[0][k]=r;
switchnum[1][k]=k;
switchnum[0][k]=0;
switchnum[1][k]=0;
J[k][k]=1/J[k][k];
for（j=0;
if（j!
J[k][j]=-J[k][j]*J[k][k];
if（i!
=k&
&
j!
J[i][j]=J[i][j]+J[i][k]*J[k][j];
}
J[i][k]=J[i][k]*J[k][k];
/*恢复成逆矩阵*/
for（k=max-1;
k>
=0;
k--）
if（switchnum[0][k]!
=0）
switchcol（switchnum[0][k],switchnum[1][k]）;
求解过程中的一些功能函数（Functions.h文件）
//=================================计算导纳矩阵=========================================
//阻抗的实部
doubleZ_Re[5][5]={{0,0.04,0,0.08,0},{0.04,0,0.06,0.06,0.02},{0,0.06,0,0.01,0.08},{0.08,0.06,0.01,0,0},{0,0.02,0.08,0,0}};
//阻抗的虚部
doubl
升级会员 