CN111900718B - 基于多级优化追赶变分迭代法的有源配电网动态仿真方法 - Google Patents

Abstract

本发明公开了一种基于多级优化追赶变分迭代法的有源配电网动态仿真方法，首先将变分迭代算法进行两个方面的改进：一方面，在校正泛函的右侧计及已经求得的近似解析解；另一方面，优先计算被耦合频率更高的微分方程。从而提出了离线收敛更快的优化追赶变分迭代法。其次，基于该法求出了有源配电网中光伏系统以及同步发电机系统中动态元件对应的近似解析解。随后，引入多级机制，即在连续时间间隔上应用近似解析解，突破了原始变分迭代法有限的收敛域，从而使其可以应用到对更大时间跨度内要求有准确解的有源配电网仿真中。该方法，比传统数值积分方法求解速度更快。

Description

基于多级优化追赶变分迭代法的有源配电网动态仿真方法

技术领域

本发明公开了一种基于多级优化追赶变分迭代法的有源配电网动态仿真方法，属于电力系统仿真与分析的技术领域。

背景技术

近年来，可再生能源正在朝着未来全球能源结构的主体迈进。像光伏这种分布式可再生能源大规模接入配电网，给配电网带来了潮流和电压分布的变化，使得单电源辐射型配电网向多电源的有源配电网方向发展。高渗透率的分布式可再生电源将会给配电网带来诸多电力系统问题，例如光伏电站是一个高度非线性的系统，其提供的电能是随机和不确定的，这会给接入的配电网在电能质量、继电保护与控制以及电压稳定性等方面带来严重的隐患，从而进一步影响配电网的规划、设计和运行。例如，像太阳辐照度这样的小扰动会导致诸如节点电压波动的电能质量问题。而诸如短路故障这样的大扰动，会导致在公共耦合点(Point of Common Coupling,PCC)处产生更加明显的电气量振荡及其恢复过程。此外，近年来，各种不同的分布式电源如分布式光伏、小型水电站等大规模化接入配电网后，其复杂和差异巨大的动态运行特性对所接入配电网遭受到各种扰动之后的动态稳定性造成了严重的影响。开展有源配电网的动态仿真工作是分析研究上述影响的基础，也是保障配电网安全稳定运行的有效工具。

但是，含高渗透率分布式电源的有源配电网的动态仿真研究工作面临着诸多的技术难题和挑战。对于有源配电网这个高维非线性系统，传统的显式和隐式积分算法，一方面在动态仿真运算中，由于迭代次数过多而导致仿真效率低下；另一方面，在多种分布式电源接入的有源配电网这种刚性微分系统的计算环境下，传统积分算法需要针对刚性特性匹配以更小的步长，从而难以兼顾动态仿真的计算效率和数值稳定性。因此，针对分布式光伏电站集群的模型维数高、非线性以及刚性强的仿真计算环境，为了克服传统数值积分算法的缺陷，开展有源配电网的动态仿真算法研究，提出兼顾稳定性的高效率仿真算法显得非常关键。

发明内容

本发明的发明目的是针对上述背景技术的不足，提供了一种基于多级优化追赶变分迭代法(Multilevel Optimization Catch-up Variational Iteration Method,MOCVIM)的有源配电网动态仿真方法，既具备隐式积分算法的数值稳定性和数值精度，在在线仿真计算时又无需进行迭代计算，从而加速了有源配电网的动态仿真计算。

本发明为实现上述发明目的采用如下技术方案：

一种基于多级优化追赶变分迭代法的有源配电网动态仿真计算方法。该方法分为两个阶段，即初始化阶段和动态过程计算阶段。

初始化阶段为前4个步骤：

第1步，建立有源配电网仿真计算所需的元件模型，并输入有源配电网的系统参数。

包括光伏发电系统模型，同步发电机模型，电网络模型等等；形成整个系统的微分方程组和代数方程组，构成全系统动态仿真数学模型，由式(1)所示的微分-代数方程组进行描述：

其中，x和y分别是代数变量和状态变量，f和h表示函数，分别为n维和m维。

输入的参数包括母线和线路的参数，光伏电站和同步发电机的参数以及所在母线的位置，系统频率f_grid，以及系统的基准电压和基准容量等。

第2步，潮流计算。

紧接着，进行潮流计算。获得整个系统的稳态值矩阵，包括母线电压v_xy＝v_x+j·v_y，负荷的节点功率s_l＝p_l+j·q_l以及各个分布式电源的节点注入功率s_g＝p_g+j·q_g等等。随后基于潮流结果，按照下面的公式计算负荷等值并联导纳y_load，并形成计及负荷的系统导纳矩阵Y。其中v是负荷节点电压的幅值。

第3步，计算所有微分方程组状态变量的初值。

根据发电机和光伏系统的物理特性计算所有微分方程组状态变量的初值。以同步发电机为例进行说明，光伏系统的计算与之类似。有源配电网受到扰动之前，根据潮流计算可以知道各个同步发电机的节点注入功率和节点电压，于是可以计算同步发电机的注入电流：

所以同步发电机的转子角初值δ₀计算如下。

δ₀＝angle{v+(r_a+j·x_q)·i_xy} (4)

此外，所有同步发电机的转速初值ω₀设置为1。

第4步，输入仿真时序数据。

包括仿真总时间t_sum，仿真步长h，第i个扰动事件的开始时刻t_start和结束时刻t_end，其中i＝1,2,3,…。

下面8个步骤(第5步-第12步)为动态过程计算阶段。

动态过程计算的本质是计算仿真时间内每一个时刻的有源配电网中微分方程组的状态变量值，除去初始时刻，共需要进行K_sum＝t_sum/h个仿真时步的计算。对于某一具体的仿真时步k＝1,2,…,K_sum，其计算可以分为进行。

第5步，计算k时刻的分布式电源的节点注入电流和节点电压。

首先需要求取同步发电机的虚拟注入电流i_xy和等值导纳矩阵Y_G。为方便后续节点电压的计算，Y_G要被并入Y形成新的系统导纳矩阵Y’，从而可以将虚拟电流i_xy看作新的节点注入电流，计算如下。

其中，E′_q是暂态电势，导纳b_x,g_y的计算如下。

E′_q＝v_q+r_a·i_q+x′_d·i_d (7)

x′_d是次同步电抗，x_q是q轴同步电抗，r_a是定子绕组电阻，i_d和i_q为和节点注入电流i_xy经过d-q坐标变换后的d轴和q轴分量。v_q是稳态时节点电压v的q轴分量。

Y_G的计算如下。

对于光伏系统而言，直接将上一时刻的i_d和i_q作为光伏接入节点的注入电流i_xy，这是因为i_d和i_q均为光伏系统微分方程组的状态变量，而状态变量是不会突变的。

如系统遇到扰动事件，就在本步骤中更新整个系统的导纳矩阵Y’即可。最后，综合各个同步发电机的注入电流i_xy以及计及Y_G和扰动事件的系统导纳矩阵Y’，使用Gauss迭代计算出k时刻的节点电压v_xy。

第6步，计算k时刻微分方程组与代数方程组之间的接口量。

该步骤的目的是将微分方程组中所有的接口量计算出来，这样就可以将这些量作为已知的常数，以便于在后面的步骤来专门来处理微分方程组的计算。根据上一步计算的k时刻的节点电压和电流，可以计算同步发电机的接口量：

p_m＝p_g+(i_x ²+i_y ²)·r_a (9)

类似的可以求出光伏系统的接口量包括，逆变器直流侧电流i_dc，PV阵列输出电流i_pv，逆变器调制电压V_dq，逆变器交流侧电压U_idq。

第7步，将微分方程组中的高阶微分方程转化为等值的1阶微分方程组。

考虑式(1)中如下的k阶微分系统：

y^(k)+f(y^k-1,y^k-2,...,y′,y)＝g(t) (11)

在绝大多数应用环境下，基于变分迭代法求取该式的拉格朗日乘子较为困难，为此，令y₁(t)＝y将该式转化为等值的k维1阶微分方程组如下所示。

对于该系统使用MOCVIM可较为容易地求出拉格朗日乘子。

对于有源配电网而言，可以推导出光伏系统和同步发电机系统的微分方程如下。其中，光伏系统选取应用广泛的双级式光伏电站的数学模型，其动态元件为9阶系统，将所有微分方程组整理为一阶微分方程的形式：

di_d/dt＝(ω_s/l_f)(-u_id+u_{sd_system})+ω_s·i_q (13)

di_q/dt＝(ω_s/l_f)(-u_iq+u_{sq_system})-ω_s·i_d (14)

di_dref1/dt＝u_dcref-u_dc (15)

du_rd1/dt＝i_dref-i_{d_ctrl} (16)

du_rq1/dt＝i_qref-i_{q_ctrl} (17)

du_pv/dt＝(i_pv-i_L)/C_pv (18)

di_L/dt＝[u_pv-(1-D)u_dc]/L_dc (19)

du_dc/dt＝[(1-D)i_L-i_dc]/C_dc (20)

dD₁/dt＝u_pv-u_m (21)

而同步发电机选取应用广泛的2阶同步发电机数学模型进行表述，其中微分方程组如下所示。

dσ/dt＝ω_s(ω-1) (22)

dω/dt＝[p_m0-p_e-D_g(ω-1)]/t_j (23)

其中i_d,i_q分别为光伏系统的实际输出电流的d轴和q轴分量，ω_s是电网角频率，l_f为滤波器电感，u_id和u_iq分别为逆变器输出电压的d轴和q轴分量，u_{sd_system},u_{sq_system}分别为并网电压的d轴和q轴分量，i_dref1是，u_dc和u_dcref分别是直流侧电压及其参考值，u_rd1,u_rq1分别为调制波的d轴和q轴分量，L_dc和C_dc分别为直流斩波器的电感和电容，D为开关的占空比，i_dc和u_dc分别为直流斩波器的输出电流和输出电压。σ和ω分别为转子角和转子角速度，p_m0和p_e分别是机械功率和电磁功率，D_g是阻尼系数，t_j是惯性常数。

第8步，给出MOCVIM的校正泛函格式。

以上个步骤的微分方程组式中最后一个微分方程为例，可得该微分方程求解第n+1步解析解的校正泛函如下：

与原始变分迭代法不同，该校正泛函将第n+1步已经计算出的所有其他状态变量的解析解作为已知量参与y_k(n+1)(t)的计算，该技术可以加速迭代解向真实解的收敛速度。

第9步，优化方程组的求解顺序。

应用第二步对微分方程组求解的重要一点就是安排方程组的求解顺序，优良的求解顺序可以进一步加速近似解析解的收敛速度。改进的校正泛函求解顺序如下所示。

这种安排顺序将会充分提高第n+1步的已知量的使用频率。上式中除了第一个校正泛函外，后面的每一个校正泛函的右侧均计及了n+1步的已知量。因此，在实际仿真应用中，仅需较少次数的迭代计算即可获得与传统变分迭代法更高阶时才具备的同等级精度的近似解析解。

第10步，求取y_k对应的拉格朗日乘子λ_k(s)。

在式(13)-(23)中，只有式(16)和(17)不与其他状态变量相互耦合，即可以直接求出解析解，而其他式子均需要进行MOCVIM推导。考虑到在有源配电网仿真运算中，在微分方程组中把代数方程组传过来的变量作为已知的常数，因此可将用下面的A₁-A₁₀来简化上述微分方程组：

A₁＝ω_s(-u_id+u_{sd_system})/l_f,A₂＝ω_s(-u_iq+u_{sq_system})/l_f,

A₃＝i_pv/C_pv,A₄＝-1/C_pv,A₅＝1/L_dc,A₆＝(D-1)/L_dc,

A₇＝(1-D)/C_dc,A₈＝-i_dc/C_dc,A₉＝(p_m0-p_e+D_g)/t_j,A₁₀＝D_g/t_j。

根据第8步可列出上述微分方程组的校正泛函，这里以i_d的校正泛函为例如下式所示，其拉格朗日乘子为λ₁。其余微分方程组的拉格朗日乘子为λ₂-λ₉形式：

为求取λ₁对上述方程两边求变分可得：

其中上方有横线的变量表示该变量为限制变分，对限制变分进行求变分等于零。

若随着迭代的进行，i_d(n+1)(t)的变分应该等于0，因此令δi_d(n+1)(t)＝0可以有限制条件：

从而可得λ₁＝-1。类似得也可以计算λ₂到λ₈均等于-1，而ω对应的

第11步，列出光伏电站和同步发电机中微分方程组相应的校正泛函，并求出各阶近似解析解。

根据第8步到第11步，即可有MOCVIM校正泛函及其求导次序如下所示。

利用该校正泛函，设状态变量对应的初值为I_d,I_q,I_dref1,U_pv,I_L,U_dc,D₁₍₀₎,σ₀,ω₀，令初值为其初始的迭代函数，代入上述校正泛函中，不断迭代即可得到不同阶的近似解析解。

第12步，基于计算获得的各微分方程组的近似解析解获取相应状态变量在第k个时段的动态轨迹。

获取相应状态变量在第k个时段的动态轨迹，也即可以获得k+1时刻的状态变量值y(k+1)。此外，状态变量u_rd1和u_rq1可以直接求出解析解，实际应用时使用下面的公式进行仿真计算：

u_rd1(k+1)＝(i_dref-i_{d_ctrl})h+u_rd1(k) (38)

u_rq1(k+1)＝(i_qref-i_{q_ctrl})h+u_rq1(k) (39)

重复第5步到第12步，直至完成K_sum个仿真时步的计算。

有益效果：

本发明针对电力系统动态仿真数值积分算法的缺点：难以同时兼顾仿真效率和数值稳定性，难适用于刚性有源配电网的动态仿真，仿真效率低、耗时长的问题，提出了一种基于多级优化追赶变分迭代法的有源配电网动态仿真方法。该方法将传统变分迭代法校正泛函中的已知量高效利用，优化了校正泛函的求解次序，将求出的近似解析解分段地应用在整个仿真时域上。在高分布式电源渗透率的有源配电网这种复杂刚性非线性系统这种应用环境下：一方面，在离线的近似解析解求解过程中，其获得同样精度近似解析解，所需的迭代次数比传统变分迭代法更少，计算量更小；另一方面，在在线仿真计算阶段，其在保证数值精度和数值稳定性的前提下，计算效率比数值积分法更快，仿真耗时更短。

附图说明

图1是基于多级优化追赶变分迭代法的有源配电网动态仿真计算流程图。

图2是金寨实际系统算例中节点7的母线电压幅值图。

具体实施方式

下面结合附图对发明的技术方案进行详细说明。本发明选择光伏电站的六个电气量波形作为聚类指标可以清晰且准确地描述光伏电站的动态特性，为了捕捉不同电站聚类指标之间的动态相似度，采用动态时间弯曲距离对基于欧式距离的聚类算法进行改进，从而提出了一种新型的动态聚类算法。

本发明公开的一种基于多级优化追赶变分迭代法的有源配电网动态仿真计算方法如图1所示，包括以下步骤：

步骤10)初始化阶段；

步骤20)动态过程计算阶段；

作为微分方程状态变量的初值计算阶段，步骤10)初始化过程如下：

步骤101)建立有源配电网仿真计算所需的元件模型，并输入有源配电网的系统参数。

包括光伏发电系统模型，同步发电机模型，电网络模型等等；形成整个系统的微分方程组和代数方程组，构成全系统动态仿真数学模型，由式(1)所示的微分-代数方程组进行描述：

其中，x和y分别是代数变量和状态变量，f和h表示函数，分别为n维和m维。

输入的参数包括母线和线路的参数，光伏电站和同步发电机的参数以及所在母线的位置，系统频率f_grid，以及系统的基准电压和基准容量等。

步骤102)潮流计算。

进行潮流计算。获得整个系统的稳态值矩阵，包括母线电压v_xy＝v_x+j·v_y，负荷的节点功率s_l＝p_l+j·q_l以及各个分布式电源的节点注入功率s_g＝p_g+j·q_g等等。随后基于潮流结果，按照下面的公式计算负荷等值并联导纳y_load，并形成计及负荷的系统导纳矩阵Y。其中v是负荷节点电压的幅值。

步骤103)计算所有微分方程组状态变量的初值。

根据发电机和光伏系统的物理特性计算所有微分方程组状态变量的初值。以同步发电机为例进行说明，光伏系统的计算与之类似。有源配电网受到扰动之前，根据潮流计算可以知道各个同步发电机的节点注入功率和节点电压，于是可以计算同步发电机的注入电流：

所以同步发电机的转子角初值δ₀计算如下。

δ₀＝angle{v+(r_a+j·x_q)·i_xy} (4)

此外，所有同步发电机的转速初值ω₀设置为1。

步骤104)输入仿真时序数据。

包括仿真总时间t_sum，仿真步长h，第i个扰动事件的开始时刻t_start和结束时刻t_end，其中i＝1,2,3,…。

作为动态仿真计算状态变量轨迹值的动态过程计算阶段，动态过程计算的本质是计算仿真时间内每一个时刻的有源配电网中微分方程组的状态变量值，除去初始时刻，共需要进行K_sum＝t_sum/h个仿真时步的计算。对于某一具体的仿真时步k＝1,2,…,K_sum。步骤20)计算过程如下：

步骤201)，计算k时刻的分布式电源的节点注入电流和节点电压。

首先需要求取同步发电机的虚拟注入电流i_xy和等值导纳矩阵Y_G。为方便后续节点电压的计算，Y_G要被并入Y形成新的系统导纳矩阵Y’，从而可以将虚拟电流i_xy看作新的节点注入电流，计算如下。

其中，E′_q是暂态电势，导纳b_x,g_y的计算如下。

E′_q＝v_q+r_a·i_q+x′_d·i_d (7)

x′_d是次同步电抗，x_q是q轴同步电抗，r_a是定子绕组电阻，i_d和i_q为和节点注入电流i_xy经过d-q坐标变换后的d轴和q轴分量。v_q是稳态时节点电压v的q轴分量。

Y_G的计算如下。

对于光伏系统而言，直接将上一时刻的i_d和i_q作为光伏接入节点的注入电流i_xy，这是因为i_d和i_q均为光伏系统微分方程组的状态变量，而状态变量是不会突变的。

如系统遇到扰动事件，就在本步骤中更新整个系统的导纳矩阵Y’即可。最后，综合各个同步发电机的注入电流i_xy以及计及Y_G和扰动事件的系统导纳矩阵Y’，使用Gauss迭代计算出k时刻的节点电压v_xy。

步骤202)，计算k时刻微分方程组与代数方程组之间的接口量。

该步骤的目的是将微分方程组中所有的接口量计算出来，这样就可以将这些量作为已知的常数，以便于在后面的步骤来专门来处理微分方程组的计算。根据上一步计算的k时刻的节点电压和电流，可以计算同步发电机的接口量：

p_m＝p_g+(i_x ²+i_y ²)·r_a (9)

类似的可以求出光伏系统的接口量包括，逆变器直流侧电流i_dc，PV阵列输出电流i_pv，逆变器调制电压V_dq，逆变器交流侧电压U_idq。

步骤203)，使用MOCVIM计算有源配电网微分方程的近似解析解。

步骤2031)，将微分方程组中的高阶微分方程转化为等值的1阶微分方程组。

考虑式(1)中如下的k阶微分系统：

y^(k)+f(y^k-1,y^k-2,...,y′,y)＝g(t) (11)

在绝大多数应用环境下，基于变分迭代法求取该式的拉格朗日乘子较为困难，为此，令y₁(t)＝y将该式转化为等值的k维1阶微分方程组如下所示。

对于该系统使用MOCVIM可较为容易地求出拉格朗日乘子。

对于有源配电网而言，可以推导出光伏系统和同步发电机系统的微分方程如下。其中，光伏系统选取应用广泛的双级式光伏电站的数学模型，其动态元件为9阶系统，将所有微分方程组整理为一阶微分方程的形式：

di_d/dt＝(ω_s/l_f)(-u_id+u_{sd_system})+ω_s·i_q (13)

di_q/dt＝(ω_s/l_f)(-u_iq+u_{sq_system})-ω_s·i_d (14)

di_dref1/dt＝u_dcref-u_dc (15)

du_rd1/dt＝i_dref-i_{d_ctrl} (16)

du_rq1/dt＝i_qref-i_{q_ctrl} (17)

du_pv/dt＝(i_pv-i_L)/C_pv (18)

di_L/dt＝[u_pv-(1-D)u_dc]/L_dc (19)

du_dc/dt＝[(1-D)i_L-i_dc]/C_dc (20)

dD₁/dt＝u_pv-u_m (21)

而同步发电机选取应用广泛的2阶同步发电机数学模型进行表述，其中微分方程组如下所示。

dσ/dt＝ω_s(ω-1) (22)

dω/dt＝[p_m0-p_e-D_g(ω-1)]/t_j (23)

其中i_d,i_q分别为光伏系统的实际输出电流的d轴和q轴分量，ω_s是电网角频率，l_f为滤波器电感，u_id和u_iq分别为逆变器输出电压的d轴和q轴分量，u_{sd_system},u_{sq_system}分别为并网电压的d轴和q轴分量，i_dref1是，u_dc和u_dcref分别是直流侧电压及其参考值，u_rd1,u_rq1分别为调制波的d轴和q轴分量，L_dc和C_dc分别为直流斩波器的电感和电容，D为开关的占空比，i_dc和u_dc分别为直流斩波器的输出电流和输出电压。σ和ω分别为转子角和转子角速度，p_m0和p_e分别是机械功率和电磁功率，D_g是阻尼系数，t_j是惯性常数。

步骤2032)，给出MOCVIM的校正泛函格式。

以上个步骤的微分方程组式中最后一个微分方程为例，可得该微分方程求解第n+1步解析解的校正泛函如下：

与原始变分迭代法不同，该校正泛函将第n+1步已经计算出的所有其他状态变量的解析解作为已知量参与y_k(n+1)(t)的计算，该技术可以加速迭代解向真实解的收敛速度。

步骤2033)，优化方程组的求解顺序。

应用第二步对微分方程组求解的重要一点就是安排方程组的求解顺序，优良的求解顺序可以进一步加速近似解析解的收敛速度。改进的校正泛函求解顺序如下所示。

这种安排顺序将会充分提高第n+1步的已知量的使用频率。上式中除了第一个校正泛函外，后面的每一个校正泛函的右侧均计及了n+1步的已知量。因此，在实际仿真应用中，仅需较少次数的迭代计算即可获得与传统变分迭代法更高阶时才具备的同等级精度的近似解析解。

步骤2034)，求取y_k对应的拉格朗日乘子λ_k(s)。

在式(13)-(23)中，只有式(16)和(17)不与其他状态变量相互耦合，即可以直接求出解析解，而其他式子均需要进行MOCVIM推导。考虑到在有源配电网仿真运算中，在微分方程组中把代数方程组传过来的变量作为已知的常数，因此可将用下面的A₁-A₁₀来简化上述微分方程组：

A₁＝ω_s(-u_id+u_{sd_system})/l_f,A₂＝ω_s(-u_iq+u_{sq_system})/l_f,

A₃＝i_pv/C_pv,A₄＝-1/C_pv,A₅＝1/L_dc,A₆＝(D-1)/L_dc,

A₇＝(1-D)/C_dc,A₈＝-i_dc/C_dc,A₉＝(p_m0-p_e+D_g)/t_j,A₁₀＝D_g/t_j。

根据步骤2032可列出上述微分方程组的校正泛函，这里以i_d的校正泛函为例如下式所示，其拉格朗日乘子为λ₁。其余微分方程组的拉格朗日乘子为λ₂-λ₉形式：

为求取λ₁对上述方程两边求变分可得：

其中上方有横线的变量表示该变量为限制变分，对限制变分进行求变分等于零。

若随着迭代的进行，i_d(n+1)(t)的变分应该等于0，因此令δi_d(n+1)(t)＝0可以有限制条件：

从而可得λ₁＝-1。类似得也可以计算λ₂到λ₈均等于-1，而ω对应的

步骤2035)，列出光伏电站和同步发电机中微分方程组相应的校正泛函，并求出各阶近似解析解。

根据第2031步到第2035步，即可有MOCVIM校正泛函及其求导次序如下所示。

利用该校正泛函，设状态变量对应的初值为I_d,I_q,I_dref1,U_pv,I_L,U_dc,D₁₍₀₎,σ₀,ω₀，令初值为其初始的迭代函数，代入上述校正泛函中，不断迭代即可得到不同阶的近似解析解。

步骤204)，基于计算获得的各微分方程组的近似解析解获取相应状态变量在第k个时段的动态轨迹。

获取相应状态变量在第k个时段的动态轨迹，也即可以获得k+1时刻的状态变量值y(k+1)。此外，状态变量u_rd1和u_rq1可以直接求出解析解，实际应用时使用下面的公式进行仿真计算：

u_rd1(k+1)＝(i_dref-i_{d_ctrl})h+u_rd1(k) (38)

u_rq1(k+1)＝(i_qref-i_{q_ctrl})h+u_rq1(k) (39)

步骤205)重复步骤201)到步骤204)，直至完成K_sum个仿真时步的计算。

下面列举一实施例以辅助说明本发明公开的动态方法的技术优势。以金寨县配网系统来验证所提算法在大规模实际系统的连续事件算例中的效率。仿真总时间设置为5s。在该系统中进行下列连续事件的仿真：首先，整体光伏电站的辐照度在1.5s时从1000W/m²上升到1300W/m²然后在2s时恢复到1000W/m²；随后，节点7处的负荷在2.5s时切除，并在3s重新接入；最后，在4s设置线路6-7靠近节点7处发生三相短路故障，并在4.1s清除。在该算例中，分别进行了步长为0.1ms、0.2ms以及0.4ms的算例，并选用迭代至第3阶和第5阶的MOCVIM的近似解析解(简写为MOCVIM-3、MOCVIM-5)与传统变分迭代法的第5阶(MVIM-5)以及隐式梯形积分法(Trapz)进行对比仿真。选取母线7的电压幅值的仿真结果展示如图2所示。根据该图可知，本章所提出的3阶MOCVIM的追踪性能已经与5阶的MVIM接近，而5阶MOCVIM的追踪性能相比前两者有大幅的提升。对于该地区的配网，应用不同的步长对连续事件算例进行仿真，获得各个方法的仿真时长如表1所示。

表1各算法在金寨实际配网中的仿真时间对比

根据该表结果，MOCVIM-3与MVIM-5在不同仿真步长下的仿真时间比较接近，而MOCVIM-5由于模型复杂度的上升，则相比于MVIM-5的仿真时间略有增加。但相比于Trapz，MOCVIM-5的仿真时间仍有明显缩短，其在步长0.1ms、0.2ms以及0.4ms下分别缩短了32.73％、22.37％和18.46％。此外，获得MVIM-5、MOCVIM-3以及MOCVIM-5三种模型的离线计算时间分别为18.35s、11.43s和31.76s。

Claims (1)

1.一种基于多级优化追赶变分迭代法的有源配电网动态仿真计算方法，其特征在于，包括以下步骤：

第1步，建立有源配电网仿真计算所需的元件模型，并输入有源配电网的系统参数，元件模型包括光伏发电系统模型，同步发电机模型，电网络模型；形成整个系统的微分方程组和代数方程组，构成全系统动态仿真数学模型，由式(1)所示的微分-代数方程组进行描述：

其中，x和y分别是代数变量和状态变量，f和h表示函数，分别为n维和m维；

输入的参数包括母线和线路的参数，光伏电站和同步发电机的参数以及所在母线的位置，系统频率f_grid，以及系统的基准电压和基准容量；

第2步，潮流计算，获得整个系统的稳态值矩阵，包括母线电压v_xy＝v_x+j·v_y，负荷的节点功率s_l＝p_l+j·q_l以及各个分布式电源的节点注入功率s_g＝p_g+j·q_g，随后基于潮流结果，按照下面的公式计算负荷等值并联导纳y_load，并形成计及负荷的系统导纳矩阵Y，其中v是负荷节点电压的幅值

第3步，计算所有微分方程组状态变量的初值：根据发电机和光伏系统的物理特性计算所有微分方程组状态变量的初值；

第4步，输入仿真时序数据，包括仿真总时间t_sum，仿真步长h，第i个扰动事件的开始时刻t_start和结束时刻t_end，其中i＝1,2,3,…；

第5步，计算k时刻的分布式电源的节点注入电流和节点电压，

首先求取同步发电机的虚拟注入电流i_xy和等值导纳矩阵Y_G，Y_G并入Y形成新的系统导纳矩阵Y’，将虚拟电流i_xy看作新的节点注入电流，计算如下

其中，E′_q是暂态电势，导纳b_x,g_y的计算如下

E′_q＝v_q+r_a·i_q+x′_d·i_d 式(7)

x′_d是次同步电抗，x_q是q轴同步电抗，r_a是定子绕组电阻，i_d和i_q为和节点注入电流i_xy经过d-q坐标变换后的d轴和q轴分量，v_q是稳态时节点电压v的q轴分量

Y_G的计算如下

综合各个同步发电机的注入电流i_xy以及计及Y_G和扰动事件的系统导纳矩阵Y’，使用Gauss迭代计算出k时刻的节点电压v_xy；

第6步，计算k时刻微分方程组与代数方程组之间的接口量，根据上一步计算的k时刻的节点电压和电流，计算同步发电机的接口量：

p_m＝p_g+(i_x ²+i_y ²)·r_a 式(9)

求出光伏系统的接口量包括，逆变器直流侧电流i_dc，PV阵列输出电流i_pv，逆变器调制电压V_dq，逆变器交流侧电压U_idq；

第7步，将微分方程组中的高阶微分方程转化为等值的1阶微分方程组，

考虑式(1)中如下的k阶微分系统：

y^(k)+f(y^k-1,y^k-2,...,y′,y)＝g(t) 式(11)

令y₁(t)＝y将该式转化为等值的k维1阶微分方程组如下所示

对于有源配电网，推导出光伏系统和同步发电机系统的微分方程如下，其中，光伏系统选取双级式光伏电站的数学模型，其动态元件为9阶系统，将所有微分方程组整理为一阶微分方程的形式：

di_d/dt＝(ω_s/l_f)(-u_id+u_{sd_system})+ω_s·i_q 式(13)

di_q/dt＝(ω_s/l_f)(-u_iq+u_{sq_system})-ω_s·i_d 式(14)

di_dref1/dt＝u_dcref-u_dc 式(15)

du_rd1/dt＝i_dref-i_{d_ctrl} 式(16)

du_rq1/dt＝i_qref-i_{q_ctrl} 式(17)

du_pv/dt＝(i_pv-i_L)/C_pv 式(18)

di_L/dt＝[u_pv-(1-D)u_dc]/L_dc 式(19)

du_dc/dt＝[(1-D)i_L-i_dc]/C_dc 式(20)

dD₁/dt＝u_pv-u_m 式(21)

同步发电机选取2阶同步发电机数学模型进行表述，其中微分方程组如下所示

dσ/dt＝ω_s(ω-1) 式(22)

dω/dt＝[p_m0-p_e-D_g(ω-1)]/t_j 式(23)

其中i_d,i_q分别为光伏系统的实际输出电流的d轴和q轴分量，ω_s是电网角频率，l_f为滤波器电感，u_id和u_iq分别为逆变器输出电压的d轴和q轴分量，u_{sd_system},u_{sq_system}分别为并网电压的d轴和q轴分量，u_dc和u_dcref分别是直流侧电压及其参考值，u_rd1,u_rq1分别为调制波的d轴和q轴分量，L_dc和C_dc分别为直流斩波器的电感和电容，D为开关的占空比，i_dc和u_dc分别为直流斩波器的输出电流和输出电压；σ和ω分别为转子角和转子角速度，p_m0和p_e分别是机械功率和电磁功率，D_g是阻尼系数，t_j是惯性常数；

第8步，给出MOCVIM的校正泛函格式；

根据式(23)得该微分方程求解第n+1步解析解的校正泛函如下：

与原始变分迭代法不同，该校正泛函将第n+1步已经计算出的所有其他状态变量的解析解作为已知量参与y_k(n+1)(t)的计算；

第9步，优化方程组的求解顺序，校正泛函求解顺序如下所示

…

第10步，求取y_k对应的拉格朗日乘子λ_k(s)，将用下面的A₁-A₁₀来简化上述微分方程组：

A₁＝ω_s(-u_id+u_{sd_system})/l_f,A₂＝ω_s(-u_iq+u_{sq_system})/l_f,

A₃＝i_pv/C_pv,A₄＝-1/C_pv,A₅＝1/L_dc,A₆＝(D-1)/L_dc,

A₇＝(1-D)/C_dc,A₈＝-i_dc/C_dc,A₉＝(p_m0-p_e+D_g)/t_j,A₁₀＝D_g/t_j；

根据第8步列出上述微分方程组的校正泛函，i_d的校正泛函如下式所示，其拉格朗日乘子为λ₁，其余微分方程组的拉格朗日乘子为λ₂-λ₉形式：

为求取λ₁对上述方程两边求变分可得：

其中上方有横线的变量表示该变量为限制变分，对限制变分进行求变分等于零，

随着迭代的进行，i_d(n+1)(t)的变分应该等于0，令δi_d(n+1)(t)＝0可以有限制条件：

从而可得λ₁＝-1；类似得也可以计算λ₂到λ₈均等于-1，而ω对应的

第11步，列出光伏电站和同步发电机中微分方程组相应的校正泛函，并求出各阶近似解析解，根据第8步到第11步，即可有MOCVIM校正泛函及其求导次序如下所示

利用该校正泛函，设状态变量对应的初值为I_d,I_q,I_dref1,U_pv,I_L,U_dc,D₁₍₀₎,σ₀,ω₀，令初值为其初始的迭代函数，代入上述校正泛函中，不断迭代即可得到不同阶的近似解析解；

第12步，基于计算获得的各微分方程组的近似解析解获取相应状态变量在第k个时段的动态轨迹，

获取相应状态变量在第k个时段的动态轨迹，也即可以获得k+1时刻的状态变量值y(k+1)，此外，状态变量u_rd1和u_rq1可以直接求出解析解，实际应用时使用下面的公式进行仿真计算：

u_rd1(k+1)＝(i_dref-i_{d_ctrl})h+u_rd1(k) 式(38)

u_rq1(k+1)＝(i_qref-i_{q_ctrl})h+u_rq1(k) 式(39)；

第13步，重复第5步到第12步，直至完成Ksum个仿真时步的计算。

Images