matlab潮流计算Word文件下载.docx
《matlab潮流计算Word文件下载.docx》由会员分享，可在线阅读，更多相关《matlab潮流计算Word文件下载.docx（20页珍藏版）》请在冰豆网上搜索。
%------------------------------------------------------------------------
clearall;
formatlong;
n=input（'
请输入节点数:
nodes='
）;
nl=input（'
请输入支路数:
lines='
isb=input（'
请输入平衡母线节点号:
balance='
pr=input（'
请输入误差精度:
precision='
B1=input（'
请输入由各支路参数形成的矩阵:
B1='
B2=input（'
请输入各节点参数形成的矩阵:
B2='
Y=zeros（n）;
e=zeros（1,n）;
f=zeros（1,n）;
V=zeros（1,n）;
sida=zeros（1,n）;
S1=zeros（nl）;
%------------------------------------------------------------------
fori=1:
nl%支路数
ifB1（i,6）==0%左节点处于1侧
p=B1（i,1）;
q=B1（i,2）;
else%左节点处于K侧
p=B1（i,2）;
q=B1（i,1）;
end
Y（p,q）=Y（p,q）-1./（B1（i,3）*B1（i,5））;
%非对角元
Y（q,p）=Y（p,q）;
Y（q,q）=Y（q,q）+1./（B1（i,3）*B1（i,5）^2）+B1（i,4）;
%对角元K侧
Y（p,p）=Y（p,p）+1./B1（i,3）+B1（i,4）;
%对角元1侧
%求导纳矩阵
disp（'
导纳矩阵Y='
disp（Y）
%-------------------------------------------------------------------
G=real（Y）;
B=imag（Y）;
%分解出导纳阵的实部和虚部
n%给定各节点初始电压的实部和虚部
e（i）=real（B2（i,3））;
f（i）=imag（B2（i,3））;
V（i）=B2（i,4）;
%PV节点电压给定模值
n%给定各节点注入功率
S（i）=B2（i,1）-B2（i,2）;
%i节点注入功率SG-SL
B（i,i）=B（i,i）+B2（i,5）;
%i节点无功补偿量
%---------------------------------------------------------------------
P=real（S）;
Q=imag（S）;
%分解出各节点注入的有功和无功功率
ICT1=0;
IT2=1;
N0=2*n;
N=N0+1;
a=0;
%迭代次数ICT1、a；
不满足收敛要求的节点数IT2
whileIT2~=0%N0=2*n雅可比矩阵的阶数；
N=N0+1扩展列
IT2=0;
a=a+1;
n
ifi~=isb%非平衡节点
C（i）=0;
D（i）=0;
forj1=1:
C（i）=C（i）+G（i,j1）*e（j1）-B（i,j1）*f（j1）;
%Σ（Gij*ej-Bij*fj）
D（i）=D（i）+G（i,j1）*f（j1）+B（i,j1）*e（j1）;
%Σ（Gij*fj+Bij*ej）
P1=C（i）*e（i）+f（i）*D（i）;
%节点功率P计算eiΣ（Gij*ej-Bij*fj）+fiΣ（Gij*fj+Bij*ej）
Q1=C（i）*f（i）-e（i）*D（i）;
%节点功率Q计算fiΣ（Gij*ej-Bij*fj）-eiΣ（Gij*fj+Bij*ej）
%求i节点有功和无功功率P'
Q'
的计算值
V2=e（i）^2+f（i）^2;
%电压模平方
%以下针对非PV节点来求取功率差及Jacobi矩阵元素-----------------------------
ifB2（i,6）~=3%非PV节点
DP=P（i）-P1;
%节点有功功率差
DQ=Q（i）-Q1;
%节点无功功率差
%以上为除平衡节点外其它节点的功率计算--------------------------------------
%求取Jacobi矩阵----------------------------------------------------------
ifj1~=isb&
j1~=i%非平衡节点&
非对角元
X1=-G（i,j1）*e（i）-B（i,j1）*f（i）;
%dP/de=-dQ/df
X2=B（i,j1）*e（i）-G（i,j1）*f（i）;
%dP/df=dQ/de
X3=X2;
%X2=dp/dfX3=dQ/de
X4=-X1;
%X1=dP/deX4=dQ/df
p=2*i-1;
q=2*j1-1;
J（p,q）=X3;
J（p,N）=DQ;
m=p+1;
%X3=dQ/deJ（p,N）=DQ节点无功功率差
J（m,q）=X1;
J（m,N）=DP;
q=q+1;
%X1=dP/deJ（m,N）=DP节点有功功率差
J（p,q）=X4;
J（m,q）=X2;
%X4=dQ/dfX2=dp/df
elseifj1==i&
j1~=isb%非平衡节点&
对角元
X1=-C（i）-G（i,i）*e（i）-B（i,i）*f（i）;
%dP/de
X2=-D（i）+B（i,i）*e（i）-G（i,i）*f（i）;
%dP/df
X3=D（i）+B（i,i）*e（i）-G（i,i）*f（i）;
%dQ/de
X4=-C（i）+G（i,i）*e（i）+B（i,i）*f（i）;
%dQ/df
J（p,q）=X3;
%扩展列△Q
m=p+1;
%扩展列△P
else
%下面是针对PV节点来求取Jacobi矩阵的元素-----------------------------------------
%PV节点有功误差
DV=V（i）^2-V2;
%PV节点电压误差
%dP/de
%dP/df
X5=0;
X6=0;
J（p,q）=X5;
J（p,N）=DV;
%PV节点电压误差
J（p,q）=X6;
%PV节点有功误差
X5=-2*e（i）;
X6=-2*f（i）;
%以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差---------------------------------------
fork=3:
N0%N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;
N1=N;
%N=N0+1即N=2*n+1扩展列△P、△Q或△U
fork2=k1:
N1%从k+1列的Jacobi元素到扩展列的△P、△Q或△U
J（k,k2）=J（k,k2）./J（k,k）;
%用K行K列对角元素去除K行K列后的非对角元素进行规格化
J（k,k）=1;
%对角元规格化K行K列对角元素赋1
%回代运算-------------------------------------------------------------------
ifk~=3%不是第三行k>
3
k4=k-1;
fork3=3:
k4%用k3行从第三行开始到当前行的前一行k4行消去
N1%k3行后各行上三角元素
J（k3,k2）=J（k3,k2）-J（k3,k）*J（k,k2）;
%消去运算（当前行k列元素消为0）
end%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J（k3,k）=0;
%当前行第k列元素已消为0
ifk==N0%若已到最后一行
break;
%前代运算------------------------------------------------------------
fork3=k1:
N0%从k+1行到2*n最后一行
N1%从k+1列到扩展列消去k+1行后各行下三角元素
%消去运算
else%是第三行k=3
%第三行k=3的前代运算----------------------------------------------------
N0%从第四行到2n行（最后一行）
N1%从第四列到2n+1列（即扩展列）
%消去运算（当前行3列元素消为0）
end%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
%当前行第3列元素已消为0
%上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵
%-----------------------------------------------------------------------------------
fork=3:
2:
N0-1
L=（k+1）./2;
e（L）=e（L）-J（k,N）;
升级会员 