2024年1月21日发(作者：)
电力系统潮流计算
Company number：【0089WT-8898YT-W8CCB-BUUT-202108】
电力系统
课程设计
题 目: 电力系统潮流计算
院系名称： 电气工程学院
专业班级： 电气F1206班
学生姓名：
学 号：
指导教师： 张孝远
成绩：
目录
日期：
原始资料 ................................................. 2
指导老师签名：
1 概述 ................................................... 4
2 潮流计算节点介绍 ....................................... 4
变量的分类 .......................................... 5
节点的分类 .......................................... 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之间。
附件一：
A2x8课程设计地理接线示意图20+8BC3x40火电厂1x50水电站130水电站54x7.5GE水电站2722x101x31.5D2x31.5F12.5+31.5L2x150水电站3牵引站火电厂水电站110kV线路220kV线路110kV变电站220kV变电站500kV变电站24H90+120M
附件二：
1、变压器：两个220kV变电站均采用参数一致的三绕组变压器，具体参数如下。
220kV变电站参数表
高压侧
中压侧绕组
低压侧绕组
60（120）
绕组
容量
电压
120
220
120
110(121)
110kV及以下的变电站的变压器省略，即可将负荷直接挂在110kV母线上。而110kV升压变只计及以下参数。
110kV变电站参数表
序号
1
2
3
5
6
7
8
9
10
11
12
13
变电站名
A
B
C
D
E
G
F
水电站1
水电站2
水电站3
水电站5
水电站4
容量
16
28
120
63
20
+
30MW、2*20MVA
72MW、3*
24MW、
4*、2*20MVA
18MW、30MVA
X1
X0
14
火电厂
50MW、+40MVA
2、线路：具体参数如下。
220kV线路参数表
序号
1
2
3
M
M
M
线路名称
L
L
H
导线牌号
2x240
2x240
2x300
线路长度km
10
10
5
110kV线路参数表
序号
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
水电站1
A
B
水电站4
B
C
水电站2
E
水电站3
水电站5
T节点
T节点
H
H
水电站2
L
线路名称
水电站2
B
C
C
D
D
C
D
D
T节点
D
G
F
D
L
G
导线牌号
150
95
95
150
95
240
240
240
150
70
70
240
150
185
240/2
240
线路长度km
30
80
1
65
63
60
75
5
5
7
130
10
17
18
D
火电厂
F
D
150
240
4
1
3、发电机
各发电机的参数如下：
水电站1
水电站2
水电站3
水电站5
水电站4
华鑫电厂
Xd’’
Xq’’
装机容量
30
72
24
30
18
5
功率因数
出力情况：
水力发电机丰大出力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节点（注入有功功率及电压的大小），平衡节点（电压的大小及相角）。
变量的分类
负荷消耗的有功、无功功率——PL1、QL1、PL2、QL2
电源发出的有功、无功功率——PG1、QG1、PG2、QG2
母线或节点的电压大小和相位——U1、U2 、1、2
在这十二个变量中，负荷消耗的有功和无功功率无法控制，因它们取决于用户，它们就称为不可控变量或是扰动变量。电源发出的有功无功功率是可以控制的自变量，因此它们就称为控制变量。母线或节点电压的大小和相位角——是受控制变量控制的因变量。其中,
U1、U2主要受QG1、QG2的控制,
1、2主要受PG1、PG2的控制。这四个变量就是简单系统的状态变量。
为了保证系统的正常运行必须满足以下的约束条件：
对控制变量
对没有电源的节点则为
对状态变量Ui的约束条件则是
对某些状态变量i还有如下的约束条件
节点的分类
⑴ 第一类称PQ节点。等值负荷功率PLi、QLi和等值电源功率PGi、QGi是给定的，从而注入功率Pi、Qi是给定的，待求的则是节点电压的大小Ui和相位角i。属于这类节点的有按给定有功、无功率发电的发电厂母线和没有其他电源的变电所母线。
⑵ 第二类称PV节点。等值负荷和等值电源的有功功率PLi、PGi是给定的，从而注入有功功率Pi是给定的。等值负荷的无功功率QLi和节点电压的大小Ui 也是给定的。待求的则是等值电源的无功功率QGi，从而注入无功功率Qi和节点电压的相位角i。有一定无功功率储备的发电厂和有一定无功功率电源的变电所母线都可以作为PV节点；
⑶ 第三类平衡节点。潮流计算时一般只设一个平衡节点。等值负荷功率PLs、QLs是给定的，节点电压的大小和相位也是给定的。担负调整系统频率任务的发电厂母线往往被选作为平衡节点。
3 计算方法简介
牛顿—拉夫逊法原理
牛顿—拉夫逊法概要
首先对一般的牛顿—拉夫逊法作一简单的说明。已知一个变量X函数为：
(0)到此方程时，由适当的近似值X(n)出发，根据：
反复进行计算，当X满足适当的收敛条件就是上面方程的根。这样的方法就是所谓的牛顿—拉夫逊法。
这一方法还可以做下面的解释，设第n次迭代得到的解语真值之差，即(n)X的误差为时，则：
把f(X(n))在X(n)附近对用泰勒级数展开
上式省略去2以后部分
(n)X的误差可以近似由上式计算出来。
(n)比较两式，可以看出牛顿—拉夫逊法的休整量和X等。
用同样的方法考虑，给出n个变量的n个方程：
的误差的一次项相可以通过解下边的方程来确定：
对其近似解X1得修正量X1式中等号右边的矩阵fn,,Xn的值。这一矩阵称为雅可都是对于X1,X2xn,,Xn后，得到如下关比（JACOBI）矩阵。按上述得到的修正向量X1,X2系
,,Xn更接近真实值。这一步在收敛到希望的值以前重复进这比X1,X2行，一般要反复计算满足
为预先规定的小正数，Xnn1是第n次迭代Xn的近似值。
牛顿法的框图及求解过程
1、用牛顿法计算潮流时，有以下的步骤：
（1）给这各节点电压初始值e(0),f(0)；
（2）将以上电压初始值代入公式，求修正方程的常数项向量
P(0),Q(0),(V2)(0)；
（3）将电压初始值在带入上述公式，求出修正方程中系数矩阵的各元素。
（4）解修正方程式e(0),f(0)；
（5）修正各节点电压e(1)e(0)e(0)，f(1)f(0)f(0)；
（6）将e(1)，f(1)在带入方程式，求出P(1),Q(1),(V2)(1)；
（7）检验是否收敛，即maxPi,Qi(k)(k)
（8）如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。如果不收敛，转回（2）进行下次迭代计算，直到收敛为止。
2、程序框图如下
MATLAB简介
MATLAB是用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
MATLAB是一种交互式、面向对象的程序设计语言术界主要用于矩阵运算广泛应用于工业界与学同时在数值分析、自动控制模拟、数字信号处理、 动态分析、绘图等方面也具有强大的功能。
MATLAB程序设计语言结构完整且具有优良的移植性它的基本数据元素特别是关于矩阵和是不需要定义的数组。它可以高效率地解决工业计算问题矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言可以用类似数学公式的方式来编写算法大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。
目前电子计算机已广泛应用于电力系统的分析计算潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算方法有五方面的要求（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信（4）适应性好亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强。
4 潮流分布计算
系统的一次接线图
图 系统的一次连接图
参数计算
设定基准值SB100MVA，Ub=，则各参数如下。
（1）发电机的次暂态电抗：X=X*Sb/Sn，ZB=UB2/SN
发电机参数 单位（MW）
装机容量
30
72
枯水出力比例
丰水出力比例
丰大有功
丰大无功
枯大有功
枯大无功
短路X*''
电厂
水电站1
水电
站2
水电站3
水电站5
华鑫电厂
24
30
50
（2）110KV升压变压器的参数：
电阻：R=PK*UN2/ 1000SN2;
电抗：X=UK(%)*UN2/ 100 SN;
电导：G=P0/ 1000 UN2;
电纳：B=I0(%)*SN/ 100 UN2;
式中UN以KV为单位，SN以MVA为单位，P0、PK以KW为单位。
110KV变压器参数
变电站容量
名
A 16
B 28
C 120
D 63
E 20
G
F +
水电站30MW、1
2*20MVA
水电站72MW、2
3*
水电站24MW、
3
水电站4*、5
2*20MVA
50MW、火电厂
+40MVA
X1
X0
100
100
100
100
100
100
Sb
R*
X*
G*
B*
（3）220KV三绕组变压器的参数：
电阻：R=PK*UN2/ 1000SN2;
电抗：X=UK(%)*UN2/ 100 SN;
电导：G=P0/ 1000 UN2;
电纳：B=I0(%)*SN/ 100 UN2
220kV变电站参数表
高压侧绕组
中压侧绕组
低压侧绕组
高压侧绕组（1）
中压侧绕组（2）
低压侧绕组（3）
SB
R*
X*
G*
B*
100
容量
120 120 120
100
电压
220 121
100
（4）110KV线路参数：
r=ρ/S； lg ； Xl*=x1*l*SB/。
110KV线路参数标么值
序号
1
2
3
4
线路名称
水电站1
A
B
水电站4
水电站2
B
C
C
导线牌号
150
95
95
150
线路长度km
30
80
1
65
UB(KV)
115
115
115
115
SB（MW)
100
100
100
100
R*
X*
B*
5
6
7
8
9
10
11
12
13
14
15
16
17
18
B
C
水电站2
E
水电站3
水电站5
T节点
T节点
H
H
水电站2
L
D
火电厂
D
D
C
D
D
T节点
D
G
F
D
L
G
F
D
95
240
240
240
150
70
70
240
150
185
240
240
150
240
63
60
75
5
5
7
130
10
4
1
115
115
115
115
115
115
115
115
115
115
115
115
115
115
100
100
100
100
100
100
100
100
100
100
100
100
100
100
变电站
A
B
C
D
E
F
G
总容量
16
28
120
63
20
44
110KV变电站负荷参数
各变电站补偿电容参数
丰变枯站丰大无枯大有站比
变比
丰大有功
功
功
枯大无功
站别
运行方式
A
3
3
B
C
0
0
D
0
E
0
0
F
G
0
0
H低压侧
L低压侧
丰大
枯大
丰大及枯大下地潮流分布情况
电压是衡量电力系统电能质量的标准之一。电压过高或过低，都将对人身及其用电设备产生重大的影响。保证用户的电压接近额定值是电力系统调度的基本任务之一。当系统的电压偏离允许值时，电力系统必须应用电压调节技术调节系统电压的大小，使其维持在允许值范围内。本文经过手算形成了等值电路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调整范围。
我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。
（1）为了保证整个系统潮流计算的完整性，我们把凡具有母线及发电机处均选作节点，这样，我们确定电厂一母线上的发电机作为平衡节点，节点号为①，其它机组作为PV节点，节点号为②③④，其余节点均为PQ节点，节点号见等值电路图。
（2）确定完节点及编号后，各条支路也相应确定了，我们对各支路参数进行了计算。根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。
该地区变压器的有功潮流分布数据
(1)该图为丰大潮流模型图：
丰大潮流模型图
该运行方式的电压合理，负荷分配也均匀，但是有些线路的负载率偏低。比如水电厂2-----L站的负载率仅仅为% ，M----L站的线路的负载率也只是%。
(2)该图为枯大潮流模型图：
) 枯大潮流模型图
该运行方式的电压合理，负荷分配也均匀，有些线路的负载率偏低但是有些线路的负载率则偏高。比如水电厂2-----L站的负载率仅仅为%，M----L站的线路的负载率也只是 % D---H站的线路却重载到84%。
重、过载负荷元件统计表
类型
线路
负载率
实际是在功率
额定容量
建议
限D站负荷到4MW以下
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 array
for 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 array
for 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 Y
Y=zeros(number,number);
YY=zeros(number,number);
yy=zeros(number,number);
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_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(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,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(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_V
for 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 correcting
end;
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 U
n=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(ii,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);
%disp('extra_Uf');
%disp(extra_UF);
Jacobian=[PF,PE;QF,QE;extra_UF,extra_UE];
%disp('Jacobian=');
%disp(Jacobian);
%creat substract result
substract_result=Initial_p_q_v-sum;
%disp('substract_result');
%disp(substract_result);
%calculate delta_f_e
delta_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)+delta_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=sum4+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:number
Q_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);
发布评论