查看原文
其他

【精彩论文】基于全状态模型的自同步电压源并网系统频率稳定分析

中国电力 中国电力 2023-12-18

基于全状态模型的自同步电压源并网系统频率稳定分析


李振垚1, 甘德强1, 栾某德1, 何国庆2

(1. 浙江大学 电气工程学院, 浙江 杭州 310027; 2. 中国电力科学研究院有限公司, 北京 100192)


摘要:使用自同步电压源控制方式的新能源是提升高比例新能源电力系统频率稳定的重要措施。针对自同步电压源并网系统,基于近似全状态模型,通过状态变量时域响应的解析解得到系统频率相关特征量与系统参数的量化关系。首先,将频率稳定问题转化为二次特征值问题,证明了系统惯量减小会降低系统频率惯量响应阶段的抗扰动性。其次,将自同步电压源模型与同步机模型参数进行对应,证明了调速器的作用相当于增加系统阻尼,能够减小受扰动后系统频率的稳态频率误差。然后,提出了一种满足预设频率动态安全要求的自同步电压源新能源虚拟惯量和下垂系数估计方法。最后,在10机39节点算例中验证了所提方法的正确性和有效性。


引文信息

李振垚, 甘德强, 栾某德, 等. 基于全状态模型的自同步电压源并网系统频率稳定分析[J]. 中国电力, 2023, 56(5): 182-192.

LI Zhenyao, GAN Deqiang, LUAN Moude, et al. Frequency stability analysis based on full state model in autonomous-synchronization voltage source interfaced power system[J]. Electric Power, 2023, 56(5): 182-192.


引言


随着新能源大容量并入电网,电力系统惯量相对减小,电力系统频率动态行为愈发复杂,频率稳定问题突出,亟需有效的频率稳定分析方法和控制手段[1]。电力系统有功功率波动导致同步机组转速发生变化,进而导致系统频率发生偏差,因此准确刻画频率动态响应过程是保证电力系统安全稳定运行的基础[2]。新能源的低惯性和低抗扰性是引发英国“8·9”大停电事故的原因之一,单一故障导致了新能源机组脱网,功率不平衡进一步扩大,使得系统频率快速跌落,触发了低频切负荷系统的动作[3]。国内也存在类似风险,2015年9月19日,锦屏—苏南直流发生双极闭锁事故,华东电网频率最低跌至49.56 Hz,并长时间处于越限状态。目前缺少新能源系统参数对电力系统频率稳定性影响的系统化、量化评估方法,亟须开展适用于高比例新能源电力系统的频率稳定评估方法研究。采用虚拟同步控制策略的自同步电压源能够模拟同步机的运行方式,在电网频率波动时调整输出功率,一定程度上有助于缓解电网频率调节压力[4-5],是提升高比例新能源电力系统频率稳定的重要措施。因此,后文所提自同步电压源均指采用虚拟同步控制策略的自同步电压源。频率动态响应基本特征包括惯量响应阶段的频率变化率(rate of change of frequency,RoCoF)、一次调频阶段的最低频率点(nadir frequency, NF)和稳态频率误差[6]。惯性是系统频率抵抗外界扰动的重要表现,在频率动态响应过程中体现为频率变化率,大规模新能源接入电网,替代了同步机出力的同时,转动惯量相对减小,从而增大了惯量响应阶段的频率变化率[7]。目前电力系统频率稳定分析主要分为3种。1)时域仿真法[8],即通过数值积分方法求解描述系统的微分代数方程组,逐步积分计算系统运动轨迹,判断系统暂态稳定性。虽然时域仿真法计算精度高,但难以解释内在稳定机理。2)简化模型方法(low order system frequency response, SFR),即基于系统平均频率的单机带集中负荷模型[9]。近年来学者基于传统模型扩展出了更加具有适应性的SFR模型。文献[10]推导出新能源渗透下电力系统的简化模型,能够给出不参与系统调频的新能源机组的渗透率上限。文献[11]提出了一种聚合系统频率响应模型,在机组参数不同时提高了等效模型的精度。文献[12]建立了一种风电场参与频率控制的系统频率响应模型,很好地模拟了风电波动后的系统频率响应。3)基于人工智能方法。文献[13]提出了一种基于深度置信网络的在线预测电力系统频率在受扰后的响应曲线的方法。文献[14]提出了一种基于多层极限学习机的在线频率安全评估方法。文献[15]结合电力系统安全风险评估的一般流程,提出了基于机器学习的频率安全评估系统的总体框架和频率动态安全风险统一评估模型结构。SFR在训练数据大的条件下能够给出精度较高的结果,但此类方法对内在物理机理解释不清晰。同时,大量自同步电压源型新能源的接入,给系统频率稳定分析带来了新的技术挑战。1)电力电子器件数量多,模型阶数高,使得时域仿真建模工作量更加繁琐;2)新能源频率控制手段、响应速度异于同步发电机,传统简化模型方法受限于其灵活性,模型聚合后精确度下降严重,影响了分析结果的实用性。

本文基于新能源并网系统近似全状态模型,借助二次特征值(the quadratic eigenvalue problem,QEP)分析理论,通过状态变量时域响应的解析解得到系统频率相关特征量与系统参数的变化关系。多机电力系统可以用一组微分方程描述,在研究频率稳定问题时其状态变量时域响应解可以用线性化模型的状态矩阵特征值及其相应的左、右特征向量求得,从而可以从系统参数对特征值的影响得到系统参数与系统频率动态响应之间的关系。基于上述结论提出了一种满足预设频率动态安全要求的自同步电压源新能源虚拟惯量和下垂系数估计方法。最后在10机39节点算例中验证上述结论和方法的正确性和有效性。


1  考虑频率问题的发电机模型


1.1  同步机模型在考虑频率稳定问题时,同步机可采用暂态电抗后的内电势 E′ 恒定的经典模型,原动机-调速器采用二阶模型,并忽略调频死区、功率饱和等因素影响,在平衡点处线性化后得到如下模型。式中: δ 为发电机转子角; ωω0 分别为发电机转速和同步转速; MD 分别为发电机惯性时间常数和阻尼系数; Pe 为电磁功率; FHP 高压缸功率占比; TCH 为主进气容积和汽室时间常数;R为调速器下垂系数; y z 分别为原动机-调速器的状态变量; TRH 为再热器时间常数。1.2  自同步电压源模型研究频率稳定问题时可忽略电流控制环,只考虑功率控制环的自同步电压源模型[16]如图1所示。

图1  自同步电压源控制框图

Fig.1  Block diagram of autonomous-synchronization voltage source control


假设新能源机组的后备储能足够充分,并认为自同步电压源内电压Evir幅值恒定,得到只考虑有功功率部分的模型为

式中: Pref 为自同步电压源的有功功率参考值; Td 为考虑频率测量和变流器控制延迟的时间常数。


2  惯量响应阶段的时域解分析


电力系统发生频率扰动后,系统频率动态响应过程主要有扰动功率分配、惯量响应和频率调节阶段3个阶段[17]。惯量响应阶段主要研究系统惯量对频率变化率的影响,新能源接入对系统惯量的影响直接,尽管部分跟网型光伏发电系统具备紧急功率控制方案,但无法像自同步电压源一样提供惯量支撑。本章采用二次特征值分析理论对惯量变化所带来的影响进行解析分析。

2.1  状态空间模型响应时域解

对于线性系统,有

式中:x为状态变量;A为状态矩阵。式(3)的零输入表达式[18]式中:viui分别为矩阵A的右、左特征值向量,相应的矩阵满足 V−1=U V是由v按列组成的矩阵,Uu按行组成的矩阵; λi 为矩阵A的特征值。引入输入矩阵后系统模型变为式中: B 为输入矩阵; r 为输入向量。引入坐标变换式中:z为坐标变换后的状态变量。将式(6)代入式(5)后状态方程变为式中: Λλ A矩阵的特征值形成的对角阵;P=UB。式(7)可表示为n个解耦的一阶方程,第i个方程为式中: pij P矩阵的i行、j列元素。假设 rj 阶跃响应,其拉氏变换 Ri(s)=aj/s ,对式(8)两边取拉氏变换,有求得时域解为其初值满足由于发电机转子的绝对角度不是唯一的,因此根据式(1)和式(2)形成的电力系统A矩阵存在一个零特征值,假设零特征值为第一个特征值,得到初始状态变量的时域解为在频率稳定问题中,系统状态变量初值通常为平衡点的值,因此只需要研究其零状态响应解的特性,即2.2  二次特征值问题2.2.1  问题描述QEP已经在声学和流体力学中的线性稳定性等领域得到了广泛应用。QEP可找到满足以下条件的标量 λ 和非零向量vu[19-20]式中:K0K1K2为方阵;vu分别为 λ 对应的右、左特征值向量。假设系统中发电机只包含同步机和自同步电压源,在功率扰动发生后的惯量响应阶段,同步机调速器和变流器下垂控制还未生效,机械功率和功率参考值基本不变,通常存在频率基本呈线性变化的时段[21]。此时电网可以用二阶模型近似表示[22],即式中: Δδ 为发电机转子角形成的向量; Δω 为发电机转速形成的向量; MD分别为发电机惯性时间常数和阻尼系数形成的对角矩阵。 ΛΓ矩阵表达式推导如下。将负荷视为恒阻抗并入网络导纳矩阵,并引入内电势节点,得到系统的增广导纳矩阵为式中:E为同步机暂态电抗后的电势或变流器内电压; VN 为联络节点和负荷节点电压; YEE 为同步机暂态电纳或者自同步电压源内电势与端电压之间的电纳所形成的导纳阵; YEV YVE YVV 分别为式中: 为原网络导纳矩阵; Yload 为负荷等值导纳所形成的导纳阵。Kron化简后得到只包含发电机内电压的网络方程为式中: Yred Kvon 化简后的网络方程。令式中:real(·)、imag(·)分别为括号内数值的实部和虚部函数。由此可得 ΛΓ矩阵表达式分别为式中: Ei0 Ej别为发电机ij在平衡点处的内电势; δij0 为发电机ij在平衡点处的相角差; Gij Bij 分别为GBi行、j列元素。并且 Λ 为对称矩阵。令 L=Λ+Γ ,得到系统的雅克比矩阵为雅克比矩阵特征值和特征向量满足展开后得到将式(23)中第一个式子代入第二个式子,得到根据特征值的定义,其满足二次特征值问题,即2.2.2  基本性质二次特征值问题的解 λ 与矩阵L性质密切相关。考虑网络损耗的L矩阵为有向图的拉普拉斯矩阵,不能保证具有实特征值,但可以根据电网特性估计其数值范围[23]L矩阵可分解为Hermitian部分 LH 和斜Hermitian部分 LS ,即式中: LL矩阵的共轭转置,并且根据电网特性可知 Gijsinδij 小于 Bijcosδij ,因此斜对称部分元素数量级远小于对称部分。令 xCn ,则有式中:上划线表示共轭。于是, xLHx 为纯实数,而 xLSx 为纯虚数,并且实部大于虚部。L矩阵的特征值和特征向量满足 Lx=γx ,并假设 xx=1 ,则

式中: γ 为一个实部大、虚部小的复数,根据Gersgorin圆盘定理可知实部为正[24](如图2所示,L矩阵行和为零,因此特征值均位于右半平面)。


图2  L矩阵的一个圆盘

Fig.2  A disk of matrix L


计算经验表明, γ 几乎均为实数,虚数部分往往也很小[25],因此以下分析假设 γ 为实数。

考虑均匀阻尼 M−1D=bI 的情况,其中b为惯量阻尼比,则有

式中: γ 为由 M−1L 特征值组成的对角矩阵,因此式(29)等价于如下所示N个独立方程的解。进而可以得到N对解为对特征值分类,并进行小干扰稳定分析,步骤如下。1)对于零特征值 γ1=0 ,得到解 λ1=0 , λ2=−;2)对于实特征值 γi>0 ,根据劳斯判据可知其对应的解 λi1,2 具有负的实数部分。因此可知 γi⩾0 是系统小干扰稳定的一个充分条件。在满足小干扰稳定条件后,进一步要明确惯性对频率动态响应的影响,需要分析特征值随惯量矩阵变化的情况。实际电力系统中参数b通常很小,因此可见决定有功子系统动态特性的关键是广义特征值问题命题1:降低任意发电机组的惯性时间常数M,系统所有特征值的虚部近似单调上升。证明:注意到正定,L近似对称,那么令 ΔM⩾0 ,再令 ΔX 为某对角元全正矩阵,那么因为 Λ 正定,那么 ΔΔX 正定,根据特征值单调定理[26]因此根据可知,随着惯量矩阵元素的减小,系统特征值频率单调上升。2.3  惯量响应阶段的时域解性质得到惯量对系统特征值的作用后,就可以通过状态变量时域响应的解析解得到系统频率相关特征量与系统参数的变化关系。减小部分同步机出力以模拟发生功率缺额的情况,则可令式中: ΔPm 为同步机机械功率变化量矩阵。得到 ΔPm 发生阶跃变化时的模型为其零状态响应为结合命题1可知系统惯量的降低使得 λ 增大,惯量响应阶段的频率动态响应变快,从而恶化系统频率稳定性。另外需要考虑惯量与零时刻频率变化率的关系,令一台同步机为参考机可消去零特征值,可以得到状态变量的变化率为式中: V′U′ 分别为消去零特征值后的右、左特征值向量矩阵; Eλ 形成的对角矩阵。系统惯性中心频率变化 Δωcoi(t) 为对时间求导得到系统惯性中心频率变化率为式中: Mcol 为发电机惯性时间常数形成的列向量; Vω 表示右特征值向量矩阵中与转速相关的部分。零时刻的频率变化率为从式(41)可以看出惯量越低的系统在零时刻的频率变化率越大,结合式(37)的结果可知低惯量系统在惯量响应阶段的频率稳定性更差。

与同步机天然具有的惯量不同,目前自同步新能源需要通过加装储能装置实现较大的虚拟惯量,因此有必要对满足频率动态安全要求的最低虚拟惯量进行估计。基于式(41),在预设的频率动态安全条件下,可以估计出新能源机组所需的最低虚拟惯量,以达到有效利用电力资源的目的。


3  频率调节阶段的时域解分析


频率调节阶段的一次调频阶段主要研究系统调频作用对稳态频率误差的影响,目前以自同步电压源接入电网的新能源机组通常采取有功-频率下垂控制方式模拟同步机的调频作用。为分析调频系统对一次调频阶段响应特征的影响,本章对比考虑调频作用前后的稳态频率误差时域解表达式,进而确定调频系统对稳态频率误差的影响。3.1  不考虑调频作用的稳态频率误差表达式

假设不考虑调速器作用时,状态变量稳态误差可写成矩阵形式为

可知状态变量稳态误差可以用消去零特征值的雅克比矩阵的逆来表示,其逆表达式为其中因此稳态频率误差可表示为首先证明所有发电机组最后收敛至同一转速。已知矩阵左乘 A12 的意义为矩阵中除第一行元素以外的每一行减去第一行元素,因此将式(45)左乘 A12 ,得到可知 (−D−1+D−1L′(A12D−1L′)−1A12D−1) 为一个行相等矩阵,证明了转速终值向量各元素均相同,因此所有发电机组最后收敛至同一转速。从式(45)可以看出,稳态频率误差与惯量大小无关,但和系统阻尼成反比,系统阻尼越小,频率跌落程度越大。综上所述,在没有调速器的情况下,当出现有功功率缺额时,系统频率会单调下降,频率跌落速度随着系统惯量的下降而加快,频率跌落程度随着系统阻尼的减小而加深。3.2  考虑调频作用的稳态频率误差表达式考虑调速器作用后的系统统一模型为式中:I为单位矩阵。上述的调速器包含了自同步电压源的功率-频率下垂控制,其参数与上述模型中的参数对应关系为对应的雅克比矩阵为按照第2章推导步骤得到一个4次特征值问题,相应的特征多项式系数为按照2.2.2节做法,则特征值多项式同样可近似等价于n个方程组,即各系数为这里需要说明的是,即使自同步电压源与同步机参数存在差异,但各系数矩阵仍然是对角占优矩阵,上述的解耦方法是近似有效的。加入调速器后的小干扰稳定分析步骤如下。1)对于 γi≠0 ,可忽略一次项中的和二次项中的方程可以近似解耦为式(54)左边已做过讨论,右边可因式分解得到2)对于 γi=0 ,不可忽略,特征值多项式变为则稳定条件为对于火电机组,有以下近似关系。因此当γ 正定时,系统总是稳定的。命题2:加入调速器后,相当于系统阻尼增加至(D + R−1 ),减小了受扰动后系统的稳态频率误差。证明:加装调速器后,功率缺额时的模型为零状态响应为同样,状态变量稳态误差可以用消去零特征值的雅克比矩阵的逆来表示,将消去零特征值的雅克比矩阵分块可得根据分块矩阵求逆公式可知在左上子矩阵中对比式(43)后可知调速器的调频作用相当于增加了系统阻尼(DR−1 ),因此减小了受扰动后系统的稳态频率误差。自同步电压源机组包含大量电力电子器件,因此其电流耐受性不强,过大的功率输出会导致电力电子器件的损坏。因此有必要对自同步电压源机组的功率下垂系数做出估计,在满足频率动态安全要求的情况下,保证电力电子器件工作在安全范围以内。基于上述结论,在预设的频率动态安全条件下,提出一种估计自同步电压源机组下垂系数的方法,步骤如下。1)在给定自同步电压源机组初值的情况下,计算系统的等效阻尼为式中: Deq0 为估计前的系统等效阻尼; RSG 为同步机的调差系数; RVSG0 为估计前的自同步电压源机组的下垂系数。2)基于式(45)得到估计前的稳态频率误差 Δω0(∞) 。3)计算 Δω0(∞) 与给定的稳态频率误差 Δωcmd(∞) 的比值,乘以 Deq0 后得到满足稳态频率误差要求的 Deq_cmd

4)根据式(65)反推得到满足频率稳态值要求的再根据每台自同步电压源机组的备用有功功率和稳态频率误差 Δωcmd(∞) 计算出每台机组能提供的最大最后按照实际需求和功率限制分配到自同步电压源机组。


4  算例分析


4.1  模型准确性验证

本节将所提出的全状态模型、SFR模型和时域仿真结果进行对比,验证其对系统频率动态描述的准确性,SFR模型聚合参数方法见文献[9]。仿真在10机39节点系统算例中进行,同步机部分参数如表1所示,其余参数TCH= 0.3 s,TRH= 7 s,FHP = 0.3。计算得到SFR模型中的聚合参数MSFR = 782 s,DSFR = 7.82,RSFR = 0.0047。


表1  同步机部分参数

Table 1  Partial parameters of synchronous generators


t = 0 s时刻减小1号机出力100 MW,利用3种方法得到30 s内系统惯性中心频率在发生功率扰动后的变化。全状态模型、SFR模型与时域仿真法计算结果的偏差如图3示。


图3  两种模型与时域仿真法的误差

Fig.3  Errors between the two models and the time-domain simulation method


由图3可以看出,全状态模型的计算精确度高于SFR模型,这是由于在SFR模型中忽略了部分环节,而全状态模型保留了包括网络在内的模型,因此具有更高的精度。4.2  广义特征值变分

为验证发电机惯量对广义特征值的影响,在仿真算例中将同步机逐台替换为低惯量的自同步电压源,其参数如表2所示,并计算其对应的广义特征值大小,结果如图4所示。


表2  自同步电压源参数

Table 2  Parameters of autonomous-synchronization voltage source


图4  惯量-广义特征值频率单调关系

Fig.4  Monotonic relationship between inertia and generalized eigenvalues


由图4可以看出,随着同步机被替换为低惯量的新能源机组,广义特征值单调上升,也意味着系统特征值频率单调上升,说明惯量减小降低了系统频率的抗干扰性,在发生扰动时系统频率振荡更快。其他仿真算例也具有类似的结果。4.3  基于全状态模型的最低虚拟惯量估计

目前对电网频率变化率没有明确要求[27],因此假设在发生100 MW的功率阶跃变化时,频率变化率0.06 Hz/s,即0.0012 (p.u.)/s作为限制,在表1中将7~10号同步机替换为自同步电压源机组,根据式(41)估计满足限制要求的系统最低惯量为833.33 s,因此4组自同步电压源机组虚拟惯量总和的最低要求为178.33 s。在满足系统最低惯量要求时,发生100 MW功率阶跃1 s内的惯性中心频率动态曲线如图5所示。


图5  1 s内的惯性中心频率动态响应

Fig.5  Inertia center frequency dynamic response in 1 s


由图5可知,1 s内系统惯性中心频率偏差小于0.06 Hz,即根据所提方法安排自同步电压源机组的虚拟惯量大小可以满足频率动态安全要求。4.4  基于全状态模型的自同步电压源功率下垂系数估计在4.3节的电网中,假设在发生100 MW的功率阶跃变化时,频率稳态误差不超过0.3 Hz。基于3.2节所提的估计自同步电压源机组下垂系数的方法,其中假设自同步电压源机组初始下垂系数均为0.5。1)计算系统的等效阻尼 Deq0=136.33 ;2)计算估计前的稳态频率误差 Δω0(∞)=−0.3180Hz ;3)计算 Δω0(∞) 与给定的系统稳态频率误差 Δωcmd(∞)=−0.3Hz 的比值,乘以 Deq0 后得到满足稳态频率误差要求的 Deq_cmd=144.53 ;4)根据式(65)反推得到满足频率稳态值要求的假设平均分配到4台自同步电压源机组,则每台自同步电压源机组的下垂系数应调整至0.247。假设每台自同步电压源机组的备用有功功率为0.05 p.u.,则每台自同步电压源机组最多能提供的因此上述下垂系数分配满足了功率限制要求。

根据上述参数估计结果,在发生100 MW功率阶跃变化时,仿真计算调整前后的系统,其惯性中心频率50 s内动态响应结果如图6所示。


图6  调整参数前后的惯性中心频率动态响应

Fig.6  Inertia center frequency dynamic response before and after adjusting parameters


由图6可知,调整自同步电压源机组下垂系数后,系统惯性中心频率的一次调频稳态频率误差达到了预设的要求,证明了所提方法的有效性。


5  结论


系统的惯性时间常数变小是新能源并入电网后带来的最显著的变化之一,其渗透率不断上升将导致电力系统暂态特性从量变到质变。为解决高比例新能源电力系统频率稳定分析技术问题,本文提出了一种基于电力系统近似全状态模型的电力系统频率分析方法,并得出以下结论。1)基于二次特征值分析证明了惯量直接影响系统受功率扰动后惯量响应阶段的频率变化率,惯量越小,系统频率变化越快,抗扰动性越差。2)基于电力系统近似全状态模型的时域解分析,证明了自同步电压源新能源的下垂特性减小了系统频率稳态误差,并且系统频率稳态误差与自同步电压源新能源的下垂系数成正比。3)基于上述结论,根据所提出的估计满足预设的频率动态安全要求的自同步电压源新能源虚拟惯量和下垂系数的方法,在充分利用电力资源和保证电力电子器件安全的条件下,使得新能源并网系统能够满足频率变化率和频率稳态值的要求。相较于时域仿真法需要不断调整参数和多次仿真,论文所提方法只通过一次计算即可达到目的,提高了计算效率。4)相较于传统分析方法,本文方法的建模工作量小,计算量少,能够更快速地获得频率动态响应特征;减小了模型聚合带来的误差,一定程度上能够弥补SFR模型不能完全满足电力电子化的电力系统频率稳定分析要求的缺陷。后续将开展含跟网型新能源电力系统的频率稳定分析工作。

(责任编辑 许晓艳)



作者介绍

李振垚(1995—),男,博士研究生,从事电力系统稳定分析与控制技术研究,E-mail:12110085@zju.edu.cn;


甘德强(1966—),男,通信作者,教授,博士生导师,从事电力系统稳定性分析与控制技术研究,E-mail:deqiang.gan@ieee.org.


欢迎点击文后“阅读原文”跳转期刊官网,获取更多信息!




 往期回顾 


《中国电力》2023年第5期目录
【精彩论文】基于t-SNE降维和放射传播聚类算法的低压配电网相位识别【精彩论文】高可靠性三相多功能并网变换器研究【精彩论文】改性氧化镁/环氧树脂复合材料的介电及亲疏水性能【精彩论文】抑制母线电压跌落的热泵温度轨迹规划算法【征稿启事】“分布式智能电网的规划、运行和电力交易”专栏征稿启事【征稿启事】“新型能源体系下电碳协同市场机制及优化运行”专栏征稿启事【征稿启事】“面向碳达峰碳中和目标的清洁高效发电技术”专题征稿启事【征稿启事】“新型电力系统低碳规划与运行”专栏征稿启事

编辑:杨彪
校对:于静茹

审核:方彤

声明

根据国家版权局最新规定,纸媒、网站、微博、微信公众号转载、摘编《中国电力》编辑部的作品,转载时要包含本微信号名称、二维码等关键信息,在文首注明《中国电力》原创。个人请按本微信原文转发、分享。欢迎大家转载分享。

继续滑动看下一个

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存