电力系统潮流计算及Matlab编程实现

您所在的位置：网站首页 › 环形网络潮流计算实验报告 › 电力系统潮流计算及Matlab编程实现

电力系统潮流计算及Matlab编程实现

2024-07-07 23:19| 来源: 网络整理| 查看: 265

1. 潮流计算：

2. 潮流计算常用算法：

2.1 牛顿-拉夫逊算法

2.1.1 牛顿-拉夫逊法的基本原理

2.1.2 潮流计算的修正方程

2.1.3 节点电压用极坐标表示时的牛顿-拉夫逊法潮流计算

2.1.4 潮流计算程序框图

2.2 PQ分解法

3. MATLAB实例计算

1. 潮流计算：

潮流计算是电力系统分析中的一种最基本的计算，对给定系统进行潮流计算可以得到各母线上的电压、网络中的功率分布及功率损耗等。

复杂电力系统分析计算的一般方法是对整个电力系统建立数学模型，并通过计算机编程求出个节点的电压及电力系统中的功率分布。

2. 潮流计算常用算法： 2.1 牛顿-拉夫逊算法 2.1.1 牛顿-拉夫逊法的基本原理

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

对于非线性代数方程组：

$f(x)=0$ 即 $f_{i}(x_{1},x_{2},x_{3}...x_{n})=0(i=1,2,3...n)$

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

$f(x^{(0)})+\dot{f}(x^{(0)})\Delta x^{(0)}=0$

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

$\Delta x^{(0)}=-[\dot{f}(x^{(0)})]^{-1}f(x^{(0)})$

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

$\dot{f}(x^{(k)})\Delta x^{(k)}=-f(x^{(k)})$

$x^{(k+1)=x^{(k)}}+\Delta x^{(k)}$

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

2.1.2 潮流计算的修正方程

n个独立节点的电力系统，其功率方程为：

$P_{i}-jQ_{i}=U_{i}\sum_{j=1}^{n}Y_{ij}U_{j}(i=1,2,3...n)$

将节点电压用极坐标表示，导纳矩阵用复数的形式表示并展开时，可以得到：

令

若各节点的电压相量 U 和输入功率（Pi ， Qi ）均为已知（正值）时，则有▲P=0 和▲Q=0；若各节点的电压相量 U 为代求量的假设初值时，而输入功率（Pi ， Qi ）均为已知，则有▲P≠0和▲Q≠0。

设电力系统有 n 个独立节点，其中有一个平衡节点，编号为 n，有 m 个 PQ 节点，编号为 1，2，⋯，m，剩下的全为 PV 节点编号为：m+1,m+2,⋯,n-1。则有 1）对 m 个 PQ 节点，P1 ,P2 ,⋯Pm 和Q1 ,Q2 ,⋯Qm 是已知的，代求量为 U 和δ。 2）对 n-m-1 个 PV 节点， $P_{1+m},P_{2+m},...P_{n-1}$ 和 $U_{1+m},U_{2+m},...U_{n-1}$ 为已知的，待求量为δ。 3）Un ，δn 是已知的。

近似得到n+m-1个修正方程式：

2.1.3 节点电压用极坐标表示时的牛顿-拉夫逊法潮流计算

修正方程如下：

2.1.4 潮流计算程序框图

2.2 PQ分解法

略详情可见《电力系统分析》课本中。

3. MATLAB实例计算

系统框图：

原始数据：

（1）电力线路数据

I侧母线

J侧母线

电阻

电抗

电纳的1/2

0.0192

0.0575

0.0264

0.0452

0.1852

0.0204

0.057

0.1737

0.0184

0.0132

0.0379

0.0042

0.0472

0.1983

0.0209

0.0581

0.1763

0.0187

0.0119

0.0414

0.0045

0.046

0.116

0.0102

0.0267

0.082

0.0085

0.012

0.042

0.0045

0.208

0.11

0.14

0.1231

0.2559

0.0662

0.1304

0.0945

0.1987

0.221

0.1997

0.0824

0.1932

0.107

0.2185

0.0639

0.1292

0.034

0.068

0.0936

0.209

0.0324

0.0845

0.0348

0.0749

0.0727

0.1499

0.0116

0.0236

0.1

0.202

0.115

0.179

0.132

0.27

0.1885

0.3292

0.2554

0.38

0.1093

0.2087

0.2198

0.4153

0.3202

0.6027

0.2399

0.4533

0.0636

0.2

0.0214

0.0169

0.0599

0.0065

表2-1 电力线路数据

（2）变压器数据

I侧母线

J侧母线

电抗

变比

0.208

1.0155

0.556

0.9629

0.256

1.0129

表2-2 变压器数据

（3）发电机数据

母线名

母线类型

有功功率

无功功率

电压幅值

电压相角

0.5756

1.0338

0.2456

1.0058

0.35

1.023

0.1793

1.0913

0.1691

1.0883

STACK

1.05

表2-3 发电机数据

（4）负荷数据

母线名

母线类型

有功功率

无功功率

电压幅值

电压相角

-0.024

-0.012

-0.076

-0.016

-0.228

-0.109

-0.058

-0.02

-0.112

-0.075

-0.062

-0.016

-0.082

-0.025

-0.035

-0.018

-0.09

-0.058

-0.032

-0.009

-0.095

-0.034

-0.022

-0.007

-0.175

-0.112

-0.032

-0.016

-0.087

-0.067

-0.035

-0.023

-0.024

-0.009

-0.106

-0.019

-0.217

1.0338

-0.942

1.0058

-0.3

1.023

表2-4 负荷数据

代码实现：

chaoliu_freedom.m %复杂潮流计算 %编辑时间：2022/11/17 11：30 %编辑人：葛仪菲 %说明： %节点类型标号 %PQ节点 1 %PV节点 2 %slack节点 3 %编号任意不强制节点编号顺序 clear; clc %总结点数 IEEE30 n=30; %计数位用于节点重编 Gen_PV = 1; Load_PQ = 0; Load_PV = 1; Line_PQ = 0; Line_PV = 0; %Change-----记录前编号 new_Change-----记录新编号（按照PQ节点---PV节点---平衡节点） Change = zeros(1,n); new_Change = zeros(1,n); %导入发电机负荷电力线路变压器补偿电容原始数据（编号任意） load 'Gen_free.mat'; load 'Load_free.mat'; load 'Line_data_free.mat'; load 'Trans_data_free.mat'; load 'Cap_data_free.mat' %%导入发电机负荷电力线路变压器补偿电容原始数据（编号按照PQ节点---PV节点---平衡节点顺序） % load 'Gen.mat'; % load 'Load.mat'; % load 'Line_data.mat'; % load 'Trans_data.mat'; % load 'Cap_data.mat' %重新编号 PQ节点 0--m PV节点 m+1----n-1 平衡节点 n for num=1:size(Gen,1) if(Gen(num,2)==2) Change(num) = Gen(num,1); Gen(num,1)=n-size(Gen,1)+Gen_PV; Gen_PV = Gen_PV+1; new_Change(num) = Gen(num,1); elseif(Gen(num,2)==3) Change(num) = Gen(num,1); Gen(num,1)=n; new_Change(num) = Gen(num,1); end end for num=1:size(Load,1) if(Load(num,2)==1) Change(num+size(Gen,1)) = Load(num,1); Load(num,1)=1+Load_PQ; Load_PQ=Load_PQ+1; new_Change(num+size(Gen,1)) = Load(num,1); end end %编号列表对比 yes = [Change;new_Change]; %负荷更新编号 for num = 1:size(Load,1) if(Load(num,2)==2) for j = 1:length(Change) if(Load(num,1)==Change(j)) Load(num,1) = new_Change(j); break; end end end end %电力线路更新编号 for num=1:size(Line_data,1) for j = 1:length(Change) if(Line_data(num,1)==Change(j)) Line_data(num,1) = new_Change(j); break; end end for k = 1:length(Change) if(Line_data(num,2)==Change(k)) Line_data(num,2) = new_Change(k); break; end end end %变压器更新编号 for num=1:size(Trans_data,1) for j = 1:length(Change) if(Trans_data(num,1)==Change(j)) Trans_data(num,1) = new_Change(j); end if(Trans_data(num,2)==Change(j)) Trans_data(num,2) = new_Change(j); end end end %补偿电容更新编号 for num=1:size(Cap_data,1) for j = 1:length(Change) if(Cap_data(num,1)==Change(j)) Cap_data(num,1) = new_Change(j); end if(Cap_data(num,2)==Change(j)) Cap_data(num,2) = new_Change(j); end end end %电力线路参数 : I侧母线 J侧母线阻抗 1/2接地导纳 Line = [Line_data(:,1:2) Line_data(:,3)+i*Line_data(:,4) Line_data(:,5)*i]; % %变压器参数： I侧母线 J侧母线电抗变比 Trans = [Trans_data(:,1:2) Trans_data(:,3)*i Trans_data(:,4)]; % %补偿电容参数：母线电纳 Cap = [Cap_data(:,1),Cap_data(:,2)*i]; %母线计数 Bus=ones(1,n); %1-牛顿拉夫逊法， 2-PQ分解法 mode=2; Tmax=10; %最大迭代次数 limit=1.0e-12; %要求精度 %变压器π型等效阻抗参数 Zt=zeros(size(Trans,1),3); Zt(:,1)=Trans(:,3)./Trans(:,4); Zt(:,2)=Trans(:,3)./(1-Trans(:,4)); Zt(:,3)=Trans(:,3)./(Trans(:,4).^2-Trans(:,4)); Trans_pi=[Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3)]; %提取同时为发电机和负荷的编号 num2 = 1; simple = zeros(1,n); for num=1:size(Gen,1) for num1=1:size(Load,1) if(Gen(num,1)==Load(num1,1)) simple(num2) = Gen(num,1); num2 = num2+1; end end end %计算PQ,PV节点数 n=numel(Bus); %总节点数 m=n-1; %PQ节点数 for num=1:size(Gen,1)%数组行数 if Gen(num,2)==2 %除去PV节点就是PQ节点 m=m-1; end end flag = 0; for num=1:size(Load,1) if Load(num,2)==2 for num1=1:length(simple) if(Load(num,1)==simple(num1)) flag = 0; break; else flag = 1; end end if(flag==1) m=m-1; end end end %PQ节点包含浮游节点，其PQ=0 %提取P,Q,U向量 P,Q为原始数据，Pi,Qi为计算结果 P=zeros(1,n); Q=zeros(1,n); U=ones(1,n); cita=zeros(1,n); %该处为角度制 for num=1:size(Gen,1) if Gen(num,2)==1 %PQ节点 P(Gen(num,1))=Gen(num,3); Q(Gen(num,1))=Gen(num,4); end if Gen(num,2)==2 %PV节点 P(Gen(num,1))=Gen(num,3); U(Gen(num,1))=Gen(num,4); end if Gen(num,2)==3 %slack节点 U(Gen(num,1))=Gen(num,3); cita(Gen(num,1))=Gen(num,4); end end for num=1:size(Load,1) if Load(num,2)==1 %PQ节点 P(Load(num,1))=P(Load(num,1))+Load(num,3); %考虑到节点可能同时存在发电机和负荷进行累加 Q(Load(num,1))=Q(Load(num,1))+Load(num,4); end if Load(num,2)==2 %PV节点 P(Load(num,1))=P(Load(num,1))+Load(num,3); U(Load(num,1))=Load(num,4); end if Load(num,2)==3 %slack节点 U(Load(num,1))=Load(num,3); cita(Load(num,1))=Load(num,4); end end disp('初始条件：') disp('各节点有功：') disp(P); disp('各节点无功：') disp(Q); disp('各节点电压幅值：') disp(U); cita=(deg2rad(cita)); %角度转换成弧度 disp('各节点电压相角(度)：') disp(rad2deg(cita)); %节点导纳矩阵的计算 Y=zeros(n); %新建节点导纳矩阵 y=zeros(n); %网络中的真实导纳 %计算y(i,j) for num=1:size(Line,1) ii=Line(num,1); jj=Line(num,2); y(ii,jj)=1/Line(num,3); y(jj,ii)=y(ii,jj); end for num=1:size(Trans_pi,1) ii=Trans_pi(num,1); jj=Trans_pi(num,2); y(ii,jj)=1/Trans_pi(num,3); y(jj,ii)=y(ii,jj); end %计算y(i,i) for num=1:size(Line,1) ii=Line(num,1); jj=Line(num,2); y(ii,ii)=y(ii,ii)+Line(num,4); y(jj,jj)=y(jj,jj)+Line(num,4); end for num=1:size(Trans_pi,1) ii=Trans_pi(num,1); jj=Trans_pi(num,2); y(ii,ii)=y(ii,ii)+Trans_pi(num,4); y(jj,jj)=y(jj,jj)+Trans_pi(num,5); end %算上补偿电容 for num=1:size(Cap,1) ii=Cap(num,1); y(ii,ii)=y(ii,ii)+Cap(num,2); end %由y计算Y ysum=sum(y,1); for num=1:n for j=1:n if num==j Y(num,j)=ysum(num); %自导纳 else Y(num,j)=-y(num,j); %互导纳 end end end disp('节点导纳矩阵：'); disp(Y); G=real(Y); %电导矩阵 B=imag(Y); %电纳矩阵 %计算▲P ▲Q Pi Qi [dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita ); disp('有功不平衡量：'); disp(dP); disp('无功不平衡量：'); disp(dQ); for num=1:Tmax fprintf('第%d次迭代\n',num); %雅可比矩阵的计算 if(mode==1) J=Jacobi( n,m,U,cita,B,G,Pi,Qi ); disp('雅可比矩阵'); disp(J); end %求解节点电压修正量 if(mode==1) [dU,dcita]=Correct( n,m,U,dP,dQ,J ); else [dU,dcita]=PQ_LJ( n,m,dP,dQ,U,B ); end disp('电压、相角修正量：'); disp(dU); disp(rad2deg(dcita)); %修正节点电压 U=U+dU; cita=cita+dcita; disp('节点电压幅值：'); disp(U); disp('节点电压相角：'); disp(rad2deg(cita)); %计算功率不平衡量 [dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita ); disp('有功不平衡量：'); disp(dP); disp('无功不平衡量：'); disp(dQ); %收敛判断 if (max(abs(dP))

【本文地址】

公司简介

联系我们