Deep Reinforcement Learning for Voltage Control and Renewable Accommodation Using Spatial-Temporal Graph Information

A DN can be modeled as an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\bm{A})$, as illustrated in Fig. 2, where $\mathcal{V}$, $\mathcal{E}$, and $\bm{A}$ represent the node set, edge set, and the adjacency matrix, respectively. Each load node $v_{n}$ generates a feature vector $\bm{x}_{t,n}$ for gathering local information, e.g., load voltage ($\mathbb{V}_{t,n}^{\text{L}}=v_{t,n}^{\text{L}}\angle\delta_{t,n}^{\text{L}}$), load power ($p_{t,n}^{\text{L}},q_{t,n}^{\text{L}}$), and connected branch power ($\mathbb{S}_{t,b}^{\text{B}}=p_{t,b}^{\text{B}}+jq_{t,b}^{\text{B}}$), at each time step, which can be formulated as

where $F_{n}$ is the dimension of node features. The feature matrix of graph $\mathcal{G}$ can be aggregated as

where $F^{\prime}$ is the highest node feature dimension defined as $F^{\prime}=\max_{n}F_{n}$.

Considering periodic patterns in the power flow, especially the daily and weekly ones [2], we develop a multi-grained vector constructor to better capture temporal correlations of the DN, where the recent, daily, and weekly graph segments can be formulated as

where $T^{\text{r}}$, $T^{\text{d}}$, and $T^{\text{w}}$ indicate the length of recent, daily, and weekly segments, respectively, and $N^{\text{d}}$ represents the number of resolving the optimization problem per day.

The framework of the proposed MG-ASTGCN is illustrated in Fig. 3, consisting of graph conversion operation, multi-grained vector constructor, and the most essential ASTGCN which takes multi-grained segments as inputs and employs stacked ST components to extract ST information. The structure of one ST component is depicted in Fig. 4, including ST attention mechanism and ST convolution, which are presented in detail in Section III-A2 and III-A3, respectively.

The key idea of the ST attention mechanism is to pay more attention to valuable graphical information in both spatial and temporal perspectives, assisting its sequential ST convolution for extracting more useful features for the DRL.

Spatial Attention: Mutual influences of each neighboring node pair in the DN dynamically vary due to the changes in power flow. To explore such mutual influences, an attention mechanism in spatial dimension is thus developed to capture the dynamic correlations [31], which can be formulated as (we omit the notations of multi-grained segments in the rest of this section for brevity)

where $\odot$ represents element-wise multiplication, $\sigma(\cdot)$ is the sigmoid activation function, $\bm{V}_{s}$, $\bm{W}_{T}$, $\bm{W}_{FT}$, $\bm{W}_{F}$, and $\bm{b}_{s}$ are all learnable parameters, and $\bm{S}$ represents the spatial attention matrix, whose element $s_{i,j}$, namely the attention weight, semantically describes the correlation strength between the $i$-th and $j$-th nodes. The derived spatial attention matrix can be adopted for spatial graph convolution to adjust spatial correlation strengths of node pairs, as shown in Fig. 4.

Temporal Attention: Similar to the spatial attention, the temporal attention mechanism [32] aims to track temporal correlations of the changing node features, which can be formulated as

where $\bm{V}_{e}$, $\bm{U}_{N}$, $\bm{U}_{FN}$, $\bm{U}_{F}$, and $\bm{b}_{e}$ are learnable parameters, $\bm{\mathcal{X}}_{\text{tatt}}$ is the transposed form of the input segment $\bm{\mathcal{X}}$, and $\bm{E}$ is referred to as the temporal attention matrix whose attention weight $e_{i,j}$ represents the temporal dependencies between two graph feature vectors $\bm{X}_{t-i}$ and $\bm{X}_{t-j}$. The temporal attention matrix is used to add temporal correlation information to the original input segment $\bm{\mathcal{X}}$ as shown in Fig. 4, which can be expressed as $\tilde{\bm{\mathcal{X}}}=\bm{\mathcal{X}}\bm{E}$.

The ST convolution consists of spatial graph convolution and temporal convolution, aiming to compress the three-dimensional input segments and extract ST features that can be integrated with the DRL algorithm.

Spatial Graph Convolution: Graph convolution is defined as a convolution operation implemented by using linear operators diagonalizing in the Fourier domain to replace the classical convolution operator [33], which can be expressed as

where $*_{G}$ represents the graph convolution operator, $g_{\theta}$ is a convolution filter, $\bm{L}$ is the Laplacian matrix of the derived undirected graph $\mathcal{G}$, $\bm{\Lambda}$ is the result of eigenvalue decomposition, and the rectified linear unit (ReLU) is adopted as the activation function. Providing that conducting eigenvalue decomposition on the Laplacian matrix is computationally expensive, Chebyshev polynomials are often used for approximating the solution of eigenvalue decomposition in practice [34]. With the addition of spatial correlation information from the spatial attention matrix $\bm{S}$ as shown in Fig. 4, the spatial graph convolution can be rewritten as

where $\tilde{\bm{L}}$ is the normalized Laplacian matrix, $\theta_{k}$ is the coefficient of Chebyshev polynomials, $K$ is the highest order of Chebyshev polynomials, and $T_{k}$ is the $k$-th order Chebyshev polynomial.

Temporal Convolution and Feature Compression: The temporal convolution operation takes the result of the spatial graph convolution as input and performs convolution along the temporal dimension, i.e., $T^{\text{r}}$, $T^{\text{d}}$, and $T^{\text{w}}$, which can be formulated as

where $h_{\theta}$ is the temporal convolution filter and $\bm{\mathcal{X}}^{\text{next}}$ is the input for the following ST component.

To pass the extracted ST information to the DRL algorithm, the multi-time-scale outputs are fused and then fed into a fully-connected neural network layer (FCNNL) [35] for feature compression as shown in Fig. 3, which can be formulated as

where concat represents the concatenating operation for the multi-time-scale outputs.

In summary, ST graphical information can be fully exploited by our MG-ASTGCN. The spatial attention mechanism explores spatial correlations of neighboring node pairs, while the temporal attention focuses on mining self-correlations of node features in the temporal view. Furthermore, based on ST correlations provided by the preceding ST attention, ST convolution extracts ST features of the underlying DN. The detailed process of the recent graph segment $\bm{\mathcal{X}}^{\text{r}}$ passing through one ST component is illustrated in Fig. 5.

Balancing the trade-off between voltage control and renewable accommodation in the DN can be considered as a consecutive decision-making process, which can be further modeled as an MDP [36] consisting of four critical parts: state space $\mathcal{S}$, action space $\mathcal{A}$, probability space $\mathcal{P}$, and reward space $\mathcal{R}$.

State Space $\mathcal{S}$: The state of the $n$-th load node, denoted by $\bm{s}_{t,n}$, is its feature vector $\bm{x}_{t,n}$ defined in Eq. (26). Moreover, the extracted ST features $\bm{y}_{t}$ is aggregated with the DN’s state. Thus, the state of the DN can be expressed as

Action Space $\mathcal{A}$: For each power generator, only its active power $p_{t,i}$ and voltage magnitude $v_{t,i}$ can be manipulated. Thus, the actions of the $i$-th power generator can be expressed as $\bm{a}_{t,i}=\left[p_{t,i},v_{t,i}\right]$. Actions of all generators in the DN can be defined as

where $\bm{a}_{t,i}^{\text{T}},\bm{a}_{t,j}^{\text{W}},\bm{a}_{t,k}^{\text{S}}$ are actions of thermoelectric, wind, and solar PV generators. Note that, our proposed strategy is also applicable to generators operating in the P/Q control mode by changing the generator’s action space into active power $p_{t,i}$ and reactive power $q_{t,i}$.

We would like to clarify that our proposed DRL-based strategy can be easily applicable to generators operating in the P/Q control mode by slightly modifying the generator’s action space. In our original manuscript, the action space of each generator is defined as $\bm{a}_{t,i}=[p_{t,i},v_{t,i}]$, where $p_{t,i}$ and $v_{t,i}$ are its active power and voltage magnitude, respectively. The generator working in the P/Q control mode can be defined as $\bm{a}_{t,i}=[p_{t,i},q_{t,i}]$, where $q_{t,i}$ is the reactive power of the generator.

Probability Space $\mathcal{P}$: $\mathcal{P}$ is the probability set of transitioning to the next state $\bm{s}_{t+1,n}$ from the current state $\bm{s}_{t}$ after taking a deterministic action $\bm{a}_{t}$.

Reward Space $\mathcal{R}$: A reward $r_{t}$ is obtained after taking action $\bm{a}_{t}$ at state $\bm{s}_{t}$, which indicates the effectiveness of the selected action. The goal of DRL is to learn an optimal action strategy $\pi(\bm{a}_{t}|\bm{s}_{t})$ to maximize the expected cumulative rewards. Hence, it is essential to encode the objective of the optimization problem into a reward function to facilitate the DRL training. The reward function of the MDP can be formulated as

where the above three reward terms correspond to the objective functions of the optimization problem, i.e., voltage control $J^{\text{vol}}$, renewable accommodation $J^{\text{RER}}$, and generation cost minimization $J^{\text{gen}}$.

We then introduce how we handle the optimization constraints defined from Eq. (10) to (25) and encode them as rewards into the DRL training process. For the power balance constraints in Eq. (10) and (11), if these two constraints cannot be satisfied during the DRL training process, a constant penalty term ($-10$ in our simulation by default) will be added to the reward for constraint violation. Moreover, the current training episode will be terminated. Subsequently, the DN environment will be reset and the input MDP state for the next episode will be randomly initialized. For the load voltage constraint, generator reactive power constraints, and branch flow constraint presented in Eq. (14), (22)-(24), and (25), respectively, if these constraints are violated, the same constant penalty term ($-10$ by default) is added to the reward for violation feedback. For the generator’s voltage and power constraints defined in Eq. (15)-(17) and (18)-(21), respectively, their corresponding values will not violate these inequality constraints. This is because the actions of each generator are its power and voltage defined as $\bm{a}_{t,i}=[p_{t,i},v_{t,i}]$, indicating that the lower and upper limits of these constraints are encoded as the bounds of the MDP’s action space. As a result, the action value of each generator will always meet these constraints.

The objective of DRL is to maximize the expected cumulative rewards, denoted by $\bar{R}_{\theta}$, which can be formulated as

where $\tau$ represents the trajectory of the MDP transitions, recording all 4-tuple transitions denoted by $\{\bm{s}_{t},\bm{a}_{t},r_{t},\bm{s}_{t+1}\}$ from the beginning of $\tau$ to its end, $R(\tau)$ represents the cumulative rewards of the trajectory, $\theta$ represents the parameters of our action strategy $\pi$, $P\left(\tau|\theta\right)$ is the occurrence probability of the trajectory $\tau$.

We then introduce DDPG [37] to maximize $\bar{R}_{\theta}$. DDPG is the most representative actor-critic DRL algorithm to optimize the derived MDP. The major difference between DDPG and other DRL algorithms is that the action policy $\pi(\bm{a}_{t}|\bm{s}_{t})$ deterministically outputs the values of actions instead of the probability distribution of actions, which dramatically decreases the computation cost and makes it much easier to implement. Specifically, policy gradient method is applied in DDPG to update our action policy, which can be formulated as

with the gradient of the DRL objective defined as

where $\eta_{\theta}$ is the learning rate, $N_{\tau}$ is the number of trajectories, $T_{\tau}$ is the length of each trajectory, and $\mathcal{A}_{\theta}\left(\bm{s},\bm{a}\right)$ represents the advantage function assessing the effectiveness of the state-action pair compared to a certain baseline, which can be formulated as

where $\gamma$ is a discounting factor, $b$ represents the baseline of reward considering all possible actions, and $\mathbb{E}_{s\sim\mathbb{P}}\left[R(\bm{s})\right]$ is the expected cumulative rewards considering all possible states at the current time step. Due to the uncertainty of the MDP transitions, both $G_{t}$ and $b$ are random variables. To accurately estimate the advantage function, DDPG introduces a critic network $Q_{\phi}(\bm{s}_{t},\bm{a}_{t})$ formulated as

where $\phi$ represents parameters of the critic networks. With the critic network, the advantage function in Eq. (45) can be rewritten as

Providing that calculating the output value of the critic network depends on environments rather than our action policy $\pi$, it is applicable to learn the critic network $Q_{\phi}$ in an off-policy manner taking advantage of transitions generated from a different action policy denoted by $\pi_{\theta^{\prime}}$. The critic network can be updated by minimizing the root mean square error formulated as

with the estimated critic network via the target action and critic networks $\pi_{\theta^{\prime}},Q_{\phi^{\prime}}$ defined as

Similarly, the critic network is updated via gradient descent formulated as

where $\eta_{\phi}$ is the learning rate for updating the critic network.

The proposed two target networks are updated in an exponential moving average manner using parameters of their corresponding action and critic networks, which can be formulated as

where $\rho$ is the smoothing parameter.

In summary, DDPG follows the policy-gradient way to optimize the derived MDP, where a critic network is introduced to assess the learned action policy. The workflow illustration and algorithmic procedure of the DDPG are presented in Fig. 6 and Algorithm 1, respectively.

The proposed DRL-based strategy is tested on the modified IEEE 33-bus, 69-bus, and 118-bus radial distribution systems (RDSs) [38] with their detailed statistics provided in Table I. The time interval between two consecutively operational time steps is one hour, aligning with previous studies in DN’s voltage control [22, 39, 25, 40]. The simulation lasts for a one-year duration to ensure sufficient DRL training. One Nvidia TITAN RTX graphics processing unit is used for DRL training. The batch size of the DDPG algorithm is set as $256$. The random noise employed in the DDPG is a Gaussian noise, denoted by $\mathcal{N}\sim(\mu,\sigma)$ with its mean $\mu$ and standard deviation $\sigma$ setting as $0$ and $0.1$, respectively. The number of hidden layers for the actor and critic networks is $3$, each of which consists of $400$ neurons. We adopt the Adam optimizer to update our actor, i.e., $\pi_{\theta}$ and critic, i.e., $Q_{\phi}$, networks. The number of ST components inside our MG-ASTGCN is $3$. The initialized parameters of our DRL-based control strategy are provided in Table II.

Note that we use the same reserve cost coefficients and penalty cost coefficients for wind and solar generators provided by [41], i.e., $c_{j}^{\text{W,r}}=c_{k}^{\text{S,r}}=1.5$ and $c_{j}^{\text{W,p}}=c_{k}^{\text{S,p}}=3$. For the cost coefficients of the thermoelectric generator, we use the same coefficients provided by [28], i.e., $a_{i}=0.0175$, $b_{i}=1.75$, and $c_{i}=0$. Additionally, we provide the power output profile of one wind power generator and one solar generator throughout one day in a part of the IEEE $118$-bus test system in Fig. 7.

In our experiments, the reward function defined in Eq. (41) is adopted to measure the performance of both DDPG and benchmark algorithms. Specifically, we introduce three metrics, namely SCORE, voltage fluctuation rate denoted by $\alpha^{\text{vol}}$, and RER accommodation rate denoted by $\alpha^{\text{RER}}$, to examine the effectiveness of the learned control strategy, voltage control, and renewable accommodation, respectively, which are defined as

where $N^{\text{eval}}$ is the number of episodes for evaluation and $T^{\text{end}}$ is the length of each episode. Both $N_{\text{eval}}$ and $T_{\text{end}}$ are initialized as $100$.

Two representative optimization-based algorithms—HHO [19] and GWO [20], are adopted to compare with our DRL-based strategy. Benefiting from their meta-heuristic characteristic, both the HHO and GWO algorithms can approach the optimal control strategy without modifying or relaxing the optimization formulation. Moreover, the meta-heuristic-based methods also present good convergence and robustness abilities in the previous studies [41, 42, 43]. Moreover, we also introduced two relaxation-based optimization methods, namely interior-point (IP)-based method [44] and linear/quadratic programming (LQP)-based method [45] as benchmarks. The evaluation results derived by these methods on IEEE 33, 69, and 118-bus RDSs are illustrated in Fig. 8 with associated statistics presented in Table III. The results reveal that our proposed DRL-based strategy outperforms all benchmarks by significant margins in terms of faster testing time, higher SCOREs, lower voltage fluctuation rates, and more renewable accommodation. The reason for such performance gaps is twofold:

• •

  The optimization-based benchmarks do not need training, resulting in longer computation time for evaluation at each time step. Meanwhile, since they tend to get stuck in local optimums, their output actions are less likely to be the optimal ones, resulting in smaller rewards.

• •

  Tremendous data collected by the DRL-based strategy contributes to its more effective and efficient searching in the action space. Therefore, the DRL reacts much faster at the beginning of evaluation and gradually obtains a higher reward, especially in large-scale DNs.