含电动汽车集群调频的信息能源系统谱特征和稳定性评估

张一媚，董朝宇，董晓红，肖 迁，张明东，贾宏杰

（1. 智能电网教育部重点实验室（天津大学）,天津市300072；2. 河北工业大学人工智能与数据科学学院,天津市300130；3. 天津大学智能与计算学部,天津市300072）

0 引言

随着环境污染加剧和化石资源逐渐枯竭以及间歇性可再生能源的大规模利用,经聚合后形成的电动汽车（electric vehicle,EV）集群,凭借快速的功率输出与响应特性,成为一种很好的调控资源,可用于平滑可再生能源出力、提供辅助服务等领域［1-5］。利用EV 集群提供调频服务已成为近期的一个研究热点。

EV 集群参与调频时,首先需收集和发送EV 状态到控制中心,然后根据控制信号,对集群内EV 充放电行为实施控制［6］,这一过程涉及对大量EV 的信息收集和协调控制。毫无疑问,此时的电力系统成为一个典型的信息能源系统,即信息与能源技术深度融合并集成能源网、通信和计算的复杂系统［7-11］,评估其稳定性需要同时考虑能量流和信息流的影响。在对大量分散EV 进行调控时,出于经济性考虑,应用异构通信手段甚至直接借助公共通信媒介成为必然选择,由此导致信息传输环节的数据拥塞、丢包、延时等情况必然时有发生。在数学模型中可用时滞对其统一加以刻画［12］。而时滞往往会影响系统动态性能甚至引发失稳［13-16］。频率小扰动稳定过程既可能是秒级的短期现象［17］,也可能是从几十秒到几分钟的长期现象［18-20］。研究EV 集群参与系统调频时,需借助时滞稳定性分析方法来科学评估时滞的影响［21-26］。

时滞稳定性分析方法包含时域法和频域法［16］：时域法通过对Lyapunov 泛函进行放缩得到充分不必要条件,尽管可以通过等效变换［21］、新泛函构造［22-23］及优化判据推导［24］等手段对其加以改进,但依然无法消除保守性；而频域法的基础是频域变换技术和特征值理论,尽管没有保守性,但时滞引入难以处理的超越项［16］。在文献［26］中,利用指数延迟项的Pade 有理多项式逼近,避免了直接求解时滞特征方程的困难,但难以应用于时滞较大或规模较大的系统中。文献［15］采用解算子离散化方案求解系统关键特征值,但是存在时滞研究范围较小的局限性。

为分析时滞对EV 聚合调控效果的影响,本文将基于频域法,推导一种信息能源系统谱特征和稳定性评估方法。首先,建立EV 集群调频的系统模型。进一步,将含物理-信息环节的系统谱转换成无穷小算子的谱,消除超越项；然后,基于用户准确度与速度需求选择离散化方案,利用切比雪夫离散化生成有限维近似矩阵,进而利用牛顿迭代校正近似矩阵的特征值,得到系统准确特征值；最后,刻画系统稳定域,可多角度评估时滞对系统稳定性的影响。提出的方法解决了时域法存在保守性的问题,也消除了频域法中难以处理的超越项,通过将谱变换、谱离散、牛顿迭代等关键技术与用户需求相结合,可以实现高效准确评估规模或时滞较大系统。

1 系统模型

1.1 EV 集群动态特性

考虑大量EV 集群行为时,一般可用EV 集群来表示。假设系统中存在k个EV 集群,即

式中：GA为EV 集群构成的集合；GAi为第i个EV 集群；i=1,2,…,k。

对于任意用于频率调节研究的EV 集群GAi,其中的EV 模型由电池模型和电池充电器模型（包括DC/AC 逆变器）组成。电池充电器控制交流电气系统和电池之间的电力交换,电池储存驱动EV 所需的能量。本文采用基于戴维南的模型作为电池的等效电路模型［27-28］。

为简化分析,忽略了EV 集群中EV 的个体差异。因此,在分析时,假设集群内EV 的荷电状态（SOC）为常数,认为集群具有等效集群电池系统增益KEVA,i,且集群内EV 的时间常数相同,表示为TEVA,i。可用如下一阶惯性环节表示EV 集群的动态特性［29-32］：

式中：GEVA,i和KEVA,i分别为集群GAi的动态特性函数和增益；TEVA,i=Lc/Rc为EV 集群电池系统时间常数,即EV 充电器的时间常数,其中Lc为电池充电器中与电网相连的滤波电抗器,其阻抗通常设计为系统总阻抗的5%,电阻Rc为电抗器电阻及逆变器损耗,取决于电缆尺寸、安装媒介及导体材料［4,33］。

KEVA,i由频率调节效果和电池SOC 约束之间的相互限制来决定,以表示EV 输出功率对频率偏差的敏感性,其数值主要取决于集群内EV 数量、EV电池的开路电压、单个EV 采用的控制策略及控制增益最大值［3,28,34］。而控制增益最大值可由EV 入网（V2G）最大功率（或电流）对系统允许频率偏差的敏感度计算［35］。综合考虑上述因素,可以对EV调频等值模型中的参数KEVA,i和TEVA,i进行确定。

需要指出的是,式（2）模型主要用于模拟EV 集群在参与调频时的总体动态特征,它忽略了EV 个体之间的细微差异,在EV 数量较大时,这一处理误差较小［36］。

1.2 考虑时滞时的调频模型

多个EV 集群参与系统调频的实现原理如图1所示。当EV 接入电网时,可通过充放电装置并利用V2G 技术,实现EV 与电网间能量和信息的双向流动。因EV 单体功率较小,大量分散接入电网后,需引入EV 集成商来负责EV 的充放电管理,以及与上级调度中心的交互。图1 中的每一个集成商负责监视、收集、分析和计算EV 集群的调控能力,然后上传至调度中心,同时负责执行调度中心下达的控制指令。

图1 EV 集群参与调频的信息能源系统控制框架Fig.1 Framework of frequency regulation for cyber-energy system with EV cluster

信息在EV 和调度中心之间的双向传输过程中必然存在一定时滞。为简单起见,假设每个EV 集群内所有的EV 拥有相同数据传输时滞,则由时滞构成的集合如下：

式中：τi为集群GAi的时滞。

本文在研究EV 集群参与系统调频时,将采用文献［6］的多区域电力系统模型。任何一个区域w的情况如图2 所示（假设系统中共包含g个如图2 所示区域）,它包含一个等效发电单元和多个EV 集群。发电单元包含调速器和再热式涡轮机,EV 集群考虑了信息传输时滞,它们通过参与因子α0,α1,…,αk实现功率分担。不同区域间通过区域控制误差（area control error,ACE）实现有机协调［6］。

图2 包括多个EV 集群的电力系统调频模型Fig.2 Frequency regulation model of power system with multiple EV clusters

单个区域的状态空间方程为［37-38］：

式中：xw,u,yw分别为系统的状态、控制、输出变量；A,B0,Bi,F,C分别为状态矩阵、非时滞控制矩阵、时滞控制矩阵、扰动矩阵、输出矩阵；ΔPL为负荷偏差。式（4）的详细推导及各矩阵组成见附录A。

经简单变换,式（4）可表述为时滞微分方程：

式 中：φw( =xw0)为 区 域w的 初 始 状 态；τw,max=max {τi},其中i=1,2,…,k；K=［KP,KI］,其中KP和KI分别为区域w中比例-积分（PI）控制器的比例增益和积分增益。

需要指出的是,当系统存在g（g＞1）个区域时,每一个区域内的系统方程均可通过上述过程推导。最后,再将各区域方程进一步整合,得到整个系统的动态方程［6］,即

其中

式中：x为系统的状态；φ为系统的初始状态；x0为系统初始状态。

进而可以得到系统的特征方程为：

式中：λ为特征值；υ为对应特征向量。

2 基于谱特征提取的稳定域求解方法

由式（8）和式（16）可知,系统的渐进稳定性将由φ、系数矩阵A͂和A͂i（i=1,2,…,m）以及信息环节时滞τ1,τ2,…,τm唯一决定。理论上,上述参数给定后,时滞系统稳定运行区域（稳定域）可唯一确定,而稳定域可用于多角度评估系统稳定性。本章将推导一种稳定域求解方法,其中间过程的一些结果还有助于分析系统关键特征值信息。其主要步骤如下。

1）谱变换：利用无穷小算子,将反映信息能源系统模型动态的时滞微分方程映射为常微分方程,从而消除难以处理的超越项。

2）切比雪夫离散化：对谱变换得到的无穷小算子进行切比雪夫离散化,得到近似离散化矩阵,将无穷维特征值问题转化为有限维特征值问题。

3）牛顿迭代校正：计算离散化矩阵的特征值及特征向量,将其作为牛顿法初值进行迭代校正,得到系统准确特征值。

4）稳定域求解：基于用户需求,快速推导出信息能源系统的临界特征值,进而通过搜索法求解系统稳定域,多角度评估系统稳定性。

2.1 谱变换

对于式（8）系统,令X：=C([-τmax,0],Cn)为由将区间[-τmax,0]映射到n维复数向量空间Cn的连续函数构成的巴拿赫空间。因为式（8）在［0,∞）内连续,假设初始状态φ存在于巴拿赫空间上,系统x（t）轨迹将在时域内连续且唯一,因此可将式（8）重新表述为巴拿赫空间X上的抽象常微分方程［39］：

式中：H：D(Hφ)⊆X→X是由式（18）定义的无穷小算子。

利用上式将t处系统状态映射到t+θ处,根据谱映射原理,原系统特征值将等于H的特征值［40］。

2.2 切比雪夫离散化

谱变换后得到的式（17）仍为无穷维系统,为简化求解,利用切比雪夫多项式插值近似初始函数,从而构造离散矩阵来逼近H。

首先,给定正整数N,考虑由N+1 个在区间[-τmax,0]内的离散点组成的网格ΩN：

离散化矩阵HN的详细推导及表达形式如附录B 所示。

经上述处理后,无限维的H被简化成有限维的HN。进而可直接求解系统的一组临界特征值,而不必求解其全部（无穷多）特征值。

2.3 牛顿迭代

在计算得到系统临界特征值的近似结果λ̂后（即HN最右侧特征值）,将其和对应Krylov 向量的前n个元素v̂,一起作为式（16）的初值{λ̂,v̂},然后利用牛顿法迭代校正,以得到准确值。

首先,将特征方程式（16）写成如下格式：

式中：In为n阶单位矩阵。

在迭代过程中v需要有一些标准化处理,以保证解的唯一性,迭代过程基于式（23）或式（24）。

式中：λp和vp分别为第p次迭代得到的特征值及特征向量。

式（24）可写为如下形式：

式中：fp为第p次迭代的函数值f；Δz包括Δv的所有元素,其中Δv中的单位元素被替换为Δλ；J为雅可比矩阵。

将式（25）转换为迭代公式,求得第p+1 次的迭代结果为：

经迭代达到收敛精度后,可得到系统准确特征值和特征向量{λ,v}。

2.4 系统稳定域求解

根据牛顿迭代得到的准确特征值,可以评估系统稳定性,进而可根据用户需求刻画稳定域。

通常用户需求体现在2 个方面：①稳定域准确度；②计算时间。基于此可以对离散化参数N进行选择,以寻求计算准确度与速度指标的平衡。

利用选择的特征值计算方法对某一参数下的系统特征值进行求解,进而可以使用成熟的特征值分析方法,对系统稳定性进行分析。

基于此,设置一定的搜索步长,沿时滞增大方向前进,在时滞空间中搜索系统的负实部临界特征值,直至到达研究区域边界。可以在不同的控制器参数设置下,对系统的稳定域进行刻画。

本文提出的谱特征和稳定性评估方法流程如附录C 图C1 所示：包括系统初始化和参数输入、建模、谱变换、切比雪夫离散化、牛顿校正、输出特征值。而后续进行稳定域刻画则是在时滞空间中对负实部临界特征值进行搜索,搜索过程对上述算法进行重复以求解各搜索点特征值。

3 算例分析

以含单个EV 集群调频的信息能源系统为例,具体参数列在文献［25］中。基于提出的谱特征提取法对引入通信环节时滞后系统的稳定性进行综合评估,并将其与以不考虑通信环节时滞的系统为初始状态的牛顿迭代法及线性矩阵不等式（LMI）方法进行对比算例分析。

为了进一步验证所提出方法,研究EV 时滞与电力系统调度网络时滞的区别,在3.3 节设置考虑发电单元引入时滞的信息能源系统,对其稳定域进行刻画。

3.1 单时滞含EV 调频的信息能源系统算例

对含单EV 集群的调频信息能源系统进行分析,其中控制器增益KP设置为0.2,KI设置为0.2,时滞τ设置为100 s,其余系统参数如文献［25］所示。对此典型算例,利用本文提出方法对系统稳定性进行分析,实现过程如下。

步骤1：系统初始化和参数输入。以设置N=30 为例,即离散化矩阵维数（N+1）n=186；设置牛顿校正收敛精度为10-6,最大迭代次数为100。

步骤2：建模。基于1.2 节方法对系统进行建模,推导得到系统非时滞参数矩阵和时滞参数矩阵。

步骤3：谱变换。为了避免直接求解含超越项的时滞微分方程,利用无穷小算子将其转化为抽象常微分方程。这样,时滞系统的谱转化为无穷小算子H的谱,进而可以通过计算H的特征值（σ(H)）来获得原系统的特征值。

步骤4：切比雪夫离散化。基于切比雪夫多项式对H进行离散,得到有限维离散矩阵HN。进而可以计算HN的特征值（σ(HN)）,对H的特征值进行近似求解,如附录D 图D1（a）所示。

步骤5：牛顿迭代。由于离散化步骤是近似逼近,因此需要进行校正。以HN的n（n=6）维临界特征 值 及 特 征 向 量{λ̂,v̂}为 初 值,应 用 牛 顿 法 求 解式（16）。第p次迭代的J和函数值fp可以由式（24）计算。

步骤6：输出特征值并进行稳定性分析。得到准确特征值后,可以有效使用特征值分析方法分析当前参数下系统稳定性。在此参数下,临界特征值实部为0.002 1＞0,意味着系统为负阻尼状态,会发生增幅振荡失稳。

步骤7：稳定域求解。基于上述分析过程,可以对研究范围内其他搜索点进行特征值计算与稳定性分析。设置时滞搜索步长为0.1 s,控制器比例增益步长为0.002,以刻画在不同的控制器参数下系统的稳定域,如图3 中红色区域所示。

3.2 对比算例分析

3.2.1 与以不考虑通信环节时滞的系统为初始状态的牛顿迭代法比较

为了验证所提方法的有效性,同时定量分析谱离散化对结果的影响及必要性,设置以不考虑通信环节时滞的系统为初始状态的牛顿迭代法进行对比分析。

同时,为对离散化维数设置对所提方法影响进行分析,设置不同离散化维数进行对比。

基于上述方法,分别对系统进行稳定性评估,得到了指定系统参数下能获取准确特征值允许的时滞研究范围（在本文研究范围内,采用所提出方法N=200 后再增大离散化维数,特征值结果不变,如附录E 表E1 所示,因此将该方法得到的特征值,作为验证方法精确性的基准）。

图3 本文方法与文献[25]方法结果比较Fig.3 Comparison between the proposed method and that in reference [25]

首先采用以不考虑通信环节时滞的系统为初始状态的牛顿迭代法,验证直接进行牛顿迭代的结果,对无时滞的系统(KP=0.2,KI=0.2)求取特征值,将其作为牛顿法初值,经过迭代得到最终值,结果显示：当时滞τ≤0.5 s 时,可以准确计算系统全部特征值；当0.5 s ＜τ≤36 s 时,不能准确计算全部特征值但仍可准确捕捉临界（最右侧）特征值；而当τ＞36 s时,全部特征值计算不准确。

随后基于谱特征方法,即对牛顿法初值进行预处理,在相同控制器参数设置下,改变N（N=10,N=20,N=30）,研究全部特征值准确、临界特征值准确、全部特征值不准确这3 种准确度标准对应的时滞范围,如附录E 图E1 所示。

总结上述结果,可以得出如下结论。

1）牛顿法实现准确校正的条件是迭代初值接近准确值：由谱离散化获得的初值影响牛顿校正结果,进而影响稳定性分析准确性。

2）时滞较小时,不考虑时滞的系统特征值与时滞系统近似,即牛顿迭代初值在成功校正范围内,经校正可以获得准确的全部特征值。

3）随时滞增大,计算特征值不再全部准确,只有临界特征值准确,时滞继续增大,迭代后全部特征值均存在误差。

4）所提出方法会显著增大对信息网络的时滞分析范围。

5）随离散化矩阵维数增大,特征值计算准确度提高,时滞分析范围增大,但计算时间增加,如附录E 表E2 所 示。

根据附录E 图E1 所示各方法的适用范围,针对不同的时滞研究范围和特征值提取要求,选择不同的离散化方案。

对于全部特征值准确计算的目标,由附录E 图E1 中蓝色折线可以看出：①当τ≤0.5 s 时,选择不考虑通信环节时滞系统为初始状态的牛顿迭代法；②当时滞范围扩大到0.5～30 s 时,采用N=10 的谱特征方法；③继续增大研究时滞范围,相应采用的N也变大。而当特征值提取要求降低时,如进行稳定性分析时,只需保证临界特征值计算准确,如附录E图E1 中红色折线所示。

基于此,可以针对用户需求（时滞研究范围和特征值提取要求（全部/临界特征值准确））对方法进行选择。

随后在不同离散化维数设置下,对系统稳定域进行求解。首先改变离散化矩阵维数,在N=10,N=20,N≥30 的 设 置 下,在N=10 和N=20 离 散化方案下刻画3.1 节算例中稳定域,如附录E 图E2和图E3 所示。

在N≥30 方案下,获得的稳定域与3.1 节一致,如图3 中红色区域所示。经验证该结果与仿真（MATLAB,Simulink）一致,如附录E 图E4 所示。

结合附录E 图E2、图E3 和图3 分析可知,在N较小时,随时滞增大,会出现稳定域刻画不准确的现象。但是在N≥30 参数设置下,均可以获得准确稳定域,且计算效率较高。

附录E 表E2 给出了不同离散化维数设置下的稳定域求解方法的迭代次数和计算时间等信息。通过对比分析可知：随离散化维数增大,谱离散所得特征值准确性提高,导致牛顿校正迭代次数变小。但由于矩阵维数变大,单次特征值计算时间增加,刻画稳定域效率下降,但时滞分析范围显著上升。

对附录E 图E2、图E3 和图3 呈现的规律进一步细化,以N=2 为步长,研究稳定域计算准确度（判断准确区域与总区域点数之比）与N的对应关系,结果显示：N增大到24 后,已经可以实现准确刻画稳定域,随后增大N,准确度不会再提升但计算效率下降。基于此,可依据用户需求,选择离散化方案,以寻求计算准确度与计算速度之间的平衡。

总结上述结果,可得出如下结论。

1）本文提出方法,与以不考虑通信环节时滞的系统为初始状态进行牛顿迭代相比,会显著增大对信息网络的时滞分析范围,实现在时滞较大时,仍能准确高效捕捉系统特征值。

2）采用谱特征提取法评估系统稳定性时,随离散化维数增大,时滞分析范围增大,稳定域准确度提高但速度下降。在N=30 内可以获取准确稳定域且实现较高效率。

3.2.2 与LMI 方法的比较

一般来说,利用李雅普诺夫定理引入保守性。一些研究不断对时域法进行优化试图减少判据保守性。参考文献［25］基于Lyaponuv 理论和LMI 方法,利用基于Writinger 的改进积分不等式提出了时滞稳定性准则。其判据被公认为保守性很小。为此,将本文方法与其对比,并基于2 种方法分别对3.1 节算例中的系统稳定域进行刻画。结果如图3所示。从图中可以看出,该算法仍存在不可避免的保守性,求得的稳定裕度小于系统实际稳定边界值,并且不能计算出第1 个稳定域之后随时滞增大再次出现的稳定域。而本文提出方法解决了传统LMI方法保守性的问题,可以准确刻画系统全部稳定域。

3.3 考虑EV 时滞与电力系统调度网络时滞的双时滞系统算例

为进一步验证本文所提方法,并研究EV 时滞与系统调度网络时滞的区别。本节构建考虑EV 时滞与系统调度网络时滞的信息能源系统,综合考虑发电单元与EV 引入时滞,考虑时滞后发电单元模型如图4 所示,整个系统动态及状态空间方程如附录F 所示。

图4 考虑时滞的发电单元模型Fig.4 Power generation unit model considering delays

采用本文提出方法,在KP=0.2 和KI=0.2 的参数设置下,首先固定发电单元时滞不变（τ0=5 s）,研究不同离散化维数下特征值准确度与EV 时滞的对应关系、结果如附录G 表G1 所示。随后固定EV时滞（τ1=5 s）,特征值准确度与发电单元时滞的关系如附录G 表G2 所示。

对比2 个表格,可以观察到固定τ0时该方法对τ1的时滞分析范围比固定τ1时对τ0的分析范围更大,附录G 图G1 解释了该结果：固定τ0时,随着τ1增大,临界特征值会逐渐收敛于固定值不再变化,因而在较小的离散化维数下可以获得更大的时滞分析范围。而固定τ1增大τ0,临界特征值不会收敛,因而时滞分析范围相对较小。

随后采用提出方法（N=30）,在τ0∈[0,100]s,τ1∈[0,100]s 的研究范围内,求解系统稳定域,如图5 所示。

图5 信息能源系统稳定域Fig.5 Stability region of cyber-energy system

1）EV 时滞的稳定域呈现3 个区间（τ0=0 s,τ1=30 s 处存在一个小区间）,而对于发电单元时滞,则存在一个确定的稳定裕度,当时滞小于裕度系统稳定,超过裕度系统失稳。

2）EV 时滞比发电单元时滞对系统稳定域影响大,如图5 所示。当τ0在[0,10]s 范围内取值时,系统存在稳定域（在τ0∈[0,10]s 范围画竖线,均有稳定域与直线相交）；当τ1∈[4,10]s 时,无论τ0取值为多少,系统均不稳定。

以N=2 为步长,研究稳定域计算准确度与N的关系,结果绘制于附录G 图G2。在N增大到14后,可以实现准确刻画稳定域。

进一步,对EV 与发电单元对系统稳定性产生不同影响的原因进行探究,结果如附录H 所示。

4 结语

本文提出了一种含EV 集群调频的信息能源系统谱特征和稳定性评估方法。利用谱变换、谱离散、牛顿迭代等关键技术,计算得到系统特征值,对信息能源系统进行稳定性综合评估。

1）含EV 集群调频的信息能源系统的稳定性分析结果表明,该方法能够较准确且高效地求解信息能源系统特征值。

2）基于谱特征提取方法分析稳定性时,准确求解特征值的时滞范围与离散化维数有关,当N在30以内,可以获得全部准确稳定域；再增加系统维数,计算效率下降。

3）基于用户需求,可以对稳定性评估方法进行选择,以寻求计算准确度与速度间的平衡。

4）综合考虑EV 时滞与电力系统调度网络时滞的算例结果揭示了EV 时滞与调度网络时滞对系统稳定域影响的区别：前者的稳定域呈现若干个区间,而后者则是存在一个确定的稳定裕度。同时,EV时滞对系统稳定域的影响更大。

本文提出方法存在当时滞很大时为了保证准确性使得离散化维数增大从而导致计算效率下降的局限性。今后的工作中将进一步对系统进行降维以实现更高的计算效率。

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