光伏发电单元电磁暂态解耦与快速仿真方法

姚蜀军，张春强，刘 刚，马嘉昊，汪 燕

（华北电力大学电气与电子工程学院，北京市 102206）

0 引言

近年来，随着能源危机与环境问题的加剧，新能源发电不断受到重视，光伏新增装机容量不断攀升，电站规模向着大型化、集中化、光伏基地方向发展［1-2］。但是，光伏电站中所含非线性电池组件和电力电子换流器等设备数目较多，基于详细模型的电磁暂态仿真的仿真速度较慢［3-4］。

电力电子换流器作为光伏电站的重要设备，其详细模型精度较高，但仿真速度慢。为了提高仿真速度，文献［5-6］提出基于电感/电容等效的恒导纳开关模型，计算过程中开关切换时无须重新形成节点导纳矩阵，提高了仿真速度。然而，由于该模型存在暂态过程，会在开关动作后产生大于实际的虚拟功率损耗［7］。文献［8］进一步提出了基于响应匹配的恒导纳换流器建模方法，该方法可以消除暂态过程导致的误差，但在开关状态切换时会存在数值振荡问题。外特性建模仿真精度不如详细模型，但仿真速度快。 文献［9］提出了电压源型换流器（voltage source converter，VSC）的平均值模型，文献［10］提出了换流器的分段平均化模型，文献［11］提出了分段广义状态空间平均模型。这类模型的仿真步长会受到窗口长度（开关周期）的限定。文献［12-13］建立了VSC 的动态相量模型，文献［14］提出动态/扩展谐波域法，以上两种方法通过傅里叶变换将实数高频信号变换为复数低频信号，从而增大了仿真步长。当计及足够阶数时，动态相量法可接近详细模型的精度。但在实际应用中，为了保证仿真速度会忽略高频谐波和高阶向量，导致模型的谐波截断误差较大。

考虑光伏发电系统的建模问题，文献［15］采用详细模型建模，该模型仿真精度高、内部故障及变量易于观察，但模型复杂、仿真需要的时间长，不能满足大电网仿真的需要。文献［16-18］使用容量加权法对光伏发电单元进行单机等效，等效模型简单，可以适用于一些较为严重的运行工况，但只能用于精度要求不高的场合。文献［19-20］使用参数辨识法进行单机等效，该方法提高了仿真精度，但无法顾及逆变器类型及光伏电站内部结构。针对光伏电站采用单机等值误差较大的情况，文献［21-23］分别采用不同的方法以逆变器为核心对光伏组件进行分群聚类，建立了光伏电站的多机等效模型，该类方法考虑了不同群组的运行参数，也能较好地模拟光伏电站的动态响应过程，但如何基于实际工程易于观测的量来提取分群指标值，以及寻找合适的分群方法仍有待进一步研究。

文献［24］提出电磁暂态半隐式延迟解耦并行仿真方法（semi-implicit latency decoupling and paralleling，SILDP），并验证了其应用于交流电网的高效性和准确性，附录A 中给出了该方法的原理。本文将其应用于光伏电站的电磁暂态解耦与快速仿真方法研究，通过建立光伏发电单元各设备的解耦模型，既实现了各设备间的并行计算，又保证了开关动作时节点导纳矩阵不变，提高了计算效率。最后，通过光伏电站的算例验证了本文所建立模型的有效性和准确性。

1 光伏发电单元的解耦

文献［24］提出了半隐式延迟解耦法，并给出了其在交流系统分网中的应用。本文利用该方法实现了光伏发电单元的解耦。

1.1 光伏电站的结构

光伏发电系统由光伏阵列将太阳能转化为直流电能，直流电能经逆变器转换为交流电能并网发电。如附录B 图B1 所示，光伏阵列、逆变器、滤波器、变压器共同构成一个光伏发电单元，大量光伏发电单元并联汇流组成光伏电站。

1.2 光伏阵列的数学模型

光伏电池模型如图1 所示。该等效电路包含电流源Isc、反并联的二极管D、分流电阻Rsh和串联电阻Rsr。当光伏电池暴露于光照下时，二极管两端产生电压Vd，流过分流电阻Rsh的电流为Ish，流过二极管的电流Id产生光伏电池的非线性特性，输出直流电流I随光照强度而变化。

由文献［25］可知，光伏电池的输出电流与输出电压之间有如下关系：

式中：q为电子电荷；k为玻尔兹曼常数；η为二极管理想因子；V为光伏电池输出电压；Isc为光电流，是太阳能电池平面上的太阳辐射G和电池表面温度Tc的函数，即

式中：IoR为参考温度下的暗电流；eg为太阳能电池材料的带隙能量。

1.3 Boost 换流器解耦模型

图2（a）为Boost 换流器拓扑，该电路由电感L、二极管VD、输出电容C，以及反并联二极管D 与绝缘栅双极型晶体管T 构成的开关组VT 共同组成。电路中vi为输入电压，iL为输入电流，vc为输出电压，io为输出电流，二极管VD 及开关组VT 均采用二值电阻模型，导通时电阻Ron=0.01 Ω，关断时电阻Roff=106Ω。Boost 换流器的二值电阻等效电路如图2（b）所示，RVD和RVT分别表示二极管及开关组的等值电阻。

以图2 中的电感电流和电容电压为状态变量，列写Boost 连续工作模式下的状态方程：

二极管与开关组通常互斥导通，令Rsum=RVT+RVD=Ron+Roff，Rmul=RVTRVD=RonRoff，则式（4）的状态空间表达式为：

式中：Req为解耦电路中输入侧的等效电阻，Req=Rmul/Rsum；Geq为解耦电路输出侧的等效并联导纳，Geq=1/Rsum；kv和ki分别为受控电压源和受控电流源的 系 数，kv=ki=RVT/Rsum。由 于Roff≫Ron，若 将Roff看作无穷大，则Req≈Ron，Geq≈0。

根据半隐式延迟解耦法，对式（5）进行矩阵分裂并计及Geq≈0，则可得

根据半隐式延迟解耦法，对式（6）进行半步延时，得到其差分方程为：

式中：上标n-1/2、n、n+1/2和n+1分别表示变量在第n-1/2、n、n+1/2 和n+1 时刻的值；Δt为仿真步长。由式（6）和式（7）可得图2（c）所示的Boost 换流器半隐式延迟解耦等效电路图。图中：Veq为解耦电路输入侧的受控电压源，Veq=kvvc；Jeq为解耦电路输出侧的受控电流源，Jeq=kiiL。开关的切换表现为kv、ki的改变，而开关损耗表现为电阻Req的损耗。各状态变量的迭代式见附录C 式（C1）。

1.4 VSC 的解耦模型

VSC 结构如图3（a）所示。图中：T1～T6为各桥臂绝缘栅双极型晶体管；D1～D6为各桥臂反并联二极管；iR1～iR6为 各桥臂电流；La、Lb和Lc为交流侧等效电感；va、vb、vc为交流侧相电压；ia、ib、ic为交流侧流入电流；O点为接地点；C1和C2为直流侧正负极电容；vc1和vc2为直流侧电容电压；ic1和ic2为直流侧电容电流；id1和id2为直流侧输入电流。图3（b）为逆变器的等效电路，反并联的绝缘栅双极型晶体管和二极管构成的开关组以二值电阻R1～R6代替。

以电感电流和电容电压为状态变量，列写图3电路的状态方程可得：

式中：Geq1和Geq2为解耦电路中直流侧并联导纳，Kvc1和Kvc2为解耦电路中直流侧并联电流源的系数，Geq1=Geq2=Kvc1=Kvc2=1/（R1+R4）+1/（R3+R6）+1/（R5+R2）。

正常工作时，由于VSC 中每相的上、下桥臂互斥 导 通，令R1+R4=R3+R6=R5+R2=Ron+Roff，则可得式（8）的状态空间方程为：

ReqA、ReqB、ReqC分别为VSC 解耦电路交流侧A、B、C相的等效电阻；Kv1～Kv6为解耦电路中交流侧受控电压源系数；Ki1～Ki6为解耦电路中交流侧受控电流源系数。类似的，由于Roff≫Ron，如将Roff看作无穷大，则ReqA=ReqB=ReqC≈Ron，Geq1、Geq2、Kvc1、Kvc2均 可 近似为0。

根据半隐式延迟解耦方法，状态变量以电感电流和电容电压分组，对式（9）进行矩阵分裂可得：

根据半隐式延迟解耦方法，可得式（10）的半步延时差分方程为：

与式（10）和式（11）对应的解耦电路如图3（c）所示。图中：VeqA、VeqB、VeqC分别为解耦电路交流侧A、B、C 相的受控电压源；Jeq1和Jeq2分别为VSC 解耦电路直流侧的两个受控电流源。开关的切换表现为Kv1～Kv4、Ki1～Ki4的 改 变，而 开 关 损 耗 表 现 为 电 阻ReqA、ReqB和ReqC的损耗。附录C 式（C2）给出了相关变量的表达式和各状态变量的迭代式。

1.5 LC/LCL 滤波器的解耦模型

光伏发电单元中，LC 滤波器和与其相连的电感（或变压器）一起构成LCL 电路。为了实现细粒度解耦，进一步提高计算的并行性，可按文献［24］中的T 形电路进行解耦，其解耦模型见附录B 图B2。

1.6 光伏发电单元解耦电路

如图4 所示，分别对LC 滤波器、Boost 换流器和VSC 进行解耦，得到各设备的半隐式延迟解耦等效电路，前一设备的输出与后一设备的输入在计算上属于同一时刻，共同构成子电路，可以将光伏发电单元解耦成以受控源为激励的多个独立子系统。图中：Req为Boost 解耦电路输入侧的等效电阻；V1～V6为解耦电路中的受控电压源；I1～I6为解耦电路中的受控电流源；Veq1～Veq6为传递到受控电压源的电压值；Jeq1～Jeq6为传递到受控电流源的电流值。如果整个光伏电站内的发电单元都按这种方式解耦，将极大地提高光伏电站的仿真效率。

1.7 解耦模型的稳定性分析

关于本文方法的数值稳定性，文献［26］已进行了分析，这里以Boost 电路为例，简述如下。

Boost 换流器的差分方程式（7）可重写为：

B1为矩阵u1对应的系数矩阵。

由线性系统的稳定性判据可知，若要保持系统稳定，需要满足矩阵A1的所有特征值的幅值应小于1。当仿真步长Δt取50 μs 时，可得矩阵A1的特征值幅值|λ（A1（t））|=0.995＜1，该参数可以保证稳定性。对于VSC 采用同样的方法，其证明过程见附录D。可见，在Boost 电路及VSC 均满足稳定性判据的情况下，采用本文参数的串联电路仍然满足稳定性。因此，当光伏电站运行在算例所示的常用参数附近时，本文方法可满足稳定性。

1.8 计算时序及流程

将图4 解耦电路中的各子系统按受控电流源和受控电压源分成两组。不失一般性，令x1包含所有受控电流源子系统；x2包含所有受控电压源子系统。按照半隐式延迟解耦法的并行计算原理，x1和x2两组变量之间使用1/2 步长顺序求解，而x1和x2内部的变量之间可并行求解。该方法的时序图如附录E图E1 所示，将各独立子电路依次编号为子系统1～N，并行计算流程如附录E 图E2 所示。

2 模型特点分析

本文所建立的光伏电站半隐式延迟解耦模型的特点如下：

1）传统的平均值和动态相量模型难以计及设备的内部特性，开关函数模型可计及内部动作过程，但忽略了换流器损耗。此外，这3 种模型交直流侧均耦合，没有解耦。而本文的解耦模型在实现交直流侧解耦的同时，还考虑了解耦模型内部的动态过程及开关的损耗，计算结果接近于详细模型。

2）发电单元各部分解耦，在实现大系统的分解从而减小计算规模的同时，子系统间又可并行，进一步提高了计算效率。此外，解耦模型中，等效电阻Req≈Ron为定值，电导Geq≈0，开关动作导致的系统拓扑变化只体现在受控源的系数上。

3）PSCAD 为了避免在开关动作时非状态变量突变引发的数值振荡，采用了半步插值法。每次开关动作都需进行插值回退，EMTP-RV 以及PSModel 采用了半步长的后退欧拉法，而半隐式延迟解耦法对系统的状态变量（电感电流和电容电压）进行分组和延时，解耦的状态变量间不会由于开关动作而突变，数值稳定性好，当换流器开关状态改变时无须切换积分格式而改变解耦形式，也无须进行后向的半步插值。因此，能始终保持计算过程中解耦的一致性和并行性。

4）解耦时采用的中心积分形式，其面积与隐式梯形积分的面积等效（电磁暂态仿真尺度下），也使得解耦模型与详细模型的精度相仿。

5）与LC 恒导纳模型相比，本文方法在保证节点导纳矩阵恒定的同时不会引入小电感及电容。因此，在开关切换过程中不存在虚拟损耗的问题，在相同精度下可以采用相对较大的步长进行仿真。

3 算例验证

以图5 所示拓扑作为算例验证系统。图中：LD1为直流侧滤波电路电感；C0为直流侧滤波电路电容；LD2为Boost 电 路 电 感；C1为Boost 电 路 电 容；Ra1、Rb1、Rc1为交流侧LC 滤波电路电阻；La1、Lb1、Lc1为交流侧LC 滤波电路电感；Ca、Cb、Cc为交流侧LC滤波电路电容；La2、Lb2、Lc2为变压器及线路的等效总电感；Ra2、Rb2、Rc2为变压器及线路的等效总电阻；Rlim为限流电阻；LS为外电路电感；RS为外电路电阻；AC 为外部无穷大电源。采用C++编程实现本文解耦模型，系统参数如表1 所示。计算结果与详细模型的仿真结果进行对比。计算平台为PC机，CPU 为Intel Core i5-6300HQ（4 核心、4 线程），RAM 为16 GB。

表1 并网光伏发电单元参数Table 1 Parameters of grid-connected photovoltaic power generation unit

仿真过程设置如下：仿真总时长为1.5 s。开始时刻，接入限流电阻Rlim（即断开断路器B）；0.05 s 时切断限流电阻Rlim（即闭合断路器B）；0.6 s 时交流侧负载处发生a 相接地短路故障，持续时间为20 ms；1 s 时交流侧负载处发生三相接地短路故障，持续时间为20 ms；以上故障发生后，系统中均未采取任何保护措施。光伏电源及Boost 电路采用最大功率点跟踪（maximum power point tracking，MPPT）算法。

3.1 仿真精度验证

以国际通用的电磁暂态仿真程序PSCAD/EMTDC 的计算结果为基准，计算本文方法与PSCAD 仿真波形的均方根相对误差，其中均方根相对误差的计算公式详见附录F。

开关频率设置为20 kHz，仿真波形对比见附录G 图G1。无论是在暂态还是稳态情况下，各处的电压、电流以及有功功率误差均在0.6%以内，可以满足仿真精度的要求。为了验证大步长情况下解耦模型的准确性，以PSCAD 中1 μs 步长的仿真结果为基准，列出解耦模型在1 μs 及50 μs 步长下的最大误差，作为对比列出PSCAD 在50 μs 步长下的仿真结果，如表2 所示。

表2 不同仿真步长下的最大误差Table 2 Maximum error with different simulation steps

由表2 可见，在50 μs 的仿真步长下，负载处有功功率解耦模型误差略大于详细模型，其他变量的误差二者接近。因此，该解耦模型在较大步长的情况下仍然可以保证与详细模型接近的仿真精度。

相比于传统的开关函数模型、平均值模型、动态相量模型等方法，本文模型保留了内部开关损耗信息。附录G 图G2 为开关频率为20 kHz、仿真步长为1 μs 的情况下，详细模型、半隐式延迟解耦模型以及开关函数模型的交流侧与直流侧功率比较。可见，本文提出的半隐式延迟解耦模型，其交流侧的有功和无功功率以及直流侧有功功率误差相比于开关函数模型具有显著优势。

3.2 仿真效率验证

分别在不同仿真步长、不同开关频率以及不同发电单元个数情况下进行仿真效率对比。在解耦模型并行求解过程中，使用OpenMP 将解耦后形成的多个小系统的求解任务分配至4 个核心。为了减小控制部分求解对结果的影响，在3 个算例中采用同样的方法将控制系统中光伏电源及不同换流器的PQ控制等自然解耦部分分配给4 个核心并行计算。在每一步的计算过程中，先进行控制系统的求解，控制系统求解完成后求解一次系统。本算例中总仿真时间为1 s，表3 所示为不同仿真步长下两种模型的CPU 用时和加速比。表4 所示为不同开关频率下两种模型的CPU 用时和加速比。表5 所示为不同光伏发电单元个数下的CPU 用时和加速比。表中，SSR（serial speedup ratio）为C++实现的详细模型串行计算用时与C++实现的解耦模型串行方式的CPU 用 时 之 比，PSR（parallel speedup ratio）为C++实现的详细模型串行计算用时与C++实现的解耦模型并行仿真的CPU 用时之比。

表3 不同仿真步长下的用时和加速比Table 3 Time consumption and acceleration ratio with different simulation steps

表4 不同开关频率下的用时和加速比Table 4 Time consumption and acceleration ratio at different switching frequencies

表5 不同发电单元个数下的用时和加速比Table 5 Time consumption and acceleration ratio with different number of power generation units

由表3 至表5 可以看出，本文所提方法在不同仿真步长、不同开关频率及不同发电单元个数的情况下都可以较大地提高仿真速度。当发电单元个数增加时，解耦并行计算相对于解耦串行计算有较大优势。解耦电路在求解过程中无须计算大矩阵求逆，可以预存LU 分解的结果，每一步只需多个前代回代运算，故解耦串行计算得以加速。在并行计算中将任务同时分配至4 个核心，相比于解耦串行计算理论加速比可以接近4。实际仿真结果并没有达到理论值，其主要原因是因为控制系统没有进行细粒度解耦，当一次系统进行细粒度解耦后，控制系统的求解用时成为限制进一步加速的主要因素。如果对控制系统进行细粒度解耦，使用核心数更高的CPU、GPU 计算，仿真速度提升将会更加显著。

4 结语

本文基于电磁暂态半隐式延迟解耦并行仿真方法，建立了光伏发电单元各换流环节和整体的半隐式延迟解耦模型，并进行了算例验证。所提出的解耦模型，一方面通过解耦，将发电单元分解以减小系统规模，并实现子系统间并行计算以提高仿真效率；另一方面，解耦电路的导纳矩阵恒定，在含有大量高频电力电子开关的电路仿真中，极大地简化了计算过程。当光伏电站进行扩展时，通过对电站不同部分进行划分，只要可以列写各部分设备及拓扑的状态方程，仍可用本文方法进行等效。

算例结果表明，本文提出的解耦模型能够兼顾仿真的精度与效率。当光伏发电系统规模增大时，该模型计算效率高、系统导纳矩阵恒定的优势更加显著。此外，本文提出的解耦方法实现了光伏发电单元的细粒度解耦，既可应用于离线仿真，也为实时仿真提供了算法支持，可基于该方法研究相应的实时仿真器。另外，本文模型及算法难以在已有的仿真软件中直接实现，兼容本文方法的通用电磁暂态仿真软件将会在后续工作中进行开发。

附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx），扫英文摘要后二维码可以阅读网络全文。