潮流计算

信息工程学院

课程设计报告书

题目: 潮流计算的手算以及matlab计算

专业：电气工程及其自动化

班级：班

学号：

学生姓名：

指导教师：

20##年6月1日

信息工程学院课程设计任务书

年月日

信息工程学院课程设计成绩评定表

摘要

潮流计算是电力系统非常重要的分析计算，用以研究系统规划和运行中提出的各种问题。对规划中的电力系统，通过潮流计算可以检验所提出的电力系统规划方案能否满足各种运行方式的要求；对运行中的电力系统，通过潮流计算可以预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以内，系统中各种元件(线路、变压器等)是否会出现过负荷，以及可能出现过负荷时应事先采取哪些预防措施等。

潮流计算是电力系统分析最基本的计算。除它自身的重要作用之外，在《电力系统分析综合程序》(PSASP)中，潮流计算还是网损计算、静态安全分析、暂态稳定计算、小干扰静态稳定计算、短路计算、静态和动态等值计算的基础。

传统的潮流计算程序缺乏图形用户界面，结果显示不直接难与其他分析功能集成。网络原始数据输入工作大量且易于出错。本文采用MATLAB语言运行WINDOWS操作系统的潮流计算软件。而采用MATLAB界面直观，运行稳定，计算准确。

关键词：电力系统潮流计算牛顿-拉夫逊算法 MATLAB

目录

1 需求分析.................................................................................................................... 1

2 原理概述.................................................................................................................... 2

2.1 手算潮流原理................................................... 2

2.2牛顿-拉夫逊算法................................................. 2

3 详细设计.................................................................................................................... 4

3.1等值电路........................................................ 4

3.2手算潮流（手算过程在附页上）.................................... 4

3.3牛顿拉夫逊算法的MATLAB实现..................................... 4

3.4程序结果....................................................... 13

4 心得体会.................................................................................................................. 16

参考文献...................................................................................................................... 17

1.需求分析

电力系统潮流计算是电力系统分析中的一种最基本的计算，是对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标是求取电力系统在给定运行状态的计算，即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有的电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和稳态分析都是以潮流计算为基础。实际电力系统的潮流技术那主要采用牛顿—拉夫逊法。

运行方式管理中，潮流是确定电网运行方式的基本出发点；在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了多个在预想操作情况下电网的潮流分布以及校验运行可靠性。在电力系统调度运行的多个领域问题是研究电力系统稳态问题的基础和前提。

在用数字见算计算机解电力系统潮流问题的开始阶段，普遍采取以节点导纳矩阵为基础的逐次代入法。这个方法的原理比较简单，要求的数字计算机内存量比较差下，适应50年代电子计算机制造水平和当时电力系统理论水平，但它的收敛性较差，当系统规模变大时，迭代次数急剧上升，在计算中往往出现迭代不收敛的情况。这就迫使电力系统的计算人员转向以阻抗矩阵为基础的逐次代入法。阻抗法改善了系统潮流计算问题的收敛性，解决了导纳无法求解的一些系统的潮流计算，在60年代获得了广泛的应用，阻抗法德主要缺点是占用计算机内存大，每次迭代的计算量大。当系统不断扩大时，这些缺点就更加突出，为了克服这些缺点，60年代中期发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需要存储各个地区系统的阻抗矩阵及它们之间联络的阻抗，这样不仅大幅度的节省了内存容量，同时也提高了计算速度。

克服阻抗法缺点是另一个途径是采用牛顿-拉夫逊法。这是数学中解决非线性方程式的典型方法，有较好的收敛性。在解决电力系统潮流计算问题时，是以导纳矩阵为基础的，因此，只要我们能在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿法潮流程序的效率。自从60年代中期，牛顿法中利用了最佳顺序消去法以后，牛顿法在收敛性。内存要求。速度方面都超过了阻抗法，成为了60年代末期以后广泛采用的优秀方法。

2．原理概述

2．1 手算潮流原理

过程1解法：总结为从已知功率、电压端，用齐头并进逐段求解功率和电压。

过程2解法：总结为“一来、二去”共两步，一来即：设所有未知电压节点的电压为线路额定电压，用第十章中适当的公式从已知功率端开始逐段求功率，直到推得已知电压点的功率；二去即：从已知电压点开始，用推得的功率和已知电压点的电压，选用第十章中适当的公式，往回逐段向未知电压点求电压。

2．2牛顿拉-夫逊算法

牛顿--拉夫逊法(简称牛顿法)在数学上是求解非线性代数方程式的有效方法。其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求解的过程。即通常所称的逐次线性化过程。

对于非线性代数方程组：

即 (2-2-1)

在待求量x的某一个初始估计值附近，将上式展开成泰勒级数并略去二阶及以上的高阶项，得到如下的经线性化的方程组：

(2-2-2)

上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量

(2-2-3)

将和相加，得到变量的第一次改进值。接着就从出发，重复上述计算过程。因此从一定的初值出发，应用牛顿法求解的迭代格式为：

(2-2-4)

(2-2-5)

上两式中：是函数对于变量x的一阶偏导数矩阵，即雅可比矩阵J;k为迭代次数。

有上式可见，牛顿法的核心便是反复形式并求解修正方程式。牛顿法当初始估计值和方程的精确解足够接近时，收敛速度非常快，具有平方收敛特性。

牛顿潮流算法突出的优点是收敛速度快，若选择到一个较好的初值，算法将具有平方收敛特性，一般迭代4~5次便可以收敛到一个非常精确的解。而且其迭代次数与所计算网络的规模基本无关。牛顿法也具有良好的收敛可靠性，对于对以节点导纳矩阵为基础的高斯法呈病态的系统，牛顿法也能可靠收敛。牛顿法所需的内存量及每次迭代所需时间均较高斯法多。

牛顿法的可靠收敛取决于有一个良好的启动初值。如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的节点上。对于正常运行的系统，各节点电压一般均在额定值附近，偏移不会太大，并且各节点间的相位角差也不大，所以对各节点可以采用统一的电压初值(也称为平直电压)，如假定：

或 (2-2-6)

这样一般能得到满意的结果。但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时，仍用上述初始电压就有可能出现问题。解决这个问题的办法可以用高斯法迭代1~2次，以此迭代结果作为牛顿法的初值。也可以先用直流法潮流求解一次以求得一个较好的角度初值，然后转入牛顿法迭代。

3.详细设计

3.1等值电路

潮流计算题目如下：

(a) 潮流计算用的电网结构图

化作等值电路

(b) 潮流计算等值网络

3.2手算潮流（手算过程在附页上）

3.3牛顿拉夫逊算法的MATLAB实现

（1）主程序

clear

basemva = 100; %基准容量

accuracy = 0.0001; %给定的收敛标准

maxiter = 10; %最大的迭代次数

% 节点数据

% IEEE 30-BUS TEST SYSTEM (American Electric Power)

% Bus Bus Voltage Angle ---Load---- -------Generator----- Injected

% No code Mag. Degree MW Mvar MW Mvar Qmin Qmax Mvar

busdata=[1 2 1.05 0.0 0.2 0.0 0.6 0.45 -1 1 0

2 0 1.00 0.0 0.25 0.0 0.0 0.0 0 0 0

3 0 1.00 0.0 0.0 0.0 0.0 0.0 0 0 0

4 0 1.00 0.0 0.0 0.0 0.0 0.0 0 0 0

5 2 1.05 0.0 0.3 0.0 0.36 0.27 -1 1 0

6 0 1.00 0.0 0.25 0.0 0.0 0.0 0 0 0

7 2 1.05 0.0 0.0 0.0 0.5 0.31 -1 1 0

8 0 1.00 0.0 0.8 0.0 0.0 0.0 0 0 0

9 2 1.05 0.0 0.35 0.0 0.5 0.375 -1 1 0

10 0 1.00 0.0 0.05 0.0 0.0 0.0 0 0 0

11 0 1.00 0.0 0.0 0.0 0.0 0.0 0 0 0

12 0 1.00 0.15 0.0 0.0 0.0 0.0 0 0 0

13 1 1.05 0.0 0.0 0.0 2.52 1.56 -10 10 0

14 0 1.00 0.0 1.2 0.0 0.0 0.0 0 0 0];

% 节点节点电压电压 ---负荷---- -------发电机----- 补偿无功

% 编号类型模值相角有功无功有功无功无功上限无功下限

%在节点类型中，1表示平衡节点，2表示PV节点，0表示PQ节点

% 线路数据

% Line code

% Bus bus R X 1/2 B = 1 for lines

% nl nr p.u. p.u. p.u. > 1 or < 1 tr. tap at bus nl

linedata=[1 2 0.011 0.263 0 0.909

2 3 0.069 0.138 0.052 1

2 14 0.243 0.481 0.023 1

3 4 0.029 0.656 0 0.909

5 6 0.014 0.328 0 0.909

6 8 0.121 0.241 0.012 1

6 14 0.069 0.138 0.027 1

7 8 0.003 0.167 0 0.909

9 8 0.003 0.054 0 0.909

8 10 0.03 0.034 0 0.909

10 11 0.026 0.052 0.001 1

11 12 0.003 0.053 0 0.909

13 14 0 0.042 0 0.909 ];

lfybus % form the bus admittance matrix 形成节点导纳矩阵

lfnewton % Load flow solution by Newton-Raphson method 用牛顿拉夫逊法计算潮流

busout % Prints the power flow solution on the screen 显示计算结果

lineflow % computes and displays the line flow and losses 计算并显示线路潮流和损耗

（2）% This program prints the power flow solution in a tabulated form

% on the screen.

%clc

disp(tech)

fprintf(' Maximum Power Mismatch = %g \n', maxerror)

fprintf(' No. of Iterations = %g \n\n', iter)

head =[' Bus Voltage Angle ------Load------ ---Generation--- Injected'

' No. Mag. Degree MW Mvar MW Mvar Mvar '

' '];

disp(head)

for n=1:nbus

fprintf(' %5g', n), fprintf(' %7.3f', Vm(n)),

fprintf(' %8.3f', deltad(n)), fprintf(' %9.3f', Pd(n)),

fprintf(' %9.3f', Qd(n)), fprintf(' %9.3f', Pg(n)),

fprintf(' %9.3f ', Qg(n)), fprintf(' %8.3f\n', Qsh(n))

end

fprintf(' \n'), fprintf(' Total ')

fprintf(' %9.3f', Pdt), fprintf(' %9.3f', Qdt),

fprintf(' %9.3f', Pgt), fprintf(' %9.3f', Qgt), fprintf(' %9.3f\n\n', Qsht)

（3）% Power flow solution by Newton-Raphson method

% 牛顿拉夫逊法计算潮流

ns=0; ng=0; Vm=0; delta=0; yload=0; deltad=0;%平衡节点数目 PV节点数目电压幅值相角线路幅值相角

nbus = length(busdata(:,1)); %节点数

for k=1:nbus %节点数据的代入，存入矩阵

n=busdata(k,1);

kb(n)=busdata(k,2);

Vm(n)=busdata(k,3);

delta(n)=busdata(k,4);

Pd(n)=busdata(k,5);

Qd(n)=busdata(k,6);

Pg(n)=busdata(k,7);

Qg(n) = busdata(k,8);

Qmin(n)=busdata(k, 9);

Qmax(n)=busdata(k, 10);

Qsh(n)=busdata(k, 11); %分别获得节点矩阵各列的数据，存入相应的列矩阵

if Vm(n) <= 0 Vm(n) = 1.0; V(n) = 1 + j*0; %当节点电压模值小于0时就赋值为1并把相角赋为0

else delta(n) = pi/180*delta(n); %把角度转换成弧度

V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));%形成复数型电压

E(n)=Vm(n)*cos(delta(n));

F(n)=Vm(n)*sin(delta(n)); %获得电压的实部和虚部

P(n)=(Pg(n)-Pd(n))/basemva;

Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;

S(n) = P(n) + j*Q(n); %获得复数型功率

end

for k=1:nbus

if kb(k) == 1 ns = ns+1; else, end %平衡节点数

if kb(k) == 2 ng = ng+1; else, end %PV节点数

ngs(k) = ng; %已有的PV节点数

nss(k) = ns; %已有的平衡节点数

end

Ym=abs(Ybus); t = angle(Ybus); %将导纳矩阵转换成相应的模值和相角矩阵

g=Ym.*cos(t);%得到电导矩阵

b=Ym.*sin(t);%得到电纳矩阵

ss=2;sr=nbus-ng+1;

for k=1:nbus

if kb(k)==1

kc(1)=1;Vm1(1)=Vm(k);delta1(1)=delta(k);Pd1(1)=Pd(k);Qd1(1)=Qd(k);Pg1(1)=Pg(k);Qg1(1)=Qg(k);Qmin1(1)=Qmin(k);Qmax1(1)=Qmax(k);

Qsh1(1)=Qsh(k);V1(1)=V(k);E1(1)=E(k);F1(1)=F(k);P1(1)=P(k);Q1(1)=Q(k);S1(1)=S(k);

end

if kb(k)==0

kc(ss)=0;Vm1(ss)=Vm(k);delta1(ss)=delta(k);Pd1(ss)=Pd(k);Qd1(ss)=Qd(k);Pg1(ss)=Pg(k);Qg1(ss)=Qg(k);Qmin1(ss)=Qmin(k);Qmax1(ss)=Qmax(k);

Qsh1(ss)=Qsh(k);V1(ss)=V(k);E1(ss)=E(k);F1(ss)=F(k);P1(ss)=P(k);Q1(ss)=Q(k);S1(ss)=S(k);ss=ss+1;

end

if kb(k)==2

kc(sr)=2;Vm1(sr)=Vm(k);delta1(sr)=delta(k);Pd1(sr)=Pd(k);Qd1(sr)=Qd(k);Pg1(sr)=Pg(k);Qg1(sr)=Qg(k);Qmin1(sr)=Qmin(k);Qmax1(sr)=Qmax(k);

Qsh1(sr)=Qsh(k);V1(sr)=V(k);E1(sr)=E(k);F1(sr)=F(k);P1(sr)=P(k);Q1(sr)=Q(k);S1(sr)=S(k);sr=sr+1;

end

end%对节点进行排序【平衡节点 PQ节点 PV节点】

ss=2;sr=nbus-ng+1;%【PV节点数=总结点数-PQ节点数+1】

for n=1:nbus

if kb(n)==1

ss1=2;sr1=nbus-ng+1;

for k=1:nbus

if kb(k)==1

Ybus1(1,1)=Ybus(n,k);end

if kb(k)==0

Ybus1(1,ss1)=Ybus(n,k);ss1=ss1+1;end

if kb(k)==2

Ybus1(1,sr1)=Ybus(n,k);sr1=sr1+1;end

end

end %平衡节点导纳数据代入导纳矩阵

if kb(n)==0

ss1=2;sr1=nbus-ng+1;

for k=1:nbus

if kb(k)==1

Ybus1(ss,1)=Ybus(n,k);end

if kb(k)==0

Ybus1(ss,ss1)=Ybus(n,k);ss1=ss1+1;end

if kb(k)==2

Ybus1(ss,sr1)=Ybus(n,k);sr1=sr1+1;end

end

ss=ss+1;

end %PQ节点导纳数据代入导纳矩阵

if kb(n)==2

ss1=2;sr1=nbus-ng+1;

for k=1:nbus

if kb(k)==1

Ybus1(sr,1)=Ybus(n,k);end

if kb(k)==0

Ybus1(sr,ss1)=Ybus(n,k);ss1=ss1+1;end

if kb(k)==2

Ybus1(sr,sr1)=Ybus(n,k);sr1=sr1+1;end

end

sr=sr+1; %PV节点导纳数据代入导纳矩阵

end

end%重新形成节点导纳矩阵

Ym1=abs(Ybus1); t1 = angle(Ybus1); %将导纳矩阵转换成相应的模值和相角矩阵

g1=Ym1.*cos(t1);%得到电导矩阵

b1=Ym1.*sin(t1);%得到电纳矩阵

m=2*nbus-2; %雅可比矩阵的行数即雅可比矩阵的列数

maxerror = 1; %最大误差

converge=1; %是否收敛

iter = 0; %迭代次数

% Start of iterations %迭代开始

clear A DC J DX

while maxerror >= accuracy & iter <= maxiter % Test for max. power mismatch

for i=1:m

for k=1:m

A(i,k)=0; %Initializing Jacobian matrix初始化雅可比矩阵为全0阵

end, end

iter = iter+1;

for n=1:nbus

nn=n-ns;%n之前的节点数

%lm=nbus+n-ngs(n)-nss(n)-ns;%全部PQ节点数+全部PV节点数+n之前的PQ节点数

ln=nbus+n-ns-ns;%n之前的节点数+全部PQ节点数+全部PV节点数

A11=0; B22=0; P11=0;Q22=0;%J33=0; L44=0;

for l=1:nbus

if l~=n

%if nl(i) == n | nr(i) == n

% if nl(i) == n, l = nr(i); end

% if nr(i) == n, l = nl(i); end%找到与n节点相连的支路

B22=B22+ F1(l)*g1(n,l)+E1(l)*b1(n,l);

A11=A11+ E1(l)*g1(n,l)-F1(l)*b1(n,l);

%if kb(n)~=1

%else, end

if kc(n) ~= 1 & kc(l) ~=1%n和l都不是平衡节点

ll = l -ns;%l之前的节点数

lk = nbus+l-ns-ns;%全部PQ节点数+全部PV节点数+l之前节点数

%H阵的非对角线元素

A(nn, ll) =-E1(n)*b1(n,l)+g1(n,l)*F1(n);

%N阵的非对角线元素

A(nn, lk) =E1(n)*g1(n,l)+b1(n,l)*F1(n);

% if kb(n) == 0

% end

if kc(n) == 0

%J阵的非对角线元素

A(ln, ll)=-A(nn, lk);

%L阵的非对角线元素

A(ln, lk)=A(nn, ll);end

%if kb(n) == 0 & kb(l) == 0

%A(lm, lk) =-Vm(n)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));end

else end

end

P11=A11+E1(n)*g1(n,n)-F1(n)*b1(n,n);Q22=B22+F1(n)*g1(n,n)+E1(n)*b1(n,n);

Pk =P11*E1(n)+Q22*F1(n);

Qk =P11*F1(n)-Q22*E1(n);

if kc(n) == 1 P1(n)=Pk; Q1(n) = Qk; end % Swing bus P平衡节点功率

if kc(n) == 2 Q1(n)=Qk;

if Qmax1(n) ~= 0

Qgc1 = Q1(n)*basemva + Qd1(n) - Qsh1(n);

if iter <= 7 % Between the 2th & 6th iterations

if iter > 2 % the Mvar of generator buses are

if Qgc1 < Qmin1(n), % tested. If not within limits Vm(n)

Vm1(n) = Vm1(n) + 0.01;

E1(n)=Vm1(n)*cos(delta1(n));

F1(n)=Vm1(n)*sin(delta1(n)); % is changed in steps of 0.01 pu to

elseif Qgc1 > Qmax1(n), % bring the generator Mvar within

Vm1(n) = Vm1(n) - 0.01;

E1(n)=Vm1(n)*cos(delta1(n));

F1(n)=Vm1(n)*sin(delta1(n));end % the specified limits.

else, end

else,end

end

if kc(n) ~= 1

A(nn,nn)=-b1(n,n)*E1(n)+g1(n,n)*F1(n)+Q22; %H阵的对角线元素

A(nn,ln)=g1(n,n)*E1(n)+b1(n,n)*F1(n)+P11; %N阵的对角线元素

DC(nn) = P1(n)-Pk;%求取有功不平衡量

end

if kc(n)==0

A(ln,nn)=P11-g1(n,n)*E1(n)-b1(n,n)*F1(n); %J阵的对角线元素

A(ln,ln)=-b1(n,n)*E1(n)+g1(n,n)*F1(n)-Q22; %L阵的对角线元素

DC(ln) = Q1(n)-Qk;

end

if kc(n)==2

A(ln,nn)=2*F1(n);%R阵的对角线元素

A(ln,ln)=2*E1(n);%S阵的对角线元素

DC(ln) = Vm1(n)^2-(E1(n)^2+F1(n)^2);

end

DX=A\DC';%求修正量

%计算新的迭代解

for n=1:nbus

nn=n-ns;

ln=nbus+n-ns-ns;

if kc(n)==0

F1(n)=F1(n)+DX(nn);E1(n)=E1(n)+DX(ln);end

if kc(n)==2

F1(n)=F1(n)+DX(nn);E1(n)=E1(n)+DX(ln);end

end

maxerror=max(abs(DC));%求最大误差

if iter == maxiter & maxerror > accuracy %当迭代达到最大次数仍不能满足要求时输出警告

fprintf('\nWARNING: Iterative solution did not converged after ')

fprintf('%g', iter), fprintf(' iterations.\n\n')

fprintf('Press Enter to terminate the iterations and print the results \n')

converge = 0; pause, else, end

end

if converge ~= 1%如果不收敛

tech= (' ITERATIVE SOLUTION DID NOT CONVERGE'); else,

tech=(' Power Flow Solution by Newton-Raphson Method');

end

V1 = E1+j*F1;%复数电压矩阵

delta1=angle(V1);

deltad=180/pi*delta1;

i=sqrt(-1);

k=0;

for n = 1:nbus

if kc(n) == 1%平衡节点

k=k+1;

S1(n)= P1(n)+j*Q1(n);

Pg1(n) = P1(n)*basemva + Pd1(n);

Qg1(n) = Q1(n)*basemva + Qd1(n) - Qsh1(n);

Pgg1(k)=Pg1(n);

Qgg1(k)=Qg1(n); %june 97

elseif kc(n) ==2%PV节点

k=k+1;

S1(n)=P1(n)+j*Q1(n);

Qg1(n) = Q1(n)*basemva + Qd1(n) - Qsh1(n);

Pgg1(k)=Pg1(n);

Qgg1(k)=Qg1(n); % June 1997

end

yload1(n) = (Pd1(n)- j*Qd1(n)+j*Qsh1(n))/(basemva*Vm1(n)^2);

end

ss=2;sr=nbus-ng+1;

for k=1:nbus

if kb(k)==1

Vm(k)=Vm1(1);delta(k)=delta1(1);Pd(k)=Pd1(1);Qd(k)=Qd1(1);Pg(k)=Pg1(1);Qg(k)=Qg1(1);Qmin(k)=Qmin1(1);Qmax(k)=Qmax1(1);

Qsh(k)=Qsh1(1);V(k)=V1(1);E(k)=E1(1);F(k)=F1(1);P(k)=P1(1);Q(k)=Q1(1);S(k)=S1(1);

end

if kb(k)==0

Vm(k)=Vm1(ss);delta(k)=delta1(ss);Pd(k)=Pd1(ss);Qd(k)=Qd1(ss);Pg(k)=Pg1(ss);Qg(k)=Qg1(ss);Qmin(k)=Qmin1(ss);Qmax(k)=Qmax1(ss);

Qsh(k)=Qsh1(ss);V(k)=V1(ss);E(k)=E1(ss);F(k)=F1(ss);P(k)=P1(ss);Q(k)=Q1(ss);S(k)=S1(ss);ss=ss+1;

end

if kb(k)==2

Vm(k)=Vm1(sr);delta(k)=delta1(sr);Pd(k)=Pd1(sr);Qd(k)=Qd1(sr);Pg(k)=Pg1(sr);Qg(k)=Qg1(sr);Qmin(k)=Qmin1(sr);Qmax(k)=Qmax1(sr);

Qsh(k)=Qsh1(sr);V(k)=V1(sr);E(k)=E1(sr);F(k)=F1(sr);P(k)=P1(sr);Q(k)=Q1(sr);S(k)=S1(sr);sr=sr+1;

end

deltad=180/pi*delta;

busdata(:,3)=Vm'; busdata(:,4)=deltad';

Pgt = sum(Pg); Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);%求取各项之和

%clear A DC DX J11 J22 J33 J44 Qk delta lk ll lm

%clear A DC DX J11 J22 J33 Qk delta lk ll lm

（4）% This program obtains th Bus Admittance Matrix for power flow solution

% 【潮流的节点导纳矩阵】

j=sqrt(-1); i = sqrt(-1); %引入虚数单位

nl = linedata(:,1);% GGG(:,k)表示取数组GGG中第k列的所有元

nr = linedata(:,2); R = linedata(:,3);

X = linedata(:,4); Bc = j*linedata(:,5); a = linedata(:, 6);

nbr=length(linedata(:,1));%支路数目

nbus = max(max(nl), max(nr));%第1、nr列的最大值构成的数列中再次取出最大值

Z = R + j*X; y= ones(nbr,1)./Z; %branch admittance【支路导纳】

%Y = ones(m,n) 或 Y = ones([m n])返回元素都为1的m*n矩阵，m和n都为标量。

for n = 1:nbr %支路1至nbr的循环

if a(n) <= 0 a(n) = 1; else end

Ybus=zeros(nbus,nbus); % initialize Ybus to zero %初始化为零矩阵

% formation of the off diagonal elements【形成非对角元素】

for k=1:nbr;

Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);%第1至nbr节点的互导纳

Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));

end

% formation of the diagonal elements【形成对角元素】

for n=1:nbus

for k=1:nbr

if nl(k)==n

Ybus(n,n) = Ybus(n,n)+y(k)/(a(k)^2) + Bc(k);%y(k)/(a(k)^2)为线路导纳经过变压器的归算

elseif nr(k)==n

Ybus(n,n) = Ybus(n,n)+y(k) +Bc(k);%Bc为线路电纳 y为n条支路的电导

else, end %对角元素的修改

end

clear Pgg

（5）% This program is used in conjunction with lfgauss or lf Newton

% for the computation of line flow and line losses.

SLT = 0;

fprintf('\n')

fprintf(' Line Flow and Losses \n\n')

fprintf(' --Line-- Power at bus & line flow --Line loss-- Transformer\n')

fprintf(' from to MW Mvar MVA MW Mvar tap\n')

for n = 1:nbus

busprt = 0;

for L = 1:nbr;

if busprt == 0

fprintf(' \n'), fprintf('%6g', n), fprintf(' %9.3f', P(n)*basemva)

fprintf('%9.3f', Q(n)*basemva), fprintf('%9.3f\n', abs(S(n)*basemva))

busprt = 1;

else, end

if nl(L)==n k = nr(L);

In = (V(n) - a(L)*V(k))*y(L)/a(L)^2 + Bc(L)/a(L)^2*V(n);

Ik = (V(k) - V(n)/a(L))*y(L) + Bc(L)*V(k);

Snk = V(n)*conj(In)*basemva;

Skn = V(k)*conj(Ik)*basemva;

SL = Snk + Skn;

SLT = SLT + SL;

elseif nr(L)==n k = nl(L);

In = (V(n) - V(k)/a(L))*y(L) + Bc(L)*V(n);

Ik = (V(k) - a(L)*V(n))*y(L)/a(L)^2 + Bc(L)/a(L)^2*V(k);

Snk = V(n)*conj(In)*basemva;

Skn = V(k)*conj(Ik)*basemva;

SL = Snk + Skn;

SLT = SLT + SL;

else, end

if nl(L)==n | nr(L)==n

fprintf('%12g', k),

fprintf('%9.3f', real(Snk)), fprintf('%9.3f', imag(Snk))

fprintf('%9.3f', abs(Snk)),

fprintf('%9.3f', real(SL)),

if nl(L) ==n & a(L) ~= 1

fprintf('%9.3f', imag(SL)), fprintf('%9.3f\n', a(L))

else, fprintf('%9.3f\n', imag(SL))

end

else, end

end

SLT = SLT/2;

fprintf(' \n'), fprintf(' Total loss ')

fprintf('%9.3f', real(SLT)), fprintf('%9.3f\n', imag(SLT))

clear Ik In SL SLT Skn Snk

3.4程序结果

Power Flow Solution by Newton-Raphson Method

Maximum Power Mismatch = 6.05634e-006

No. of Iterations = 10

Bus Voltage Angle ------Load------ ---Generation--- Injected

No. Mag. Degree MW Mvar MW Mvar Mvar

1 1.100 -1.350 0.200 0.000 0.600 -4.632 0.000

2 1.000 -1.410 0.250 0.000 0.000 0.000 0.000

3 1.000 -1.618 0.000 0.000 0.000 0.000 0.000

4 1.000 -1.618 0.000 0.000 0.000 0.000 0.000

5 1.060 -0.157 0.300 0.000 0.360 -0.798 0.000

6 1.000 -0.170 0.250 0.000 0.000 0.000 0.000

7 1.060 -0.126 0.000 0.000 0.500 -0.621 0.000

8 1.000 -0.162 0.800 0.000 0.000 0.000 0.000

9 1.060 -0.156 0.350 0.000 0.500 -1.900 0.000

10 1.000 -0.166 0.050 0.000 0.000 0.000 0.000

11 1.000 -0.168 0.000 0.000 0.000 0.000 0.000

12 1.000 -0.168 0.000 0.000 0.000 0.000 0.000

13 1.050 0.000 0.000 0.000 1.812 -24.170 0.000

14 1.000 -0.032 1.200 0.000 0.000 0.000 0.000

Total 3.400 0.000 3.772 -32.120 0.000

Line Flow and Losses

--Line-- Power at bus & line flow --Line loss-- Transformer

from to MW Mvar MVA MW Mvar tap

1 0.400 -4.632 4.649

2 0.400 -4.632 4.649 0.002 0.039 0.909

2 -0.250 0.000 0.250

1 -0.398 4.670 4.687 0.002 0.039

3 0.028 -15.539 15.539 0.028 -15.539

14 0.120 10.869 10.869 0.333 -5.880

3 0.000 0.000 0.000

2 -0.000 -0.000 0.000 0.028 -15.539

4 -0.000 0.000 0.000 0.000 -0.000 0.909

4 0.000 0.000 0.000

3 0.000 -0.000 0.000 0.000 -0.000

5 0.060 -0.798 0.800

6 0.060 -0.798 0.800 0.000 0.002 0.909

6 -0.250 0.000 0.250

5 -0.060 0.799 0.802 0.000 0.002

8 0.200 -1.079 1.097 0.000 -3.272

14 -0.390 0.280 0.480 0.008 -7.327

7 0.500 -0.621 0.797

8 0.500 -0.620 0.797 0.000 0.001 0.909

8 -0.800 0.000 0.800

6 -0.200 -2.193 2.202 0.000 -3.272

7 -0.500 0.621 0.797 0.000 0.001

9 -0.150 1.901 1.907 0.000 0.001

10 0.050 -0.330 0.333 0.000 0.000 0.909

9 0.150 -1.900 1.906

8 0.150 -1.900 1.906 0.000 0.001 0.909

10 -0.050 0.000 0.050

8 -0.050 0.330 0.333 0.000 0.000

11 0.000 -0.330 0.330 0.000 -0.330

11 0.000 0.000 0.000

10 0.000 0.000 0.000 0.000 -0.330

12 -0.000 0.000 0.000 -0.000 0.000 0.909

12 0.000 0.000 0.000

11 0.000 -0.000 0.000 -0.000 0.000

13 1.812 -24.170 24.238

14 1.812 -24.170 24.238 0.000 0.185 0.909

14 -1.200 0.000 1.200

2 0.213 -16.749 16.750 0.333 -5.880

6 0.399 -7.607 7.617 0.008 -7.327

13 -1.812 24.355 24.423 0.000 0.185

Total loss 0.372 -32.120

4.心得体会

本次课程设计，我选择了“电力系统潮流计算”这一课题。鉴于以往单独MATLAB编程实现复杂系统的潮流计算的经验几乎没有，所以刚开始进行还是困难重重的。最大的困难莫过于对于潮流计算的理解以及大系统中的MATLAB潮流具体实现中的细节调试。

在完成课程设计期间，无论是手算潮流对于计算能力的要求还是MATLAB计算对于我的编程都是一大挑战。手算潮流的时候，因为节点较多，计算的繁复。在有限的时间里，我也只算了部分，其余的参考了一些资料。另外，MATLAB实现潮流计算的时候，借鉴了老师的程序，自己多读多想，由整体到部分分析了程序。期间遇到了各种不易理解的细节。通过网上找资料和翻阅相关图书，网上提问中网友的帮忙解决了大部分，其中目前还有一些需要向老师咨询的地方。总体上，通过自己的不断尝试，最后终于艰难地完成了潮流手算以及MATLAB编程计算。

这次电力系统及其自动化的课程设计，我开始进一步认识到了自己在专业课程上的巨大不足——实践太少，理论基础也并不扎实。的确，实践是检验真理的唯一标准，只有通过实践的考验来检验自己的能力，才会是真正学有所用的的能力，才会是作为一个合格的电气专业生所应具备的能力。另一方面，在实践中，在不断地对于实际问题的探究中，随着新问题的不断出现，如何自主的去思考与解决未曾涉及的领域，实现学习上的自我突破，这又是一个我相当欠缺的能力。

虽然距离毕业上有一段时间，但这次的惨痛教训算是给我敲响了一个绕梁三日的警钟吧。下一步，更要找好方向，给自己加油！

参考文献

[1] 何仰赞.温增银.《电力系统分析》（第三版）：华中科技大学，2002

[2] 于群.曹娜.《MATLAB/SIMULINK电力系统建模与仿真》（第一版）：机械工业出版社，2013

相关推荐

潮流计算

专栏推荐