电力系统潮流计算

电力系统潮流计算The final revision was on November 23, 2020
电力系统
课程设计题目: 电力系统潮流计算
院系名称：电气工程学院
专业班级：电气F1206班
学生姓名：
学号：
指导教师：张孝远
1
2
节点的分类 (5)
3 计算方法简介 (6)
牛顿—拉夫逊法原理 (6)
牛顿—拉夫逊法概要 (6)
牛顿法的框图及求解过程 (8)
MATLAB简介 (9)
4 潮流分布计算 (10)
系统的一次接线图 (10)
参数计算 (10)
丰大及枯大下地潮流分布情况 (14)
该地区变压器的有功潮流分布数据 (15)
重、过载负荷元件统计表 (17)
5 设计心得 (17)
参考文献 (18)
附录:程序 (19)
原始资料
一、系统接线图见附件1。

二、系统中包含发电厂、变电站、及其间的联络线路。

500kV变电站以外的系统以一个等值发电机代替。

各元件的参数见附件2。

设计任务
1、手动画出该系统的电气一次接线图，建立实际网络和模拟网络之间的联系。

2、根据已有资料，先手算出各元件的参数，后再用Matlab表格核算出各元件的参数。

3、潮流计算
1）对两种不同运行方式进行潮流计算，注意110kV电网开环运行。

2）注意将电压调整到合理的范围
110kV母线电压控制在106kV～117kV之间;
220kV母线电压控制在220 kV～242kV之间。

附件一：
72
水电站2
水电站1
30
3x40
C
20+8
B 2x8
A
2x31.5
D
4x7.5
水电站5
E
2x10
90+120
H
12.5+31.5
F
G
1x31.5
水电站3
24
L
2x150
火电厂
1x50
M
110kV线路220kV线路课程设计地理接线示意图
110kV变电站220kV变电站牵引站火电厂水电站500kV变电站
附件二：
1、变压器：两个220kV变电站均采用参数一致的三绕组变压器，具体参数如下。

220kV变电站参数表
110kV及以下的变电站的变压器省略，即可将负荷直接挂在110kV母线上。

而110kV升压变只计及以下参数。

2、线路：具体参数如下。

3、发电机
各发电机的参数如下：
出力情况：
水力发电机丰大出力70%,枯大出力20%。

火力发电机丰大出力80%,枯大出力80%。

4、负荷
各110kV变电站丰大负荷按该站变电容量的50%估算，枯大负荷按该站变电容量的60%估算。

两个220kV变电站的低压侧上各挂10MW的负荷，中压侧各挂20MW负荷。

功率因素均为。

5、并联电容器
两个220kV变电站的低压侧上均装设并联补偿。

补偿总量按该站变电容量的20%装设，分组原则以每组电容器的容量不超过10MVar且经济性较好为准。

1 概述
潮流计算是电力系统最基本最常用的计算。

根据系统给定的运行条件，网络接线及元件参数，通过潮流计算可以确定各母线的电压，包括电压的幅值和相角，各元件流过的功率，整个系统的功率损耗等一系列系统中的潮流数据。

近几年，对潮流算法的研究仍然是如何改善传统的潮流算法，即高斯-塞德尔法、牛顿法和快速解耦法。

牛顿法，由于其在求解非线性潮流方程时采用的是逐次线性化的方法，为了进一步提高算法的收敛性和计算速度，人们考虑采用将泰勒级数的高阶项或非线性项也考虑进来，于是产生了二阶潮流算法。

后来又提出了根据直角坐标形式的潮流方程是一个二次代数方程的特点，提出了采用直角坐标的保留非线性快速潮流算法。

潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种，牛顿—拉夫逊Newton-Raphson法是数学上解非线性方程组的有效方法，有较好的收敛性。

将N-R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性，稀疏性及节点编号顺序优划等技巧，使N-R法在收敛性，占用内存，计算速度等方面的优点都超过了阻抗法
总结为在电力系统运行方式和规划方案的研究中都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。

同时为了实时监控电力系统的运行状态也需要进行大量而快速的潮流计算。

因此潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。

在系统规划设计和安排系统的运行方式时采用离线潮流计算在电力系统运行状态的实时监控中则采用在线潮流计算。

2 潮流计算节点介绍
常规的电力系统潮流计算中一般具有三种类型的节点：PQ 、PV 及平衡节点。

一个节点有四个变量，即注入有功功率、注入无功功率，电压大小及相角。

常规的潮流计算一般给定其中的二个变量：PQ 节点（注入有功功率及无功功率），PV 节点（注入有功功率及电压的大小），平衡节点（电压的大小及相角）。

变量的分类
负荷消耗的有功、无功功率——1L P 、1L Q 、2L P 、2L Q 电源发出的有功、无功功率——1G P 、1G Q 、2G P 、2G Q 母线或节点的电压大小和相位——1U 、2U 、1δ、2δ
在这十二个变量中，负荷消耗的有功和无功功率无法控制，因它们取决于用户，它们就称为不可控变量或是扰动变量。

电源发出的有功无功功率是可以控制的自变量，因此它们就称为控制变量。

母线或节点电压的大小和相位角——是受控制变量控制的因变量。

其中, 1U 、2U 主要受1G Q 、2G Q 的控制, 1δ、
2δ主要受1G P 、2G P 的控制。

这四个变量就是简单系统的状态变量。

为了保证系统的正常运行必须满足以下的约束条件： 对控制变量 对没有电源的节点则为 对状态变量i U 的约束条件则是 对某些状态变量i δ还有如下的约束条件
节点的分类
⑴ 第一类称PQ 节点。

等值负荷功率Li P 、Li Q 和等值电源功率Gi P 、Gi Q 是给定的，从而注入功率i P 、i Q 是给定的，待求的则是节点电压的大小i U 和相位角i δ。

属于这类节点的有按给定有功、无功率发电的发电厂母线和没有其他电源的变电所母线。

⑵ 第二类称PV 节点。

等值负荷和等值电源的有功功率Li P 、Gi P 是给定的，从而注入有功功率i P 是给定的。

等值负荷的无功功率Li Q 和节点电压的大小
i U 也是给定的。

待求的则是等值电源的无功功率Gi Q ，从而注入无功功率i
Q 和节点电压的相位角i δ。

有一定无功功率储备的发电厂和有一定无功功率电源的变电所母线都可以作为PV 节点；
⑶ 第三类平衡节点。

潮流计算时一般只设一个平衡节点。

等值负荷功率
Ls P 、Ls Q 是给定的，节点电压的大小和相位也是给定的。

担负调整系统频率任
务的发电厂母线往往被选作为平衡节点。

3 计算方法简介
牛顿—拉夫逊法原理
牛顿—拉夫逊法概要
首先对一般的牛顿—拉夫逊法作一简单的说明。

已知一个变量X 函数为：
到此方程时，由适当的近似值)
0(X
出发，根据：
反复进行计算，当)
(n X 满足适当的收敛条件就是上面方程的根。

这样的方
法就是所谓的牛顿—拉夫逊法。

这一方法还可以做下面的解释，设第n 次迭代得到的解语真值之差，即
)
(n X 的误差为ε时，则： 把)()
(ε+n X
f 在)
(n X
附近对ε用泰勒级数展开
上式省略去2ε以后部分
)
(n X 的误差可以近似由上式计算出来。

比较两式，可以看出牛顿—拉夫逊法的休整量和)
(n X 的误差的一次项相
等。

用同样的方法考虑，给出n 个变量的n 个方程：
对其近似解1X '得修正量1
X '∆可以通过解下边的方程来确定： 式中等号右边的矩阵
n
n
x f ∂∂都是对于n X X X ''',,,2
1 的值。

这一矩阵称为雅可比（JACOBI ）矩阵。

按上述得到的修正向量n X X X '∆'∆'∆,,,2
1 后，得到如下关系
这比n X X X ''',,,2
1 更接近真实值。

这一步在收敛到希望的值以前重复进行，一般要反复计算满足
ε为预先规定的小正数，1+n n X 是第n 次迭代n X 的近似值。

牛顿法的框图及求解过程
1、用牛顿法计算潮流时，有以下的步骤： （1）给这各节点电压初始值)
0()0(,f
e ；
（2）将以上电压初始值代入公式，求修正方程的常数项向量
)0(2)0()0()(,,V Q P ∆∆∆；
（3）将电压初始值在带入上述公式，求出修正方程中系数矩阵的各元素。

（4）解修正方程式)0()0(,f e ∆∆；
（5）修正各节点电压)0()0()1(e e e ∆+=，)0()0()1(f f f ∆+=； （6）将)1(e ，)1(f 在带入方程式，求出)1(2)1()1()(,,V Q P ∆∆∆； （7）检验是否收敛，即{
}ε<∆∆)
()(,max k i
k i Q P
（8）如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。

如果不收敛，转回（2）进行下次迭代计算，直到收敛为止。

2、程序框图如下
MATLAB 简介
MATLAB 是用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB 和Simulink 两大部分。

是由美国mathworks 公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。

它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C 、Fortran ）的编辑模式，代表了当今国际科学计算软件的先进水平。

MATLAB 是一种交互式、面向对象的程序设计语言广泛应用于工业界与学
术界主要用于矩阵运算
同时在数值分析、自动控制模拟、数字信号处理、
动态分析、绘图等方面也具有强大的功能。

MATLAB 程序设计语言结构完整且具有优良的移植性它的基本数据元素是不需要定义的数组。

它可以高效率地解决工业计算问题
特别是关于矩阵和
矢量的计算。

MATLAB 与C 语言和FORTRAN 语言相比更容易被掌握。

通过M 语言可以用类似数学公式的方式来编写算法
大大降低了程序所需的难度并节省了
时间，从而可把主要的精力集中在算法的构思而不是编程上。

目前电子计算机已广泛应用于电力系统的分析计算潮流计算是其基本应用软件之一。

现有很多潮流计算方法。

对潮流计算方法有五方面的要求（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信（4）适应性好亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强。

4 潮流分布计算
系统的一次接线图
图 系统的一次连接图
参数计算
设定基准值MVA S B 100 ，Ub=，则各参数如下。

（1）发电机的次暂态电抗：X=X*Sb/Sn ，Z B =U B 2/S N
发电机参数 单位（MW ）
电厂 装机容量 枯水出力比例 丰水出力比例 丰大有功 丰大无功 枯大有功 枯大无功 短路X*'' 水电站1 30 水电
72
（2）110KV升压变压器的参数：
电阻：R=P K*U N2/ 1000S N2;
电抗：X=U K(%)*U N2/ 100 S N;
电导：G=P0/ 1000 U N2;
电纳：B=I0(%)*S N/ 100 U N2;
式中U N以KV为单位，S N以MVA为单位，P0、P K以KW为单位。

110KV变压器参数
（3）220KV三绕组变压器的参数：
电阻：R=P K*U N2/ 1000S N2;
电抗：X=U K(%)*U N2/ 100 S N;
电导：G=P0/ 1000 U N2;
电纳：B=I0(%)*S N/ 100 U N2
220kV 变电站参数表SB R* X* G*B*
高压侧绕组中压
侧绕
组
低压
侧绕
组
高压
侧绕
组
（1
）
100
容量120120120中压
侧绕
组
（2
）
100
电压220121低压
侧绕
组
（3
）
100
（4）110KV线路参数：
r=ρ/S；lg； Xl*=x1*l*SB/。

110KV线路参数标么值
序号线路名称导线
牌号
线路
长度
km
UB(K
V)
SB
（MW)
R*X*B*
1水电
站1
水电
站2
150********
2A B9580115100 3B C951115100
4水电
站4
C150********
110KV变电站负荷参数
丰大及枯大下地潮流分布情况
电压是衡量电力系统电能质量的标准之一。

电压过高或过低，都将对人身及其用电设备产生重大的影响。

保证用户的电压接近额定值是电力系统调度的基本任务之一。

当系统的电压偏离允许值时，电力系统必须应用电压调节技术调节系统电压的大小，使其维持在允许值范围内。

本文经过手算形成了等值电路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调整范围。

我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。

（1）为了保证整个系统潮流计算的完整性，我们把凡具有母线及发电机处均选作节点，这样，我们确定电厂一母线上的发电机作为平衡节点，节点号为①，其它机组作为PV节点，节点号为②③④，其余节点均为PQ节点，节点号见等值电路图。

（2）确定完节点及编号后，各条支路也相应确定了，我们对各支路参数进行了计算。

根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。

该地区变压器的有功潮流分布数据
(1)该图为丰大潮流模型图：
丰大潮流模型图
该运行方式的电压合理，负荷分配也均匀，但是有些线路的负载率偏低。

比如水电厂2-----L站的负载率仅仅为% ，M----L站的线路的负载率也只是%。

(2)该图为枯大潮流模型图：
) 枯大潮流模型图
该运行方式的电压合理，负荷分配也均匀，有些线路的负载率偏低但是有些线路的负载率则偏高。

比如水电厂2-----L站的负载率仅仅为%，M----L站的线路的负载率也只是 % D---H站的线路却重载到84%。

重、过载负荷元件统计表
5 设计心得
通过这次课程设计,我发现自己有很多不足的地方,如基础知识掌握不牢固,很多知识点都忘记了,计算速度慢及准确性低,分析问题能力不够全面等等。

同时，在设计的过程中遇到很多问题，如怎样使用WORD的工具,计算公式输入，画图等。

明白了有些东西看起来很简单，但一旦做起来却需要很多心思，要注意到很多细节问题。

要做到能好好理解课本的内容，一定要认认真真做一次计算。

因此，完成课程设计使我对课本的内容加深了理解。

总体来说，这次的课
程设计不单在专业基础方面反映了我的学习还要加倍努力，还在对一些软件的应用需要加强。

计算在各种情况下的潮流分布，对于丰大和枯大情况下的潮流分布有了明确的认识。

在本次课程设计过程当中，锻炼了自己实际操作分析能力，理论联系实际，对运行中的电力系统，通过潮流计算可以预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以内，系统中各种元件(线路、变压器等)是否会出现过负荷，以及可能出现过负荷时应事先采取哪些预防措施等。

同时采用MATLAB进行计算机的计算，在计算时采用特殊算法使得潮流计算的过程更快，效率更高。

而采用计算机的运算应该是未来的一种趋势，所以我会学习一定的编辑语言如C,C++等，以提高运算准确性和快速性。

总体而言，这次的课程设计对我们运用所学知识，发现、提出、分析和解决实际问题、锻炼实践能力的考察，使我们更清楚地知道不足之出，从而提高我们。

参考文献
[1]于永源主编.电力系统分析湖南师范大学出版社[M].1992年7月
[2]陈珩编.电力系统稳态分析水利电力出版社[M].1995年1月第二版
[3]邱晓燕.刘天琪电力系统分析的计算机算法北京中国电力出版 200
[4]李光琦.电力系统暂态分析[M].北京：水利电力出版社，
[5]陆敏政主编.电力系统习题集水利电力出版社[M].1990年
附录:程序
%潮流计算MATLAB 粗略程序
%creat a new_data
t=0;
s=0;
r=0;
w=0;
number=input('How many node are there=');
% Convert Pq to a new array
for ii=1:number
if data(ii,4)==1
t=t+1;
for jj=1:14
new_data1(t,jj)=data(ii,jj);
end;
a(1,t)=ii;
s=s+1;
%record the number of the PQ node
end;
end;
%Convert pv to a new arrayfor ii=1:number if
data(ii,4)==2 t=t+1; for
jj=1:14 new_data1(t,jj)=data(ii,jj);
end; a(1,t)=ii;
r=r+1;
%record the number of the PV
node end;end;%Convert set_v to a new arrayfor
ii=1:number if
data(ii,4)==3 t=t+1; for
jj=1:14 new_data1(t,jj)=data(ii,jj);
end; a(1,t)=ii;
w=w+1; end;end;
%creat a new_data2[x,y]=size(data2) for
ii=1:x for jj=1:2 for
mm=1:number if
data2(ii,jj)==a(1,mm) new_data2(ii,jj) =mm; end; end; end;end;
for ii=1:x for
jj=3:14 new_data2(ii,jj)=data2(ii,jj); end;end; %creat a
YY=zeros(number,number);YY=zeros(number,number);yy=zeros(number,numb er);
for ii=1:x
% for jj=1:14
iii=new_data2(ii,1);
jjj=new_data2(ii,2);
if new_data2(ii,5)==2
sub=new_data2(ii,6)./(new_data2(ii,7).*new_da ta2(ii,7)+new_data2(ii,6).*new_data2(ii,6))-
new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).* new_data2(ii,6))*i;
Y(iii,jjj)=-sub./new_data2(ii,14);
YY(iii,jjj)=sub./new_data2(ii,14);
Y(jjj,iii)=-sub/new_data2(ii,14);
YY(jjj,iii)=sub./new_data2(ii,14);
yy(iii,jjj)=(ii,14))./(new_data2(ii,14).*new_ data2(ii,14)).*sub; yy(jjj,iii)=(new_data2(ii, 14)-
1)./(new_data2(ii,14)).*sub; else
Y(iii,jjj)=-
new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).* new_data2(ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+ new_data2(ii,6).*new_data2(ii,6))*i; YY(iii, jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii ,6).*new_data2(ii,6))-
new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).* new_data2(ii,6))*i; Y(jjj,iii)=-
new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).* new_data2(ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+ new_data2(ii,6).*new_data2(ii,6))*i; YY(jjj,ii i)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6 ).*new_data2(ii,6))-
new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).* new_data2(ii,6))*i; yy(iii,jjj)=new_data2(ii,8 )./2.*i; yy(jjj,iii)=new_data2(ii,8)./2.*i;
end; %end;end;
for iii=1:number Y(iii,iii)=0;end;
%for ii=1:x % for
jj=1:14 for iii=1:number for
jj=1:number % if
iii~=jj Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);
% end; end;end;
%creat B, G
for ii=1:number for
jj=1:number G(ii,jj)=
real(Y(ii,jj)); B(ii,jj)=
imag(Y(ii,jj)); end;end;
%creat Initial_P Initial_Q Initial_Vfor
ii=1:(s+r) set_P(ii,1)=(new_data1(ii,9)-
new_data1(ii,7))./100;end;for
ii=1:s; set_Q(ii,1)=(new_data1(ii,10)-
new_data1(ii,8))./100;end;
for
ii=1:r set_V(ii,1)=new_data1(ii+s,12).*new_data1(ii+s,12);%try to modify for sike of correctingend;
Initial_p_q_v=[set_P;set_Q;set_V];
disp(Initial_p_q_v);
%creat Initial_e,Initial_f
for ii=1:number-1
e(ii,1)=1;
f(ii,1)=;%change f to test used to be end;
e(number,1)=new_data1(number,12);
f(number,1)=0;
% e(64,1)=;%test 118ieee
% f(14,1)=0;
% e(10,1)=;
%e(11,1)=;
%e(12,1)=;
%e(13,1)=;
% Start NEWTOWN CALULATION
for try_time=1:25%Creat every node consume P Q and Un=s;m=r;for
ii=1:(n+m) sum1=0; for
jj=1:(n+m+1) sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-
B(ii,jj).*f(jj,1))+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));
end; p(ii,1)=sum1;end;
for ii=1:n sum2=0; for
jj=1:(n+m+1) sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-
B(ii,jj).*f(jj,1))-
e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1)); end; q(i
i,1)=sum2;end;
disp('q=');disp(q);
u=zeros((n+m),1);for
ii=(n+1):(n+m) u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;
for ii=n+1:(n+m) extra_u((ii-
n),1)=u(ii,1);end;disp('extra_u=');disp(extra_u);
sum=[p;q;extra_u];
disp(sum)
disp(s);
disp(p);
%creat Jacobian
disp(n);
disp(m);
for ii=1:(n+m)
for jj=1:(n+m)
if (ii~=jj)
PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);
PE(ii,jj)=-G(ii,jj).*e(ii,1)-
B(ii,jj).*f(ii,1); else ss=0;
qq=0; for
num=1:(n+m+1) ss=ss+G(ii,num).*f(num,1 )+B(ii,num).*e(num,1); qq=qq+G(ii,num) .*e(num,1)-
B(ii,num).*f(num,1); end;
PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-
G(ii,jj).*f(ii,1);%TEST+1 PE(ii,jj)=-qq-
G(ii,jj).*e(ii,1)-
B(ii,jj).*f(ii,1);%TEST+1 end; end;end;
copy=;
disp('================copy================')for ii=1:n for
jj=1:m+n if
(ii~=jj) QE(ii,jj)=B(ii,jj). *e(ii,1)-
G(ii,jj).*f(ii,1);%TEST+1
QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1
else ss=0;
qq=0; for num=1:(n+m+1) ss=ss+G(ii,num).*f(num,1 )+B(ii,num).*e(num,1); qq=qq+G(ii,num) .*e(num,1)-
B(ii,num).*f(num,1); end;
QF(ii,jj)=-
qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1
QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-
G(ii,jj).*f(ii,1);%TEST+1 end;
end;end;%disp('QF');%disp(QF);%disp('QE');%disp(QE);UE=zeros((n+m ),(n+m));UF=zeros((n+m),(n+m));
for ii=n+1:n+m
for jj=1:(n+m)
if (ii~=jj)
UE(ii,jj)=0;
UF(ii,jj)=0;
else
ss=0;
qq=0;
for num=1:(n+m+1)
ss=ss+G(ii,num).*f(num,1)+B(ii,num).* e(num,1);
qq=qq+G(ii,num).*e(num,1)-
B(ii,num).*f(num,1);
end; UF(ii,jj)=-
2.*f(ii,1); UE(ii,jj)=-
2.*e(ii,1); end; end;end;
for ii=(n+1):(n+m) for jj=1:(n+m) extra_UE((ii-
n),jj)=UE(ii,jj); extra_UF((ii-
n),jj)=UF(ii,jj); end;end;%disp('extra_UE');%disp(extra_UE);%d isp('extra_Uf');%disp(extra_UF);
Jacobian=[PF,PE;QF,QE;extra_UF,extra_UE];
%disp('Jacobian=');%disp(Jacobian);
%creat substract resultsubstract_result=Initial_p_q_v-
sum;%disp('substract_result');%disp(substract_result);%calculate delta_f_edelta_f_e=-
inv(Jacobian)*substract_result;%disp(delta_f_e);for ii=1:number-
1; f(ii,1)=f(ii,1)+delta_f_e(ii,1); e(ii,1)=e(ii,1)+delt a_f_e(ii+number-1,1);end;
if max(substract_result)<1e-4 break;end ;end;
%disp('substract_result');%disp(substract_result);%disp('e=');
%disp(e);%disp('f=');%disp(f);
for ii=1:number uuu(ii,1)=
e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);U_RESULT(ii,1)=sqrt(uuu(ii,1));end; for ii=1:number for jj=1:number if
ii==a(1,jj) Old_Uresult(ii,1)=U_RESULT(jj,1) end ;end;end;
for
ii=1:number Old_Uresult(ii,2)=ii;end;%disp('U_result'); %disp(U_RESULT);disp('=====================================');disp(' The last result is :')
disp('===========U===================BUS-NO.');
disp('U=')
disp(Old_Uresult);
%calculate the angle
PI=for ii=1:number Angle(ii,1)=atan(f(ii,1)./e(ii,1))./PI*180; end;
for ii=1:number for jj=1:number if
ii==a(1,jj) Old_Angle(ii,1)=Angle(jj,1);
Old_Angle(ii,2)=ii; end;end;end;disp('=======Angle============ =======BUS-NO.');disp('Angle=');disp(Old_Angle);disp('=====Try-
times=======================')disp('Try-times=')disp(try_time);
%disp('p====================');%disp(p);
% for jj=1:number % if
a(1,jj)==118 % man=jj % end;
%end;%disp('man=========');%disp(man)
sum4=0;for
jj=1:number Y_conj(number,jj)=conj(Y(number,jj)); sum4=s
um4+Y_conj(number,jj).*(e(jj,1)-
f(jj,1)*i);end;%sum4=sum4*e(number,1);disp('===============Balance P Q=========BUS-NO');%disp(sum4);
Blance_Q(1,1)=imag(sum4)*100;Blance_Q(1,2)=a(1,number);Blance_P(1,1) =real(sum4)*100;Blance_P(1,2)=a(1,number);disp('Q Of the Balance node= ');
disp(Blance_Q);disp('P Of the Balance node=
')disp(Blance_P);disp('=================================BUS-
NO');%calculate the Q of the P-V node
Q_PV_node=zeros(number,2);Y_conj=conj(Y);for
ii=(s+1):(s+r) for
jj=1:number Q_PV_node(ii,1)=Q_PV_node(ii,1)+(e(ii,1)+f(ii,1)*i ).*(Y_conj(ii,jj).*(e(jj,1)-f(jj,1)*i)); end;end;for
ii=(s+1):(s+r);Q_PV_node(ii,1)=Q_PV_node(ii,1).*100+new_data1(ii,8)* i;end;
for
ii=1:number old_number=a(1,ii); Q_PV_node_old(old_number ,1)=Q_PV_node(ii,1);end;for
ii=1:numberQ_PV_node_old(ii,1)=imag(Q_PV_node_old(ii,1));end;for
ii=1:number Q_PV_node_old(ii,2)=ii;end;
disp('Q gen=');
disp(Q_PV_node_old);。