电力系统分析潮流计算的计算机算法

姓名： 学号： 班级：

一、

实验目的

掌握潮流计算的计算机算法。

熟悉MATLAB,并掌握MATLAB程序的基本调试方法。 二、

实验准备

根据课程内容，熟悉MATLAB软件的使用方法，自行学习MATLAB程序的基础语法，并根据所学知识编写潮流计算牛顿拉夫逊法（或PQ分解法） 的计算程序，用相应的算例在MATLAB上进行计算、调试和验证。 三、

实验要求

每人一组，在实验课时内，用MATLAB调试和修改运行程序，用算例计算输出潮流结果。 四、

实验程序

clear； %清空内存 n=input('请输入节点数:n=’）； n1=input(’请输入支路数：n1=’); isb=input（’请输入平衡节点号:isb=’）； pr=input（'请输入误差精度：pr=’）; B1=input(’请输入支路参数:B1=’); B2=input('请输入节点参数：B2=’)； X=input（'节点号和对地参数：X=’）; Y=zeros(n)； Times=1;

%一：创建节点导纳矩阵

for i=1:n1

if B1(i，6)==0 %不含变压器的支路 p=B1（i,1）; q=B1（i,2）；

Y(p，q)=Y（p，q)-1/B1（i,3）； Y（q,p)=Y（p，q)；

Y（p,p)=Y（p,p)+1/B1(i,3)+0。5＊B1（i,4）; Y(q,q)=Y(q,q）+1/B1(i，3)+0。5＊B1（i,4）; else ％含有变压器的支路 p=B1（i,1）； q=B1(i,2)；

Y（p,q）=Y(p,q）—1/（B1(i,3）*B1(i,5)）； Y（q,p)=Y（p,q)；

Y（p，p）=Y(p,p)+1/B1(i,3）;

Y(q，q）=Y（q，q）+1/(B1(i，5)^2*B1（i，3)）; end end Y;

％将OrgS、DetaS初始化 OrgS=zeros(2＊n-2,1); DetaS=zeros（2＊n-2，1)；

%二：创建OrgS，用于存储初始功率参数

h=0; j=0;

for i=1:n %对PQ节点的处理 if i~=isb&B2（i，6)==2

h=h+1；

for j=1：n OrgS(2＊h—1，1)=OrgS(2＊h—1,1)+real(B2（i,3))＊（real（Y（i，j))＊real(B2（j,3))—imag(Y（i，j））*I mag（B2（j，3)）)+imag(B2(i，3))＊（real（Y（i，j）)*imag（B2（j,3))+imag(Y(i,j）)＊real(B2（j，3）)）; OrgS(2*h，1)=OrgS(2*h,1）+imag(B2(i，3))＊(real（Y(i，j）)＊real（B2(j，3))-imag(Y（i,j））*imag（B2（j，3)））—real(B2（i，3））＊（real（Y（i，j））*imag(B2(j,3））+imag（Y（i，j））*real（B2（j,3)));

end end end

%三:对PV节点的处理，注意这时不可再将h初始化为0

for i=1:n if i~=isb&B2（i,6)==3 h=h+1;

for j=1：n OrgS（2*h—1，1）=OrgS(2*h—1，1)+real

（B2(i,3）)*(real（Y(i,j））*real(B2（j,3))—imag(Y（i，j）)＊imag（B2(j，3)))+imag(B2（i，3））*（real（Y（i，j）)*imag（B2(j，3))+imag(Y（i,j)）*real(B2（j，3））); OrgS(2＊h，1）=OrgS（2＊h,1）+imag（B2(i，3））*(real(Y（i，j））*real（B2(j,3）)-imag（Y(i,j)）＊imag(B2(j，3）))—real(B2（i,3））＊（real（Y（i,j))＊imag(B2（j,3)）+imag（Y（i，j）)*real（B2（j,3）））；

end end end OrgS；

％四：创建PVU 用于存储PV节点的初始电压

PVU=zeros(n-h—1,1）; t=0; for i=1:n

if B2（i,6)==3 t=t+1；

PVU（t,1）=B2(i,3)； end end PVU;

％五:创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量

h=0;

for i=1：n %对PQ节点的处理 if i～=isb＆B2（i,6)==2 h=h+1；

DetaS（2*h-1，1）=real（B2(i,2)）—OrgS(2*h-1,1); DetaS（2*h,1）=imag（B2(i,2)）-OrgS（2＊h,1); end end t=0；

for i=1:n

%六：对PV节点的处理，注意这时不可再将h初始化为0

if i～=isb&B2（i，6)==3 h=h+1； t=t+1；

DetaS(2＊h—1,1)=real（B2（i,2））—OrgS（2＊h—1,1)； DetaS(2

＊h，1）=real(PVU（t，1））^2+imag(PVU(t,1）)^2-real(B2(i,3)）^2-imag（ B2（i，3）)^2；

end end

DetaS；

％七：创建I,用于存储节点电流参数

i=zeros（n—1,1）; h=0; for i=1：n if i～=isb h=h+1；

I(h,1)=(OrgS（2*h-1，1)-OrgS(2＊h,1)*sqrt(—1）)/conj(B2(i，3))； end end I;

%八：创建Jacbi（雅可比矩阵)

Jacbi=zeros（2*n-2)； h=0; k=0;

for i=1：n ％对PQ节点的处理 if B2(i，6)==2 h=h+1; for j=1:n if j~=isb k=k+1；

if i==j ％对角元素的处理

Jacbi(2*h—1,2*k—1)=—imag(Y（i,j））＊real（B2（i，3））+real(Y(i，j))＊imag（B2(i,3))+imag(I（h,1））; Jacbi（2＊h—1,2＊k）=real(Y(i，j））＊real（B2（i,3)）+imag(Y（i，j)）＊imag(B2(i，3)）+real（I（h，1）);

Jacbi(2＊h，2*k-1）=—Jacbi(2＊h-1,2*k）+2＊real(I(h,1));

Jacbi（2*h，2*k）=Jacbi（2*h—1,2＊k-1）-2＊imag（I(h,1））；

else ％非对角元素的处理

Jacbi（2＊h-1,2*k-1)=-imag（Y（i，j））＊real（B2（i，3)）+real(Y(i，j))*imag(B2(i,3）);

Jacbi（2＊h-1,2*k)=real(Y（i,j））*real(B2(i,3）)+imag（Y(i，j))＊imag（B2（i，3)）;

Jacbi（2＊h,2＊k—1)=-Jacbi(2＊h-1,2*k）； Jacbi(2*h，2＊k）=Jacbi（2＊h-1,2＊k—1）; end

if k==(n-1） ％将用于内循环的指针置于初始值，以确保雅可比矩阵换行

k=0; end end end end end k=0；

for i=1:n ％对PV节点的处理 if B2(i，6）==3 h=h+1； for j=1：n if j～=isb k=k+1；

if i==j %对角元素的处理

Jacbi(2*h—1，2＊k-1）=—imag（Y（i，j)）＊real(B2（i，3）)+real(Y(i，j））＊imag(B2(i,3)）+imag(I（h,1))； Jacbi（2＊h-1,2＊k）=real(Y（i，j)）*real（B2（i,3））+imag(Y（i，j)）*imag(B2(i,3))+real（I（h,1))；

Jacbi(2＊h，2*k—1)=2*imag（B2(i，3)); Jacbi（2＊h，2＊k）=2＊real(B2（i,3））； else ％非对角元素的处理

Jacbi（2＊h—1,2*k-1）=—imag（Y（i，j）)*real(B2(i,3)）+real（Y（i,j)）＊imag（B2（i，3)）； Jacbi(2＊h-1,2*k）=real(Y（i,j））*real(B2(i，3））+imag(Y(i，j））*imag（B2(i,3));

Jacbi(2＊h，2*k-1)=0; Jacbi（2＊h,2*k）=0； end

if k==（n-1) ％将用于内循环的指针置于初始值，以确保雅可比矩阵换行

k=0; end end end end end Jacbi;

%九：求解修正方程，获取节点电压的不平衡量

DetaU=zeros(2＊n-2，1）; DetaU=inv（Jacbi)＊DetaS；

DetaU; %修正节点电压 j=0;

for i=1:n ％对PQ节点处理 if B2（i，6）==2 j=j+1;

B2（i，3）=B2（i，3)+DetaU(2＊j,1）+DetaU（2＊j—1，1）*sqrt(—1）; end end

for i=1:n %对PV节点的处理 if B2（i，6)==3 j=j+1；

B2(i，3)=B2(i,3)+DetaU（2*j，1）+DetaU（2＊j-1，1)*sqrt（—1）； end end B2;

％十：开始循环＊＊＊**＊＊***＊*＊*＊*＊＊*＊＊＊＊＊*****＊*******

＊**＊＊＊*＊＊*＊＊**＊＊＊********＊***＊＊＊*

while abs（max(DetaU)）〉pr

OrgS=zeros(2*n-2,1); %!！!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵

h=0； j=0； for i=1：n

if i~=isb&B2（i,6）==2 h=h+1；

for j=1:n

OrgS（2*h-1,1)=OrgS(2*h-1，1）+real(B2（i，3））*（real(Y(i，j））*real（B2（j,3））-imag(Y(i,j))＊imag(B2(j，3）））+imag(B2（i，3))＊(real(Y(i，j））＊imag(B2（j,3)）+imag(Y（i,j））*real（B2（j,3）))； OrgS(2＊h，1）=OrgS（2*h，1)+imag（B2(i，3）)＊（real(Y（i，j)）*real（B2（j，3)）-imag(Y（i,j)）*imag（B2(j，3)）)—real（B2(i，3)）＊（real（Y(i,j）)＊imag(B2（j，3）)+imag(Y（i,j）)*real(B2（j，3))）;

end end end for i=1:n

if i～=isb＆B2(i,6）==3 h=h+1; for j=1：n

OrgS(2＊h—1,1）=OrgS（2＊h-1,1)+real（B2（i,3)）*(real（Y（i,j))＊real(B2(j,3))-imag(Y(i,j))*imag（B2（j,3)))+imag(B2（i,3)）＊（real(Y（i，j)）*imag(B2（j,3)）+imag(Y(i,j))*real(B2(j，3）））; OrgS（2*h，1)=OrgS（2*h,1）+imag（B2（i,3）)＊（real(Y（i，j）)*real（B2(j，3）)-imag（Y(i，j）)*imag（B2（j,3))）—real(B2(i,3）)＊（real(Y（i,j）)*imag(B2（j，3)）+imag(Y（i,j)）＊real（B2（j,3)));

end end end OrgS; %创建DetaS h=0； for i=1:n

if i~=isb＆B2(i，6)==2

h=h+1；

DetaS（2*h—1，1)=real（B2（i,2））—OrgS（2＊h-1，1); DetaS（2＊h，1）=imag（B2(i,2))—OrgS(2*h，1）; end end t=0; for i=1:n

if i~=isb&B2（i,6)==3 h=h+1; t=t+1；

DetaS（2*h-1，1)=real(B2（i,2））-OrgS（2＊h-1，1)； DetaS（2＊h，1)=real(PVU(t，1））^2+imag（PVU（t，1）)^2—real(B2（i,3））^2—imag（B2(i,3）)^2；

end end DetaS; ％创建I i=zeros（n-1,1）; h=0; for i=1：n if i~=isb h=h+1；

I（h，1）=(OrgS（2＊h—1,1)—OrgS（2*h，1)*sqrt（—1））/conj（B2（i,3))；

end end

I； %创建Jacbi Jacbi=zeros(2＊n-2); h=0； k=0; for i=1:n

if B2(i，6)==2 h=h+1； for j=1:n if j~=isb k=k+1； if i==j

Jacbi(2*h-1,2*k-1)=—imag（Y(i,j））*real(B2(i,3））+real(Y（i,j））*imag(B2(i,3））+imag(I(h,1））； Jacbi(2＊h-1，2*k)=real(Y（i,j))*real（B2（i，3））+imag(Y(i，j）)*imag(B2(i,3））+real（I（h，1））; Jacbi（2＊h,2*k-1）=-Jacbi（2＊h-1，2＊k)+2*real(I(h,1)）;

Jacbi(2*h,2＊k）=Jacbi（2＊h-1,2*k-1)-2*imag（I(h,1)）；

else

Jacbi（2＊h-1,2＊k-1）=—imag(Y(i,j)）＊real(B2(i，3))+real(Y（i，j））＊imag（B2(i，3)); Jacbi（2*h-1,2＊k)=real（Y（i,j))＊real（B2（i,3））+imag（Y(i,j))＊imag（B2(i,3））；

Jacbi（2＊h,2*k-1）=-Jacbi(2*h—1，2＊k）; Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);

end if k==(n—1） k=0; end end end end end k=0; for i=1:n

if B2(i,6）==3 h=h+1; for j=1:n if j~=isb k=k+1； if i==j

Jacbi(2＊h—1,2＊k—1）=—imag（Y（i,j）)＊real(B2(i,3）)+real（Y（i，j)）＊imag(B2（i，3）)+imag（I（h,1)）； Jacbi（2*h-1,2＊k）=real(Y(i,j)）＊real(B2（i,3））+imag（Y(i，j）)*imag(B2（i，3）)+real（I（h,1）);

Jacbi（2＊h,2＊k—1)=2＊imag(B2(i,3）)； Jacbi（2*h,2*k)=2＊real（B2（i，3)）； else

Jacbi(2*h—1，2*k-1)=—imag(Y(i,j)）*real(B2(i,3））+real（Y(i,j)）*imag(B2(i，3））;

Jacbi（2*h-1，2＊k）=real(Y(i,j))*real(B2(i，3））

+imag(Y（i,j)）*imag（B2（i，3）);

Jacbi（2＊h,2＊k-1)=0； Jacbi（2*h,2*k）=0; end if k==(n—1) k=0； end end end end end Jacbi；

DetaU=zeros（2＊n—2，1); DetaU=inv（Jacbi）*DetaS; DetaU； %修正节点电压 j=0; for i=1：n if B2(i,6)==2 j=j+1;

B2（i,3)=B2(i，3）+DetaU(2＊j，1）+DetaU(2*j-1,1）*sqrt（—1）; end end for i=1:n

if B2（i，6）==3

j=j+1；

B2（i，3)=B2(i，3）+DetaU(2*j，1)+DetaU（2＊j-1,1）＊sqrt（—1)；

end end B2；

Times=Times+1; ％迭代次数加1 end Times;

五、

实验流程

六、实验结果

参数输入： 运行结果：

七、实验体会

通过这次实验，让我第一次接触到了MATLAB，并深切体会到了它的强大之处；潮流计算的计算机算法的实现不仅巩固了我的学过的知识，还让我学到一些MATLAB的编程,虽然在实验的过程中出现了很多的错误，但在老师的细心指导下,问题都解决啦；计算机为我们省去了大量的人工计算，希望在以后的学习中能接触到更多的软件，学习到更多的知识。