matlab潮流计算备课讲稿Word下载.docx

matlab潮流计算备课讲稿Word下载.docx

《matlab潮流计算备课讲稿Word下载.docx》由会员分享，可在线阅读，更多相关《matlab潮流计算备课讲稿Word下载.docx（20页珍藏版）》请在冰豆网上搜索。

2为PQ节点；

3为PV节点；

%------------------------------------------------------------------------

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;