电力系统运行方式分析和计算
设计报告
】
专
业：
电气工程及其自动化
班 级： 11级电气1班 学 号： 3166 2176
姓 名： 杨玉豪 潘鸣
·
华南理工大学电力学院
2015-01-05
0、课程设计题目A3：电力系统运行方式分析和计算
姓名：指导教师：
一、一个220kV分网结构和参数如下：
变电站
变电站
#3 #5
500kV站（#1）的220kV母线视为无穷大母线，电压恒定在230kV。
*
#6
,
220kV站
220+j30
各变电站负荷曲线基本一致。
日负荷曲线主要参数为：
日负荷率：，日最小负荷系数：
各线路长度如图所示。
所有线路型号均为LGJ-2*300，基本电气参数为：
正序参数：r = Ω/km, x = Ω/km, C = µF/km；
零序参数：r0 = Ω/km, x0 = Ω/km, C0 = µF/km；
40ºC长期运行允许的最大电流：1190A。
|
燃煤发电厂G有三台机组，均采用单元接线。
电厂220kV侧采用双母接线。
发电机组主要参数如下表（在PowerWorld中选择GENTRA模型）：
机组台数
单台
容量
（MW
）
额定电
压
（EV）
功
率
因
数
升
压
变
容
量
MVA
Xd Xd’Xq。
Td0’
TJ=2
H a i,2
t/(MW2
h)
a i,1
t/(MW
h)
a i,
t/h
Pmax
(MW)
@
Pmin
(MW)
1300350；
8
7300120
1300…35087】
300
120
1250300`
7
6250100
稳定仿真中发电机采用无阻尼绕组的凸极机模型。
不考虑调速器和原动机模型。
不考虑电力系统稳定器模型。
励磁系统模型为：
！
该模型在PowerWorld 中为BPA_EG 模型，主要参数如下：
KA=40 TA= TA1= KF= TF= VRmax= VRmin= 发电厂按PV 方式运行，高压母线电压定值为。
考虑两种有功出力安排方式： 满发方式： 开机三台，所有发电机保留10%的功率裕度； 轻载方式： 仅开250MW 机组，且保留10%的功率裕度； 发电厂厂用电均按出力的7%考虑。
二、 设计的主要内容：
1、根据负荷变化和机组出力变化，拟定至少两种典型运行方式；（完成）
2、进行参数计算和标幺化，形成两种典型运行方式的潮流计算参数；（完成） —
3、用Matlab 编制潮流计算程序，可任选一种潮流计算方法；（完成）
4、用所编制的潮流程序完成典型运行方式的潮流计算，进行电压和网损分析；（完成）
5、用PowerWorld 软件进行潮流计算并与自己编制的软件计算结果进行校核和分析；（完成）
6、用所编制的潮流程序完成大方式的“N-1”潮流校核，进行线路载流能力和电压水平分析；（完成）
7、用Matlab 编制三相短路的短路容量计算程序；（完成）
8、对主要220kV 母线进行三相短路容量测算，并与PowerWorld 的计算结果进行校核；（完成）
9、自行选择2-3种故障方案，用PowerWorld 进行稳定计算，给出摇摆曲线，并计算故障的极限切除时间。
10、假定电网公司下发给燃煤发电厂G 的日发电计划曲线如下图，按照等微增率准则对三台机组进行经济负荷分配，同时采用matlab 中的quadprog 函数对三台机组进行负荷优化分配，并对两种分配结果进行分析比较。
要求给出三台机组的日发电计划曲线。
11、编制课程设计报告
/M W
P /t h
92124
6400
600795
三、 -
四、
设计要求和设计成果：
1、每两位同学为一组，自行分工，但任务不能重复；
2、每位同学对自己的设计任务编写课程设计说明书一份；
3、一组同学共同完成一份完整的设计报告；
4、设计说明和报告应包含：
以上设计任务每一部分的计算过程和结果分析；
所编制的潮流、短路和机组经济负荷分配源程序（主要语句应加注释）；
潮流计算结果（潮流图）
稳定计算的功角曲线等；
1、电力系统参数计算及其标幺化电力系统等效电路图
π型等效电路运行方式拟定
1．满载发电负荷最大运行方式：
发电厂：满发。
取发电机容量的10%为裕量，再按已知保留出力的7%作为厂用电，即发出功率为总容量的%。
负荷：采用最大负荷计算。
2．满载发电负荷最小运行方式：
发电厂：满发。
取发电机容量的10%为裕量，再按已知保留出力的7%作为厂用电，即发出功率为总容量的%。
负荷：将最大负荷与日最小负荷系数相乘，得负荷最小值。
线路参数计算
线路参数给定如下：
正序参数：r = Ω/km x = Ω/km C = µF/km；
零序参数：r0 = Ω/km x0 = Ω/km C0 = µF/km；
线路长度：L12:30km; L23:20km; L24:11km; L36:9km;
L45:11km;L6G:16km; L5G:25km;
1.线路参数有名值计算：
按照双回路线路参数考虑，应用如下公式进行有名值计算：
R= X= =ωcl
线路L12L23L24L36L45L6G L5G R/Ω
X/Ω
Ω
线路L12L23L24L36L45L6G L5G R/Ω
X/Ω
Ω
所选基准电压：230KV；基准容量：100MW。
即：
应用如下公式进行标幺值计算：
线路L12L23L24L36L45L6G L5G R
X
线路L12L23L24L36L45L6G L5G
R X
发电机参数计算
采用作为发电机端的基准电压，230kV 为系统侧的基准电压。
将三台机组分别赋予编号，两个300MW 机组为1、2号，剩余一台250MW 机组为3号。
1.功率输出
P G1=P G2=300×%=，P G3=250×%= Q G1=Q G2= = MVar Q G3=
=
相应的标幺值： PG1*=PG2*=, PG3*= QG1*=QG2*=, QG3*= 2.机端电压
VG=, 取基准电压，VG*=1 3.相关电抗值归算
应用如下公式进行归算：2
2
*X N B N B
V S S V X =⨯⨯机组 发电容量 X d *
X d ’*
X q *
1 300MW
2 300MW
3 250MW
变压器参数计算
分别将与1、2、3号机组相连的变压器分别编号为1、2、3号。
根据题目计算X 2
2
%
*100S N B N B V V S T S V X =⨯⨯；k *=(242*/(230*=
变压器 容量/MVA 变比k *
X T /Ω
X T * 1 350 2 350 3 300
负荷参数计算
已知的日负荷率和日最小负荷系数，可得到以下数据：
节点 类型
P max /MW Q max /MW P min /MW
Q min /MW
#1 平衡节点
- - - -
2、潮流计算的Matlab编程及PowerWorld仿真
牛顿拉夫逊法Matlab计算程序
牛顿拉夫逊法是计算潮流时常见的方法，该方法具有广泛的应用。
该方法是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。
多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。
方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。
牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。
相应程序详见附件1.
满发满载运行方式下计算
在这种运行方式下，三台机组全部开启，各自保留10%作为功率裕度，重用负荷占发电量7%，负荷为最大负荷。
此方式下，节点1为平衡节点，节点2、3、4、5、6、7为PQ节点，节点8、9、10为PV节点。
经过Matlab的运算，我们得到如下结果：
满发轻载运行方式下计算经Matlab运算，我们得到：
PowerWorld仿真电路图
满发满载方式运行下电力系统的潮流计算
满发满载方式运行下，发电厂三台机组均满载，取发电机容量10%为裕量，并保留出力的7%作厂用电，即发出功率为总容量的%，负荷采用最大负荷计算。
满发满载电路图：
编辑模式下，#1为平衡节点，电压恒定为。
为了反映这一特性，#1应接一台容量无穷大的发电机。
具体参数设置如图3所示（机组出力暂不确定，可填为100MW。
）
接在平衡节点处的发电机参数设置
#7为PV节点，P=，V=，正常运行时各机组发出的有功功率P1=P2 =，P3=。
具体参数设置如图4、图5所示。
接在#7处的发电机#8、#9参数设置
接在#7处的发电机#10参数设置
其余节点为PQ节点，将基准电压设置为230kV，带上相应大小的负荷。
最后，用输电线路将各节点连接起来。
潮流仿真结果如下：
各节点电压
节点编号基准电压kV标幺电压实际电压(kV)相角(度) 12300
2230
3230
4230
5230
6 7230 230
各发电机运行状态
节点编号发电机编号有功出力（MW）无功出力（MVar）电压标幺值11
88
99
1010
各支路状态
首端节点编号末端节
点编号
首端有功
功率（MW）
首端无功功
率（Mvar）
首端视在功
率（MVA）
有功损
耗（MW）
无功损耗
（Mvar）
12
23
24
36
45-47
5 6
7 7
电网线损率计算：
线损率按以下公式计算：
-
=100%
⨯
线路首端有功功率线路末端有功功率
线损率
线路首端有功功率
支路#1#2#2#3#2#4#3#6#4#5#5#7#6#7线损
率
%%%%%%%
网损率按以下公式计算：
-
=100%
⨯
D L
D
发电机输出功率P负荷P
网损率
发电机输出功率P
其中，发电机输出功率包括发电厂和平衡节点的功率，发电机输出功率和负荷均只取有功分量。
则发电机输出功率P D=+664MW=
总负荷P L =230MW+210MW+300MW+410MW+220MW =1370MW
网损率=(P D-P L)/P D×100%=%
满发轻载方式运行下电力系统的潮流计算
满发轻载方式运行下，发电厂三台机组均满载，取发电机容量10%为裕量，并保留出力的7%作为厂用电，即发出功率为总容量的%，负荷采用最小负荷计算。
满发轻载电路图
潮流仿真结果如下：
各节点电压
节点编号基准电压kV标幺电压实际电压(kV)相角(度) 12300
2230
3230
4230
5230
6 7230 230
各发电机运行状态
节点编号发电机编号有功出力（MW）无功出力（MVar）电压标幺值11
88
99
1010
各支路状态
首端节点编号末端节
点编号
首端有功
功率（MW）
首端无功功
率（Mvar）
首端视在功
率（MVA）
有功损
耗（MW）
无功损耗
（Mvar）
12 23 24 36 45
5 67 7
电网线损率计算：
线损率按以下公式计算：
-
=100%
⨯
线路首端有功功率线路末端有功功率
线损率
线路首端有功功率
网损率按以下公式计算：
-
=100%
⨯
D L
D
发电机输出功率P负荷P
网损率
发电机输出功率P
其中，发电机输出功率包括发电厂和平衡节点的功率，发电机输出功率和负荷均只取有功分量。
则发电机输出功率P D=+=
总负荷P L =++192MW++ =
网损率=(P D-P L)/P D×100%=%
Matlab和PowerWorld计算的结果比较
以大方式的结果为例：电压标幺值很接近；电角度有一定差别，但
一般小数点第一位都相同；线损有一定差别，但也都不大。
造成结果差别的原因，首先，两种程序输入参数上存在误差，Matlab
的参数会保留一定位小数，不能达到最精确值，而PowerWorld会将线
路参数强制以美制单位保存，每次查看和修改参数自动转换成公制单位
造成误差，而且这种误差还会随着查看与修改的次数累积；使用Matlab
的算法较为理论和简单，系统中各个元件和部分的反馈和相互作用没什
么体现，而PowerWorld仿真的算法更为精确和复杂，如发电机的模型
更加详细。
3、大方式下“N-1”潮流校核
使用满发大负荷运行方式，即：三台机组保留10%容量备用，发出功率的7%为厂用电，出力为机组总容量的%；负荷使用提供的220kV最大负荷。
N-1情况下，将3发电机组及其变压器统一以#7节点代表进行计算。
依次分别断开各条支路中双回线的一条线，计算节点电压和各个功率及损耗；根据设计要求：40°C长期运行允许的最大电流：1190A。
Matlab“N-1”潮流校核
N-1故障可理解为双回线中的一回被切除，只剩一回在维持运行。
这种情况会造成电力系统潮流状态的变更，需要进行潮流校核来检验系统的稳定性。
此方式中，我们依次分别将其中一条双回线参数改成相应的单回线参数带入Matlab计算。
经过Matlab运算，得到以下结果：
L67
标幺值为。
PowerWorld“N-1”潮流校核
断线电路图（以L12为例）
L12断开一路：
节点1234567
电压基准值（kV）230230230230230230
230
标幺值1
有名值（kV）230
功角（°）0
线路12232436455767电流值（kA）
分析：由于题目限定40°C长期运行允许最大电流是1190A，而表中12
回线路的断线电流为，已经超过极限，说明该线路载流能力较差，支路
12不能承受断一回线路的电流，而对于各母线的电压水平均无大变化。
L23断开一路：
节点1234567
L24断开一路：
L36断开一路：
L45断开一路：
L57断开一路：
L67断开一路：
Matlab与PowerWorld“N-1”潮流校核对比：
通过就“N-1”潮流校核分别进行Matlab和PowerWorld计算、仿真，我们发现结果差距较大，其中PowerWorld仿真下只有L12电流值超过给定值，而Matlab计算下过流线路较多。
经过反复计算、仿真，我们大胆认为中间存在的较大差异并非本身操作失误，而是软件在运算中机理不同而生：节点导纳矩阵形成过程不同。
4、220kV母线的三相短路容量测算
选用满发重荷的运行方式，平衡机基准容量为100MVA，发电机基准容量为各自的额定容量，变压器改为原来的标幺变比和标幺阻抗，对非平衡母线的所有220kv的母线进行三相短路测算。
短路容量的Matlab计算
相应Matlab计算程序详见附件2.
满发满载运行方式下短路容量计算
经过Matlab运算，分别得到短路电流和短路容量：
节点#1#2#3#4#5#6#7电流
容量1390150010501510151015101530满发轻载运行方式下短路容量计算
经过Matlab运算，分别得到短路电流和短路容量：
节点#1#2#3#4#5#6#7电流
容量970103010401040104010401060短路容量的PowerWorld仿真
短路容量仿真操作
满发满载运行方式下短路容量计算
经过PowerWorld仿真，分别得到短路电流和短路容量：
节点#1#2#3#4#5#6#7容量
满发轻载运行方式下短路容量计算
经过PowerWorld仿真，分别得到短路电流和短路容量：
节点#1#2#3#4#5#6#7容量
两种计算方法结果比较
使用Matlab和PowerWorld进行计算和仿真的结果有一定相关性，
但是存在一定差别。
我们认为除了两者输入各种参数的误差之外，还可
能是导纳矩阵算法之间的差距。
PowerWorld形成的节点导纳矩阵与我
们处理的不同，短路电流的计算准确。
Matlab的算法是我们根据书本
概念按自己的理解写成的程序代码，核心算法也遵照书本所述，可能是
我们对负荷模型处理不够准确造成了这种差别。
5、PowerWorld的网络稳定扫描
如果短路点越靠近发电机节点，则短路故障对电力系统扰动越大。
因此在以下稳定扫描仿真中，均把短路点尽可能靠近母线7上。
若在此情况下，发电机摇摆曲线仍能保持暂态稳定，则电力系统稳定性符合要求。
暂态参数设置
先按照要求添加发电机的暂态模型，设置好发电机的暂态参数。
线路的零序阻抗参数。
具体参数设置如下图：
发电机#8、#9 GENTRA 模型参数
发电机#10 GENTRA 模型参数
发电机#8、#9、#10 BPA_EG模型参数
发电机升压变压器零序参数设定
线路零序参数设定（以L67为例）故障方案
故障均设置在t=发生，仿真时间为20s。
故障方案一：母线7三相短路，切除故障线路。
故障方案一设计
故障方案一运行
从故障切除时间为的各发电机转子角摇摆曲线，我们可以看出个发电机的转子角由于故障而产生扰动，而在切除故障之后逐渐趋于稳定，各转子角的差（功角）趋于稳定，所以系统暂态稳定。
故障极限切除时间求解：
用二分法通过不断地设置切除时间，知道系统在故障切除之后不能保持暂态稳定性，到达临界值，从而确定大致的系统故障的极限切除时间。
故障切除时：
功角图
功角图
故障方案二：线路6-7三相短路，切除故障线路。
从故障切除时间为的各发电机转子角摇摆曲线，我们可以看出个发电机的转子角由于故障而产生扰动，而在切除故障之后逐渐趋于稳定，各转子角的差（功角）趋于稳定，所以系统暂态稳定。
故障极限切除时间求解：
用二分法通过不断地设置切除时间，知道系统在故障切除之后不能保持暂态稳定性，到达临界值，从而确定大致的系统故障的极限切除时间。
故障切除时：
功角图
功角图
方案总结
由此可以得出故障的极限切除时间的大致区间为[，]，我们可以取故障的极限切除时间为：。
6、发电机组的负荷经济分配
遵守等微增率原则的经济分配方案
F(P)=a i2p2+a i1p+a i0
F1(P1)=++ 150≤P1≤300
F2(P2)=++ 150≤P2≤300
F3(P3)=++ 100≤P3≤250
由题目可知，三个发电厂的微增率表达式如下所示：
Λ1=+
Λ2=+
Λ3=+
令λ1=λ2=λ3，
①P1+P2+P3=400,
P1=150; P2=150; P3=100.
②P1+P2+P3=600,
P1=300; P2=2000; P3=100.
③P1+P2+P3=795,
P1=300; P2=300; P3=195.
使用Matlab函数进行优化的分配方案
使用Matlab优化的过程中，我们按要求使用qua函数进行优化分配，这个函数专门适用于进行二次优化，使用广泛，准确度高。
详情请见附件3.
得出以下结果：
①P1=×102 P2=×102 P3=×102
②P1=300 P2=200 P3=100
③P1=×102 P2=×102 P3=×102
7、课程设计总结
关于电力系统的计算，很重要的一个步骤就是形成节点导纳矩阵。
该矩阵的正确与否能很大程度上决定运算结果的正确性。
在计算和仿真的过程中，PowerWorld的导纳矩阵形成的算法更精确，跟Matlab上常用的代码算出来的结果有所不同，很可能这些不同造成了两种方法的结果差异，比如对“N-1”潮流计算，结果的大相径庭使我们对此产生了深刻的认识。
关于PowerWorld，参数的正确设置非常重要，标幺值和有名值的区分很关键，只要这几点做好，得出来的结果比Matlab运算的更加可靠。
经过这次课程设计，我们学会了常用的电力系统仿真软件的使用。
提前熟悉了很多行业知识，为将来就业后可能遇到的跟专业知识相关的简单情况模拟和处理，提高了我们处理问题的能力。
经过这次课程设计，我们初步了解了进行电力系统分析的大致步骤和方法，将电力系统课程所学知识联系起来解决实际问题。
虽然在完成设计任务的过程中所涉及到的工具在之前的学习中没有学到，但应用到的理论、原理却是跟我们书本，甚至是跟实际的电力系统研究息息相关。
这加深了我们对以前所学知识的理解。
总的来说，这次课程设计的经历是是一笔宝贵的财富，让我们获益匪浅。
8、附件
附件一：牛顿拉夫逊法Matlab计算程序
% 牛拉法计算潮流程序
%-----------------------------------------------------------------------
% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳
% 5、支路的变比；6、支路首端处于K侧为1，1侧为0
% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；3为PV节点；
%------------------------------------------------------------------------
clear all;
format long;
n=input('请输入节点数:nodes=');
nl=input('请输入支路数:lines=');
isb=input('请输入平衡母线节点号:balance=');
pr=input('请输入误差精度:precision=');
B1=input('请输入由各支路参数形成的矩阵:B1=');
B2=input('请输入各节点参数形成的矩阵:B2=');
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
%------------------------------------------------------------------
for i=1:nl %支路数
if B1(i,6)==0 %左节点处于1侧
p=B1(i,1);q=B1(i,2);
else %左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元
Y(q,p)=Y(p,q); %非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4); %对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); %对角元1侧
end
%求导纳矩阵
disp('导纳矩阵Y=');
disp(Y)
%-------------------------------------------------------------------
G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部
for i=1:n %给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4); %PV节点电压给定模值
end
for i=1:n %给定各节点注入功率
S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量
end
%---------------------------------------------------------------------
P=real(S);Q=imag(S); %分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2
while IT2~=0 % N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb %非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)
%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2; %电压模平方
%以下针对非PV节点来求取功率差及Jacobi矩阵元素-----------------------------
if B2(i,6)~=3 %非PV节点
DP=P(i)-P1; %节点有功功率差
DQ=Q(i)-Q1; %节点无功功率差
%以上为除平衡节点外其它节点的功率计算--------------------------------------
%求取Jacobi矩阵----------------------------------------------------------
for j1=1:n
if j1~=isb&j1~=i %非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df=dQ/de
X3=X2; % X2=dp/df X3=dQ/de
X4=-X1; % X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1; % X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1; % X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%下面是针对PV节点来求取Jacobi矩阵的元素-----------------------------------------
DP=P(i)-P1; % PV节点有功误差
DV=V(i)^2-V2; % PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i %非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; % PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; % PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; % PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; % PV节点有功误差
J(m,q)=X2;
end
end
end
end
end
%以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差---------------------------------------
for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Q 或△U
for k2=k1:N1 % 从k+1列的Jacobi元素到扩展列的△P、△Q 或△U
J(k,k2)=J(k,k2)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1; % 对角元规格化K行K列对角元素赋1
%回代运算-------------------------------------------------------------------
if k~=3 % 不是第三行k > 3
k4=k-1;
for k3=3:k4 % 用k3行从第三行开始到当前行的前一行k4行消去
for k2=k1:N1 % k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end %用当前行K2列元素减去当前行k列元素乘以第k 行K2列元素
J(k3,k)=0; %当前行第k列元素已消为0
end
if k==N0 %若已到最后一行
break;
end
%前代运算------------------------------------------------------------
for k3=k1:N0 % 从k+1行到2*n最后一行
for k2=k1:N1 % 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end %用当前行K2列元素减去当前行k列元素乘以第k 行K2列元素
J(k3,k)=0; %当前行第k列元素已消为0
end
else %是第三行k=3
%第三行k=3的前代运算----------------------------------------------------
for k3=k1:N0 %从第四行到2n行（最后一行）
for k2=k1:N1 %从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end %用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0; %当前行第3列元素已消为0
end
end
%上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵
%-----------------------------------------------------------------------------------
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N); %修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N); %修改节点电压虚部
end
%修改节点电压---------------------------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr %电压偏差量是否满足要求
IT2=IT2+1; %不满足要求的节点数加1
end
end
ICT2(a)=IT2; %不满足要求的节点数
ICT1=ICT1+1; %迭代次数
end
%用高斯消去法解"w=-J*V"
disp('迭代次数：');
disp(ICT1);
disp('没有达到精度要求的个数：');
disp(ICT2);
for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2); %计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi; %计算各节点电压的角度
E(k)=e(k)+f(k)*1i; %将各节点电压用复数表示
end
%计算各输出量------------------------------------------------------
disp('各节点的实际电压标幺值E为：');
disp(E); %显示各节点的实际电压标幺值E用复数表示
disp('-----------------------------------------------------');
disp('各节点的电压大小V为：');
disp(V); %显示各节点的电压大小V的模值
disp('-----------------------------------------------------');
disp('各节点的电压相角deg为：');
disp(sida); %显示各节点的电压相角
for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点的注入电流的共轭值end
S(p)=E(p)*C(p); %计算各节点的功率S = 电压X 注入电流的共轭值
disp('各节点的功率S为：');
disp(S); %显示各节点的注入功率
disp('-----------------------------------------------------');
disp('各条支路的首端功率Si为：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...
-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...
-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('-----------------------------------------------------');
end
disp('各条支路的末端功率Sj为：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...
-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...
-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('-----------------------------------------------------');
end
disp('各条支路的功率损耗DS为：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('-----------------------------------------------------');
end
附件二：三相短路容量计算的Matlab程序代码
%短路容量计算程序
%---------------------------------------------------------------------
% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳
% 5、支路的变比；6、支路首端处于K侧为1，1侧为0
% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量
% 6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；
%Yd为修改后的节点导纳矩阵
%-----------------------------------------------------------------------
clear all;
format long;
g1=input('300MW发电机数：g1=');
g2=input('250MW发电机数：g2=');
n=input('请输入节点数:n=');
nl=input('请输入支路数:nl=');
B1=input('请输入由各支路参数形成的矩阵:B1=');
B2=input('请输入各节点参数形成的矩阵:B2=');
Y=zeros(n);
% Y为修改前节点导纳矩阵
for i=1:nl %支路数
if B1(i,6)==0 %左节点处于1侧
p=B1(i,1);q=B1(i,2);
else %左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元
Y(q,p)=Y(p,q); %非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4); %对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); %对角元1侧
end
%----------------------------------------------------------
%Y2-Y6为各PQ节点负荷的导纳
Y1=0;Y2=conj(B2(2,2));Y3=conj(B2(3,2));Y4=conj(B2(4,2));Y5=conj(B2(5,2)); Y6=conj(B 2(6,2));
Xd300=^2;
XT300=^2;
Xd250=^2;
XT250=^2;
Y7=g1/(XT300+Xd300)+g2/(XT250+Xd250);
%处理相应的负荷及机组部分的导纳
%------------------------------------------------------------
C=[Y1,Y2,Y3,Y4,Y5,Y6,Y7];
Yd=Y;
for i=1:n
Yd(i,i)=Yd(i,i)+C(i); %修改各节点自导纳end
disp(Yd);
Z = inv(Yd); %求节点阻抗矩阵
for j=1:n
I(j) = 1/Z(j,j); %电压故障前电压标幺值为1 S(j)=abs(I(j));
Sn(j)=S(j)*100;
%短路电流有名值
end
%计算完毕----------------------------------------------------
disp('各节点短路时的短路电流幅值标幺值')
disp(abs(I))
disp('短路容量有名值Sn=');
disp(Sn);
附件三：使用Matlab中quadprog函数进行优化的代码function qua( )
%主程序
clear all;
format long;
H=input('请输入由发电机参数a2除以2形成的对角矩阵：a2='); %输入a2 f=input('请输入由发电机参数a1形成的行矩阵：a1='); %输入a1
UB=input('请输入由发电机最大输出功率形成的列矩阵：UB='); %输入UB LB=input('请输入由发电机最小输出功率形成的列矩阵：LB='); %输入LB Beq=input('请输入负荷功率之和：Pd='); %输入Beq
Aeq=[1,1,1] %线性方程的系数约束矩阵
[P,F]=quadprog(H,f,[],[],Aeq,Beq,LB,UB,[],[]);
disp('各发电机的负荷功率P为：');
disp(P);
disp('3台发电机的耗量总和F为：');
disp(F+++; %加上a0系数
—。
