智能科学与工程学院
《电力系统仿真实训》报告书
题目: 计算机编程实现系统潮流计算
专 业： 电气工程及其自动化
班 级： 0322405
学 号： 032240536
学生姓名： 黄卓
指导教师： 钟建伟
2024年 11 月 1 日
《电力系统仿真实训》任务书
| 学生姓名 | 黄卓 | 成绩 | |
| 设计题目 | 计算机编程实现系统潮流计算 | ||
| 设 计 技 术 参 数 | 线路每公里参数参阅《电力系统分析》书中给出的典型数据 线路假设 L1=40KM, L5=30KM r=0.45Ω/km,x=0.44Ω/km,b1=2.58*10^-6 S/km L2=140KM r=0.17Ω/km,x=0.409Ω/km,b=2.82*10^-6 S/km L3=80KM,L4=70KM r=0.27Ω/km,x=0.423Ω/km,b=2.69*10^-6 S/km 负荷假设 已知负荷的视在功率，负荷有功，无功量自行给定 | ||
| 设 计 要 求 | 输入数据后可以通过计算机实现系统潮流计算 | ||
| 工 作 量 | 字数6527 包含潮流计算的基本概念，潮流计算的程序实现，总结，参考文献和附录五个部分 | ||
| 工 作 计 划 |
| ||
| 参 考 资 料 | 参考文献
[3]《电力系统分析的计算机算法》邱晓燕 刘天琪 黄媛，中国电力出版社，2015. | ||
潮流计算是电力系统分析中的基本任务，用于确定网络中的电压分布和功率流。本文介绍了潮流计算的基本原理，并根据潮流计算的电路模型和数学模型进行分析，使用matlab平台得到一种潮流计算算法。本文阐述了牛顿-拉夫逊法进行潮流计算的基本原理并建立了数学模型，使用matlab实现该算法。该算法可以计算给定电力网络的潮流，包括电压幅值和相位角、线路功率流和节点功率平衡。
关键词：电力系统 潮流计算 牛顿-拉夫逊法 MATLAB
总结 16
附录 18
1.1 潮流计算的概述
电力系统潮流计算是电力系统最基本的计算，也是最重要的计算。所谓潮流计算，是指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电网中的分布。即根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。
电力系统潮流计算的结果是电力系统稳定计算和故障分析的基础。在电网规划阶段,通过潮流计算,合理规划电源容量及接入点,合理规划网架,选择无功补偿方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。对于正在运行的电力系统，通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提供原始数据。
总结为在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。
1.2 潮流计算的发展状况
电力系统潮流计算方法有两类，即手算潮流和计算机潮流计算。手算潮流主要借助于形成简化的等值电路来实现，这种方法只适用于规模不大的辐射型电力网潮流计算。计算机潮流计算的实现有两种途径：其一是编程实现网络方程的迭代求解；其二是借助于电力系统分析仿真软件，搭建系统模型来完成潮流计算，如MATLAB、PSASP、POWERWORLD、SIMULATOR等软件，这种计算机辅助潮流计算方法在较大规模的电力系统分析中应用广泛。
计算机潮流计算也分为离线计算和在线计算两种， 前者主要用于系统规划设计和安排系统的运行方式， 后者则用于正在运行系统的经常监视及实时控制。在线计算要求计算实时性号，计算速度快。
利用电子数字计算机进行电力系统潮流计算从 50 年代中期就已经开始。 在这 20 年内， 潮流计算曾采用了各种不同的方法， 这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。 对潮流计算的要求可以归纳为下面几点：
（1） 计算方法的可靠性或收敛性；
（2） 对计算机内存量的要求；
（3） 计算速度；
（4） 计算的方便性和灵活性。
电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题， 其解法都离不开迭代。 因此， 对潮流计算方法， 首先要求它能可靠地收敛， 并给出正确答案。 由于电力系统结构及参数的一些特点， 并且随着电力系统不断扩大， 潮流问题的方程式阶数越来越高， 对这样的方程式并不是任何数学方法都能保证给出正确答案的。这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。
在用数字计算机解电力系统潮流问题的开始阶段， 普遍采取以节点导纳矩阵为基础的逐次代入法。 这个方法的原理比较简单， 要求的数字计算机内存量比较小， 适应 50 年代电子计算机制造水平和当时电力系统理论水平。 但它的收敛性较差， 当系统规模变大时， 迭代次数急剧上升， 在计算中往往出现迭代不收敛的情况。 这就迫使电力系统计算人员转向以阻抗矩阵为基础的逐次代入法。阻抗法改善了系统潮流计算问题的收敛性， 解决了导纳法无法求解的一些系统的潮流计算， 在 60 年代获得了广泛的应用。
阻抗法的主要缺点是占用计算机内存大， 每次迭代的计算量大。克服阻抗法缺点的一条途径是采用牛顿-拉夫逊法。 这是数学中解决非线性方程式的典型方法， 有较好的收敛性。 在解决电力系统潮流计算问题时， 是以导纳矩阵为基础的， 因此， 只要我们能在迭代过程中尽可能保持方程式系数矩阵的稀疏性， 就可以大大提高牛顿法潮流程序的效率。
2.1 复杂电力系统潮流计算的数学模型
2.1.1 潮流计算的基本方程
潮流计算所用的电力网络由变压器、输电线路、电容器、电抗器等静止线性元件所构成，并用集中参数表示的串联或并联等支路来模拟。对这样的线性网络进行分析，一般采用的是节点法。以节点导纳矩阵表示的网络方程式为
其展开式为
(2-1)
式中：
为节点注入电流相量；
为节点电压相量；
为节点导纳矩阵；n为网络节点数。
节点电流可以用节点功率和电压表示
(2-2)
对于n个节点的网络，可列写2n个方程，但却有6n个变量。通常把负荷功率作为已知量，并把节点功率
和
引入网络方程。这样，n个节点电力系统的潮流方程的一般形式是
(2-3)
2.1.2 节点的分类
将上述方程的实部和虚部分开，对每一节点可得两个实数方程，但是变量仍有4个，即P、Q、U、δ。我们必须给定其中2个，只留下两个作为待求变量，方程组才可以求解。根据电力系统的实际运行条件，按给定变量的不同，一般将节点分为以下三种类型。
PQ节点：
这类节点的有功功率P和无功功率Q是已知的，节点电压幅值U和电压相角δ是待求量。通常变电所和发电厂是这一类型的节点。电力系统忠的绝大多数节点属于这一类型。浮游节点是一种特殊的PQ节点，其P、Q给定值为零。
PU节点：
这类节点的有功功率P和电压幅值U是已知的，待求量是无功功率Q和电压相角δ。这类节点需要有足够的可调无功容量，用来维持电压幅值稳定。这类节点在电力系统中数目很少。
平衡节点：
平衡节点又称slack节点，其已知量为电压幅值U和相角δ，待求量为有功功率P和无功功率Q。平衡节点只有一个，这个节点承担了系统的有功功率平衡。平衡节点也作基准节点，指定其相角为零，作为计算各节点电压相位的参考。
2.2 牛顿-拉夫逊法潮流计算
2.2.1 牛顿-拉夫逊法的基本原理
假设有n个联立的非线性代数方程
(2-4)
假设以给出各变量的初值
，令
分别为各变量的修正量， 使满足以上方程， 所以：
(2-5)
将上式中的 n 个多元函数在初始值附近分别展开成泰勒级数， 并略去含有
的二次及以上阶次的各项， 便得：
(2-6)
以上方程是对于修正量
的线性方程组，称为牛顿-拉夫逊法的修正方程， 可解出
。对初始近似解进行修正：
(2-7)
反复迭代， 在进行 k+1 次迭代时， 从求解修正方程式：
(2-8)
上式可以简写成
(2-9)
其中
称为雅克比矩阵。
得到修正量
，对各量进行修正
(i=1,2,…,n) 迭代过程一直进行到满足收敛判据
(2-10)
2.2.2 节点电压用极坐标表示时的牛顿-拉夫逊法潮流计算
采用极坐标时，节点电压表示为
(2-11)
节点功率方程(2-3)将写成
(2-12)
式中，
，是i、j两节点电压的相角差。
在有n个节点的系统中，假定第1~m号为PQ节点，第m+1~n-1号节点为PU节点，第n号节点为平衡节点。
和
是给定的，PU节点的电压幅值
~
也是给定的。因此，只剩下n-1个节点的电压相角
和m个节点的电压幅值
是未知量。
由PQ节点和PU节点共可以列写n-1个有功功率不平衡方程，式中
为已知功率。
(2-13)
由PQ节点共可以列写m个无功功率不平衡方程，式中
(2-14)
对于方程式(2-12)和(2-13)可以写出修正方程如下：
(2-15)
式中
；
；
；
；
；
为雅可比矩阵的分块形式。
使用牛顿-拉夫逊法计算潮流的程序框图如图(2-1)所示。首先要输入网络的原始数据，由此形成节点导纳矩阵。确定节点电压初值，迭代次数限制，迭代结果精度要求。迭代步骤如下：
(1)由上次迭代出的节点电压计算各类节点的不平衡量
、
和
。
(2)判断是否达到要求精度，达到则跳出循环。
(3)计算雅可比矩阵的各元素。
(4)解修正方程，得到修正量。
(5)修正各节点电压。
(6)迭代次数加1，返回第(1)步继续迭代。
若计算收敛，计算平衡节点功率及全部线路功率，否则返回“计算不收敛”。
牛顿-拉夫逊法潮流计算的程序框图如下：
图2-1
2.3 初始数据的输入
基准容量默认选取100MVA，要求输入的数据有母线基准电压、发电机参数、交流线参数、变压器参数、负荷参数、最大迭代次数、要求精度。除母线电压外，其余数据都采用标幺值。
发电机的参数包括：母线名、节点类型（节点类型PQ、PV、slack分别用1、2、3表示），若是PQ节点，则标明P、Q值，若是PV节点，则标明P、V值，若是slack节点，则标明U、
。
交流线参数包括：I侧母线、J侧母线、阻抗
、接地导纳
。交流线模型如下：
图2-2
变压器参数包括：I侧母线、J侧母线、阻抗
、变比k。
负荷参数包括：母线名、节点类型，若是PQ节点，则标明P、Q值，若是PV节点，则标明P、V值，若是slack节点，则标明U、
。
最后需要输入的是迭代次数限制Tmax和要求精度limit。
2.4 变压器π型等效
依照输入参数，变压器的模型为：
图2-3
由于理想变压器的存在，有两个电压等级，需要进行参数折算，不便计算，通常将变压器转成π型等效电路，使其模型与交流线相似。π型等效模型如下：
图2-4
由模型可知，需要根据参数
和k转成
、
、
。
进一步，
、
。
2.5 节点导纳矩阵的形成
有了交流线、变压器π型电路参数，就能进行节点导纳矩阵Y的计算。它的对角线元素
称为节点i的自导纳，其值等于接于节点i的所有支路导纳之和。
(2-16)
非对角线元素
称为节点i、j间的互导纳，它等于直接联接在节点i、j间的支路导纳的负值，
。节点导纳矩阵是稀疏对称矩阵，
，若节点i、j间不存在直接支路，则有
。
2.6 功率不平衡量的计算
由于功率不平衡量的计算需要调用多次，故写成子程序的形式。其输入参数有：n、m、P、Q、U、cita即
、G、B，输出参数有dP即
、dQ即
、Pi、Qi。计算公式为
(2-17)
2.7 雅可比矩阵的计算
为使表达式更加整齐，计算更加方便清晰，雅可比矩阵J采用矩阵分块的形式
。式中H为(n-1)×(n-1)阶方阵，其元素
；N是(n-1)×m阶矩阵，其元素
；K是m×(n-1)阶矩阵，其元素
；L是m×m阶矩阵，其元素
。
在这里把节点不平衡功率对节点电压幅值的偏导数都乘以该节点电压，相应地把节点电压的修正量都除以该节点的电压幅值，这样，雅克比矩阵元素的表达式就具有比较整齐的形式。
当
时，有
(2-18)
当
时，
(2-19)
2.8 求解节点电压修正量
根据修正方程(2-15)，可解出修正量
(2-20)
进一步
。
节点有n个，电压相角修正量有n-1个，电压幅值修正量有m个，PV节点无电压幅值修正量，slack节点无电压幅值和相角修正量，因此这些修正量均以0补上。
由此，可以对节点电压幅值和电压相角进行修正，得到新的节点电压，用于下一轮迭代。
2.9 迭代过程
根据第2.2.2节所述，迭代过程共有6步：
(1)由上次迭代出的节点电压计算各类节点的不平衡量
、
。
(2)判断是否达到要求精度，达到则跳出循环。
(3)计算雅可比矩阵的各元素。
(4)解修正方程，得到修正量。
(5)修正各节点电压。
(6)迭代次数加1，返回第(1)步继续迭代。
在程序中采用for循环，循环次数为设定的迭代次数上限Tmax，满足精度要求则可跳出循环，计算功率不平衡量、雅可比矩阵、节点电压修正量均调用子程序。迭代结束后，再次判断是否达到精度要求来确定迭代的收敛性。
是否达到精度的判断条件为
。
2.10 平衡节点及各支路功率计算
节点功率在计算功率不平衡量的时候已经计算过了，公式详见第2.6节。有功功率储存在
，无功功率储存在
，可直接输出。支路功率损耗计算公式为
(2-21)
总结
经过两周的课程设计，我对电力系统潮流计算有了非常深刻的认识。我通过手算就能明显的感觉到潮流计算的复杂性，因为节点越多，方程的数量也越多，对应矩阵的阶数也越多，手算只能解决简单的网络，对于一般化的网络，使用程序计算是必然的选择，因此编写一个潮流计算程序是非常实用的。在这两周里，我一边查阅书籍复习有关潮流计算的知识，一边借助网络资源，突击matlab编程，通过一个项目来学习一门语言是非常高效的，两周时间足以掌握那些基础的操作，编写出实现潮流计算的程序。
我编写的这一程序可以在给定各个标幺值参数的情况下，计算潮流，其参数输入条理清晰，仿照PSASP程序设计，计算收敛速度快，具有普遍适用性。这个程序需要事先计算出标幺值；迭代过程中未对电压、功率越限进行分析。
参考文献
- 《电力系统分析（第四版）》 刘天琪，邱晓燕编著．—北京：科学出版社，2005
- 《MATLAB/Simulink电力系统建模与仿真（第2版）》 于群 曹娜，机械工业出版社，2017.
- 《电力系统分析的计算机算法》邱晓燕 刘天琪 黄媛，中国电力出版社，2015.
附录
附录1 Y矩阵
| % //**** 参数数据 ******* // parameter data % ********发电机参数******* G parameter %发电机1 g1=4*15; x_g1_d=0.136; x_g1_x2=0.16; x_g1_0=0.073; x_g1_cos=0.8; %发电机2 g2=3*12; x_g2_d=0.136; x_g2_x2=0.161; x_g2_0=0.075; x_g2_cos=0.8; %发电机3 g3=4*63; x_g3_d=0.134; x_g3_x2=0.161; x_g3_0=0.06; x_g3_cos=0.85; %发电机4 g4=1*50; x_g4_d=0.138; x_g4_x2=0.154; x_g4_0=0.054; x_g4_cos=0.85; %发电机5 g5=4*63; x_g5_d=0.128; x_g5_x2=0.157; x_g5_0=0.0591; x_g5_cos=0.8; %************潮流参数************ power flow parameter %第一段潮流 s1=20; %第二段潮流 s2=25; %第三段潮流 s3=10; %第四段潮流 s4=30; %第五段潮流 s5=25; %第六段潮流 s6=120; %第七段潮流 s7=35; %第八段潮流 s8=80; %第九段潮流 s9=5; %第十段潮流 s10=15; %**************变压器参数********** transformer parameter %变压器1 10/110 deltap0_q1 =18.6; deltaps_q1 =89; i0_q1=0.5;% I0代表IO% vs_q1=10.5;% vs代表vs% sn_q1=20; vn_q1=110; %变压器2 110/10 deltap0_q2 =15.7; deltaps_q2 =73; i0_q2=0.5;% I0代表IO% vs_q2=10.5;% vs代表vs% sn_q2=16; vn_q2=110; %变压器3 10/110 deltap0_q3 =15.7; deltaps_q3 =73; i0_q3=0.5;% I0代表IO% vs_q3=10.5;% vs代表vs% sn_q3=16; vn_q3=110; %变压器4 10/110 deltap0_q4 =44; deltaps_q4 =7121; i0_q4=0.35;% I0代表IO% vs_q4=10.5;% vs代表vs% sn_q4=63; vn_q4=110; %变压器5 三绕组变压器 升压 10/110 10-s8 10-s9 deltap0_q5 =13.2; deltaps_q5 =63; i0_q5=0.55;% I0代表IO% vs_1_q5=10.5; vs_2_q5=6.5; vs_3_q5=17.5; sn_q5=10; vn_q5=110; %变压器6 deltap0_q6 =11; deltaps_q6 =50; i0_q6=0.55;% I0代表IO% vs_q6=10.5;% vs代表vs% sn_q6=10; vn_q6=110; %************线路参数************* line data %线路1 2 3 4 5 line1=40;%这里是L1 line2=140; line3=80; line4=70; line5=30; %/ end end end end /// 参数数据 end of parameter data %///***** 线路假设条件 *****/// % L1=40KM, L5=30KM % r=0.45 ,x=0.44,b1=2.58*10^-6 % L2=140KM % r=0.17Ω/km,x=0.409Ω/km,b=2.82*10^-6 S/km % L3=80KM,L4=70KM % R=0.27 X=0.423 b=2.69*10^-6 %/// 计算后得到的 R X G B Y 如下 %变压器 Q1=[calc_R_compute(deltaps_q1, vn_q1, sn_q1),calc_X_compute(vs_q1 , vn_q1, sn_q1),calc_G_compute(deltap0_q1,vn_q1),calc_B_compute(i0_q1,sn_q1,vn_q1)]; Q2=[calc_R_compute(deltaps_q2, vn_q2, sn_q2),calc_X_compute(vs_q2 , vn_q2, sn_q2),calc_G_compute(deltap0_q2,vn_q2),calc_B_compute(i0_q2,sn_q2,vn_q2)]; Q3=[calc_R_compute(deltaps_q3, vn_q3, sn_q3),calc_X_compute(vs_q3 , vn_q3, sn_q3),calc_G_compute(deltap0_q3,vn_q3),calc_B_compute(i0_q3,sn_q3,vn_q3)]; Q4=[calc_R_compute(deltaps_q4, vn_q4, sn_q4),calc_X_compute(vs_q4 , vn_q4, sn_q4),calc_G_compute(deltap0_q4,vn_q4),calc_B_compute(i0_q4,sn_q4,vn_q4)]; Q5=[calc_R_compute(deltaps_q5, vn_q5, sn_q5),calc_X_compute(vs_1_q5 , vn_q5, sn_q5),calc_G_compute(deltap0_q5,vn_q5),calc_B_compute(i0_q5,sn_q5,vn_q5);calc_R_compute(deltaps_q5, vn_q5, sn_q5),calc_X_compute(vs_2_q5 , vn_q5, sn_q5),calc_G_compute(deltap0_q5,vn_q5),calc_B_compute(i0_q5,sn_q5,vn_q5);calc_R_compute(deltaps_q5, vn_q5, sn_q5),calc_X_compute(vs_3_q5 , vn_q5, sn_q5),calc_G_compute(deltap0_q5,vn_q5),calc_B_compute(i0_q5,sn_q5,vn_q5)]; Q6=[calc_R_compute(deltaps_q6, vn_q6, sn_q6),calc_X_compute(vs_q6 , vn_q6, sn_q6),calc_G_compute(deltap0_q6,vn_q6),calc_B_compute(i0_q6,sn_q6,vn_q6)]; %线路的 R,X,B 用阻抗Z表示 L1=complex(0.45*line1,0.44*line1); L3=complex(0.17*line3,0.409*line3); L4=complex(0.17*line4,0.409*line4); L5=complex(0.45*line5,0.44*line5); % 线路L2采用 Z和Y的表达形式 L2=[complex(0.17*line2,0.409*line2),2.82*10^-6*line2]; % ***** 根据上述数据选定节点后生成 Y矩阵 **** %按需要选取九个结点 Z1=complex(Q1(1,1),Q1(1,2)) Z2=complex(Q2(1,1),Q2(1,2))+L1 Z3=L2(1,1) Z4=L3 Z5=complex(Q3(1,1),Q3(1,2)) Z6=L4 Z7=complex(Q4(1,1),Q4(1,2)) Z8=complex(Q5(1,1)+Q5(2,1)+Q5(3,1),Q5(1,2)+Q5(2,2)+Q5(3,2)) Z9=complex(Q6(1,1),Q6(1,2))+L5 %转化为Y矩阵 Z_single=[Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8,Z9] Z=[Z1+Z2+Z3 -Z2 -Z3 0 0 0 0 0 ; -Z2 Z2 0 0 0 0 0 0 ; -Z3 0 Z3+Z4+Z5 -Z5 0 Z4 0 0 ; 0 0 -Z5 Z5 0 0 0 0 ; 0 0 0 0 Z7 -Z7 0 0 ; 0 0 -Z4 0 -Z7 Z4+Z7+Z8+Z9 -Z8 -Z9; 0 0 0 0 0 -Z8 Z8 0; 0 0 0 0 0 -Z9 0 Z9]; Y=1./Z %****** /// 结点的选取 /// ****** %PQ 1 5 7 %PV 2 3 6 8 %平衡 4 %******/// 结点数据 /// ***** %结点理想电压 ulx=[110 10 110 10 10 110 10 110] PQ_P=[calc_P(g1,s1,x_g1_cos,ulx(1,1),Y(1,1)) calc_P(0,s2,x_g2_cos,ulx(1,2),Y(2,2)) calc_P(g3,s3,x_g3_cos,ulx(1,3),Y(3,3)) calc_P(g2,s4,x_g4_cos,ulx(1,4),Y(4,4)) calc_P(g4,s5,x_g5_cos,ulx(1,5),Y(5,5)) calc_P(0,s6,0,ulx(1,6),Y(6,6)) calc_P(g5,s7,x_g5_cos,ulx(1,7),Y(7,7)) calc_P(0,s8,0,ulx(1,8),Y(8,8)) ] PQ_Q=PQ_P' function P_compute =calc_P(g,s,cos,ulxdy,ydy) P_compute =g.*cos-s+sqrt(1-cos.*cos) -ulxdy.*ulxdy.*ydy; end function deltaP=calc_deltaP(g,s,cos,ulxdy,ydy,utzdy) deltaP=calc_p(g,s,cos,ulxdy,ydy)-utzdy.*utzdy.*ydy; end %///***** 等效电路参数计算 *****/// Equivalent parameters compute %******** R ******* function R_compute = calc_R_compute(deltaps, vn, sn) R_compute = deltaps * vn^2 ./ (1000 * sn^2); end %******* X ********* function X_compute = calc_X_compute(vs, vn, sn) X_compute =vs*vn*vn./(100*sn); end %******* G ********* function G_compute = calc_G_compute(deltap0,vn) G_compute =deltap0./(1000*vn*vn); end %******* B ********* function B_compute =calc_B_compute(i0, sn ,vn) B_compute = i0*sn./(100*vn*vn); end |
附录2 雅各比矩阵
| % function: 计算雅可比矩阵 function [ J ] = Jacobi( n,m,U,cita,B,G,Pi,Qi ) %雅可比矩阵的计算 %分块 H N K L %i!=j时 for i=1:n-1 for j=1:n-1 H(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j))); end end for i=1:n-1 for j=1:m N(i,j)=-U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j))); end end for i=1:m for j=1:n-1 K(i,j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j))); end end for i=1:m for j=1:m L(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j))); end end %i==j时 for i=1:n-1 H(i,i)=U(i).^2*B(i,i)+Qi(i); end for i=1:m N(i,i)=-U(i).^2*G(i,i)-Pi(i); end for i=1:m K(i,i)=U(i).^2*G(i,i)-Pi(i); end for i=1:m L(i,i)=U(i).^2*B(i,i)-Qi(i); end %合成雅可比矩阵 J=[H N;K L]; end |
附录3 功率不平衡量
| % function: 计算功率不平衡量 function [ dP,dQ,Pi,Qi ] = Unbalanced( n,m,P,Q,U,G,B,cita ) %计算ΔPi有功的不平衡量 for i=1:n for j=1:n Pn(j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j))); end Pi(i)=sum(Pn); end dP=P(1:n-1)-Pi(1:n-1); %dP有n-1个 %计算ΔQi无功的不平衡量 for i=1:n for j=1:n Qn(j)=U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j))); end Qi(i)=sum(Qn); end dQ=Q(1:m)-Qi(1:m); %dQ有m个 end%func |
附录4 PQ分解法计算修正量并求节点电压修正量
| % function:使用PQ分解法计算电压修正量 function [ dU,dcita ] = PQ_LJ( n,m,dP,dQ,U,B ) dP_U=dP./U(1:n-1); dQ_U=dQ./U(1:m); dUdcita=(-inv(B(1:n-1,1:n-1))*dP_U')'; dcita=dUdcita./U(1:n-1); dU=(-inv(B(1:m,1:m))*dQ_U')'; dU=[dU zeros(1,n-m)]; dcita=[dcita 0];%补零 End % function:修正节点电压 function [ dU,dcita ] = Correct( n,m,U,dP,dQ,J ) %求解节点电压修正量 for i=1:m Ud2(i,i)=U(i); end dPQ=[dP dQ]'; dUcita=(-inv(J)*dPQ)'; dcita=dUcita(1:n-1); dcita=[dcita 0]; dU=(Ud2*dUcita(n:n+m-1)')'; dU=[dU zeros(1,n-m)]; end |
扫一扫
7465
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
