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

潮流计算的计算机算法实验报告
姓名： 学号： 班级：
一、 实验目的
掌握潮流计算的计算机算法。
熟悉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的编程,虽然在实验的过程中出现了很多的错误，但在老师的细心指导下，问题都解决啦；计算机为我们省去了大量的人工计算，希望在以后的学习中能接触到更多的软件，学习到更多的知识.