本发明公开一种电力系统随机稳定性分析方法，包括包括步骤：S100确定系统参数；S200求取发电机节点无功功率对电压幅值的导数，建立发电机判据列表；S300计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率；S400通过步骤S200和S300，计算系统中每个发电机失去静态稳定性的概率；S500根据对每个发电机失去静态稳定性的概率的分析，得到系统失去静态稳定性的概率。本发明以无功功率对电压幅值的导数作为判据，与蒙特卡洛法相结合判断系统的静态稳定性；解决现有系统稳定性判据阈值难以直接得到，当系统中存在大量随机因素时，解析法和近似法建模困难，无法全面考虑随机因素对系统稳定性的影响。
1.一种电力系统随机稳定性分析方法，其特征在于，包括步骤：
S100确定系统参数，生成节点信息矩阵、发电机信息矩阵和支路信息矩阵；
S200求取发电机节点无功功率对电压幅值的导数，建立发电机判据列表，包括步骤：S201设置计算次数n_monte，迭代初值q＝1；
S202根据发电机停运概率生成发电机工作状态；根据每台发电机停运概率，抽取所有发电机的工作状态，组成发电机工作状态矩阵；根据发电机工作状态矩阵，修改发电机信息矩阵；
S203设置系统的负荷为正态分布；设置正态分布的期望和方差，方差值取均值的0.1，抽取以上均值和方差服从正态分布的随机数；将产生的随机数对应填入节点信息矩阵的Pd和Qd处；
S204根据支路信息矩阵生成节点导纳矩阵，结合节点电压初值、以及步骤S202和步骤S203中得到的各发电机和负荷的功率，根据牛顿拉夫逊算法计算最后一次迭代得到雅可比矩阵；
S205根据所述雅可比矩阵计算发电机节点的无功功率对电压幅值的导数；
S206判断计算次数q<＝n_monte时，则返回步骤S202重新产生发电机和负荷的功率，进行下一次牛顿拉夫逊潮流计算，且q＝q+1；若q>n_monte，则计算结束；
S207当计算结束后，将得到计算次数为n_monte的计算结果，建立发电机判据列表，其行数为n_monte，列数为发电机个数；
S300计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率；
S400通过步骤S200和S300，计算系统中每个发电机失去静态稳定性的概率；
S500根据对每个发电机失去静态稳定性的概率的分析，得到系统失去静态稳定性的概率。
2.根据权利要求1所述的一种电力系统随机稳定性分析方法，其特征在于，步骤S100中，确所述定系统参数，包括步骤：S101确定系统中发电机的额定功率、无功功率和停运概率，确定系统中负荷的有功功率、无功功率的期望和标准差；
S102确定系统网络结构以及各支路的阻抗大小，确定各节点的节点类型和节点电压初值；
S103根据步骤S101和S102中确定的参数值，生成节点信息矩阵、发电机信息矩阵和支路信息矩阵。
3.根据权利要求2所述的一种电力系统随机稳定性分析方法，其特征在于，所述节点信息矩阵格式为：[bus_i,type,Pd,Qd,Vm,Va,baseKV]；
其中，bus_i为节点编号；type为节点类型，PQ节点标记为1，PV节点标记为2，平衡节点标记为3；Pd为节点负荷有功功率，Qd为节点负荷无功功率；Vm为节点电压幅值初值，Va为节点电压相角初值；baseKV为该节点电压基准值，节点信息矩阵的行数为节点的数目。
4.根据权利要求3所述的一种电力系统随机稳定性分析方法，其特征在于，所述发电机信息矩阵格式如下：[Bus_i,Pg,Qg,Qmax,Qmin,Pmax,Pmin,Vg,mBase,status_f]其中，Bus_i为发电机接入节点的编号；Pg为发电机额定功率，Qg为发电机无功功率；
Qmin为发电机输出无功功率的最小值，Qmax为发电机输出无功功率的最大值，Pmin为发电机输出有功功率的最小值，Pmax为发电机输出有功功率最大值；Vg为发电机端节点电压初值；mBase为系统基准容量值；status_f为发电机工作状态，发电机工作标记为1，发电机停运标记为0；发电机信息矩阵行数为发电机的数目。
5.根据权利要求4所述的一种电力系统随机稳定性分析方法，其特征在于，所述支路信息矩阵格式如下：[fbus,tbus,r,x,b,ratio,status_b]
其中，fbus为支路起始节点编号，tbus为支路末端节点；r为支路电阻，x为支路电抗，b为支路对地导纳；ratio为支路中变压器变比；status_b为支路状态，支路闭合为1，支路断开为0。
6.根据权利要求5所述的一种电力系统随机稳定性分析方法，其特征在于，在步骤S202中，根据发电机停运概率生成发电机工作状态时：发电机工作标记为1，发电机停运标记为
0，发电机运行状态服从停运概率为sd的两点分布，根据sd抽取发电机工作状态；根据每台发电机停运概率，抽取所有发电机的工作状态，组成发电机工作状态矩阵；
在步骤S202中，根据发电机工作状态矩阵，修改发电机信息矩阵时：工作状态为0的发电机，Pg和Qg设置为0，Vg初值设置为1，工作状态status_f标记为0，节点信息矩阵发电机接入节点类型设置为PQ节点；工作状态为1的发电机，Pg和Qg设置功率值，Vg设置发电机电势，工作状态status_f标记为1，节点信息矩阵发电机接入节点类型设置为PV节点。
7.根据权利要求6所述的一种电力系统随机稳定性分析方法，其特征在于，步骤S300中，计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率，包括步骤：选取发电机判据列表中发电机列元素，统计发电机列元素小于0的次数为m1，则发电机失去静态稳定性的概率为：
8.根据权利要求7所述的一种电力系统随机稳定性分析方法，其特征在于，步骤S500中，选择所有发电机的失去静态稳定性的概率中最大的概率，作为系统失去静态稳定性的概率。
一种电力系统随机稳定性分析方法
技术领域
[0001] 本发明属于电力系统技术领域，特别是涉及一种电力系统随机稳定性分析方法。
背景技术
[0002] 系统的静态安全评估保证系统安全可靠运行的必不可少的工作，随着新能源的大规模接入，随机静态安全评估的概念被引入。目前的随机静态安全评估算法都是直接根据概率潮流的计算结果分析节点电压越限或支路潮流过载的概率。但是，对于不同结构网络，其节点电压和支路潮流阈值是不同的，这些判据需要经过对多种运行方式进行计算才能得到。因此，寻找一种计算量小、可靠性高、适用性强的分析电力系统随机静态安全评估方法具有非常重要的研究价值。
[0003] 电力系统静态稳定是限制输电系统输送功率的一个重要因素。目前随机静态安全评估的主要判据如下：节点电压、支路潮流和支路两端相角差。但是，在实际的电力系统中，每个节点电压越限或每条支路潮流过载的的阈值需要进行多种方式的计算，迫使运行条件不断恶化才能得到，且不同大小规模的电网其阈值存在较大的差别。当电力系统的结构不同时，对失去静态稳定的最大相角差的设定有很大影响。因此，在实际分析电力系统稳定性的问题时，以上判据是很难得到的，直接根据经验规定电力系统失去静态稳定性的阈值会对计算结果产生不可预知的影响。
[0004] 由于系统中存在大量的随机因素，因此概率潮流的研究方法被广泛应用于系统的安全性分析。电力系统概率潮流方法主要分为以下三种：模拟法、近似法和解析法。点估计法是近似法的代表，但是节点电压或支路潮流的高阶矩存在较大的误差，无法得到准确的概率分布。半不变量法是解析法的代表，但是目前半不变量法在考虑存在大规模随机因素的问题时，存在较大的误差。同时，随着系统规模的增大和结构的复杂，采用近似法和解析法很难建立准确的数学模型来模拟变量之间的复杂关系。
发明内容
[0005] 为了解决上述问题，本发明提出了一种电力系统随机稳定性分析方法，以无功功率对电压幅值的导数作为判据，与蒙特卡洛法相结合判断系统的静态稳定性。解决现有系统稳定性判据阈值难以直接得到，当系统中存在大量随机因素时，解析法和近似法建模困难，无法全面考虑随机因素对系统稳定性的影响。
[0006] 为达到上述目的，本发明采用的技术方案是：一种电力系统随机稳定性分析方法，包括步骤：
[0007] S100确定系统参数；
[0008] S200求取发电机节点无功功率对电压幅值的导数，建立发电机判据列表；
[0009] S300计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率；
[0010] S400通过步骤S200和S300，计算系统中每个发电机失去静态稳定性的概率；
[0011] S500根据对每个发电机失去静态稳定性的概率的分析，得到系统失去静态稳定性的概率。
[0012] 进一步的是，步骤S100中，确所述定系统参数，包括步骤：
[0013] S101确定系统中发电机的额定功率、无功功率和停运概率，确定系统中负荷的有功功率、无功功率的期望和标准差；
[0014] S102确定系统网络结构以及各支路的阻抗大小，确定各节点的节点类型和节点电压初值；
[0015] S103根据步骤S101和S102中确定的参数值，生成节点信息矩阵、发电机信息矩阵和支路信息矩阵。
[0016] 进一步的是，所述节点信息矩阵格式为：
[0017] [bus_i,type,Pd,Qd,Vm,Va,baseKV]；
[0018] 其中，bus_i为节点编号；type为节点类型，PQ节点标记为1，PV节点标记为2，平衡节点标记为3；Pd为节点负荷有功功率，Qd为节点负荷无功功率；Vm为节点电压幅值初值，Va为节点电压相角初值；baseKV为该节点电压基准值，节点信息矩阵的行数为节点的数目。
[0019] 进一步的是，所述发电机信息矩阵格式如下：
[0020] [Bus_i,Pg,Qg,Qmax,Qmin,Pmax,Pmin,Vg,mBase,status_f]
[0021] 其中，Bus_i为发电机接入节点的编号；Pg为发电机额定功率，Qg为发电机无功功率；Qmin为发电机输出无功功率的最小值，Qmax为发电机输出无功功率的最大值，Pmin为发电机输出有功功率的最小值，Pmax为发电机输出有功功率最大值；Vg为发电机端节点电压初值；mBase为系统基准容量值；status_f为发电机工作状态，发电机工作标记为1，发电机停运标记为0；发电机信息矩阵行数为发电机的数目。
[0022] 进一步的是，所述支路信息矩阵格式如下：
[0023] [fbus,tbus,r,x,b,ratio,status_b]
[0024] 其中，fbus为支路起始节点编号，tbus为支路末端节点；r为支路电阻，x为支路电抗，b为支路对地导纳；ratio为支路中变压器变比；status_b为支路状态，支路闭合为1，支路断开为0。
[0025] 进一步的是，步骤S200中，求取发电机节点无功功率对电压幅值的导数，建立发电机判据列表，包括步骤：
[0026] S201设置计算次数n_monte，迭代初值q＝1；
[0027] S202根据发电机停运概率生成发电机工作状态，发电机工作标记为1，发电机停运标记为0，发电机运行状态服从停运概率为sd的两点分布，根据sd抽取发电机工作状态；根据每台发电机停运概率，抽取所有发电机的工作状态，组成发电机工作状态矩阵；
[0028] 根据发电机工作状态矩阵，修改发电机信息矩阵；其中，工作状态为0的发电机，Pg和Qg设置为0，Vg初值设置为1，工作状态status_f标记为0，节点信息矩阵发电机接入节点类型设置为PQ节点；工作状态为1的发电机，Pg和Qg设置功率值，Vg设置发电机电势，工作状态status_f标记为1，节点信息矩阵发电机接入节点类型设置为PV节点；
[0029] S203设置系统的负荷为正态分布；设置正态分布的期望和方差，方差值取均值的0.1，抽取以上均值和方差服从正态分布的随机数；将产生的随机数对应填入节点信息矩阵的Pd和Qd处；
[0030] S204根据支路信息矩阵生成节点导纳矩阵，结合节点电压初值、以及步骤S202和步骤S203中得到的各发电机和负荷的功率，根据牛顿拉夫逊算法计算最后一次迭代得到雅可比矩阵；
[0031] S205根据所述雅可比矩阵计算发电机节点的无功功率对电压幅值的导数；雅可比矩阵中提取相应的元素，这些元素就是相关“发电机节点的无功功率对电压幅值的导数”；
[0032] S206判断计算次数q<＝n_monte时，则返回步骤S202重新产生发电机和负荷的功率，进行下一次牛顿拉夫逊潮流计算，且q＝q+1；若q>n_monte，则计算结束；
[0033] S207当计算结束后，将得到计算次数为n_monte的计算结果，建立发电机判据列表，其行数为n_monte，列数为发电机个数。
[0034] 进一步的是，步骤S300中，计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率，包括步骤：
[0035] 选取发电机判据列表中发电机列元素，统计发电机列元素小于0的次数为m1，则发电机失去静态稳定性的概率为：
[0036]
[0037] 进一步的是，选择所有发电机的失去静态稳定性的概率中最大的概率，作为系统失去静态稳定性的概率。
[0038] 采用本技术方案的有益效果：
[0039] 本发明提出了以无功功率对电压幅值的导数为判据，避免了传统稳定性分析中难以计算节点电压阈值、支路潮流阈值和支路最大相角差的问题，同时采用蒙特卡洛法准确计算整个系统失去静态稳定性的概率，避免了近似法和模拟法概率潮流计算的误差；静态稳定性是离线计算的，所提方法能够满足系统静态稳定性计算精度和效率的要求。
附图说明
[0040] 图1为本发明的一种电力系统随机稳定性分析方法流程示意图；
[0041] 图2为本发明实施例中标准IEEE-14节点电路图。
具体实施方式
[0042] 为了使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明作进一步阐述。
[0043] 在本实施例中，参见图1所示，本发明提出了一种电力系统随机稳定性分析方法，包括步骤：S100-S500。
[0044] S100确定系统参数；
[0045] 确所述定系统参数，包括步骤：
[0046] S101确定系统中发电机的额定功率、无功功率和停运概率，确定系统中负荷的有功功率、无功功率的期望和标准差；
[0047] S102确定系统网络结构以及各支路的阻抗大小，确定各节点的节点类型和节点电压初值；
[0048] S103根据步骤S101和S102中确定的参数值，生成节点信息矩阵、发电机信息矩阵和支路信息矩阵。
[0049] 所述节点信息矩阵格式为：
[0050] [bus_i,type,Pd,Qd,Vm,Va,baseKV]；
[0051] 其中，bus_i为节点编号；type为节点类型，PQ节点标记为1，PV节点标记为2，平衡节点标记为3；Pd为节点负荷有功功率，Qd为节点负荷无功功率；Vm为节点电压幅值初值，Va为节点电压相角初值；baseKV为该节点电压基准值，节点信息矩阵的行数为节点的数目。
[0052] 所述发电机信息矩阵格式如下：
[0053] [Bus_i,Pg,Qg,Qmax,Qmin,Pmax,Pmin,Vg,mBase,status_f]
[0054] 其中，Bus_i为发电机接入节点的编号；Pg为发电机额定功率，Qg为发电机无功功率；Qmin为发电机输出无功功率的最小值，Qmax为发电机输出无功功率的最大值，Pmin为发电机输出有功功率的最小值，Pmax为发电机输出有功功率最大值；Vg为发电机端节点电压初值；mBase为系统基准容量值；status_f为发电机工作状态，发电机工作标记为1，发电机停运标记为0；发电机信息矩阵行数为发电机的数目。
[0055] 所述支路信息矩阵格式如下：
[0056] [fbus,tbus,r,x,b,ratio,status_b]
[0057] 其中，fbus为支路起始节点编号，tbus为支路末端节点；r为支路电阻，x为支路电抗，b为支路对地导纳；ratio为支路中变压器变比；status_b为支路状态，支路闭合为1，支路断开为0。
[0058] S200求取发电机节点无功功率对电压幅值的导数，建立发电机判据列表；
[0059] 求取发电机接点无功功率对电压幅值的导数，建立发电机判据列表，包括步骤：
[0060] S201设置计算次数n_monte，迭代初值q＝1；
[0061] S202根据发电机停运概率生成发电机工作状态，发电机工作标记为1，发电机停运标记为0，发电机运行状态服从停运概率为sd的两点分布，根据sd抽取发电机工作状态；根据每台发电机停运概率，抽取所有发电机的工作状态，组成发电机工作状态矩阵；
[0062] 根据发电机工作状态矩阵，修改发电机信息矩阵；其中，工作状态为0的发电机，Pg和Qg设置为0，Vg初值设置为1，工作状态status_f标记为0，节点信息矩阵发电机接入节点类型设置为PQ节点；工作状态为1的发电机，Pg和Qg设置功率值，Vg设置发电机电势，工作状态status_f标记为1，节点信息矩阵发电机接入节点类型设置为PV节点；
[0063] S203设置系统的负荷为正态分布；设置正态分布的期望和方差，方差值取均值的0.1，抽取以上均值和方差服从正态分布的随机数；将产生的随机数对应填入节点信息矩阵的Pd和Qd处；
[0064] S204根据支路信息矩阵生成节点导纳矩阵，结合节点电压初值、以及步骤S202和步骤S203中得到的各发电机和负荷的功率，根据牛顿拉夫逊算法计算最后一次迭代得到雅可比矩阵。
[0065] 雅可比矩阵元素如下：
[0066]
[0067]
[0068]
[0069] 其中Vi、Vj分别为节点i和j的电压幅值，θi、θj为节点i和j的电压相角，Pi、Qi为节点i注入的有功和无功功率，Gij、Bij分别为节点导纳矩阵中第i行j列元素的实部和虚部，δij为从节点i到节点j间支路的相角差。
[0070] 根据以上公式可以计算得到牛顿拉夫逊潮流算法最后一次迭代的雅可比矩阵ΔJ。
[0071]
[0072] H、N、J、L分别为由Hij、Nij、Jij、Lij组成的矩阵。
[0073] S205根据所述雅可比矩阵计算发电机节点的无功功率对电压幅值的导数；雅可比矩阵中提取相应的元素，这些元素就是相关“发电机节点的无功功率对电压幅值的导数”。
[0074] 其定义如下：
[0075]
[0076] 节点i所接发电机的比同步功率因数为矩阵H中Hii。
[0077] 抽取矩阵H的对角线元素，记录所有发电机节点的比同步功率因数。
[0078] S206判断计算次数q<＝n_monte时，则返回步骤S202重新产生发电机和负荷的功率，进行下一次牛顿拉夫逊潮流计算，且q＝q+1；若q>n_monte，则计算结束；
[0079] S207当计算结束后，将得到计算次数为n_monte的计算结果，建立发电机判据列表，其行数为n_monte，列数为发电机个数。
[0080] 例如下表：
[0081]计算次数 一号发电机D1 二号发电机D2 三号发电机D3 四号发电机D4 五号发电机D5
1
2
…
n_monte
[0082] S300计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率。
[0083] S400通过步骤S200和S300，计算系统中每个发电机失去静态稳定性的概率。
[0084] 计算发电机节点无功功率对电压幅值的导数小于0的概率，作为发电机失去静态稳定性的概率，包括步骤：
[0085] 选取发电机判据列表中发电机列元素，统计发电机列元素小于0的次数为m1，则发电机失去静态稳定性的概率为：
[0086]
[0087] 按照上述方法计算所有发电机的失去静态稳定性的概率。
[0088] S500根据对每个发电机失去静态稳定性的概率的分析，得到系统失去静态稳定性的概率。
[0089] 选择所有发电机的失去静态稳定性的概率中最大的概率，作为系统失去静态稳定性的概率。
[0090] 下边举例说明本方法过程：
[0091] 在标准IEEE-14节点系统中进行测试，其中包含14个节点，20条支路，如图2负荷服从正态分布，均值取系统原始负荷峰值，标准差取均值的10％。IEEE-14节点系统中共有5台发电机，分别接在1号、2号、3号、6号和8号节点处,将其编号为1号、2号、3号、4号和5号发电机.发电机正常工作时，输出有功功率和无功功率为系统原始额定功率，发电机的强迫停运概率分别为0.1、0.2、0.05、0.15和0.25.蒙特卡洛法计算次数为10000，即n_monte＝10000，设置1号节点为平衡节点，其电压相角初值为1.06p.u,相角初值为0.由于1号发电机位于平衡节点处，不用计算其比同步功率。根据本专利方法，计算IEEE-14节点系统在如上条件下各发电机失去静态稳定性的概率和系统失去静态稳定性的概率如下表所示：
[0092]
[0093]
[0094] 以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
