电力系统运行方式分析和计算

.-实用文档-.可修编.电力系统运行方式分析和计算设计报告专业：电气工程及其自动化班级：11级电气1班学号：66XX：杨玉豪潘鸣华南理工大学电力学院2015-01-050、课程设计题目A3：电力系统运行方式分析和计算：指导教师：一个220kV分网结构和参数如下：500kV站（#1）的220kV母线视为无穷大母线，电压恒定在230kV。图中，各变电站参数如下表：编号类型220kV最大负荷，MVA#1500kV站平衡节点#2220kV站230+j40#3220kV站210+j25#4220kV站300+j85#5220kV站410+j110#6220kV站220+j30各变电站负荷曲线基本一致。日负荷曲线主要参数为：日负荷率：0.85，日最小负荷系数：0.64各线路长度如图所示。所有线路型号均为LGJ-2*300，基本电气参数为：正序参数：r=0.054Ω/km,x=0.308Ω/km,C=0.0116µF/km；零序参数：r0=0.204Ω/km,x0=0.968Ω/km,C0=0.0078µF/km；40ºC长期运行允许的最大电流：1190A。燃煤发电厂G有三台机组，均采用单元接线。电厂220kV侧采用双母接线。发电机组主要参数如下表（在PowerWorld中选择GENTRA模型）：机组台数单台容量（MW）（EV）功率因数升压变容量MVAXdXd’XqTd0’TJ=2Hai,2t/(MW2h)ai,1t/(MWh)ai,0t/hPmax(MW)Pmin(MW)130010.50.853501.80.181.2870.000040.29810.22300120130010.50.853501.80.181.2870.000030.30510.32300120125010.50.853002.10.21.5760.000030.3219.38250100升压变参数均为Vs%=10.5%，变比10.5kV/242kV。不计内阻和空载损耗。稳定仿真中发电机采用无阻尼绕组的凸极机模型。不考虑调速器和原动机模型。不考虑电力系统稳定器模型。励磁系统模型为：该模型在PowerWorld中为BPA_EG模型，主要参数如下：KA=40TA=0.1TA1=0.1KF=0.05TF=0.7VRmax=3.7VRmin=0.0发电厂按PV方式运行，高压母线电压定值为1.05VN。考虑两种有功出力安排方式：满发方式：开机三台，所有发电机保留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、编制课程设计报告设计要求和设计成果：1、每两位同学为一组，自行分工，但任务不能重复；2、每位同学对自己的设计任务编写课程设计说明书一份；3、一组同学共同完成一份完整的设计报告；4、设计说明和报告应包含：以上设计任务每一部分的计算过程和结果分析；所编制的潮流、短路和机组经济负荷分配源程序（主要语句应加注释）；潮流计算结果（潮流图）稳定计算的功角曲线等；1、电力系统参数计算及其标幺化1.1电力系统等效电路图π型等效电路1.2运行方式拟定1．满载发电负荷最大运行方式：发电厂：满发。取发电机容量的10%为裕量，再按已知保留出力的7%作为厂用电，即发出功率为总容量的83.7%。负荷：采用最大负荷计算。2．满载发电负荷最小运行方式：发电厂：满发。取发电机容量的10%为裕量，再按已知保留出力的7%作为厂用电，即发出功率为总容量的83.7%。负荷：将最大负荷与日最小负荷系数相乘，得负荷最小值。1.3线路参数计算线路参数给定如下：正序参数：r=0.054Ω/kmx=0.308Ω/kmC=0.0116µF/km；零序参数：r0=0.204Ω/kmx0=0.968Ω/kmC0=0.0078µF/km；线路长度：L12:30km;L23:20km;L24:11km;L36:9km;L45:11km;L6G:16km;L5G:25km;1.线路参数有名值计算：按照双回路线路参数考虑，应用如下公式进行有名值计算：R=0.5rlX=0.5xl0.5B=ωcl正序参数：线路L12L23L24L36L45L6GL5GR/Ω0.810.540.2970.2430.2970.4320.675X/Ω4.623.081.6941.3861.6942.4643.850.5B/Ω0.00010930.000072890.000040090.000032800.000040090.000058310.00009111零序参数：线路L12L23L24L36L45L6GL5GR/Ω3.062.041.1220.9181.1221.6322.55X/Ω14.529.685.3244.3565.3247.74412.10.5B/Ω0.000073510.000049010.000026950.000022050.000026950.000039210.000061262.线路参数标幺值计算：所选基准电压：230KV；基准容量：100MW。即：应用如下公式进行标幺值计算：正序参数：线路L12L23L24L36L45L6GL5GR0.0015310.0010210.0005610.0004590.0005610.0008170.001276X0.0087330.0058220.0032020.0026200.0032020.0046580.0072780.5B0.05781970.03855880.02120760.01735120.02120760.03084600.0481972零序参数：线路L12L23L24L36L45L6GL5GR0.0057840.0038560.0021210.0017350.0021210.0030850.004820X0.0274480.0182990.0100640.0082340.0100640.0146390.0228730.5B0.03888690.02592630.01425660.01166450.01425660.02074210.03240651.4发电机参数计算采用10.5KV作为发电机端的基准电压，230kV为系统侧的基准电压。将三台机组分别赋予编号，两个300MW机组为1、2号，剩余一台250MW机组为3号。1.功率输出PG1=PG2=300×83.7%=251.1MW，PG3=250×83.7%=209.25MWQG1=QG2=QUOTE=155.617803MVarQG3=QUOTE=129.681503MVar相应的标幺值：PG1*=PG2*=2.511000,PG3*=2.092500QG1*=QG2*=1.556178,QG3*=1.2968152.机端电压VG=10.5kV,取基准电压10.5kV，VG*=13.相关电抗值归算应用如下公式进行归算：机组发电容量Xd*Xd’*Xq*1300MW0.51000.0510000.4080002300MW0.51000.0510000.4080003250MW0.71400.0680000.5100001.5变压器参数计算分别将与1、2、3号机组相连的变压器分别编号为1、2、3号。根据题目计算XT，相关计算公式为：；k*=(242*10.5)/(230*10.5)=1.05变压器容量/MVA变比k*XT/ΩXT*13500.95041317.5692000.03321223500.95041317.5692000.03321233000.95041320.4974000.0387471.6负荷参数计算已知的日负荷率0.85和日最小负荷系数0.64，可得到以下数据：负荷有名值：节点类型Pmax/MWQmax/MWPmin/MWQmin/MW#1平衡节点----#2PQ23040147.225.6#3PQ21025134.416#4PQ3008519254.4#5PQ410110262.470.4#6PQ22030140.819.2#7PQ0000#8PV300185.9233251.1155.6178#9PV300185.9233251.1155.6178#10PV250154.9361209.25129.6815负荷标幺值：节点类型电压Pmax/MWQmax/MWPmin/MWQmin/MW#1平衡节点1----#2PQ12.30.41.4720.256#3PQ12.10.251.3440.16#4PQ130.851.920.544#5PQ14.11.12.6240.704#6PQ12.20.31.4080.192#7PQ10000#8PV1.0531.8592332.5111.556178#9PV1.0531.8592332.5111.556178#10PV1.052.51.5493612.09251.2968152、潮流计算的Matlab编程及PowerWorld仿真2.1牛顿拉夫逊法Matlab计算程序牛顿拉夫逊法是计算潮流时常见的方法，该方法具有广泛的应用。该方法是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x)=0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x)=0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。相应程序详见附件1.2.1.1满发满载运行方式下计算在这种运行方式下，三台机组全部开启，各自保留10%作为功率裕度，重用负荷占发电量7%，负荷为最大负荷。此方式下，节点1为平衡节点，节点2、3、4、5、6、7为PQ节点，节点8、9、10为PV节点。经过Matlab的运算，我们得到如下结果：2.1.2满发轻载运行方式下计算经Matlab运算，我们得到：2.2PowerWorld仿真电路图2.2.1满发满载方式运行下电力系统的潮流计算满发满载方式运行下，发电厂三台机组均满载，取发电机容量10%为裕量，并保留出力的7%作厂用电，即发出功率为总容量的83.7%，负荷采用最大负荷计算。满发满载电路图：编辑模式下，#1为平衡节点，电压恒定为1.05VN。为了反映这一特性，#1应接一台容量无穷大的发电机。具体参数设置如图3所示（机组出力暂不确定，可填为100MW。）接在平衡节点处的发电机参数设置#7为PV节点，P=711.45MW，V=1.05VN，正常运行时各机组发出的有功功率P1=P2=251.1MW，P3=209.25MW。具体参数设置如图4、图5所示。接在#7处的发电机#8、#9参数设置接在#7处的发电机#10参数设置其余节点为PQ节点，将基准电压设置为230kV，带上相应大小的负荷。最后，用输电线路将各节点连接起来。潮流仿真结果如下：各节点电压节点编号基准电压kV标幺电压实际电压(kV)相角(度)12301.00000230.0000022300.9992229.816-1.7032301.0021230.483-1.8842300.9985229.655-2.0252300.9998229.954-2.08672302301.00411.0090230.943232.07-1.80-1.39各发电机运行状态节点编号发电机编号有功出力（MW）无功出力（MVar）电压标幺值11664.18-100.661.000088251.10139.611.050099251.10139.611.05001010209.25119.201.0500各支路状态首端节点编号末端节点编号首端有功功率（MW）首端无功功率（Mvar）首端视在功率（MVA）有功损耗（MW）无功损耗（Mvar）12332.1-50.3335.91.72-1.752343.3-60.874.60.05-7.4224172.1-7.8172.30.17-3.2836-61.8-65.990.30.04-3.284521.9-4751.80.01-4.155677-183.1-171.8-97.8-77.6207.6188.50.540.28-6.65-4.63电网线损率计算：线损率按以下公式计算：支路#1#2#2#3#2#4#3#6#4#5#5#7#6#7线损率0.518%0.115%0.099%0.065%0.046%0.295%0.163%电网网损率的计算：网损率按以下公式计算：其中，发电机输出功率包括发电厂和平衡节点的功率，发电机输出功率和负荷均只取有功分量。则发电机输出功率PD=711.45MW+664MW=1375.45MW总负荷PL=230MW+210MW+300MW+410MW+220MW=1370MW网损率=(PD-PL)/PD×100%=0.396%2.2.2满发轻载方式运行下电力系统的潮流计算满发轻载方式运行下，发电厂三台机组均满载，取发电机容量10%为裕量，并保留出力的7%作为厂用电，即发出功率为总容量的83.7%，负荷采用最小负荷计算。满发轻载电路图潮流仿真结果如下：各节点电压节点编号基准电压kV标幺电压实际电压(kV)相角(度)12301.00000230.0000022301.0056231.2880-0.4932301.0093232.1390-0.3942301.0061231.4030-0.5852301.0079231.8170-0.51672302301.01131.0157232.5990233.6110-0.260.17各发电机运行状态节点编号发电机编号有功出力（MW）无功出力（MVar）电压标幺值11167.57-169.221.000088251.10118.211.050099251.10118.211.05001010209.25100.861.0500各支路状态首端节点编号末端节点编号首端有功功率（MW）首端无功功率（Mvar）首端视在功率（MVA）有功损耗（MW）无功损耗（Mvar）1283.8-84.6119.10.2-10.4723-37.6-60.871.30.05-7.562447.6-26.454.40.02-4.2036-104.9-61.0121.30.07-3.1745-48.4-49.469.20.03-4.165677-179.6-175.3-80.4-67.4196.8187.80.480.28-7.15-4.75电网线损率计算：线损率按以下公式计算：支路#1#2#2#3#2#4#3#6#4#5#5#7#6#7线损率0.239%0.133%0.042%0.067%0.062%0.267%0.160%电网网损率的计算网损率按以下公式计算：其中，发电机输出功率包括发电厂和平衡节点的功率，发电机输出功率和负荷均只取有功分量。则发电机输出功率PD=711.45MW+167.571MW=879.021MW总负荷PL=147.2MW+134.4MW+192MW+262.4MW+140.8MW=876.8MW网损率=(PD-PL)/PD×100%=0.253%2.3Matlab和PowerWorld计算的结果比较以大方式的结果为例：电压标幺值很接近；电角度有一定差别，但一般小数点第一位都相同；线损有一定差别，但也都不大。造成结果差别的原因，首先，两种程序输入参数上存在误差，Matlab的参数会保留一定位小数，不能达到最精确值，而PowerWorld会将线路参数强制以美制单位保存，每次查看和修改参数自动转换成公制单位造成误差，而且这种误差还会随着查看与修改的次数累积；使用Matlab的算法较为理论和简单，系统中各个元件和部分的反馈和相互作用没什么体现，而PowerWorld仿真的算法更为精确和复杂，如发电机的模型更加详细。3、大方式下“N-1”潮流校核使用满发大负荷运行方式，即：三台机组保留10%容量备用，发出功率的7%为厂用电，出力为机组总容量的83.7%；负荷使用提供的220kV最大负荷。N-1情况下，将3发电机组及其变压器统一以#7节点代表进行计算。依次分别断开各条支路中双回线的一条线，计算节点电压和各个功率及损耗；根据设计要求：40°C长期运行允许的最大电流：1190A。3.1Matlab“N-1”潮流校核N-1故障可理解为双回线中的一回被切除，只剩一回在维持运行。这种情况会造成电力系统潮流状态的变更，需要进行潮流校核来检验系统的稳定性。此方式中，我们依次分别将其中一条双回线参数改成相应的单回线参数带入Matlab计算。经过Matlab运算，得到以下结果：N-1故障下各节点电压故障线路#1#2#3#4#5#6#7L121.0001.0111.0201.0111.0161.0261.039L231.0001.0071.0221.0081.0141.0271.039L241.0001.0081.0181.0091.0141.0241.037L361.0001.0071.0161.0081.0141.0261.039L451.0001.0081.0181.0081.0161.0241.039L571.0001.0051.0181.0031.0051.0251.041L671.0001.0061.0131.0081.0141.0181.040N-1故障各线路电流故障线路L12L23L24L36L45L57L67L122.840.961.370.671.162.091.7L233.650.571.590.841.032.611.39L243.690.751.660.71.022.731.22L363.651.031.350.91.082.651.32L453.650.841.350.631.252.471.49L573.650.781.970.951.212.211.97L673.650.981.041.061.082.951.47根据基准值的选取，可以得到题目已知的40℃长期允许最大电流1190A的标幺值为2.737。3.2PowerWorld“N-1”潮流校核断线电路图（以L12为例）L12断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.997851.000870.997150.998541.002891.00788有名值（kV）230229.506230.200229.344229.663230.665231.813功角（°）0-3.4-3.6-3.7-3.8-3.5-3.1一回线相电流线路12232436455767电流值（kA）1.68740.19230.43370.23070.13540.52440.4743分析：由于题目限定40°C长期运行允许最大电流是1190A，而表中12回线路的断线电流为1687.4A，已经超过极限，说明该线路载流能力较差，支路12不能承受断一回线路的电流，而对于各母线的电压水平均无大变化。L23断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.998731.003230.998190.999731.004951.00942有名值（kV）230229.709230.743229.583229.938231.140232.167功角（°）0-1.7-1.99-2.03-2.11-1.90-1.45一回线相电流线路12232436455767电流值（kA）0.84110.29450.45390.22280.15810.51480.4791L24断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.999121.001890.997710.999131.003801.00858有名值（kV）230229.798230.434229.473229.801230.873231.974功角（°）0-1.7-1.94-2.27-2.29-1.89-1.53一回线相电流线路12232436455767电流值（kA）0.8430.21990.7730.20030.11880.56290.4312L36断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.998941.001310.998340.999821.004831.00938有名值（kV）230229.756230.302229.618229.958231.110232.157功角（°）0-1.7-1.89-2.01-2.06-1.75-1.35一回线相电流线路12232436455767电流值（kA）0.84190.1810.41980.40550.13520.53830.4524L45断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.998991.002080.998021.000311.004131.00917有名值（kV）230229.769230.479229.544230.071230.951232.110功角（°）0-1.7-1.89-2.01-2.11-1.82-1.41一回线相电流线路12232436455767电流值（kA）0.8420.1990.42650.22960.22460.52080.4700L57断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.998081.002140.996230.996491.004631.01047有名值（kV）230229.557230.492229.133229.192231.065232.408功角（°）0-1.69-1.75-2.08-2.20-1.62-1.11一回线相电流线路12232436455767电流值（kA）0.83970.18700.53500.32120.16280.80950.5740L67断开一路：节点电压和功角节点1234567电压基准值（kV）230230230230230230230标幺值10.998541.000280.998140.999821.001731.00985有名值（kV）230229.664230.064229.572229.959230.398232.266功角（°）0-1.70-1.95-1.97-1.98-1.91-1.19一回线相电流线路12232436455767电流值（kA）0.84090.20740.3710.15650.13590.58550.79713.3Matlab与PowerWorld“N-1”潮流校核对比：通过就“N-1”潮流校核分别进行Matlab和PowerWorld计算、仿真，我们发现结果差距较大，其中PowerWorld仿真下只有L12电流值超过给定值，而Matlab计算下过流线路较多。经过反复计算、仿真，我们大胆认为中间存在的较大差异并非本身操作失误，而是软件在运算中机理不同而生：节点导纳矩阵形成过程不同。4、220kV母线的三相短路容量测算选用满发重荷的运行方式，平衡机基准容量为100MVA，发电机基准容量为各自的额定容量，变压器改为原来的标幺变比和标幺阻抗，对非平衡母线的所有220kv的母线进行三相短路测算。4.1短路容量的Matlab计算相应Matlab计算程序详见附件2.4.1.1满发满载运行方式下短路容量计算经过Matlab运算，分别得到短路电流和短路容量：节点#1#2#3#4#5#6#7电流13.78714.99514.99815.09315.13915.06915.291容量13901500105015101510151015304.1.2满发轻载运行方式下短路容量计算经过Matlab运算，分别得到短路电流和短路容量：节点#1#2#3#4#5#6#7电流9.69310.34310.37010.38710.43210.42710.591容量9701030104010401040104010604.2短路容量的PowerWorld仿真短路容量仿真操作4.2.1满发满载运行方式下短路容量计算经过PowerWorld仿真，分别得到短路电流和短路容量：节点#1#2#3#4#5#6#7容量1433.71475.51475.51476.41478.814.7951492.64.2.2满发轻载运行方式下短路容量计算经过PowerWorld仿真，分别得到短路电流和短路容量：节点#1#2#3#4#5#6#7容量961.5986.4989.5986.5990.4992.81001.64.3两种计算方法结果比较使用Matlab和PowerWorld进行计算和仿真的结果有一定相关性，但是存在一定差别。我们认为除了两者输入各种参数的误差之外，还可能是导纳矩阵算法之间的差距。PowerWorld形成的节点导纳矩阵与我们处理的不同，短路电流的计算准确。Matlab的算法是我们根据书本概念按自己的理解写成的程序代码，核心算法也遵照书本所述，可能是我们对负荷模型处理不够准确造成了这种差别。5、PowerWorld的网络稳定扫描如果短路点越靠近发电机节点，则短路故障对电力系统扰动越大。因此在以下稳定扫描仿真中，均把短路点尽可能靠近母线7上。若在此情况下，发电机摇摆曲线仍能保持暂态稳定，则电力系统稳定性符合要求。5.1暂态参数设置先按照要求添加发电机的暂态模型，设置好发电机的暂态参数。线路的零序阻抗参数。具体参数设置如下图：发电机#8、#9GENTRA模型参数发电机#10GENTRA模型参数发电机#8、#9、#10BPA_EG模型参数发电机升压变压器零序参数设定线路零序参数设定（以L67为例）5.2故障方案故障均设置在t=1.0s发生，仿真时间为20s。故障方案一：母线7三相短路，1.10s切除故障线路。故障方案一设计故障方案一运行从故障切除时间为0.1s的各发电机转子角摇摆曲线，我们可以看出个发电机的转子角由于故障而产生扰动，而在切除故障之后逐渐趋于稳定，各转子角的差（功角）趋于稳定，所以系统暂态稳定。故障极限切除时间求解：用二分法通过不断地设置切除时间，知道系统在故障切除之后不能保持暂态稳定性，到达临界值，从而确定大致的系统故障的极限切除时间。故障切除时：1.4218s2.4218s功角图2.4219s功角图故障方案二：线路6-7三相短路，1.20s切除故障线路。从故障切除时间为0.2s的各发电机转子角摇摆曲线，我们可以看出个发电机的转子角由于故障而产生扰动，而在切除故障之后逐渐趋于稳定，各转子角的差（功角）趋于稳定，所以系统暂态稳定。故障极限切除时间求解：用二分法通过不断地设置切除时间，知道系统在故障切除之后不能保持暂态稳定性，到达临界值，从而确定大致的系统故障的极限切除时间。故障切除时：1.4217s2.4217s功角图2.4218s功角图5.3方案总结由此可以得出故障的极限切除时间的大致区间为[1.4217，1.4218]，我们可以取故障的极限切除时间为：1.4217s。6、发电机组的负荷经济分配6.1遵守等微增率原则的经济分配方案F(P)=ai2p2+ai1p+ai0F1(P1)=0.00004P2+0.298P+10.22150≤P1≤300F2(P2)=0.00034P2+0.305P+10.32150≤P2≤300F3(P3)=0.00004P2+0.321P+9.38100≤P3≤250由题目可知，三个发电厂的微增率表达式如下所示：Λ1=0.00008P1+0.298Λ2=0.00006P1+0.305Λ3=0.00006P1+0.321令λ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.6.2使用Matlab函数进行优化的分配方案使用Matlab优化的过程中，我们按要求使用qua函数进行优化分配，这个函数专门适用于进行二次优化，使用广泛，准确度高。详情请见附件3.得出以下结果：①P1=1.7999×102P2=1.2×102P3=1.00×102②P1=300P2=200P3=100③P1=2.999999×102P2=2.99999×102P3=1.95×1027、课程设计总结关于电力系统的计算，很重要的一个步骤就是形成节点导纳矩阵。该矩阵的正确与否能很大程度上决定运算结果的正确性。在计算和仿真的过程中，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节点；%------------------------------------------------------------------------clearall;formatlong;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);%------------------------------------------------------------------fori=1:nl %支路数ifB1(i,6)==0 %左节点处于1侧p=B1(i,1);q=B1(i,2);else%左节点处于K侧p=B1(i,2);q=B1(i,1);endY(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); %分解出导纳阵的实部和虚部fori=1:n %给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV节点电压给定模值endfori=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(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；不满足收敛要求的节点数IT2whileIT2~=0%N0=2*n雅可比矩阵的阶数；N=N0+1扩展列IT2=0;a=a+1;fori=1:nifi~=isb %非平衡节点C(i)=0;D(i)=0;forj1=1:nC(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)endP1=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矩阵元素-----------------------------ifB2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1;%节点无功功率差%以上为除平衡节点外其它节点的功率计算--------------------------------------%求取Jacobi矩阵----------------------------------------------------------forj1=1:nifj1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); %dP/df=dQ/deX3=X2;%X2=dp/dfX3=dQ/deX4=-X1;%X1=dP/deX4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;%X3=dQ/deJ(p,N)=DQ节点无功功率差J(m,q)=X1;J(m,N)=DP;q=q+1;%X1=dP/deJ(m,N)=DP节点有功功率差J(p,q)=X4;J(m,q)=X2;%X4=dQ/dfX2=dp/dfelseifj1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%下面是针对PV节点来求取Jacobi矩阵的元素-----------------------------------------DP=P(i)-P1; %PV节点有功误差DV=V(i)^2-V2; %PV节点电压误差forj1=1:nifj1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i);%dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i);%dP/dfX5=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;elseifj1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/dfX5=-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;endendendendend%以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差---------------------------------------fork=3:N0%N0=2*n（从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N;%N=N0+1即N=2*n+1扩展列△P、△Q或△Ufork2=k1:N1 %从k+1列的Jacobi元素到扩展列的△P、△Q或△UJ(k,k2)=J(k,k2)./J(k,k);%用K行K列对角元素去除K行K列后的非对角元素进行规格化endJ(k,k)=1;%对角元规格化K行K列对角元素赋1 %回代运算-------------------------------------------------------------------ifk~=3%不是第三行k>3k4=k-1;fork3=3:k4 %用k3行从第三行开始到当前行的前一行k4行消去fork2=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列元素已消为0endifk==N0%若已到最后一行break;end%前代运算------------------------------------------------------------fork3=k1:N0%从k+1行到2*n最后一行fork2=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列元素已消为0endelse%是第三行k=3%第三行k=3的前代运算----------------------------------------------------fork3=k1:N0%从第四行到2n行（最后一行）fork2=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列元素已消为0endendend%上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵%-----------------------------------------------------------------------------------fork=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N);%修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N);%修改节点电压虚部end%修改节点电压---------------------------fork=3:N0DET=abs(J(k,N));ifDET>=pr%电压偏差量是否满足要求IT2=IT2+1;%不满足要求的节点数加1endendICT2(a)=IT2;%不满足要求的节点数ICT1=ICT1+1;%迭代次数end%用高斯消去法解"w=-J*V"disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);fork=1:nV(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);%显示各节点的电压相角forp=1:nC(p)=0;forq=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值endS(p)=E(p)*C(p);%计算各节点的功率S=电压X注入电流的共轭值enddisp('各节点的功率S为：');disp(S);%显示各节点的注入功率disp('-----------------------------------------------------');disp('各条支路的首端功率Si为：');fori=1:nlp=B1(i,1);q=B1(i,2);ifB1(i,6)==0Si(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);elseSi(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);enddisp(Si(p,q));SSi(p,q)=Si(p,q);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为：');fori=1:nlp=B1(i,1);q=B1(i,2);ifB1(i,6)==0Sj(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);elseSj(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);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为：');fori=1:nlp=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),')=',nu