CN114139427A - 基于双网格差分法的含风电电力系统频率特征分析方法 - Google Patents

基于双网格差分法的含风电电力系统频率特征分析方法 Download PDF

Info

Publication number
CN114139427A
CN114139427A CN202111533048.5A CN202111533048A CN114139427A CN 114139427 A CN114139427 A CN 114139427A CN 202111533048 A CN202111533048 A CN 202111533048A CN 114139427 A CN114139427 A CN 114139427A
Authority
CN
China
Prior art keywords
wind power
grid
power system
frequency
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111533048.5A
Other languages
English (en)
Inventor
张东寅
许汉平
杨东俊
魏聪
李亚馨
王怡聪
胡婷
王法靖
乔立
童雅芳
颜炯
桑子夏
徐秋实
阮博
周思璇
孙建军
盛举
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Economic and Technological Research Institute of State Grid Hubei Electric Power Co Ltd
Original Assignee
Wuhan University WHU
Economic and Technological Research Institute of State Grid Hubei Electric Power Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU, Economic and Technological Research Institute of State Grid Hubei Electric Power Co Ltd filed Critical Wuhan University WHU
Priority to CN202111533048.5A priority Critical patent/CN114139427A/zh
Publication of CN114139427A publication Critical patent/CN114139427A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/067—Enterprise or organisation modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Energy or water supply
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/04—Power grid distribution networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
    • Y02E40/70—Smart grids as climate change mitigation technology in the energy generation sector
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/50—Systems or methods supporting the power network operation or management, involving a certain degree of interaction with the load-side end user applications

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Strategic Management (AREA)
  • Operations Research (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Geometry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Game Theory and Decision Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Water Supply & Treatment (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Public Health (AREA)
  • Quality & Reliability (AREA)
  • Probability & Statistics with Applications (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • Computer Hardware Design (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)

Abstract

一种基于双网格差分法的含风电电力系统频率特征分析方法，包括以下步骤：S1、建立含风电电力系统功率波动模型；S2、建立电力系统随机动力学模型；S3、建立含风电电力系统一次调频模型；S4、分析含风电电力系统一次调频模型特点，得到含风电电力系统随机动力学模型的FPK方程；S5、求解含风电电力系统频率概率分布密度函数。本设计通过建立含风电电力系统的一次调频模型，进而得到关于系统频率概率分布密度函数的FPK方程，并采用双网格差分迭代数值算法，结合有限元法和有限差分法的优点，准确快速地获得系统频率的概率分布函数，为电网调频能力分析、调频参数优化、机组经济运行等提供参考。

Description

基于双网格差分法的含风电电力系统频率特征分析方法
技术领域
本发明涉及电力系统技术领域，尤其涉及一种基于双网格差分法的含风电电力系统频率特征分析方法，主要适用于提高分析准确度与分析效率。
背景技术
电力系统频率是衡量电能质量的重要标准之一，也是描述电网中发电和用电功率平衡的重要参数。新型电力系统中的新能源出力占比不断提高，常规传统机组被大量 替代，电力系统的惯性支撑能力和频率调节能力发生变化；另外，新能源出力的随机 性增加了电力系统功率波动的复杂性，系统频率概率分布特性的分析难度增大。
一般情况下，由于负荷波动呈现正态分布，电力系统频率的概率分布也是呈现正态分布的。电网实际运行过程中，电力系统的调频环节存在大量非线性因素，如调速 器死区环节、气阀门特性、机组出力限幅环节等。上述一系列因素的影响下，系统频 率的概率分布已经不是简单的正态分布。对电力系统频率概率分布问题的常用研究方 法有两种：一是采用随机信号分析方法进行研究，但是通过分段分析获得的系统频率 概率分布曲线不连续，与电网实际情况有一定偏差；二是通过时域仿真的方法得到一 系列频率数据并对其统计分析，长时间时域仿真具有普适性，但运算量巨大且机理不 明确。
发明内容
本发明的目的是克服现有技术中存在的分析准确度低、分析效率低的缺陷与问题，提供一种分析准确度高、分析效率高的基于双网格差分法的含风电电力系统频率 特征分析方法。
为实现以上目的，本发明的技术解决方案是：一种基于双网格差分法的含风电电力系统频率特征分析方法，包括以下步骤：
S1、建立含风电电力系统功率波动模型；
S2、建立电力系统随机动力学模型；
S3、建立含风电电力系统一次调频模型；
S4、分析含风电电力系统一次调频模型特点，得到含风电电力系统随机动力学模型的FPK方程；
S5、求解含风电电力系统频率概率分布密度函数。
步骤S1具体包括以下步骤：
S11、计算修正后的风电功率波动标准差
Figure BDA0003411542280000021
Figure BDA0003411542280000022
上式中，a为风电功率波动时间尺度系数，表示风电出力变化量随时间尺度变化而变化的规律；zα1为标准正态分布α1分位点，表示风电功率波动分布特性；η为风 电出力占比；b为风电功率波动平滑效应指数，表示风电的空间平滑效应；
S12、计算含风电电力系统总功率波动X：
X＝XW+XL～N(0，σ*2) (2)
Figure BDA0003411542280000023
上式中，负荷功率波动XL、风电功率波动XW相互独立且均服从正态分布，其中
Figure BDA0003411542280000024
σ*为风电接入后的波动强度，σ1为负荷功 率波动标准差。
步骤S2中，电力系统随机动力学模型为：
Figure BDA0003411542280000025
上式中，p为功率波动偏差；γ为负荷波动的均值恢复系数，表示使负荷波动的 随机过程在长时间下接近均值μ的能力；σ为波动强度，ξ(t)为高斯白噪声过程；
设μ为0，则式(4)为：
Figure BDA0003411542280000026
步骤S3中，电力系统调频模型的伊藤随机微分方程组为：
Figure BDA0003411542280000027
上式中，ω为系统频率偏差，δ为调差系数，H为系统惯性时间常数；
风电接入后的机组调差系数δ*为：
Figure BDA0003411542280000031
风电接入后的系统惯性时间常数H*为：
H*＝(1-η)H(8)
含风电电力系统一次调频模型用伊藤随机微分方程组表示为：
Figure BDA0003411542280000032
上式中，ε(ω)为死区函数，死区分为两种，分别为普通型死区ε1(ω)和增强型死区ε2(ω)，表达式为：
Figure BDA0003411542280000033
Figure BDA0003411542280000034
步骤S4具体包括以下步骤：
S41、将含风电电力系统一次调频模型写作向量形式；
对于式(9)，令x＝[x1，x2]＝[ω，p]，以向量形式写出为：
Figure BDA0003411542280000035
其中：
Figure BDA0003411542280000036
Figure BDA0003411542280000037
上式中，f(x)为漂移项，为电力系统动态过程的确定性部分，是状态变量x的 非线性函数；g(x)ξ为扩散项，g(x)是方程的扩散系数，ξ为白噪声向量；
S42、得到含风电电力系统随机动力学模型的FPK方程：
Figure BDA0003411542280000041
上式中，Gi＝[G1，...，Gn]为：
Figure BDA0003411542280000042
上式中，a为漂移系数，b为扩散系数；
Figure BDA0003411542280000043
Figure BDA0003411542280000044
式(9)关于概率密度函数Γ(x，t)的FPK方程为：
Figure BDA0003411542280000045
令：
Figure BDA0003411542280000046
从而得到：
Figure BDA0003411542280000047
步骤S5中，通过以下步骤求解不考虑死区情况下的含风电电力系统频率概率分布密度函数：
对于式(21)，不考虑死区函数，解得Γ(ω，p)的表达式为：
Figure BDA0003411542280000048
其中：
Figure BDA0003411542280000049
Figure BDA0003411542280000051
在p上积分可得关于系统频率偏差ω的概率密度函数F(ω)：
Figure BDA0003411542280000052
通过解析法得到含风电电力系统频率概率分布密度函数；
通过以下步骤求解考虑死区情况下的含风电电力系统频率概率分布密度函数：
式(21)可用如下的偏微分方程表示：
Figure BDA0003411542280000053
上式中，x、y为变量，L是一个微分算子，Ω为边值；
将平面区域进行剖分，分别用网格步长为d和D将Ω剖分为细网格区域Ωd和粗 网格区域ΩD；
在某网格节点上，用五点差分法得到式(24)的差分方程：
Figure BDA0003411542280000054
用矩阵形式表示为：
Figure BDA0003411542280000055
Figure BDA0003411542280000056
上式中，Ld为N×N维的矩阵，LD为M×M维的矩阵，
Figure BDA0003411542280000057
为细网格节点自由度 组成的N维列向量，
Figure BDA0003411542280000058
为粗网格节点自由度组成的M维列向量，
Figure BDA0003411542280000059
为细网格节点上 的广义载荷N维列向量，
Figure BDA00034115422800000510
为粗网格节点上的广义载荷M维列向量；
当粗网格节点与细网格节点重合时，采用直接映射的方式，即粗网格节点值取该节点的细网格节点值：
Figure BDA0003411542280000061
上式中，
Figure BDA0003411542280000062
为粗网格到细网格上网格函数的转移算子；
采用线性插值方法，将粗网格节点值按加权分配到邻近的细网格节点上去：
Figure BDA0003411542280000063
上式中，
Figure BDA0003411542280000064
为细网格到粗网格上网格函数的转移算子；
用双网格有限差分迭代求解方程，按照如下步骤进行运算：
(1)设定初值
Figure BDA0003411542280000065
(2)在细网格区域进行迭代松弛运算κ1次，得到近似值
Figure BDA0003411542280000066
(3)计算细网格损耗量
Figure BDA0003411542280000067
与粗网格到细网格上网格函数的转移算 子
Figure BDA0003411542280000068
(4)在粗网格区域计算精确解
Figure BDA0003411542280000069
并修正细网格节点值；
(5)以
Figure BDA00034115422800000610
为初值在细网格区域进行迭代松弛运算κ2次，得到近似值
Figure BDA00034115422800000611
回到步骤(2)，进行n次运算，从而得到含风电电 力系统频率概率分布密度函数。
与现有技术相比，本发明的有益效果为：
本发明一种基于双网格差分法的含风电电力系统频率特征分析方法中，通过建立含风电电力系统的一次调频模型并得到关于系统频率概率分布密度函数的FPK方程， 采用双网格差分迭代数值算法，通过在细网格区域松弛和在粗网格区域修正，结合有 限元法和有限差分法的优点，可以快速准确地获得系统频率概率分布函数，有效缩短 了数值求解时间，从而为电网调频能力分析、调频参数优化、机组经济运行等提供分 析工具。
附图说明
图1是本发明基于双网格差分法的含风电电力系统频率特征分析方法的流程图。
图2是本发明中普通型和增强型调速器死区函数示意图。
图3是本发明中双网格差分法中网格节点分布情况图。
图4是本发明中双网格有限差分迭代法的算法流程图。
图5是本发明的实施例中不同死区下的系统频率概率分布密度函数曲线。
图6是本发明的实施例中不同方法得到的系统频率概率分布密度函数曲线对比图。
图7是本发明的实施例中修改的3机9节点电力系统示意图。
图8是本发明的实施例中不同死区下的电力系统频率概率分布统计图。
图9是本发明的实施例中中不同风电渗透率下的系统频率概率分布密度函数曲线。
图10是本发明的实施例中不同新能源渗透率下的电力系统频率概率分布统计图。
图11是本发明的实施例中不同调频参数下的系统频率概率分布密度函数曲线。
图12是本发明的实施例中不同调频参数下的电力系统频率概率分布统计图。
具体实施方式
以下结合附图说明和具体实施方式对本发明作进一步详细的说明。
参见图1至图4，一种基于双网格差分法的含风电电力系统频率特征分析方法， 包括以下步骤：
S1、建立含风电电力系统功率波动模型；
S2、建立电力系统随机动力学模型；
S3、建立含风电电力系统一次调频模型；
S4、分析含风电电力系统一次调频模型特点，得到含风电电力系统随机动力学模型的FPK方程；
S5、求解含风电电力系统频率概率分布密度函数。
步骤S1具体包括以下步骤：
S11、计算修正后的风电功率波动标准差
Figure BDA0003411542280000071
Figure BDA0003411542280000072
上式中，a为风电功率波动时间尺度系数，表示风电出力变化量随时间尺度变化而变化的规律；zα1为标准正态分布α1分位点，表示风电功率波动分布特性；η为风 电出力占比；b为风电功率波动平滑效应指数，表示风电的空间平滑效应；
S12、计算含风电电力系统总功率波动X：
X＝XW+XL～N(0，σ*2) (2)
Figure BDA0003411542280000081
上式中，负荷功率波动XL、风电功率波动XW相互独立且均服从正态分布，其中
Figure BDA0003411542280000082
σ*为风电接入后的波动强度，σ1为负荷功 率波动标准差。
步骤S2中，电力系统随机动力学模型为：
Figure BDA0003411542280000083
上式中，p为功率波动偏差；γ为负荷波动的均值恢复系数，表示使负荷波动的 随机过程在长时间下接近均值μ的能力；σ为波动强度，ξ(t)为高斯白噪声过程；
设μ为0，则式(4)为：
Figure BDA0003411542280000084
步骤S3中，电力系统调频模型的伊藤随机微分方程组为：
Figure BDA0003411542280000085
上式中，ω为系统频率偏差，δ为调差系数，H为系统惯性时间常数；
风电接入后的机组调差系数δ*为：
Figure BDA0003411542280000086
风电接入后的系统惯性时间常数H*为：
H*＝(1-η)H (8)
含风电电力系统一次调频模型用伊藤随机微分方程组表示为：
Figure BDA0003411542280000087
上式中，ε(ω)为死区函数，死区分为两种，分别为普通型死区ε1(ω)和增强型死区ε2(ω)，表达式为：
Figure BDA0003411542280000091
Figure BDA0003411542280000092
步骤S4具体包括以下步骤：
S41、将含风电电力系统一次调频模型写作向量形式；
对于式(9)，令x＝[x1，x2]＝[ω，p]，以向量形式写出为：
Figure BDA0003411542280000093
其中：
Figure BDA0003411542280000094
Figure BDA0003411542280000095
上式中，f(x)为漂移项，为电力系统动态过程的确定性部分，是状态变量x的 非线性函数；g(x)ξ为扩散项，g(x)是方程的扩散系数，ξ为白噪声向量；
S42、得到含风电电力系统随机动力学模型的FPK方程：
Figure BDA0003411542280000096
上式中，Gi＝[G1，...，Gn]为：
Figure BDA0003411542280000097
上式中，a为漂移系数，b为扩散系数；
Figure BDA0003411542280000098
Figure BDA0003411542280000101
式(9)关于概率密度函数Γ(x，t)的FPK方程为：
Figure BDA0003411542280000102
令：
Figure BDA0003411542280000103
从而得到：
Figure BDA0003411542280000104
步骤S5中，通过以下步骤求解不考虑死区情况下的含风电电力系统频率概率分布密度函数：
对于式(21)，不考虑死区函数，解得Γ(ω，p)的表达式为：
Figure BDA0003411542280000105
其中：
Figure BDA0003411542280000106
Figure BDA0003411542280000107
在p上积分可得关于系统频率偏差ω的概率密度函数F(ω)：
Figure BDA0003411542280000108
通过解析法得到含风电电力系统频率概率分布密度函数；
通过以下步骤求解考虑死区情况下的含风电电力系统频率概率分布密度函数：
式(21)可用如下的偏微分方程表示：
Figure BDA0003411542280000111
上式中，x、y为变量，L是一个微分算子，Ω为边值；
将平面区域进行剖分，分别用网格步长为d和D将Ω剖分为细网格区域Ωd和粗 网格区域ΩD；
在某网格节点上，用五点差分法得到式(24)的差分方程：
Figure BDA0003411542280000112
用矩阵形式表示为：
Figure BDA0003411542280000113
Figure BDA0003411542280000114
上式中，Ld为N×N维的矩阵，LD为M×M维的矩阵，
Figure BDA00034115422800001113
为细网格节点自由度 组成的N维列向量，
Figure BDA0003411542280000115
为粗网格节点自由度组成的M维列向量，
Figure BDA0003411542280000116
为细网格节点上 的广义载荷N维列向量，
Figure BDA0003411542280000117
为粗网格节点上的广义载荷M维列向量；
当粗网格节点与细网格节点重合时，采用直接映射的方式，即粗网格节点值取该节点的细网格节点值：
Figure BDA0003411542280000118
上式中，
Figure BDA0003411542280000119
为粗网格到细网格上网格函数的转移算子；
采用线性插值方法，将粗网格节点值按加权分配到邻近的细网格节点上去：
Figure BDA00034115422800001110
上式中，
Figure BDA00034115422800001111
为细网格到粗网格上网格函数的转移算子；
用双网格有限差分迭代求解方程，按照如下步骤进行运算：
(1)设定初值
Figure BDA00034115422800001112
(2)在细网格区域进行迭代松弛运算κ1次，得到近似值
Figure BDA0003411542280000121
(3)计算细网格损耗量
Figure BDA0003411542280000122
与粗网格到细网格上网格函数的转移算 子
Figure BDA0003411542280000123
(4)在粗网格区域计算精确解
Figure BDA0003411542280000124
并修正细网格节点值；
(5)以
Figure BDA0003411542280000125
为初值在细网格区域进行迭代松弛运算κ2次，得到近似值
Figure BDA0003411542280000126
回到步骤(2)，进行n次运算，从而得到含风电电 力系统频率概率分布密度函数。
本发明的原理说明如下：
目前对电网频率的研究主要侧重于系统调频及评价指标，关于电网频率分布的研究较少，若能快速准确地获得频率概率分布函数曲线，挖掘频率概率分布特性及其演 化包含的电网频率控制信息，了解系统中影响频率调节能力与频率概率分布的关键因 素，可为电网调频控制的能力分析、参数优化、经济运行等提供参考。
FPK方程是描述粒子在势能场中受到随机力后，随时间演化的位置或速度的分布函数，可以用来计算随机微分方程中某变量分布函数的解；将FPK方程用在电力系统 这一随机动力学系统中，可以研究系统频率这一变量的概率密度函数。
实施例：
参见图1，一种基于双网格差分法的含风电电力系统频率特征分析方法，包括以下步骤：
S1、建立含风电电力系统功率波动模型；具体包括以下步骤：
S11、风电功率波动标准差σ2与整个系统中风电机组出力占比密切相关，计算修正后的风电功率波动标准差
Figure BDA0003411542280000127
Figure BDA0003411542280000128
上式中，a为风电功率波动时间尺度系数，表示风电出力变化量随时间尺度变化而变化的规律；zα1为标准正态分布α1分位点，表示风电功率波动分布特性；η为风 电出力占比；b为风电功率波动平滑效应指数，表示风电的空间平滑效应；上述系数 都可通过历史数据统计而得；
S12、当考虑风电出力时，电力系统净功率波动的概率函数为负荷功率波动和风电出力波动的复合概率函数，工程实践的数据统计证实了负荷的功率波动规律符合正 态分布；计算含风电电力系统总功率波动X：
X＝XW+XL～N(0，σ*2) (2)
Figure BDA0003411542280000131
上式中，负荷功率波动XL、风电功率波动XW相互独立且均服从正态分布，其中
Figure BDA0003411542280000132
σ*为风电接入后的波动强度，σ1为负荷功 率波动标准差(由历史数据统计而得)；
S2、考虑系统功率随机波动的电力系统模型往往用随机微分方程组(stochasticdifferential equation，SDE)来表示，能体现负荷功率波动的随机过程是奥恩斯坦 -乌伦贝克(Ornstein-Uhlenbeck，OU)过程，与标准Wiener过程相比，OU过程更适 合于对动力学系统中的随机扰动建模；建立电力系统随机动力学模型为：
Figure BDA0003411542280000133
上式中，p为功率波动偏差；γ为负荷波动的均值恢复系数，表示使负荷波动的 随机过程在长时间下接近均值μ的能力；σ为波动强度，ξ(t)为高斯白噪声过程；
设μ为0，则式(4)为：
Figure BDA0003411542280000134
S3、不考虑调速器时滞特性、死区特性、原动机时间常数与二次调频的情况下， 建立含风电电力系统一次调频模型；
电力系统调频模型的伊藤随机微分方程组为：
Figure BDA0003411542280000135
上式中，ω为系统频率偏差，δ为调差系数，H为系统惯性时间常数；
由于风电机组本身只能提供有限的虚拟惯性，且风电机组参与一次调频控制的技术尚不成熟，电力系统中的调频任务还是由常规同步发电机组承担，随着系统中的风 电机组出力占比提高，电力系统的一次调频特性也随之发生改变；
电力系统发电机组总的单位调节功率降低，可得风电接入后的机组调差系数δ*为：
Figure BDA0003411542280000141
风电接入后的系统惯性时间常数H*为：
H*＝(1-η)H (8)
考虑死区环节的含风电电力系统一次调频模型用伊藤随机微分方程组表示为：
Figure BDA0003411542280000142
上式中，ε(ω)为死区函数，死区分为两种，分别为普通型死区ε1(ω)和增强型死区ε2(ω)，如图2所示，表达式为：
Figure BDA0003411542280000143
Figure BDA0003411542280000144
S4、分析含风电电力系统一次调频模型特点，得到含风电电力系统随机动力学模型的FPK方程；具体包括以下步骤：
S41、将含风电电力系统一次调频模型写作向量形式；
对于式(9)，令x＝[x1，x2]＝[ω，p]，以向量形式写出为：
Figure BDA0003411542280000145
其中：
Figure BDA0003411542280000146
Figure BDA0003411542280000147
上式中，f(x)为漂移项，为电力系统动态过程的确定性部分，是状态变量x的 非线性函数；g(x)ξ为扩散项，g(x)是方程的扩散系数，ξ为白噪声向量；
S42、得到含风电电力系统随机动力学模型的FPK方程：
Figure BDA0003411542280000151
上式中，Γ表示概率密度函数Γ(x，y)，x表示ω，t为时间；Gi＝[G1，...，Gn]为：
Figure BDA0003411542280000152
上式中，a为漂移系数，b为扩散系数；n为维数，本设计为2；
Figure BDA0003411542280000153
Figure BDA0003411542280000154
式(9)关于概率密度函数Γ(x，t)的FPK方程为：
Figure BDA0003411542280000155
为了得到长时间内的概率密度，令：
Figure BDA0003411542280000156
从而得到：
Figure BDA0003411542280000157
S5、求解含风电电力系统频率概率分布密度函数；
通过以下步骤求解不考虑死区情况下的含风电电力系统频率概率分布密度函数：
对于式(21)，不考虑死区函数，解得Γ(ω，p)的表达式为：
Figure BDA0003411542280000158
其中：
Figure BDA0003411542280000161
Figure BDA0003411542280000162
在p上积分可得关于系统频率偏差ω的概率密度函数F(ω)：
Figure BDA0003411542280000163
通过解析法得到含风电电力系统频率概率分布密度函数；
通过以下步骤求解考虑死区情况下的含风电电力系统频率概率分布密度函数：
式(21)的FPK方程是一个椭圆型变系数偏微分方程，可用如下的偏微分方程表示：
Figure BDA0003411542280000164
上式中，x、y为变量，L是一个微分算子，Ω为边值；
将平面区域进行剖分，分别用网格步长为d和D将Ω剖分为细网格区域Ωd和粗 网格区域ΩD，且D＝2d；细网格节点与粗网格节点的分布情况如图3所示；
在某网格节点上，用五点差分法得到式(24)的差分方程：
Figure BDA0003411542280000165
上式中，Δd表示离散化的微分算子，即Lu的离散化为Δdui，j；fij表示f(x，y)； gi，j表示g(x，y)；i，j表示横纵离散点；
用矩阵形式表示为：
Figure BDA0003411542280000166
Figure BDA0003411542280000167
上式中，Ld为N×N维的矩阵，LD为M×M维的矩阵，
Figure BDA0003411542280000171
为细网格节点自由度 组成的N维列向量，
Figure BDA0003411542280000172
为粗网格节点自由度组成的M维列向量，
Figure BDA0003411542280000173
为细网格节点上 的广义载荷N维列向量，
Figure BDA0003411542280000174
为粗网格节点上的广义载荷M维列向量；
当粗网格节点与细网格节点重合时，采用直接映射的方式，即粗网格节点值取该节点的细网格节点值：
Figure BDA0003411542280000175
上式中，
Figure BDA0003411542280000176
为粗网格到细网格上网格函数的转移算子，ud为向量
Figure BDA0003411542280000177
里的具体 数值；
采用线性插值方法，将粗网格节点值按加权分配到邻近的细网格节点上去：
Figure BDA0003411542280000178
上式中，
Figure BDA0003411542280000179
为细网格到粗网格上网格函数的转移算子；
参见图4，用双网格有限差分迭代求解方程，按照如下步骤进行运算：
(1)设定初值
Figure BDA00034115422800001710
(2)在细网格区域进行迭代松弛运算κ1次，得到近似值
Figure BDA00034115422800001711
(3)计算细网格损耗量
Figure BDA00034115422800001712
与粗网格到细网格上网格函数的转移算 子
Figure BDA00034115422800001713
(4)在粗网格区域计算精确解
Figure BDA00034115422800001714
并修正细网格节点值；
(5)以
Figure BDA00034115422800001715
为初值在细网格区域进行迭代松弛运算κ2次，得到近似值
Figure BDA00034115422800001716
回到步骤(2)，进行n次运算，从而得到含风电电 力系统频率概率分布密度函数。
采用的初始条件和无穷边界条件为：
Figure BDA00034115422800001717
Figure BDA0003411542280000181
运算过程中，粗网格精度为40×40，细网格精度为80×80，求得方程的数值解， 从而得到频率偏差的概率密度曲线。
在单机系统随机动力学模型中，求得不考虑死区、考虑普通型死区、考虑增强型死区三种情况下的系统频率概率分布密度(Probability distribution density，PDF) 函数曲线，如图5所示。与蒙特卡洛模拟法得到的频率概率分布统计进行对比，如图 6所示。图6(a)中，蒙特卡洛模拟法的仿真时间为1×105s，仿真步长为0.05s，计 算机实际的运行时长为12.6s；图6(b)中，仿真时间为5×105s，仿真步长为0.01s， 实际运行时长为43.2s；图6(c)中，仿真时间为1.5×106s，仿真步长为0.01s，实 际运行时长为383.3s；本设计的数值解法运算时长为40.1s。对比图6(a)、图6(b) 和图6(c)，采取蒙特卡洛法只能得到离散非光滑的电网频率概率分布密度函数曲线， 样本数量越小，系统频率概率分布的统计图体现的频率概率分布双驼峰效果也越差， 如需得到同样光滑的结果，需要将蒙特卡洛法里面的模拟次数提高数量级，从而其计 算时间也要相应数量级提高。
为了进一步证明本设计方法的快速性与精确性，表1给出了在同样的求解精度下，蒙特卡洛模拟法、有限元法和本设计的双网格有限差分迭代法的耗时对比。以有限元 法在10×10的网格中1×10-5的剖分精度计算得到的概率密度函数作为标准解，将概 率密度函数数值解与标准解的差的绝对值作为精度标准。其中，所有模拟和数值计算 都是在64位Windows 10的计算机上进行，该计算机具有2.9GHz主频的AMD R7-4800H 和16GB内存，软件平台为Mathematica。
表1不同方法耗时对比
Figure BDA0003411542280000182
对比可见，本设计的方法能更加快速准确地得到系统频率概率分布密度函数曲线。表2给出了不同精度下本设计方法相比有限元的减少时间的百分比。不同的求解 精度下，本设计都可以减少求解时间，并且随着求解精度的不断提高，本设计减少的 时间百分比也在增加。
表2不同精度下的减少时间百分比
Figure BDA0003411542280000183
Figure BDA0003411542280000191
为了验证本设计算法的合理性，在DIgSILENT/PowerFactory环境搭建如图7所 示的含风、火电的3机9节点电力系统。该系统包含3台火电机组G1、G2与G3，1 台风电机组，其中，调频机组G1-G3的总装机容量分别为1000MW、600MW、800MW；风 电机组不参于电力系统一次调频，装机容量为450MW。在负荷L1-L3和风电机组出力 端注入高斯白噪声，模拟负荷端和风电机组出力的功率波动。仿真时间为80000s，选 取节点B5的频率作为采用数据，采样间隔为1s。在仿真模型中，对G1-G3分别不设 死区、设置普通型死区、设置增强型死区，电力频率概率分布统计如图8所示。分别 采用不同的风电出力占比，系统频率概率分布密度函数曲线如图9所示。当风电出力 占比从10％逐渐增加到30％，系统等效惯性时间常数不断降低，等效调差系统不断增 大，总功率波动的标准差也随之增大。整个系统的一次调频能力不断下降，概率密度 分布曲线的驼峰顶与双驼峰间的鞍部都不断变低。在仿真模型中，减小系统中火电机 组容量，增加风电机组容量，将风电出力功率等比例增大或减小，以改变风电在系统 中的出力占比。三种不同风电出力占比情况下的电网频率概率分布统计如图10所示。
分别采用不同的系统参数，系统频率概率分布密度函数曲线如图11所示。采用 不同的惯性时间常数H，系统频率概率分布密度函数曲线变化不大，所以惯性时间常 数H的影响较小。图11(b)中，当负荷均值灰度系统从0.06变为0.04时，概率密 度分布曲线的驼峰顶与双驼峰间的鞍部都不断变低。图11(c)中，当调速器死区从 1.5r/min变为2.5r/min时，概率密度分布曲线的驼峰顶不降低，双驼峰间的鞍部宽 度增大，鞍部最低点降低。系统频率处于死区之外的时间从85.3％变为78.4％。随着 死区的增加，双驼峰现象愈加明细。图11(d)中，当调速不等率从0.06变为0.04 时，系统的一次调频能力不断增强，概率密度分布曲线的驼峰顶不断变高，双驼峰间 的鞍部则基本保持不变。图12给出了仿真模型中，不同系统参数对系统频率概率分 布的影响。
上述的对比体现了本设计算法的合理性。本设计通过建立含风电电力系统的一次调频模型，得到关于系统频率概率分布密度函数的FPK方程，采用双网格差分迭代数 值算法结合了有限元法和有限差分法的优点，可以快速准确地获得频率概率分布函 数，从而为电网调频能力分析、调频参数优化、机组经济运行等提供参考目标。

Claims (6)

1.一种基于双网格差分法的含风电电力系统频率特征分析方法，其特征在于，包括以下步骤：
S1、建立含风电电力系统功率波动模型；
S2、建立电力系统随机动力学模型；
S3、建立含风电电力系统一次调频模型；
S4、分析含风电电力系统一次调频模型特点，得到含风电电力系统随机动力学模型的FPK方程；
S5、求解含风电电力系统频率概率分布密度函数。
2.根据权利要求1所述的一种基于双网格差分法的含风电电力系统频率特征分析方法，其特征在于：步骤S1具体包括以下步骤：
S11、计算修正后的风电功率波动标准差
Figure FDA0003411542270000011
Figure FDA0003411542270000012
上式中，a为风电功率波动时间尺度系数，表示风电出力变化量随时间尺度变化而变化的规律；zα1为标准正态分布α1分位点，表示风电功率波动分布特性；η为风电出力占比；b为风电功率波动平滑效应指数，表示风电的空间平滑效应；
S12、计算含风电电力系统总功率波动X：
X＝XW+XL～N(0，σ*2) (2)
Figure FDA0003411542270000013
上式中，负荷功率波动XL、风电功率波动XW相互独立且均服从正态分布，其中
Figure FDA0003411542270000014
σ*为风电接入后的波动强度，σ1为负荷功率波动标准差。
3.根据权利要求2所述的一种基于双网格差分法的含风电电力系统频率特征分析方法，其特征在于：步骤S2中，电力系统随机动力学模型为：
Figure FDA0003411542270000015
上式中，p为功率波动偏差；γ为负荷波动的均值恢复系数，表示使负荷波动的随机过程在长时间下接近均值μ的能力；σ为波动强度，ξ(t)为高斯白噪声过程；
设μ为0，则式(4)为：
Figure FDA0003411542270000021
4.根据权利要求3所述的一种基于双网格差分法的含风电电力系统频率特征分析方法，其特征在于：
步骤S3中，电力系统调频模型的伊藤随机微分方程组为：
Figure FDA0003411542270000022
上式中，ω为系统频率偏差，δ为调差系数，H为系统惯性时间常数；
风电接入后的机组调差系数δ*为：
Figure FDA0003411542270000023
风电接入后的系统惯性时间常数H*为：
H*＝(1-η)H (8)
含风电电力系统一次调频模型用伊藤随机微分方程组表示为：
Figure FDA0003411542270000024
上式中，ε(ω)为死区函数，死区分为两种，分别为普通型死区ε1(ω)和增强型死区ε2(ω)，表达式为：
Figure FDA0003411542270000025
Figure FDA0003411542270000026
5.根据权利要求4所述的一种基于双网格差分法的含风电电力系统频率特征分析方法，其特征在于：步骤S4具体包括以下步骤：
S41、将含风电电力系统一次调频模型写作向量形式；
对于式(9)，令x＝[x1，x2]＝[ω，p]，以向量形式写出为：
Figure FDA0003411542270000031
其中：
Figure FDA0003411542270000032
Figure FDA0003411542270000033
上式中，f(x)为漂移项，为电力系统动态过程的确定性部分，是状态变量x的非线性函数；g(x)ξ为扩散项，g(x)是方程的扩散系数，ξ为白噪声向量；
S42、得到含风电电力系统随机动力学模型的FPK方程：
Figure FDA0003411542270000034
上式中，Gi＝[G1，...，Gn]为：
Figure FDA0003411542270000035
上式中，a为漂移系数，b为扩散系数；
Figure FDA0003411542270000036
Figure FDA0003411542270000037
式(9)关于概率密度函数Γ(x，t)的FPK方程为：
Figure FDA0003411542270000038
令：
Figure FDA0003411542270000039
从而得到：
Figure FDA0003411542270000041
6.根据权利要求5所述的一种基于双网格差分法的含风电电力系统频率特征分析方法，其特征在于：
步骤S5中，通过以下步骤求解不考虑死区情况下的含风电电力系统频率概率分布密度函数：
对于式(21)，不考虑死区函数，解得Γ(ω，p)的表达式为：
Figure FDA0003411542270000042
其中：
Figure FDA0003411542270000043
Figure FDA0003411542270000044
在p上积分可得关于系统频率偏差ω的概率密度函数F(ω)：
Figure FDA0003411542270000045
通过解析法得到含风电电力系统频率概率分布密度函数；
通过以下步骤求解考虑死区情况下的含风电电力系统频率概率分布密度函数：
式(21)可用如下的偏微分方程表示：
Figure FDA0003411542270000046
上式中，x、y为变量，L是一个微分算子，Ω为边值；
将平面区域进行剖分，分别用网格步长为d和D将Ω剖分为细网格区域Ωd和粗网格区域ΩD；
在某网格节点上，用五点差分法得到式(24)的差分方程：
Figure FDA0003411542270000051
用矩阵形式表示为：
Figure FDA0003411542270000052
Figure FDA0003411542270000053
上式中，Ld为N×N维的矩阵，LD为M×M维的矩阵，
Figure FDA0003411542270000054
为细网格节点自由度组成的N维列向量，
Figure FDA0003411542270000055
为粗网格节点自由度组成的M维列向量，
Figure FDA0003411542270000056
为细网格节点上的广义载荷N维列向量，
Figure FDA0003411542270000057
为粗网格节点上的广义载荷M维列向量；
当粗网格节点与细网格节点重合时，采用直接映射的方式，即粗网格节点值取该节点的细网格节点值：
Figure FDA0003411542270000058
上式中，
Figure FDA0003411542270000059
为粗网格到细网格上网格函数的转移算子；
采用线性插值方法，将粗网格节点值按加权分配到邻近的细网格节点上去：
Figure FDA00034115422700000510
上式中，
Figure FDA00034115422700000511
为细网格到粗网格上网格函数的转移算子；
用双网格有限差分迭代求解方程，按照如下步骤进行运算：
(1)设定初值
Figure FDA00034115422700000512
(2)在细网格区域进行迭代松弛运算κ1次，得到近似值
Figure FDA00034115422700000513
(3)计算细网格损耗量
Figure FDA00034115422700000514
与粗网格到细网格上网格函数的转移算子
Figure FDA0003411542270000061
(4)在粗网格区域计算精确解
Figure FDA0003411542270000062
并修正细网格节点值；
(5)以
Figure FDA0003411542270000063
为初值在细网格区域进行迭代松弛运算κ2次，得到近似值
Figure FDA0003411542270000064
回到步骤(2)，进行n次运算，从而得到含风电电力系统频率概率分布密度函数。
CN202111533048.5A 2021-12-15 2021-12-15 基于双网格差分法的含风电电力系统频率特征分析方法 Pending CN114139427A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111533048.5A CN114139427A (zh) 2021-12-15 2021-12-15 基于双网格差分法的含风电电力系统频率特征分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111533048.5A CN114139427A (zh) 2021-12-15 2021-12-15 基于双网格差分法的含风电电力系统频率特征分析方法

Publications (1)

Publication Number Publication Date
CN114139427A true CN114139427A (zh) 2022-03-04

Family

ID=80382501

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111533048.5A Pending CN114139427A (zh) 2021-12-15 2021-12-15 基于双网格差分法的含风电电力系统频率特征分析方法

Country Status (1)

Country Link
CN (1) CN114139427A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117554862A (zh) * 2024-01-11 2024-02-13 山东康吉诺技术有限公司 一种变压器智能检测预警方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117554862A (zh) * 2024-01-11 2024-02-13 山东康吉诺技术有限公司 一种变压器智能检测预警方法及系统
CN117554862B (zh) * 2024-01-11 2024-03-29 山东康吉诺技术有限公司 一种变压器智能检测预警方法及系统

Similar Documents

Publication Publication Date Title
Bemporad et al. Ultra-fast stabilizing model predictive control via canonical piecewise affine approximations
Cohen Optimization by decomposition and coordination: A unified approach
Glizer Asymptotic analysis and solution of a finite-horizon H∞ control problem for singularly-perturbed linear systems with small state delay
Hajivassiliou et al. Simulation of multivariate normal rectangle probabilities and their derivatives theoretical and computational results
Orcan et al. Elastic–plastic stresses in linearly hardening rotating solid disks of variable thickness
Thangavel et al. Dynamical analysis of T–S fuzzy financial systems: A sampled-data control approach
CN114139427A (zh) 基于双网格差分法的含风电电力系统频率特征分析方法
CN110649588B (zh) 一种柔性直流输电控制系统攻击量化评估方法
CN111158241B (zh) 具有不确定时滞的线性奇异系统的时滞相关h∞控制方法
Davis A differential equation approach to linear combinations of independent chi-squares
Tuan et al. Monotonic relaxations for robust control: New characterizations
CN111399376A (zh) 一种t-s模糊系统的二维重复控制器设计优化方法
CN113591417B (zh) 一种应用于高精度间断迦辽金流体仿真的粘性项处理方法
CN115276026A (zh) 一种时滞电力系统稳定性分析方法及系统
Xie et al. Nonuniform piecewise membership function approximation methods based robust tracking control design of T‐S fuzzy systems
Li et al. Development of an efficient wetting and drying treatment for shallow‐water modeling using the quadrature‐free Runge‐Kutta discontinuous Galerkin method
Gonçalves et al. H∞ state feedback control of discrete-time Markov jump linear systems through linear matrix inequalities
Aguila-Camacho et al. Error-Based Switched Fractional Order Model Reference Adaptive Control for MIMO Linear Time Invariant Systems
Wang et al. Optimization of transient performance for continuous-time switched linear autonomous systems
Santos et al. On the structure of the equilibrium price set of overlapping-generations economies
CN106649203A (zh) 一种提高大数据处理质量的方法
Baár et al. Robust minimum gain lemma
Tu et al. Complexity and control of a cournot duopoly game in exploitation of a renewable resource with bounded rationality players
Liu et al. A new approach to stabilization of uncertain nonlinear systems
CN116205074B (zh) 一种反应堆堆芯功率非线性控制算法的稳定域获取方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination