CN104298809B - 一种基于矩阵指数电磁暂态仿真的非线性建模求解方法 - Google Patents

Abstract

一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，首先在状态分析框架下建立线性网络和并网非线性元件的电磁暂态仿真模型，然后通过子系统互联关系形成待研究电力系统的整体电磁暂态仿真模型；设定仿真步长、收敛精度等相关仿真参数后，启动仿真程序；在每一个仿真步长内，使用不动点迭代方法求解包含有矩阵指数函数的非线性方程，取其结果作为当前时刻的状态变量；由输出方程得到输出向量yn+1，写入输出文件，仿真推进一个步长；依此迭代进行，直到仿真结束。本发明保留了矩阵指数积分方法良好的数值精度和刚性处理能力，并且对电力系统元件的非线性特性具有一般性的建模和仿真能力，扩展了矩阵指数积分方法的在电力系统电磁暂态仿真领域的适用范围。

Description

一种基于矩阵指数电磁暂态仿真的非线性建模求解方法

技术领域

本发明涉及一种仿真的非线性求解方法。特别是涉及一种适用于电力系统电磁暂态建模仿真应用的基于矩阵指数电磁暂态仿真的非线性建模求解方法。

背景技术

电力系统电磁暂态仿真主要反映系统中电场与磁场的相互影响产生的电气量的变化过程，可得到从工频到几十kHz频谱范围内的三相电压电流瞬时值波形。为了准确获取系统动态过程中的高频特性，电磁暂态的仿真步长通常是微秒级。随着直流技术的发展与广泛应用，交直流系统互联给电磁暂态仿真提出了新的挑战。近年来我国电力系统实际运行中出现的一些新的仿真场景，系统规模大，仿真时间长，使用传统电磁仿真算法会消耗较多计算资源和计算时间，需要结合问题特性从算法层面提出针对性的改进。

电力系统电磁暂态仿真本质上可归结为对动力学系统时域响应的求取，它包括系统本身的数学模型和与之相适应的数值算法。

当前，电力系统电磁暂态仿真基本框架可分为两类，包括节点分析法(NodalAnalysis)和状态变量分析法(State-Variable Analysis)。基于节点分析框架的电磁暂态仿真方法可概括为先采用某种数值积分方法(通常为梯形积分法)将系统中动态元件的特性方程差分化，得到等效的计算电导与历史项电流源并联形式的诺顿等效电路，此时联立整个电气系统的元件特性方程形成节点电导矩阵，如式(1)所示，对其求解即可得到系统中各节点电压的瞬时值。

Gu＝i (1)

式(1)所示的节点电导矩阵为线性方程组，可使用各种成熟的线性稀疏矩阵算法库进行求解。节点分析法广泛应用于EMTP、PSCAD/EMTDC等专业的电力系统电磁暂态仿真程序中，工程上也称基于节点分析框架的电磁暂态仿真工具为EMTP类程序。节点分析法的主要优势体现在程序实现难度和仿真计算效率方面，但由于式(1)的节点电导方程本身已将数值积分方法与系统模型融为一体，导致EMTP类程序在求解算法选择方面缺乏灵活性与开放性，同时式(1)已不能给出系统本身的特征信息。

与节点分析法不同，状态变量分析法属于一般性建模方法(general purposemodeling)，不仅适于电路与电力系统仿真，同样也适于其它工程领域的动力学系统的建模与仿真。Matlab/SimPowerSystems软件是状态变量分析框架下暂态仿真程序的典型代表。与节点分析框架相比，状态方程在模型的计算求解方面具有高度的开放性和灵活性，可方便地选择与问题相适应的数值积分方法，同时能够提供关于系统各种特征的丰富信息(如系统的特征值)，进而能够从全局角度了解系统的动态特性，为各种快速、准确、高效的仿真算法的开发与测试工作提供了便利条件。

应用状态变量分析的基础是形成式(2)所示标准形式的状态-输出方程，此时系统中的电源作为输入u。

y＝Cx+Du (2)

在电力系统仿真领域，式(2)可以由改进节点法Modified Nodal Analysis(MNA)通过KCL、KVL等约束关系以及元件伏安特性进行构造而得到MNA模型，再经过一定的正规化处理(regularization)转化而来。MNA模型是形如式(3)的状态-输出方程。

y＝L^Tx (3)

也可以采用一般支路类方法，如Automated State Model Generator(ASMG)方法直接构造得到。

在数值算法方面，传统数值积分方法可分为显式和隐式两类，不同积分方法所具有的数值稳定性和数值精度各不相同。一般来说，隐式方法处理仿真模型中刚性特征的能力较强。电力系统由于动态过程时间尺度差异较大，系统模型表现出一定刚性，这使得主流电磁暂态软件EMTP类程序采用隐式方法以保证数值稳定性。从计算开销方面来看，隐式方法在每一时步内需求解线性方程组，极大限制了其在大规模系统的应用能力。与之相对的，传统显式方法无需迭代，在每一时步内的运算量较小，但其有限的数值稳定域使得仿真步长受到约束，综合来看对刚性系统的仿真性能不佳。对于现代电力系统来说，系统中既存在微秒级的电力电子开关动态过程，又存在同步机组的励磁、调速等秒级的机电动态过程，时间尺度差异极大，系统刚性特征十分显著。充分利用状态方程框架在数值算法选择方面的灵活性，结合电力系统电磁暂态仿真的应用场景与特殊需求，发展合适的数值积分方法，是提高电力系统电磁暂态仿真的计算性能和应用前景的关键。

矩阵指数积分方法(Exponential Integrator)是近年来从应用数学领域发端的一种数值积分方法。它使用矩阵指数算子e^hA精确描述动态系统的线性变化规律，可以准确求解形如

的线性动态系统，并具有计算效率高、刚性处理能力强等特点。矩阵指数积分方法已经在诸如应用物理、化学工程等领域得到一定应用。然而，现代电力系统的一个重要特征在于其元件非线性环节带来的复杂动态特征，这使得矩阵指数积分方法不能直接应用于具有强非线性特征的电力系统模型。为了在电力系统电磁暂态仿真问题中使用矩阵指数积分方法，需要考虑与其匹配的非线性环节的处理方法。

发明内容

本发明所要解决的技术问题是，提供一种结合了矩阵指数对于线性动态系统的精确仿真能力和刚性处理能力，并弥补了其处理非线性元件的能力不足缺陷的基于矩阵指数电磁暂态仿真的非线性建模求解方法。

本发明所采用的技术方案是：一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，通过子系统互联形成包含非线性元件的研究系统整体仿真模型，结合矩阵指数方法的计算精度优势以及不动点迭代的计算速度优势，实现具有强非线性特征的电力系统电磁暂态仿真，具体包括如下步骤：

1)划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，选取并网点端口电压和电流作为网络与并网元件的接口变量，在状态分析框架下，分别建立待研究电力系统中线性网络子系统、以及各非线性并网元件子系统的电磁暂态仿真模型；

2)根据非线性元件与线性网络之间的子系统互联关系，由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，所述整体状态空间仿真模型形式为

其中，x是包含当前时刻电力系统内所有储能元件与记忆元件状态的状态向量，t为时间，Ax和f(x,t)分别表示电力系统动态模型中的线性和非线性部分，它们由各个子系统电磁暂态仿真模型和子系统之间互联关系确定，y是仿真使用者所要求的仿真输出向量，输出函数g(x,t)根据仿真研究关注点不同，由使用者任意指定；

3)设定：仿真时间T，仿真步长Δt，非线性迭代收敛精度∈，设定当前时刻t_n为仿真起始时刻t₀，依照仿真需要，设置仿真初值x₀，并赋值给当前时刻状态向量x_n，计算仿真起始时刻输出向量y₀＝g(x₀,t₀)，并写入输出文件；

4)计算当前时刻t_n时整体仿真模型中的非线性项f(x_n,t_n)，并计算 其中，e^ΔtA为矩阵指数函数表示的单位步长下系统线性状态转移矩阵；

5)记下一个时刻的时间值t_n+1＝t_n+Δt，采用不动点迭代方法求解非线性方程

至设定的非线性迭代收敛精度∈，得到t_n+1时刻的状态向量x_n+1；

6)由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出文件，更新当前时刻为下一时刻t_n＝t_n+1，仿真向前推进一个步长；

7)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

步骤1)所述的划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，具体划分范围是：其中，PI型等效电路建模的输配电线路、不考虑饱和效应的电力变压器、线性恒阻抗负荷、传统串补设备、线路高抗、无功补偿电容器以及调压器设备在线性网络子系统内统一建模；而同步和异步发电机组、电力电子器件、可再生能源电源以及电弧炉非线性负荷设备作为非线性并网元件子系统分别建模。

所述的在线性网络子系统内统一建模，得到的子系统电磁暂态仿真模型具体形式为：

其中，x_net是线性网络子系统的状态向量，包括独立电感电流和电容电压，A_net是线性网络子系统的状态矩阵；u_net,i是来自第i个非线性并网元件子系统的输入变量，取决于具体接口方式，视为受控的电压源或电流源，B_net,i是对应的输入矩阵；y_net,i是线性网络子系统给第i个非线性并网元件子系统的输出变量，取决于具体接口方式，选取端口电压或输出电流，C_net,i是对应的输出矩阵。

所述的对非线性并网元件子系统分别建模，具体是指：

待研究电力系统包含有多个非线性并网元件，对于第i个元件，非线性并网元件子系统电磁暂态仿真模型形式为

其中，x_i是所述非线性并网元件子系统的状态向量，A_i是所述非线性并网元件子系统仿真模型中的线性部分；u_i是所述非线性并网元件子系统来自线性网络子系统的输入变量，B_i是对应的输入矩阵；f_i(x_i,u_i,t)是所述非线性并网元件子系统的非线性特性函数；y_i是所述非线性并网元件子系统给线性网络子系统的输出变量，C_i是对应的输出矩阵，g_i(x_i,u_i,t)是输出变量的非线性特性函数。

步骤2)所述的根据非线性元件与线性网络子系统之间的互联关系，由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，具体实现方式为：

首先由输入文件读入非线性并网元件子系统与线性网络子系统之间的互联关系，所述互联关系指的是，一个子系统的输出变量对应着另一个子系统的输入变量，用数学公式表达为

联列线性网络子系统和非线性并网元件子系统的电磁暂态仿真模型，并将输入和输出变量替换为子系统状态向量的函数，即得到待研究电力系统的整体状态空间仿真模型

其中，Ax和f(x,t)分别是系统仿真模型中的线性和非线性部分，它们与子系统电磁暂态仿真模型参数的关系是：

步骤5)所述的采用不动点迭代方法求解非线性方程，包括如下步骤：

(1)设定最大迭代次数N，外插法设置迭代初值记当前已迭代步数k＝0；

(2)比较已迭代步数k与最大迭代步数N，如果k≥N，跳转至步骤(4)，否则，使用如下公式进行第k+1步迭代：

并更新已迭代步数k＝k+1；

(3)计算当前步与前一步迭代值的偏差量并比较这一偏差dx与非线性迭代收敛精度∈，如果dx>∈，跳转到步骤(2)继续进行迭代运算，否则进行到下一步骤；

(4)将当前的迭代值作为t_n+1时刻的状态量值，即完成当前时步的迭代运算。

本发明的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，在状态分析框架下，基于矩阵指数运算和不动点迭代实现了一般电力系统电磁暂态模型的仿真计算，使得矩阵指数电磁暂态仿真方法适用范围更广，可涵盖工程实际遇到的电力系统电磁暂态仿真问题。本发明保留了矩阵指数积分方法良好的数值精度和刚性处理能力，并且对电力系统元件的非线性特性具有一般性的建模和仿真能力，扩展了矩阵指数积分方法在电力系统电磁暂态仿真领域的适用范围。本发明采用不动点迭代方法进行非线性方程求解，避免了雅各比矩阵计算和线性方程组求解，有较高的计算效率。

附图说明

图1是并网型分布式电源仿真算例结构图；

图2是本发明的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法的流程图；

图3是并网型分布式电源仿真算例中，不同仿真算法得到的分布式电源向网络输送的有功功率；

图4是并网型分布式电源仿真算例中，不同仿真算法得到的分布式电源向网络输送的有功功率在区间[0.200825，0.20087]上的放大图，图中可清楚看到本发明提出的矩阵指数方法与隐式梯形法之间的数值精度比较；

图5是并网型分布式电源仿真算例中，分布式电源向网络输送的有功功率，采用不同数值积分方法时的绝对误差比较；

图6是并网型分布式电源仿真算例中，不同仿真算法得到的LCL滤波器输出电感上流过的A相电流；

图7是并网型分布式电源仿真算例中，不同仿真算法得到的LCL滤波器输出电感上流过的A相电流在区间[0.20474，0.20483]上的放大图，图中可清楚看到本发明提出的矩阵指数方法与隐式梯形法之间的数值精度比较；

图8是并网型分布式电源仿真算例中，LCL滤波器输出电感上流过的A相电流，采用不同数值积分方法时的绝对误差比较。

具体实施方式

下面结合实施例和附图对本发明的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法做出详细说明。

本发明的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，是通过子系统互联形成包含非线性元件的研究系统整体仿真模型，结合矩阵指数方法的计算精度优势以及不动点迭代的计算速度优势，实现具有强非线性特征的电力系统电磁暂态仿真。算法流程图如图2所示，具体包括如下步骤：

1)划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，选取并网点端口电压和电流作为网络与并网元件的接口变量，在状态分析框架下，分别建立待研究电力系统中线性网络子系统、以及各并网非线性元件的子系统电磁暂态仿真模型；

所述的划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，具体划分范围是：其中，PI型等效电路建模的输配电线路、不考虑饱和效应的电力变压器、线性恒阻抗负荷、传统串补设备、线路高抗、无功补偿电容器以及调压器设备在线性网络子系统内统一建模；而同步和异步发电机组、电力电子器件、可再生能源电源以及电弧炉非线性负荷设备作为非线性并网元件子系统则分别建模。其中，

所述的在线性网络子系统内统一建模，得到的子系统电磁暂态仿真模型具体形式为：

其中，x_net是线性网络子系统的状态向量，包括独立电感电流和电容电压，A_net是线性网络子系统的状态矩阵；u_net,i是来自第i个非线性并网元件子系统的输入变量，取决于具体接口方式，视为受控的电压源或电流源，B_net,i是对应的输入矩阵；y_net,i是线性网络子系统给第i个非线性并网元件子系统的输出变量，取决于具体接口方式，选取端口电压或输出电流，C_net,i是对应的输出矩阵。

所述的对非线性并网元件子系统分别建模是：

待研究电力系统包含有多个非线性并网元件，对于第i个元件，其子系统电磁暂态仿真模型形式为

其中，x_i是所述子系统的状态向量，A_i是所述子系统仿真模型中的线性部分；u_i是所述子系统来自线性网络子系统的输入变量，B_i是对应的输入矩阵；f_i(x_i,u_i,t)是所述子系统的非线性特性函数；y_i是所述子系统给线性网络子系统的输出变量，C_i是对应的输出矩阵，g_i(x_i,u_i,t)是输出变量的非线性特性函数。

2)根据非线性元件与线性网络之间的子系统互联关系，由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，所述整体状态空间仿真模型形式为

其中，x是包含当前时刻电力系统内所有储能元件与记忆元件状态的状态向量，t为时间，Ax和f(x,t)分别表示电力系统动态模型中的线性和非线性部分，它们由各个子系统电磁暂态仿真模型和子系统之间互联关系确定，y是仿真使用者所要求的仿真输出向量，输出函数g(x,t)根据仿真研究关注点不同，由使用者任意指定。

所述的根据非线性元件与线性网络子系统之间的互联关系，由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，其具体实现方式为：

首先由输入文件读入非线性并网元件子系统与线性网络子系统之间的互联关系，所述互联关系指的是，一个子系统的输出变量对应着另一个子系统的输入变量，用数学公式表达为

联列线性网络子系统和非线性并网元件子系统的电磁暂态仿真模型，并将输入和输出变量替换为子系统状态向量的函数，即得到待研究电力系统的整体状态空间仿真模型

其中，Ax和f(x,t)分别是系统仿真模型中的线性和非线性部分，它们与子系统电磁暂态仿真模型参数的关系是：

3)设定：仿真时间T，仿真步长Δt，非线性迭代收敛精度∈，设定当前时刻t_n为仿真起始时刻t₀，依照仿真需要，设置仿真初值x₀，并赋值给当前时刻状态向量x_n，计算仿真起始时刻输出向量y₀＝g(x₀,t₀)，并写入输出文件；

4)计算当前时刻t_n时，整体仿真模型中的非线性项f(x_n,t_n)，并计算 其中，e^ΔtA为矩阵指数函数表示的单位步长下系统线性状态转移矩阵；

5)记下一个时刻的时间值t_n+1＝t_n+Δt，采用不动点迭代方法求解非线性方程

至设定的非线性迭代收敛精度∈，得到t_n+1时刻的状态向量x_n+1；

所述的采用不动点迭代方法求解非线性方程，包括如下步骤：

(1)设定最大迭代次数N，外插法设置迭代初值记当前已迭代步数k＝0；

(2)比较已迭代步数k与最大迭代步数N，如果k≥N，跳转至步骤(4)，否则，使用如下公式进行第k+1步迭代：

并更新已迭代步数k＝k+1；

(3)计算当前步与前一步迭代值的偏差量并比较这一偏差dx与非线性迭代收敛精度∈，如果dx>∈，跳转到步骤(2)继续进行迭代运算，否则进行到下一步骤；

(4)将当前的迭代值作为t_n+1时刻的状态量值，即完成当前时步的迭代运算。

6)由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出文件，更新当前时刻为下一时刻t_n＝t_n+1，仿真向前推进一个步长；

7)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

以一个典型并网分布式发电系统的暂态仿真作为实施例，其系统结构如图1所示。分布式发电单元容量为84kW，经单级逆变器接入0.4kV低压配网。直流侧建模为理想直流电压源，逆变电路采用平均值模型建模，其控制采用dq坐标系下的单环功率控制，滤波电路为LCL型。其基于矩阵指数电磁暂态仿真的非线性求解如下：

1)划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，选取并网点端口电压和电流作为网络与并网元件的接口变量，在状态分析框架下，分别建立待研究电力系统中线性网络子系统、以及各并网非线性元件的子系统电磁暂态仿真模型。线性网络子系统的电磁暂态仿真模型形式为：

其中，x_net是电网络的状态向量，包含网络中独立的电感电流和电容电压，A_net是电网络子系统的状态矩阵；u_net是来自并网非线性元件的输入电流，B_net是对应的输入矩阵；y_net是电网络与并网非线性元件接口处的端口电压(相-相)，C_net是对应的输出矩阵。

本实施例包含1个并网非线性元件，分布式发电系统的电磁暂态仿真模型为：

其中，x_DG是这个分布式发电单元的状态向量，包含LCL滤波器中的电感电流和电容电压，以及控制器中的相关状态向量；u_DG是分布式发电单元与电网络的接口处的端口电压(相-相)，A_DG是该子系统动态特性中的线性部分，而f_DG(x_DG,u_DG,t)是这个元件的非线性特性；r(t)是分布式发电单元的电流指令，B_r是对应这一指令信号的输入矩阵；y_DG是分布式发电单元LCL滤波器上的输出电流，g_DG(x_DG,u_DG)是对应的输出函数；

2)根据分布式发电单元与电网之间的互联关系，即

由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，所述整体状态空间仿真模型形式为：

其中，x是包含当前时刻电力系统内所有储能元件与记忆元件状态的状态向量，t为时间，Ax和f(x,t)分别表示电力系统动态模型中的线性和非线性部分，它们与电网络子系统和分布式发电单元子系统电磁暂态仿真模型参数的关系是：

输出函数g(x,t)根据仿真结果显示需要，由使用者任意指定。本实施例中选取分布式电源向电网输出的有功功率，以及LCL滤波器输出电感上的A相电流，作为输出变量；

3)设定：仿真时间T，仿真步长Δt，非线性迭代收敛精度∈，设定当前时刻t_n为仿真起始时刻t₀，依照仿真需要，设置仿真初值x₀，并赋值给当前时刻状态向量x_n，计算仿真起始时刻输出向量y₀＝g(x₀,t₀)，并写入输出文件；

4)计算当前时刻t_n时，整体状态空间仿真模型中的非线性项f(x_n,t_n)，并计算 其中，e^ΔtA为矩阵指数函数表示的单位步长下系统线性状态转移矩阵；

5)记下一个时刻时间值t_n+1＝t_n+Δt，采用不动点迭代方法求解非线性方程

直到满足设定的非线性收敛精度∈，得到t_n+1时刻的状态向量x_n+1。具体实现方式如下：

(1)设定最大迭代次数N，外插法设置迭代初值记当前已迭代步数k＝0；

(2)比较已迭代步数k与最大迭代步数N，如果k≥N，跳转至步骤(4)，否则，使用如下公式进行第k+1步迭代：

并更新已迭代步数k＝k+1；

(3)计算当前步与前一步迭代值的偏差量并比较这一偏差dx与非线性迭代收敛精度∈，如果dx>∈，跳转到步骤(2)继续进行迭代运算，否则进行到下一步骤；

(4)将当前的迭代值作为t_n+1时刻的状态量值，即完成当前时步的迭代运算。

6)由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出文件，更新当前时刻为下一时刻t_n＝t_n+1，仿真向前推进一个步长；

7)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

本发明的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，执行仿真计算的计算机硬件环境为Intel Core2Q8400 2.66GHz CPU，内存容量2GB；软件环境为Windows7操作系统。

分别将本发明的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法(MatrixExponential)与典型隐式方法——梯形法(Trap)进行比较，并采用有误差控制的变步长Adams方法作为仿真的基准值。对于矩阵指数方法和隐式梯形法，仿真步长设置为5us。仿真场景设置为：零状态启动，依单位功率因数满功率运行；0.2s功率指令改变为Pref＝0.2p.u.，Qref＝-0.2p.u.。

附图3给出分布式发电单元向网络输出有功曲线，在该图中，使用梯形法和本专利提出的方法所得的仿真曲线与基准曲线基本一致，难以区分。为了进一步比较，附图4是附图3曲线的局部放大图，从图中可以看出，使用本专利提出的方法所得的仿真结果与基准值更接近。附图5给出两种数值方法所得结果与基准值之间的误差，从图中可以看出，本专利提出的方法与梯形法相比具有更低的数值误差。

附图6给出滤波器出口处A相电流的仿真波形，在该图中，使用梯形法和本专利提出的方法所得的仿真曲线与基准曲线基本一致，难以区分。为了进一步比较，附图7是附图6曲线的局部放大图，从图中可以看出，使用本专利提出的方法所得的仿真结果与基准值更接近。附图8给出两种数值方法所得结果与基准值之间的误差，从图中可以看出，本专利提出的方法与梯形法相比具有更低的数值误差。

从上述数值比较可以看出，虽然矩阵指数方法与梯形法同为二阶方法，但其数值误差明显更小。这是因为矩阵指数方法可以精确地处理微分方程中的线性部分，没有截断误差，其误差仅来源于非线性部分的近似。

以上算例测试结果证明，本发明提出的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法具有良好的可行性与适用性，为解决电力系统电磁暂态仿真提供了一种很好的解决思路。

Claims (6)

1.一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，其特征在于，通过子系统互联形成包含非线性元件的研究系统整体仿真模型，结合矩阵指数方法的计算精度优势以及不动点迭代的计算速度优势，实现具有强非线性特征的电力系统电磁暂态仿真，具体包括如下步骤：

1)划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，选取并网点端口电压和电流作为网络与并网元件的接口变量，在状态分析框架下，分别建立待研究电力系统中线性网络子系统、以及各非线性并网元件子系统的电磁暂态仿真模型；

2)根据非线性元件与线性网络之间的子系统互联关系，由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，所述整体状态空间仿真模型形式为

其中，x是包含当前时刻电力系统内所有储能元件与记忆元件状态的状态向量，t为时间，Ax和f(x,t)分别表示电力系统动态模型中的线性和非线性部分，它们由各个子系统电磁暂态仿真模型和子系统之间互联关系确定，y是仿真使用者所要求的仿真输出向量，输出函数g(x,t)根据仿真研究关注点不同，由使用者任意指定；

3)设定：仿真时间T，仿真步长Δt，非线性迭代收敛精度∈，设定当前时刻t_n为仿真起始时刻t₀，依照仿真需要，设置仿真初值x₀，并赋值给当前时刻状态向量x_n，计算仿真起始时刻输出向量y₀＝g(x₀,t₀)，并写入输出文件；

4)计算当前时刻t_n时整体仿真模型中的非线性项f(x_n,t_n)，并计算 其中，e^ΔtA为矩阵指数函数表示的单位步长下系统线性状态转移矩阵；

5)记下一个时刻的时间值t_n+1＝t_n+Δt，采用不动点迭代方法求解非线性方程

至设定的非线性迭代收敛精度∈，得到t_n+1时刻的状态向量x_n+1；

6)由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出文件，更新当前时刻为下一时刻t_n＝t_n+1，仿真向前推进一个步长；

7)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

2.根据权利要求1所述的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，其特征在于，步骤1)所述的划分待研究电力系统为线性网络子系统与若干非线性并网元件子系统，具体划分范围是：其中，PI型等效电路建模的输配电线路、不考虑饱和效应的电力变压器、线性恒阻抗负荷、传统串补设备、线路高抗、无功补偿电容器以及调压器设备在线性网络子系统内统一建模；而同步和异步发电机组、电力电子器件、可再生能源电源以及电弧炉非线性负荷设备作为非线性并网元件子系统分别建模。

3.根据权利要求2所述的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，其特征在于，所述的在线性网络子系统内统一建模，得到的子系统电磁暂态仿真模型具体形式为：

其中，x_net是线性网络子系统的状态向量，包括独立电感电流和电容电压，A_net是线性网络子系统的状态矩阵；u_net,i是来自第i个非线性并网元件子系统的输入变量，取决于具体接口方式，视为受控的电压源或电流源，B_net,i是对应的输入矩阵；y_net,i是线性网络子系统给第i个非线性并网元件子系统的输出变量，取决于具体接口方式，选取端口电压或输出电流，C_net,i是对应的输出矩阵。

4.根据权利要求2所述的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，其特征在于，所述的对非线性并网元件子系统分别建模，具体是指：

待研究电力系统包含有多个非线性并网元件，对于第i个元件，非线性并网元件子系统电磁暂态仿真模型形式为

其中，x_i是所述非线性并网元件子系统的状态向量，A_i是所述非线性并网元件子系统仿真模型中的线性部分；u_i是所述非线性并网元件子系统来自线性网络子系统的输入变量，B_i是对应的输入矩阵；f_i(x_i,u_i,t)是所述非线性并网元件子系统的非线性特性函数；y_i是所述非线性并网元件子系统给线性网络子系统的输出变量，C_i是对应的输出矩阵，g_i(x_i,u_i,t)是输出变量的非线性特性函数。

5.根据权利要求1所述的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，其特征在于，步骤2)所述的根据非线性元件与线性网络子系统之间的互联关系，由子系统电磁暂态仿真模型生成待研究电力系统的整体状态空间仿真模型，具体实现方式为：

首先由输入文件读入非线性并网元件子系统与线性网络子系统之间的互联关系，所述互联关系指的是，一个子系统的输出变量对应着另一个子系统的输入变量，用数学公式表达为

联列线性网络子系统和非线性并网元件子系统的电磁暂态仿真模型，并将输入和输出变量替换为子系统状态向量的函数，即得到待研究电力系统的整体状态空间仿真模型

其中，Ax和f(x,t)分别是系统仿真模型中的线性和非线性部分，它们与子系统电磁暂态仿真模型参数的关系是：

6.根据权利要求1所述的一种基于矩阵指数电磁暂态仿真的非线性建模求解方法，其特征在于，步骤5)所述的采用不动点迭代方法求解非线性方程，包括如下步骤：

(1)设定最大迭代次数N，外插法设置迭代初值记当前已迭代步数k＝0；

(2)比较已迭代步数k与最大迭代步数N，如果k≥N，跳转至步骤(4)，否则，使用如下公式进行第k+1步迭代：

并更新已迭代步数k＝k+1；

(3)计算当前步与前一步迭代值的偏差量并比较这一偏差dx与非线性迭代收敛精度∈，如果dx>∈，跳转到步骤(2)继续进行迭代运算，否则进行到下一步骤；

(4)将当前的迭代值作为t_n+1时刻的状态量值，即完成当前时步的迭代运算。