CN104361159A - 一种大规模电力***暂态稳定时间空间并行仿真方法 - Google Patents

Abstract

本发明提供一种大规模电力***暂态稳定时间空间并行仿真方法，包括以下步骤：建立交直流***模型；对交直流***进行空间解耦，并对直流***进行时间解耦；交直流***并行仿真计算。本发明采用CPU/GPU异构平台实现了大规模交直流互联电力***的暂态仿真，按照交流***和直流***空间上的独立性以及其在电气方面容易解耦的特点，将交流***对应的计算任务分配到CPU上，将直流***对应的计算任务分配到GPU上，实现空间上的并行；采用移动时间窗的技术，避免GPU计算资源的浪费，提高了GPU的利用效率。

Description

一种大规模电力***暂态稳定时间空间并行仿真方法

技术领域

本发明涉及一种仿真方法，具体涉及一种大规模电力***暂态稳定时间空间并行仿真方法。

背景技术

暂态稳定计算的任务是通过数值计算方法判断电力***遭受较大扰动之后，各个发电机是否能继续保持同步运行。通常需要考虑的电力***大扰动包括：

1、负荷的突然变化，比如大容量用户的投入和切出；

2、切出或者投入***的主要元件，如发电机、变压器或者线路；

3、发生短路故障。

随着电力工业的发展，电网之间的互联更加紧密，我国电网正在形成大规模交直流互联电力网络。互联电网的规模扩大之后，对交直流***的暂态稳定仿真速度也提出了更高的要求。因此暂态稳定计算的并行算法研究具有实际意义。目前在并行计算框架下，仿真研究主要关注三个方面：

1、直流***建模计算方法；

2、直流***和交流***之间的接口；

3、并行计算软硬件平台。

按照对直流线路和直流控制模拟的详细程度，建模方法可分为响应模型、详细模型和电磁暂态模型。

暂态稳定的计算平台有CPU、GPU、FPGA等等。其中，CPU作为传统的通用计算处理器，单核能力强，适合各种类型的计算任务，但是包含的核心数量通常不多；GPU(GraphicProcessing Unit)为一类专门用于处理图形数据，作为一种较晚进入通用计算领域的器件，较CPU具有处理器核数众多，线程分配与销毁迅速的特点。但是其逻辑处理能力相对较弱，更加适合数据密集型计算。近些年，将GPU作为一种面向通用计算的流处理器的做法也越来越普遍。GPU可以用于多种并行计算任务，比如分子动力学计算。他们非常适合数据输入输出量非常大的计算。大量的数据使得GPU可以充分地利用GPU的向量计算单元或者单指令多数据的结构。基于GPU的计算在大规模的计算中发挥了越来越大的作用，世界上最强的十台超级计算机中，有三台都利用了GPU的优势。由CPU和GPU组成的异构平台，能够融合两种处理器的优势，具有更高的计算能力，但是也要求开发人员提出适合这种新平台的并行算法。

在大规模的交直流电网***中，直流***本身的独立性较强，容易从整个***中单独划分出来。因此可以结合交直流***仿真特性，将电网***中的交流部分和直流部分划分开来，将其计算任务分配到异构平台的CPU和GPU上。使用CPU计算交流***部分，使用GPU计算直流***部分。同时，多个直流***之间的互相影响也相对较小，具有较高的天然并行性。考虑直流详细模型的特点之后，将直流***的时域仿真任务按照时步进行划分，构成流水线的形式，可以充分利用GPU的强大计算能力，有效提高直流***的仿真计算速度。

流水线指的是按照以下规律形成的一种数据处理方法，一个元素的输出等于另外一个的输入。一条流水线上的元素经常是并行执行的，按照时间划分出区域分块放置的。多条流水线可以充分暴露任务执行的并行性，提高处理器计算资源的利用效率。

暂态稳定问题指的是电力***遭受较大扰动后，各个发电机是否能继续保持同步运行的问题。暂态稳定计算即为利用数值计算方法研究遭受扰动后，***各项指标的变化。

并行计算是一种可以同时进行多项任务的计算。其基本理念在于大的问题常常可以划分为小的问题，而这些小的问题通常可以同时解决。并行计算有多种形式：位级并行，指令级并行，数据和任务并行。并行计算在高性能计算领域已经有多年的应用历史。在民用领域，由于单核处理器性能的限制，并行计算的重要性也越来越得到重视。

发明内容

为了克服上述现有技术的不足，本发明提供一种大规模电力***暂态稳定时间空间并行仿真方法，基于大规模电力***交直流***空间分布独立，电气关系易于解耦的特点，将交流***部分仿真计算任务在CPU上进行，直流***部分仿真计算任务在GPU上使用流模式实现并行；利用异构平台CPU和GPU之间的数据传递可以异步执行的特点，开启多个流控制直流交流***间的计算数据传递。

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

本发明提供一种大规模电力***暂态稳定时间空间并行仿真方法，所述方法包括以下步骤：

步骤1：建立交直流***模型；

步骤2：对交直流***进行空间解耦，并对直流***进行时间解耦；

步骤3：交直流***并行仿真计算。

所述步骤1中，交直流***模型包括交流***动态模型、直流***模型和交直流电网模型；

所述交流***动态模型包括发电机模型、励磁***模型和原动机***模型；所述直流***模型包括直流电缆模型、直流控制***模型和换流器模型。

所述直流电缆模型采用T型模型，直流电缆模型对应的动态方程如下：

${\stackrel{\·}{I}}_{\mathrm{dcr}}=(V_{\mathrm{dcr}}-V_c-0.5R_{\mathrm{dc}}{I}_{\mathrm{dcr}})/(L_{\mathrm{sr}}+L_{\mathrm{line}}/2)$

${\stackrel{\·}{I}}_{\mathrm{dcr}}=({-V}_{\mathrm{dci}}+V_c-0.5R_{\mathrm{dc}}{I}_{\mathrm{dci}})/(L_{\mathrm{si}}+L_{\mathrm{line}}/2)$

${\stackrel{\·}{V}}_c=I_c/C_{\mathrm{cap}}=(I_{\mathrm{dcr}}-I_{\mathrm{dci}})/C_{\mathrm{cap}}$

其中，V_dcr为整流侧的直流电压，V_dci为逆变侧的直流电压；I_dcr为整流侧得直流电流，I_dci为逆变侧的直流电流；V_c和I_c为电容电压与电流；R_dc、L_line和C_cap分别为直流电缆电阻、电感和电容；L_sr为整流侧平滑电抗器电感，L_si为逆变侧平滑电抗器电感。

所述直流控制***模型从顶层到底层分为主控制级和极控制级；

所述主控制级包含紧急功率控制器，当直流***所连的交流***电压或频率发生较大波动时改变直流电流或功率参考值；

所述极控制级分为以下三个环节：

1)低压限流控制器：其输入为主控制级设定的电流参考值，在直流电压过低时，限制该设定的电流参考值，并将电流参考值输入给电流增益控制器；

2)电流增益控制器：其将低压限流控制器输入的电流参考值转换为触发角参考值，输出给熄弧角控制器；

3)熄弧角控制器：将触发角参考值转换为实际触发脉冲。

所述换流器模型为换流器准稳态模型，表示为：

$V_{\mathrm{dcr}}=N[\frac{3\sqrt{2}{T}_r}{\mathrm{\π}{T}_{\mathrm{tap}}}{V}_{\mathrm{adcr}}\cos \mathrm{\α}-\frac{3X_c{I}_{\mathrm{dcr}}}{\mathrm{\π}}-2R_c{I}_{\mathrm{dcr}}]$

$\cos (\mathrm{\δ}+\mathrm{\α})=\cos \mathrm{\α}-\frac{\sqrt{2}{T}_{\mathrm{tap}}{I}_{\mathrm{dcr}}{X}_c}{T_r{V}_{\mathrm{adcr}}}$

$P_{\mathrm{acr}}=N\frac{\sqrt{18}{T}_r}{2\mathrm{\π}{T}_{\mathrm{tap}}}{I}_{\mathrm{dcr}}{V}_{\mathrm{adcr}}[\cos \mathrm{\α}+\cos (\mathrm{\α}+\mathrm{\δ})]$

$I_{\mathrm{adcr}}=\frac{P_{\mathrm{dcr}}{V}_{\mathrm{adcrx}}+Q_{\mathrm{dcr}}{V}_{\mathrm{adcry}}}{V_{\mathrm{adcrx}}^2+V_{\mathrm{adcry}}^2}+j\frac{P_{\mathrm{dcr}}{V}_{\mathrm{adcry}}-Q_{\mathrm{dcr}}{V}_{\mathrm{adcrx}}}{V_{\mathrm{adcrx}}^2+V_{\mathrm{adcry}}^2}$

其中，V_dcr和I_dcr分别为整流侧的直流电压和直流电流，N为串联整流桥个数，T_tap为整流侧变压器分接头位置，I_adcr为整流侧注入交流***电流，V_adcr为整流侧变压器交流侧电压，V_adcrx和V_adcry分别为整流侧变压器交流侧电压实部和虚部，X_c和R_c分别为整流侧变压器电抗和电阻，α为整流侧触发角，δ为换相重叠角，T_r为整流侧变压器基本变比，P_acr和Q_acr分别为整流侧吸收的有功功率和无功功率，为功率因数角。

所述步骤2中，对建立的交直流***进行空间解耦具体过程如下：

交流***动态模型表示为：

${\stackrel{\·}{x}}_{\mathrm{ac}}=f(x_{\mathrm{ac}},V_{\mathrm{ac}})$

其中，x_ac和V_ac分别为交流***状态变量和纯交流节点电压向量；

直流***模型表示为：

${\stackrel{\·}{x}}_{\mathrm{dc}}=H(x_{\mathrm{dc}},V_{\mathrm{adc}})$

其中，x_dc和V_adc分别为直流***状态变量及直流***接入交流***处的节点电压向量；

交直流电网模型表示为：

0＝I-Y_NV

其中，I为交直流***节点注入电流，其包括V为交直流***节点电压，Y_N为网络导纳矩阵；

由于交流***状态变量和直流***状态变量相互独立，且通过网络方程相关联，按照交流***和直流***空间上的位置关系进行解耦，直流***采用诺顿等效接入交直流***电网网络，在CPU部署交流***部分计算数据，即交流***动态方程和网络方程中变量初始值及参数值，GPU部署直流***部分计算数据，即直流***动态方程中变量初始值及参数值。

所述步骤2中，对直流***进行时间解耦，具体过程如下：

对直流***动态方程略去下标，可得直流微分方程，有：

$\stackrel{\·}{x}=H(x,V)$

在交流积分步长h_ac中，直流微分方程采用单时步顺序求解的方法，计算h_ac/h_dc步，其中h_dc为直流积分步长；使用隐式梯形积分法，设下标n表示t时刻，n+1表示t+h_ac时刻，直流微分方程的差分化方程为：

$x_n=x_{n-1}+\frac{h_{\mathrm{dc}}}{2}[H(x_{n-1},V_{n-1})+H(x_n,V_n)],n=1,...,h_{\mathrm{ac}}/h_{\mathrm{dc}}$

其中，x_n为直流***t时刻的状态变量，x_n-1为直流***t-h_ac时刻的状态变量，V_n为直流***t接入交流***处的节点电压向量，V_n-1为直流***t-h_ac接入交流***处的节点电压向量；

定义中间向量R_n，其表示为：

$R_n=x_n-x_{n-1}-\frac{h_{\mathrm{dc}}}{2}[H(x_{n-1},V_{n-1})+H(x_n,V_n)]$

对上式采用牛顿-拉夫逊法，第k次迭代公式为：

$R_n^{k-1}=-\mathrm{J\Δ}{x}_n^k$

其中，k为迭代次数，为两次迭代中直流***状态变量差值，J为牛顿-拉夫逊法中第k次迭代的雅克比矩阵，和J分别表示为：

$\mathrm{\Δ}{x}_n^k=x_n^k-x_n^{k-1}$

$J=\frac{\∂R_n}{\∂x_n}{|}_{x_n=x_n^{k-1}}$

其中，为第k次迭代中直流***t时刻的状态变量，为第k-1次迭代中直流***t时刻的状态变量；

于是得到其中V₀为上一交流时步所得电压值。

所述步骤3具体包括以下步骤：

步骤3-1：交直流***进行初始化，采用OpenMP多线程在CPU上开启两个线程，分别为线程A和线程B，线程A在CPU上完成交流***相关的计算，线程B用于控制GPU，向GPU分配任务、发送控制指令；然后交直流***进入时域仿真过程，按照时间的推进，计算每个时刻对应的交直流***状态变量，设置仿真时刻t＝0，设推进步长为一个交流时步h_ac；

步骤3-2：交流***和直流***间交换边界点电压和电流值；

步骤3-3：采用交替迭代的隐式梯形积分方法，在异构平台中的CPU上计算交流***状态变量；

步骤3-4：在异构平台的GPU上完成多个直流***状态变量的计算，具体包括以下步骤：

步骤3-4-1：判断是否有多条直流线路，若是，则开启多个GPU流，每个流计算一条直流线路；若否，则只使用一个GPU流；

步骤3-4-2：每个直流线路使用结合移动窗的时间流水线并行方法计算，并结合电压值，判断流水线的条数；

步骤3-4-3：计算收敛后，计算直流***注入交流***的电流；

步骤3-5：当获得全部直流线路注入电流之后，设置仿真时刻t＝t+h_ac，进入下一时刻；

步骤3-6：重复步骤3-1至步骤3-4，直到仿真结束。

所述步骤3-4-2中，将本仿真时段划分为多个时间窗口，而每个时间窗口包含n_w个时步，多个时间窗口按顺序计算，但是在每个时间窗口内，将多个时步的计算任务装配为流水线的形式，实现并行计算；

对于时间窗口1，n_w个时步对应的差分方程为：

$\left[\begin{array}{c}{R}_1\\ R_2\\ \·\\ \·\\ \·\\ R_{n_w}\end{array}\right]=\left[\begin{array}{c}{x}_1-x_0-\frac{h_{\mathrm{dc}}}{2}[H(x_1,V_1)+H(x_0,V_0)]\\ x_2-x_1-\frac{h_{\mathrm{dc}}}{2}[H(x_2,V_2)+H(x_1,V_1)]\\ \·\\ \·\\ \·\\ x_{n_w}-x_{n_w-1}-\frac{h_{\mathrm{dc}}}{2}[H(x_{n_w},V_{n_w})+H(x_{n_w-1},V_{n_w-1})]\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \·\\ \·\\ \·\\ 0\end{array}\right]$

对上式采用牛顿-拉夫逊法，第k次迭代公式为：

其中： $J_{n,n-1}^{k-1}=\frac{\∂R_n}{\∂x_{n-1}}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^{k-1}}},J_{n,n}^{k-1}=\frac{\∂R_n}{\∂x_n}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^{k-1}}},n=1,...,n_w;$

对于时间窗口中的n_w个时步，在第一次迭代中，流水线P1进行时步1的第一次迭代求解，获得时步1的状态变量在第二次迭代中，流水线P2进行时步2的第一次迭代求解以及时步1的第二次迭代求解；依次类推，实现对当前时间窗口中全部时步的并行求解；求解公式为：

$\mathrm{\Δ}{x}_n^k=\left\{\begin{array}{cc}-{\left(J_{n,n}^{*k-1}\right)}^{-1}{R}_n^{k-1}& ,n=1\\ -{\left(J_{n,n}^{*k-1}\right)}^{-1}[R_n^{k-1}+J_{n,n-1}^{*k}\mathrm{\Δ}{x}_{n-1}^k]& ,n=2,...,n_w\end{array}\right.$

其中， $J_{n,n-1}^{*k-1}=\frac{\∂R_n}{\∂x_{n-1}}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^k}},J_{n,n}^{*k-1}=\frac{\∂R_n}{\∂x_{n-1}}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^k}}.$

与现有技术相比，本发明的有益效果在于：

1、采用CPU/GPU异构平台实现了大规模交直流互联电力***的暂态仿真。

2、按照交流***和直流***空间上的独立性以及其在电气方面容易解耦的特点，将交流***对应的计算任务分配到CPU上，将直流***对应的计算任务分配到GPU上，实现空间上的并行。

3、针对多个直流***互相独立的特点，将其对应的计算任务在GPU上实现时间上的并行。

4、针对直流***的暂态仿真，将其计算任务划分为时间窗，装配成多条流水线的形式，利用GPU的多核并行处理一个时间窗内的各个时步计算任务。

5、采用移动时间窗的技术，避免GPU计算资源的浪费，提高了GPU的利用效率。

附图说明

图1是本发明实施例中大规模电力***暂态稳定时间空间并行仿真方法流程图；

图2是本发明实施例中直流电缆模型示意图；

图3是本发明实施例中直流控制***模型示意图；

图4是本发明实施例中换流器模型示意图；

图5是本发明实施例中时间窗口序列示意图；

图6是本发明实施例中流水线及移动窗示意图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

如图1，本发明提供一种大规模电力***暂态稳定时间空间并行仿真方法，所述方法包括以下步骤：

步骤1：建立交直流***模型；

步骤2：对交直流***进行空间解耦，并对直流***进行时间解耦；

步骤3：交直流***并行仿真计算。

所述步骤1中，交直流***模型包括交流***动态模型、直流***模型和交直流电网模型；

所述交流***动态模型包括发电机模型、励磁***模型和原动机***模型；所述直流***模型包括直流电缆模型、直流控制***模型和换流器模型。

如图2，所述直流电缆模型采用T型模型，直流电缆模型对应的动态方程如下：

${\stackrel{\·}{I}}_{\mathrm{dcr}}=(V_{\mathrm{dcr}}-V_c-0.5R_{\mathrm{dc}}{I}_{\mathrm{dcr}})/(L_{\mathrm{sr}}+L_{\mathrm{line}}/2)---\left(1\right)$

${\stackrel{\·}{I}}_{\mathrm{dcr}}=({-V}_{\mathrm{dci}}+V_c-0.5R_{\mathrm{dc}}{I}_{\mathrm{dci}})/(L_{\mathrm{si}}+L_{\mathrm{line}}/2)---\left(2\right)$

${\stackrel{\·}{V}}_c=I_c/C_{\mathrm{cap}}=(I_{\mathrm{dcr}}-I_{\mathrm{dci}})/C_{\mathrm{cap}}---\left(3\right)$

其中，V_dcr为整流侧的直流电压，V_dci为逆变侧的直流电压；I_dcr为整流侧得直流电流，I_dci为逆变侧的直流电流；V_c和I_c为电容电压与电流；R_dc、L_line和C_cap分别为直流电缆电阻、电感和电容；L_sr为整流侧平滑电抗器电感，L_si为逆变侧平滑电抗器电感。

如图3，所述直流控制***模型从顶层到底层分为主控制级、极控制级和阀控制级；在详细模型中，换流器采用准稳态模型，因此忽略了阀控制级的作用；

所述主控制级包含紧急功率控制器，当直流***所连的交流***电压或频率发生较大波动时改变直流电流或功率参考值；

所述极控制级分为以下三个环节：

1)低压限流控制器：其输入为主控制级设定的电流参考值，在直流电压过低时，限制该设定的电流参考值，并将电流参考值输入给电流增益控制器；

2)电流增益控制器：其将低压限流控制器输入的电流参考值转换为触发角参考值，输出给熄弧角控制器；

3)熄弧角控制器：将触发角参考值转换为实际触发脉冲。

如图4，所述换流器模型为换流器准稳态模型，表示为：

$V_{\mathrm{dcr}}=N[\frac{3\sqrt{2}{T}_r}{\mathrm{\π}{T}_{\mathrm{tap}}}{V}_{\mathrm{adcr}}\cos \mathrm{\α}-\frac{3X_c{I}_{\mathrm{dcr}}}{\mathrm{\π}}-2R_c{I}_{\mathrm{dcr}}]---\left(4\right)$

$\cos (\mathrm{\δ}+\mathrm{\α})=\cos \mathrm{\α}-\frac{\sqrt{2}{T}_{\mathrm{tap}}{I}_{\mathrm{dcr}}{X}_c}{T_r{V}_{\mathrm{adcr}}}---\left(5\right)$

$P_{\mathrm{acr}}=N\frac{\sqrt{18}{T}_r}{2\mathrm{\π}{T}_{\mathrm{tap}}}{I}_{\mathrm{dcr}}{V}_{\mathrm{adcr}}[\cos \mathrm{\α}+\cos (\mathrm{\α}+\mathrm{\δ})]---\left(6\right)$

$I_{\mathrm{adcr}}=\frac{P_{\mathrm{dcr}}{V}_{\mathrm{adcrx}}+Q_{\mathrm{dcr}}{V}_{\mathrm{adcry}}}{V_{\mathrm{adcrx}}^2+V_{\mathrm{adcry}}^2}+j\frac{P_{\mathrm{dcr}}{V}_{\mathrm{adcry}}-Q_{\mathrm{dcr}}{V}_{\mathrm{adcrx}}}{V_{\mathrm{adcrx}}^2+V_{\mathrm{adcry}}^2}---\left(8\right)$

其中，V_dcr和I_dcr分别为整流侧的直流电压和直流电流，N为串联整流桥个数，T_tap为整流侧变压器分接头位置，I_adcr为整流侧注入交流***电流，V_adcr为整流侧变压器交流侧电压，V_adcrx和V_adcry分别为整流侧变压器交流侧电压实部和虚部，X_c和R_c分别为整流侧变压器电抗和电阻，α为整流侧触发角，δ为换相重叠角，T_r为整流侧变压器基本变比，P_acr和Q_acr分别为整流侧吸收的有功功率和无功功率，为功率因数角。

所述步骤2中，对建立的交直流***进行空间解耦具体过程如下：

交直流***动态模型表示为：

${\stackrel{\·}{x}}_{\mathrm{ac}}=f(x_{\mathrm{ac}},V_{\mathrm{ac}})---\left(9\right)$

其中，x_ac和V_ac分别为交流***状态变量和纯交流节点电压向量；

直流***模型表示为：

${\stackrel{\·}{x}}_{\mathrm{dc}}=H(x_{\mathrm{dc}},V_{\mathrm{adc}})---\left(10\right)$

其中，x_dc和V_adc分别为直流***状态变量及直流***接入交流***处的节点电压向量；

交直流电网模型表示为：

0＝I-Y_NV   (11)

其中，I为交直流***节点注入电流，其包括V为交直流***节点电压，Y_N为网络导纳矩阵；

由于交流***状态变量和直流***状态变量相互独立，且通过网络方程相关联，按照交流***和直流***空间上的位置关系进行解耦，直流***采用诺顿等效接入交直流***电网网络，在CPU部署交流***部分计算数据，即交流***动态方程和网络方程中变量初始值及参数值，GPU部署直流***部分计算数据，即直流***动态方程中变量初始值及参数值。

所述步骤2中，对直流***进行时间解耦，具体过程如下：

对直流***动态方程略去下标，可得直流微分方程，有：

$\stackrel{\·}{x}=H(x,V)---\left(12\right)$

在交流积分步长h_ac中，直流微分方程采用单时步顺序求解的方法，计算h_ac/h_dc步，其中h_dc为直流积分步长；使用隐式梯形积分法，设下标n表示t时刻，n+1表示t+h_ac时刻，直流微分方程的差分化方程为：

$x_n=x_{n-1}+\frac{h_{\mathrm{dc}}}{2}[H(x_{n-1},V_{n-1})+H(x_n,V_n)],n=1,...,h_{\mathrm{ac}}/h_{\mathrm{dc}}---\left(13\right)$

其中，x_n为直流***t时刻的状态变量，x_n-1为直流***t-h_ac时刻的状态变量，V_n为直流***t接入交流***处的节点电压向量，V_n-1为直流***t-h_ac接入交流***处的节点电压向量；

定义中间向量R_n，其表示为：

$R_n=x_n-x_{n-1}-\frac{h_{\mathrm{dc}}}{2}[H(x_{n-1},V_{n-1})+H(x_n,V_n)]---\left(14\right)$

对上式采用牛顿-拉夫逊法，第k次迭代公式为：

$R_n^{k-1}=-\mathrm{J\Δ}{x}_n^k---\left(15\right)$

其中，k为迭代次数，为两次迭代中直流***状态变量差值，J为牛顿-拉夫逊法中第k次迭代的雅克比矩阵，和J分别表示为：

$\mathrm{\Δ}{x}_n^k=x_n^k-x_n^{k-1}---\left(16\right)$

$J=\frac{\∂R_n}{\∂x_n}{|}_{x_n=x_n^{k-1}}---\left(17\right)$

其中，为第k次迭代中直流***t时刻的状态变量，为第k-1次迭代中直流***t时刻的状态变量；

于是得到其中V₀为上一交流时步所得电压值。

所述步骤3具体包括以下步骤：

步骤3-1：交直流***进行初始化，采用OpenMP多线程在CPU上开启两个线程，分别为线程A和线程B，线程A在CPU上完成交流***相关的计算，线程B用于控制GPU，向GPU分配任务、发送控制指令；然后交直流***进入时域仿真过程，按照时间的推进，计算每个时刻对应的交直流***状态变量，设置仿真时刻t＝0，设推进步长为一个交流时步h_ac；

步骤3-2：交流***和直流***间交换边界点电压和电流值；

步骤3-3：采用交替迭代的隐式梯形积分方法，在异构平台中的CPU上计算交流***状态变量；

步骤3-4：在异构平台的GPU上完成多个直流***状态变量的计算，具体包括以下步骤：

步骤3-4-1：判断是否有多条直流线路，若是，则开启多个GPU流，每个流计算一条直流线路；若否，则只使用一个GPU流；

步骤3-4-2：每个直流线路使用结合移动窗的时间流水线并行方法计算，并结合电压值，判断流水线的条数；

步骤3-4-3：计算收敛后，计算直流***注入交流***的电流；

步骤3-5：当获得全部直流线路注入电流之后，设置仿真时刻t＝t+h_ac，进入下一时刻；

步骤3-6：重复步骤3-1至步骤3-4，直到仿真结束。

所述步骤3-4-2中，将本仿真时段划分为多个时间窗口，而每个时间窗口包含n_w个时步，多个时间窗口按顺序计算，但是在每个时间窗口内，将多个时步的计算任务装配为流水线的形式，实现并行计算；

如图5，对于时间窗口1，n_w个时步对应的差分方程为：

$\left[\begin{array}{c}{R}_1\\ R_2\\ \·\\ \·\\ \·\\ R_{n_w}\end{array}\right]=\left[\begin{array}{c}{x}_1-x_0-\frac{h_{\mathrm{dc}}}{2}[H(x_1,V_1)+H(x_0,V_0)]\\ x_2-x_1-\frac{h_{\mathrm{dc}}}{2}[H(x_2,V_2)+H(x_1,V_1)]\\ \·\\ \·\\ \·\\ x_{n_w}-x_{n_w-1}-\frac{h_{\mathrm{dc}}}{2}[H(x_{n_w},V_{n_w})+H(x_{n_w-1},V_{n_w-1})]\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \·\\ \·\\ \·\\ 0\end{array}\right]---\left(18\right)$

对上式采用牛顿-拉夫逊法，第k次迭代公式为：

其中： $J_{n,n-1}^{k-1}=\frac{\∂R_n}{\∂x_{n-1}}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^{k-1}}},J_{n,n}^{k-1}=\frac{\∂R_n}{\∂x_n}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^{k-1}}},n=1,...,n_w;$

式(19)中的方程系数矩阵是大型对角矩阵，如果将其完全分解之后再将解向量代入下一次迭代，依次执行迭代过程的话，并不能对整个仿真过程进行加速。为了利用GPU的并行能力，下面引入流水线技术。

由于多次迭代过程对应的计算具有前后依赖关系，同一时间窗内的各个时步之间也具有依赖关系，即计算完成之后才可计算然而的计算可以同时进行，所以这样的计算过程还具有较大的并行潜力。按照这样的思路，可以将整个过程组成流水线形式。

对于时间窗口中的n_w个时步，在第一次迭代中，流水线P1进行时步1的第一次迭代求解，获得时步1的状态变量在第二次迭代中，流水线P2进行时步2的第一次迭代求解以及时步1的第二次迭代求解；依次类推，实现对当前时间窗口中全部时步的并行求解；求解公式为：

$\mathrm{\Δ}{x}_n^k=\left\{\begin{array}{cc}-{\left(J_{n,n}^{*k-1}\right)}^{-1}{R}_n^{k-1}& ,n=1\\ -{\left(J_{n,n}^{*k-1}\right)}^{-1}[R_n^{k-1}+J_{n,n-1}^{*k}\mathrm{\Δ}{x}_{n-1}^k]& ,n=2,...,n_w\end{array}\right.---\left(20\right)$

其中， $J_{n,n-1}^{*k-1}=\frac{\∂R_n}{\∂x_{n-1}}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^k}},J_{n,n}^{*k-1}=\frac{\∂R_n}{\∂x_{n-1}}{|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^k}}.$

多条流水线由CUDA中的block实现，每个流水线对应一个block；第n次迭代过程如下：

1)取第n-1次迭代获得的直流***状态变量作为本次迭代的初值，其中m＝1,2,…,n_W；

2)使用n_W个block模拟n_W条流水线，同时生成n_W个时步对应的雅克比矩阵和右端项，完成待求解的线性方程组；

3)求解步骤2中形成的线性方程组，计算出n_W-1个时步的直流***的状态变量；

4)检查是否收敛，如果不收敛，进入下一次的迭代。

另外，为了进一步提高本方法的效率，本发明还使用了移动时间窗的方法来避免处理器处于等待的状态。如果时间窗是固定的，第一个时间窗内包含n_w个时步，那么在第n_w个时步计算完成之后，才会进行第二个时间窗的计算。这样的任务装配方法无疑还有进一步优化的余地，移动时间窗的使用正是为了解决这样的问题。移动时间窗指的是，一旦时步1的状态变量迭代计算完毕，达到收敛，立即将原属于下一时间窗的时步n_w+1补充进现在的流水线。

由于***动态变量的剧烈变化与一个时间窗内时步的增加都会使时间并行算法的整体收敛速度下降。所以，本专利根据直流电压值，将整个仿真过程划分为故障阶段和恢复阶段。如果暂态过程中直流线路电压值低于稳态值的80％，定义为故障阶段，算法开启n_w6条流水线。如果直流线路电压超过稳态值的80％，定义为恢复阶段，开启12条流水线。

如果有多条直流线路，每条直流线路使用GPU的一个流实现，图6示即为流水线及移动窗示意图。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，所属领域的普通技术人员参照上述实施例依然可以对本发明的具体实施方式进行修改或者等同替换，这些未脱离本发明精神和范围的任何修改或者等同替换，均在申请待批的本发明的权利要求保护范围之内。

Claims (9)

1.一种大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述方法包括以下步骤：

步骤1：建立交直流***模型；

步骤2：对交直流***进行空间解耦，并对直流***进行时间解耦；

步骤3：交直流***并行仿真计算。

2.根据权利要求1所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述步骤1中，交直流***模型包括交流***动态模型、直流***模型和交直流电网模型；

所述交流***动态模型包括发电机模型、励磁***模型和原动机***模型；所述直流***模型包括直流电缆模型、直流控制***模型和换流器模型。

3.根据权利要求2所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述直流电缆模型采用T型模型，直流电缆模型对应的动态方程如下：

${\stackrel{.}{I}}_{\mathrm{dcr}}=(V_{\mathrm{dcr}}-V_c-0.5R_{\mathrm{dc}}{I}_{\mathrm{dcr}})/(L_{\mathrm{sr}}+L_{\mathrm{line}}/2)$

${\stackrel{.}{I}}_{\mathrm{dcr}}=(-V_{\mathrm{dci}}+V_c-0.5R_{\mathrm{dc}}{I}_{\mathrm{dci}})/(L_{\mathrm{si}}+L_{\mathrm{line}}/2)$

${\stackrel{.}{V}}_c=I_c/C_{\mathrm{cap}}=(I_{\mathrm{dcr}}-I_{\mathrm{dci}})/C_{\mathrm{cap}}$

其中，V_dcr为整流侧的直流电压，V_dci为逆变侧的直流电压；I_dcr为整流侧得直流电流，I_dci为逆变侧的直流电流；V_c和I_c为电容电压与电流；R_dc、L_line和C_cap分别为直流电缆电阻、电感和电容；L_sr为整流侧平滑电抗器电感，L_si为逆变侧平滑电抗器电感。

4.根据权利要求2所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述直流控制***模型从顶层到底层分为主控制级和极控制级；

所述主控制级包含紧急功率控制器，当直流***所连的交流***电压或频率发生较大波动时改变直流电流或功率参考值；

所述极控制级分为以下三个环节：

1)低压限流控制器：其输入为主控制级设定的电流参考值，在直流电压过低时，限制该设定的电流参考值，并将电流参考值输入给电流增益控制器；

2)电流增益控制器：其将低压限流控制器输入的电流参考值转换为触发角参考值，输出给熄弧角控制器；

3)熄弧角控制器：将触发角参考值转换为实际触发脉冲。

5.根据权利要求2所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述换流器模型为换流器准稳态模型，表示为：

$V_{\mathrm{dcr}}=N[\frac{3\sqrt{2}{T}_r}{\mathrm{\π}{T}_{\mathrm{tap}}}{V}_{\mathrm{adcr}}\cos \mathrm{\α}-\frac{3X_c{I}_{\mathrm{dcr}}}{\mathrm{\π}}-2R_c{I}_{\mathrm{dcr}}]$

$\cos (\mathrm{\δ}+\mathrm{\α})=\cos \mathrm{\α}-\frac{\sqrt{2}{T}_{\mathrm{tap}}{I}_{\mathrm{dcr}}{X}_c}{T_r{V}_{\mathrm{adcr}}}$

$P_{\mathrm{acr}}=N\frac{\sqrt{18}{T}_r}{2\mathrm{\π}{T}_{\mathrm{tap}}}{I}_{\mathrm{dcr}}{V}_{\mathrm{adcr}}[\cos \mathrm{\α}+\cos (\mathrm{\α}+\mathrm{\δ})]$

$I_{\mathrm{adcr}}=\frac{P_{\mathrm{dcr}}{V}_{\mathrm{adcrx}}+Q_{\mathrm{dcr}}{V}_{\mathrm{adcry}}}{V_{\mathrm{adcrx}}^2+V_{\mathrm{adcry}}^2}+j\frac{P_{\mathrm{dcr}}{V}_{\mathrm{adcry}}-Q_{\mathrm{dcr}}{V}_{\mathrm{adcrx}}}{V_{\mathrm{adcrx}}^2+V_{\mathrm{adcry}}^2}$

其中，V_dcr和I_dcr分别为整流侧的直流电压和直流电流，N为串联整流桥个数，T_tap为整流侧变压器分接头位置，I_adcr为整流侧注入交流***电流，V_adcr为整流侧变压器交流侧电压，V_adcrx和V_adcry分别为整流侧变压器交流侧电压实部和虚部，X_c和R_c分别为整流侧变压器电抗和电阻，α为整流侧触发角，δ为换相重叠角，T_r为整流侧变压器基本变比，P_acr和Q_acr分别为整流侧吸收的有功功率和无功功率，为功率因数角。

6.根据权利要求1所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述步骤2中，对建立的交直流***进行空间解耦具体过程如下：

交流***动态模型表示为：

${\stackrel{.}{x}}_{\mathrm{ac}}=f(x_{\mathrm{ac}},V_{\mathrm{ac}})$

其中，x_ac和V_ac分别为交流***状态变量和纯交流节点电压向量；

直流***模型表示为：

${\stackrel{.}{x}}_{\mathrm{dc}}=H(x_{\mathrm{dc}},V_{\mathrm{adc}})$

其中，x_dc和V_adc分别为直流***状态变量及直流***接入交流***处的节点电压向量；

交直流电网模型表示为：

0＝I-Y_NV

其中，I为交直流***节点注入电流，其包括V为交直流***节点电压，Y_N为网络导纳矩阵；

由于交流***状态变量和直流***状态变量相互独立，且通过网络方程相关联，按照交流***和直流***空间上的位置关系进行解耦，直流***采用诺顿等效接入交直流***电网网络，在CPU部署交流***部分计算数据，即交流***动态方程和网络方程中变量初始值及参数值，GPU部署直流***部分计算数据，即直流***动态方程中变量初始值及参数值。

7.根据权利要求1所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述步骤2中，对直流***进行时间解耦，具体过程如下：

对直流***动态方程略去下标，可得直流微分方程，有：

$\stackrel{.}{x}=H(x,V)$

在交流积分步长h_ac中，直流微分方程采用单时步顺序求解的方法，计算h_ac/h_dc步，其中h_dc为直流积分步长；使用隐式梯形积分法，设下标n表示t时刻，n+1表示t+h_ac时刻，直流微分方程的差分化方程为：

$x_n=x_{n-1}+\frac{h_{\mathrm{dc}}}{2}[H(x_{n-1},V_{n-1})+H(x_n,V_n)],n=1,...,h_{\mathrm{ac}}/h_{\mathrm{dc}}$

其中，x_n为直流***t时刻的状态变量，x_n-1为直流***t-h_ac时刻的状态变量，V_n为直流***t接入交流***处的节点电压向量，V_n-1为直流***t-h_ac接入交流***处的节点电压向量；

定义中间向量R_n，其表示为：

$R_n=x_n-x_{n-1}-\frac{h_{\mathrm{dc}}}{2}[H(x_{n-1},V_{n-1})+H(x_n,V_n)]$

对上式采用牛顿-拉夫逊法，第k次迭代公式为：

$R_n^{k-1}=-\mathrm{J\Δ}{x}_n^k$

其中，k为迭代次数，为两次迭代中直流***状态变量差值，J为牛顿-拉夫逊法中第k次迭代的雅克比矩阵，和J分别表示为：

$\mathrm{\Δ}{x}_n^k=x_n^k-x_n^{k-1}$

$J={\frac{\∂R_n}{\∂x_n}|}_{x_n=x_n^{k-1}}$

其中，为第k次迭代中直流***t时刻的状态变量，为第k-1次迭代中直流***t时刻的状态变量；

于是得到其中V₀为上一交流时步所得电压值。

8.根据权利要求1所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述步骤3具体包括以下步骤：

步骤3-1：交直流***进行初始化，采用OpenMP多线程在CPU上开启两个线程，分别为线程A和线程B，线程A在CPU上完成交流***相关的计算，线程B用于控制GPU，向GPU分配任务、发送控制指令；然后交直流***进入时域仿真过程，按照时间的推进，计算每个时刻对应的交直流***状态变量，设置仿真时刻t＝0，设推进步长为一个交流时步h_ac；

步骤3-2：交流***和直流***间交换边界点电压和电流值；

步骤3-3：采用交替迭代的隐式梯形积分方法，在异构平台中的CPU上计算交流***状态变量；

步骤3-4：在异构平台的GPU上完成多个直流***状态变量的计算，具体包括以下步骤：

步骤3-4-1：判断是否有多条直流线路，若是，则开启多个GPU流，每个流计算一条直流线路；若否，则只使用一个GPU流；

步骤3-4-2：每个直流线路使用结合移动窗的时间流水线并行方法计算，并结合电压值，判断流水线的条数；

步骤3-4-3：计算收敛后，计算直流***注入交流***的电流；

步骤3-5：当获得全部直流线路注入电流之后，设置仿真时刻t＝t+h_ac，进入下一时刻；

步骤3-6：重复步骤3-1至步骤3-4，直到仿真结束。

9.根据权利要求8所述的大规模电力***暂态稳定时间空间并行仿真方法，其特征在于：所述步骤3-4-2中，将本仿真时段划分为多个时间窗口，而每个时间窗口包含n_w个时步，多个时间窗口按顺序计算，但是在每个时间窗口内，将多个时步的计算任务装配为流水线的形式，实现并行计算；

对于时间窗口1，n_w个时步对应的差分方程为：

$\left[\begin{array}{c}{R}_1\\ R_2\\ .\\ .\\ .\\ R_{n_w}\end{array}\right]=\left[\begin{array}{c}{x}_1-x_0-\frac{h_{\mathrm{dc}}}{2}[H(x_1,V_1)+H(x_0,V_0)]\\ x_2-x_1-\frac{h_{\mathrm{dc}}}{2}[H(x_2,V_2)+H(x_1,V_1)]\\ .\\ .\\ .\\ x_{n_w}-x_{n_w-1}-\frac{h_{\mathrm{dc}}}{2}[H(x_{n_w},V_{n_w})+H(x_{n_w-1},V_{n_w-1})]\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ .\\ .\\ .\\ 0\end{array}\right]$

对上式采用牛顿-拉夫逊法，第k次迭代公式为：

其中： $J_{n,n-1}^{k-1}={\frac{\∂R_n}{\∂x_{n-1}}|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^{k-1}}},J_{n,n}^{k-1}={\frac{\∂R_n}{\∂x_n}|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^{k-1}}},n=1,...,n_w;$

对于时间窗口中的n_w个时步，在第一次迭代中，流水线P1进行时步1的第一次迭代求解，获得时步1的状态变量在第二次迭代中，流水线P2进行时步2的第一次迭代求解以及时步1的第二次迭代求解；依次类推，实现对当前时间窗口中全部时步的并行求解；求解公式为：

$\mathrm{\Δ}{x}_n^k=\left\{\begin{array}{c}-{\left(J_{n,n}^{*k-1}\right)}^{-1}{R}_n^{k-1},n=1\\ -{\left(J_{n,n}^{*k-1}\right)}^{-1}[R_n^{k-1}+J_{n,n-1}^{*k}\mathrm{\Δ}{x}_{n-1}^k],n=2,...,n_w\end{array}\right.$

其中， $J_{n,n-1}^{*k-1}={\frac{\∂R_n}{\∂x_{n-1}}|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^k}},J_{n,n}^{*k-1}={\frac{\∂R_n}{\∂x_{n-1}}|}_{\underset{x_n=x_n^{k-1}}{x_{n-1}=x_{n-1}^k}}.$