Figure 1 shows the conventional AGC dynamic optimization scheme for optimizing the regulated power of AGC units in steps of 1 min over a duration of 15 min based on the deviation information for system frequency, ACE, and tie-line power. Generation units in a power grid are of two types, namely AGC and non-AGC units, whose power outputs are denoted as P i AG and P g NG , respectively. Here, i and g are both the indices of the AGC and non-AGC units. In this work, the AGC units participate in both primary and secondary frequency control, while the non-AGC units only provide primary frequency control. Hence, the dynamic power outputs can be calculated as follows: P g , t + 1 NG = P g , t NG − K g Δ f t + 1 − Δ f t , (1) P i , t + 1 AG = P i , t AG − K i Δ f t + 1 − Δ f t + Δ P i , t AGC , (2) where Δ f t + 1 and Δ f t are the frequency deviations at time t + 1 and t , respectively; Δ P i , t AGC is the power adjustment of the AGC unit i at time t for secondary frequency control; K i and K g are the frequency regulation constants of the AGC unit i and non-AGC unit g , respectively. Here, Equation 2 includes two parts, which are the primary frequency control power K i ( Δ f t + 1 − Δ f t ) of the AGC unit i and power increment of the AGC unit Δ P i , t AGC at each optimization time. Thus, the AGC units respond to the uncertain power fluctuations in the grid through power regulation Δ P i , t AGC at time t .

From Figure 1, we see that the AGC strategy requires three dynamic system parameters as inputs: frequency Δ f t , ACE e t ACE , and tie-line power P t tie . Hence, their system dynamics are as follows: d ( Δ f t ) d t = − D 2 H Δ f t + 1 2 H ( Δ P t wind + Δ P t tie + Δ P t gen + Δ P t d ) (3) e t ACE = Δ P t tie + 10 B ⋅ Δ f t , (4) d Δ P t tie d t = ∑ j K s j ⋅ Δ f t − Δ f j , t , (5) where H is the equivalent inertial constant; D is the equivalent damping coefficient; Δ P t wind , Δ P t gen and Δ P t d are power adjustments of wind, generation, and demands at time t; B is the frequency regulation constant of the control system in megawatts per 0.1 Hz (positive value); K s j and Δ f j , t are the tie-line synchronization coefficient and frequency deviation of the connected j -region, respectively; P t tie is the tie-line power outflow that is considered to be positive; Δ P t tie is the tie-line power deviation. Equations 1–5 describe the system dynamics of key variables.

In the present work, our control objective is to schedule the AGC units so as to satisfy both the minimum economic cost of auxiliary services as well as stability and safety of the CPS. Therefore, the objective can be expressed mathematically as follows: min f 1 = ∑ i ∈ ACG units c i | P i , t AG − P i , 0 AG | + u i , t R i , t AG , (6) where c i is the auxiliary service cost coefficient of the AGC unit i ; R i , t AG is the climbing power at time t of the AGC unit i ; u i , t ∈ − 1,0,1 indicates the change in the direction of the power output, where u i , t = − 1 denotes a decrease, u i , t = 1 denotes an increase, and u i , t = 0 denotes that there is no change.

The assessment indexes of the CPS include CPS1 and CPS2. Here, CPS1 is defined to evaluate the correlation between the system frequency deviation and ACE, while CPS2 is defined as the average ACE over 15 min, indicating that the ACE is maintained within a tolerance range to ensure that the power exchanged between the regions does not exceed the specified limits. The detailed definitions of CPS1 and CPS2 are as follows: K cps 1 = ∑ t = 1 T e t ACE | T | ⋅ ∑ t = 1 T f t − f ref | T | ⋅ 1 10 B ϵ 1 2 , (7) K cps 2 = | 1 | T | ∑ t = 1 T Δ P t tie + 10 B ⋅ Δ f t | , (8) where T is the time interval; f ref is the rated frequency; ϵ 1 represents the frequency control target, which is usually taken as the root mean-squared (RMS) value of the mean frequency deviation of 1 min in the previous year. In practice, when CPS1 satisfies the condition K cps 1 ≥ 2 , CPS2 will not be assessed (Wu et al., 2010); otherwise, CPS2 must be less than the following threshold: K cps 2 ≤ 1.65 ϵ 15 m i n 100 B B s , (9) where ϵ 15 m i n is the RMS value of the mean frequency deviation over 15 min in the previous year; B s is the equivalent frequency regulation constant for the entire interconnection power grid. If CPS1 is suitably satisfied in this work, CPS2 will not be considered an objective or a hard constraint.

The operational constraints for AGC dynamic optimization include system power balance, AGC unit regulation characteristics, and limits for the frequency and tie-line power deviations. The basic power balance constraint must satisfy the condition ∑ i P i , t AG + ∑ g P g , t NG + P t wind − P t D − P t tie − P t loss = 0 , (10) where P t wind and P t D are the forecast values of the wind power and loads for period t ; P t loss is the line transmission loss.

As shown in Figure 1, the saturation function and ramp rate limiter will take effect on the control signals before being executed. Specifically, the power output and ramp power constraints of AGC units are defined as follows: P ̲ i AG ≤ P i , t AG ≤ P ̄ i AG , (11) R ̲ i AG ≤ R i , t AG ≤ R ̄ i AG , (12) where P ̄ i AG and P ̲ i AG are the upper and lower limits of the output power of the i -th AGC unit, respectively; R ̄ i AG and R ̲ i AG are the upper and lower limits of the ramp power of the i -th AGC unit, respectively. Moreover, the limitations on the frequency deviation and tie-line power are as follows: Δ f ̲ ≤ Δ f t ≤ Δ f ̄ , (13) P ̲ tie ≤ P t tie ≤ P ̄ tie , (14) where Δ f ̄ and Δ f ̲ are the corresponding upper and lower limits of the system frequency deviation, respectively; P ̄ tie and R ̲ tie are upper and lower limits of the tie-line power, respectively.

In this work, the MDP is a mathematical framework used to model the AGC dynamic optimization problem as a sequential decision-making process, as shown in Figure 2. The MDP is defined by a tuple S,A,P,R , γ , where S is a set of states, A is a set of actions, P is a transition model that gives the probability of moving from one state to another when given an action, R is a reward function that provides the immediate reward for state transitions, and γ ∈ [ 0,1 ] is a discount factor that determines the importance of future rewards ¹ . In the MDP, the AGC strategy is considered an agent that observes the grid operating state s t ∈ S and outputs an action a t ∈ A at each time step t . Then, the agent will receive an immediate reward r t ∈ R . The sequential continuous experience τ is then recorded as s 0 , a 0 , s 1 , … , a T − 1 , s T . The objective of the agent is to maximize its cumulative reward J π by iteratively updating its policy π : S → A , which is given by max π J π = E τ ∼ π ∑ t = 0 T γ t r t . (15)

Based on the objectives and constraints defined in Equations 6–14, we further present a well-designed MDP formulation. In this work, the control variables are the regulation direction and regulation power of each AGC unit. To simplify the action space scale using smaller action dimensions, the action is defined as the power adjustment of the AGC unit: a t = Δ P i , t AG | i ∈ I ⊺ , (16) where I denotes the set of AGC units. Before execution, all actions are subjected to the saturation and ramp limits shown in Figure 1.

The design of the state space must capture necessary information based on two aspects: (1) conditions of the current operating system; (2) uncertain environments that must be forecast. For the former, we take into account four factors, including the current power outputs of the AGC units P i , t AG , frequency deviation Δ f t , tie-line power deviation Δ P t tie , and ACE value e t ACE . For uncertain environments, we only consider wind power forecasting in this work because the load fluctuations within 15 min are usually negligible. Hence, we introduce historical wind power information to the state space for better forecasting. The system state is defined as s t = P i , t AG | i ∈ I , Δ f t , Δ P t tie , e t ACE , Δ P t − k wind | k = 1,2 , … , l , ⊺ , (17) where l is the period of historical observations; Δ P t wind = P t wind − P t − 1 wind denotes the difference between successive wind power outputs.

The reward design for the MDP should consider both objectives and constraints. Here, we design the reward function based on three aspects as follows: r t = w 1 f 1 + w 2 | K cps 1 − 2 | 2 + w 3 r penalty , (18) where w 1 , w 2 , and w 3 are weight factors that balance the tradeoffs between the three subreward items; f 1 is the economy objective defined in Equation 6. Note that the presence of too many items will complicate the design of the corresponding weights and lead to convergence failure. Hence, we design the third item r penalty as a penalty for the total violations of the upper/lower limitations, such as those of the output power, ramp power, tie-line power, and frequency in Equations 11–14. The penalty term is formulated as r penalty = ∑ | k | ReLU c k − c ̄ k , (19) where k represents the index of the constraints defined in Equations 11–14, i.e., | k | = 2 | I | + 2 ; ReLU ( x ) = max ( 0 , x ) is a linear rectification function for measuring the violations; c k and c ̄ k represent the actual index value and required limit value, respectively. Taking the constraint in Equation 13 as an example, c k and c ̄ k can be defined as c k = { Δ f t , − Δ f t } and c ̄ k = { Δ f ̄ , − Δ f ̲ } , respectively; here, we separate the original Equation 13 into two inequalities as Δ f t ≤ Δ f ̄ and − Δ f t ≤ − Δ f ̲ , which are expressed using the same structure. Hence, Equations 16–19 show our design of four key fators in MDP formulations.

Behavioral cloning (BC) is a common method for implementing imitation learning (Daftry et al., 2017), where the demonstrator (i.e., expert experiences) can be imitated directly without interacting with a real-world environment. The key idea of BC is to replicate the expert policy using a classifier or regressor based on previously collected training data from the encountered states and demonstrator actions (Rajaraman et al., 2020). Therefore, BC-based imitation learning can be used in an MDP framework without defining a reward function. The learning objective for the agent here is to obtain an imitation policy as the initial policy, i.e., π 0 : S → A , which is necessary to behave like an expert. Here, we adopt the variable ξ = { ( s 0 , a 0 ) , ( s 1 , a 1 ) , … , ( s T , a T ) } to represent the expert demonstrations.

Given the precollected set of state–action pairs ( s i , a i ) , the objective of the agent is to seek an imitation policy π 0 ( ϕ ) that best matches the provided set of state–action pairs. The policy network parameter ϕ is then updated using maximum-likelihood estimation, i.e., the optimal ϕ * is defined as ϕ * = arg max ϕ Π t = 0 T π ϕ a i | s i . (20) Considering that the designed action space is continuous, we assume that the policy follows a Gaussian distribution over each action dimension. In this work, we adopt a neural network to approximate the policy π and use the same network structure as that of the actor in the RL framework. Then, the Adam stochastic gradient descent optimizer is adopted to solve for ϕ * in Equation 20, where the gradient descent approach aims to find changes in ϕ that can increase the accuracy of each imitated demonstrator action a i based on the imitation policy π ϕ ( ⋅ | s i ) . The pseudocode for this process is summarized in Algorithm 1.

Algorithm 1

Behavioral-cloning-based imitation learning.

1: Initialize π ϕ as a random policy network

To solve the optimal policy π * in Equation 15, we adopt the SAC algorithm in this work to maximize the agent’s cumulative rewards while satisfying the safety constraints. Compared with the conventional RL algorithm, SAC uses a stochastic policy that inherently encourages exploration by adding entropy to the reward. Hence, SAC is less likely to be stuck in local optima and can better explore the action space. Moreover, incorporating policy learning with entropy regularization helps the agent to become more stable during training. The entropy term prevents the policy from becoming too deterministic too early, thereby providing a more balanced and robust learning process. The SAC objective is given as J π = ∑ t = 0 T E τ ∼ π r t + α H π ⋅ | s t , (21) where α is the hyperparameter that balances the importance of the entropy term H(·) with the reward r t . The entropy is calculated as follows: H π ⋅ | s t = − E a t ∼ π ⁡ log ⁡ π a t | s t = − log ⁡ π ⋅ | s t . (22)

For a given policy π , the state–action value function Q : S × A → R is defined to evaluate the expected values of the pair ( s t , a t ) at time step t to guide policy learning and optimization. Generally, a larger Q value indicates better policy control performance. The Q function is defined using the Bellman equation as follows: Q s t , a t ≜ r s t , a t + γ E s t + 1 ∼ ρ , a t ∼ π Q s t , a t − log ⁡ π a t | s t , (23) To improve the stability and accuracy of the estimated values, the SAC algorithm incorporates two Q networks to mitigate the overestimation bias that can occur in the Q learning process. We use Q θ 1 and Q θ 2 to denote the two Q networks with parameters θ 1 and θ 2 separately. Then, each Q network is updated by minimizing the mean-squared error (MSE) between the current and target Q values. The loss for each Q network is defined as J Q θ i = E τ ∼ π Q θ i s t , a t − y t 2 , (24) y t = r t + γ E a t + 1 ∼ π min Q θ 1 ′ s t + 1 , a t + 1 , Q θ 2 ′ s t + 1 , a t + 1 − α ⁡ log ⁡ π a t + 1 | s t + 1 (25) where θ 1 ′ and θ 2 ′ are the parameters of the two target Q networks Q θ 1 and Q θ 2 , respectively. For each iteration, the parameters of each Q network are updated using gradient descents computed from the following loss functions: θ i ← θ i − λ Q ∇ θ i J Q θ i , (26) θ i ′ ← β θ i + 1 − β θ i ′ , (27) where λ Q is the learning rate of the Q network; β is a hyperparameter that controls the update rate of each target network based on the moving average value. From Equations 20–27, we can effectively achive the policy iteration update in the SAC algotithm.

For the policy network π , the newly updated policy at each iteration is improved using the information projection defined in terms of the Kullback–Leibler (KL) divergence D K L . Specifically, for each state, the policy is updated as follows: π new = arg min π ′ D K L π ′ ⋅ | s t ‖ exp 1 α Q π old s t , ⋅ Z π old s t , (28) where Z π old ( s t ) is the normalization item that does not influence the policy gradient calculation. Based on the projection, the new policy π new has a higher value than the old one (Haarnoja et al., 2018), subject to the maximum entropy objective. Furthermore, we can rewrite the gradient of the stochastic policy π θ using a noise vector ϵ t , which is added to the action as a t = f ϕ ( ϵ t ; s t ) . The expected KL divergence in Equation 28 can be rewritten as J π ϕ = E s t ∼ D , ϵ t ∼ N log π ϕ f ϕ ϵ t ; s t | s t − Q θ s t , f ϕ ϵ t ; s t , (29) where π ϕ is defined implicitly in terms of f ϕ . Then, the policy gradient π for Equation 29 is approximated as ▽ ϕ ⁢ J π ⁢ ϕ = ▽ a t ⁡ log ⁢ π ϕ ⁢ a t | s t − ▽ a t ⁢ Q ⁡ s t , a t ⁢ ▽ ϕ ⁢ f ϕ ⁢ ϵ t ; s t + ▽ ϕ ⁢ π ϕ ⁢ a t | s t , (30) where a t is evaluated after adding noise as f ϕ ( ϵ t ; s t ) . This method can be easily extended from the determined policy gradient to any tractable stochastic policy. Finally, the policy updates itself through the learning rate λ π as ϕ = ϕ − λ p i ▽ ϕ J π ϕ . (31) Note that in the SAC approach, the update rule for the temperature parameter α involves minimizing a specific objective function to ensure that the entropy of the policy remains at a desired level. Hence, the objective function for α is designed to minimize the difference between the current policy entropy and a target entropy H target . The loss function for α is given by J α = E a t ∼ π − α ⁡ log ⁡ π a t | s t − α H target . (32) The gradient of the loss function with respect to α is further calculated as ∇ α J α = E a t ∼ π − log ⁡ π a t | s t − H target , (33) α ← α − λ α ∇ α J α , (34) where λ α is the learning rate of α ; H target is usually given as a hyperparameter according to the specific task or desired level of exploration. Equations 30–34 give the gradient calculation of the hyperparameter α and policy network parameters. The pseudocode for the SAC algorithm is summarized as Algorithm 2.

Algorithm 2

Soft actor–critic algorithm.

1: Initialize the parameters for the two Q networks and policy network θ 1 , θ 2 , ϕ

The system settings consist of two aspects, which are the physical grid environment and RL agent. For the environment settings, the single-line diagram of the modified IEEE 39-bus system is shown in Figure 4, which includes three AGC units and seven non-AGC units. The three AGC units are installed at buses 31, 38, and 39; a wind farm of 300 MW capacity is installed at bus 39, and an external power grid is connected through a tie-line at bus 29. The parameter settings for the power system are listed in Table 1. Note that the control period and control step are set as 15 min and 1 min in this work, and the initialized deviations of the system frequency and tie-line power are assumed to be 0. The wind fluctuations were obtained from the New England power grid ² , and the loads were assumed to be the same over the 15 min duration because load fluctuations are usually smooth compared to wind power fluctuations. Figure 5 shows the minute-level fluctuations of wind power in three random periods; it is seen that the stochastic wind power changes heavily even in adjacent time steps, with a maximum fluctuation of over 10 MW.

For the agent settings, the discount factor γ is set as 0.95 because the current power adjustment of the AGC units significantly impacts the future operating state of the system. The neural network structures of all the actor–critic networks are the same (i.e., Q θ 1 , Q θ 2 , Q θ 1 ′ , Q θ 2 ′ , and π ϕ ) and comprise two hidden layers of size 128 × 64 . The smoothing factor for the two target networks is designed as β = 0.01. Adam optimizer was adopted for gradient optimization with the learning rate λ α = λ Q = 0.001. The three weight factors are set as w 1 = -0.05, w 2 = -1, and w 3 = -20. The replay buffer size | D | is set to 100,000. The introduced noise follows a Gaussian distribution of the form ϵ t ∼ N ( 0 , 1 2 ) . The simulations were implemented in Python using an Intel Core i7 CPU @3.0 GHz and 16 GB memory.

In this work, we evaluated three benchmarks to validate the benefits of the proposed method: 1) proposed imitation-based SAC strategy denoted as ISAC; 2) traditional SAC strategy that uses the random initialization policy; 3) a classical proportional integral (PI) strategy. The PI strategy regulates the AGC units in proportion to the system frequency deviations. The following subsections show the results of both the training and control processes.

The training process of the RL agent involves two stages: offline imitation learning and online RL training. The offline dataset is split into two subsets by random sampling, where 80% of the data are used for training and the remaining 20% are reserved for testing (i.e., validation). For the imitation network, the input data are the observed system state, and the label is the power regulation of the AGC units based on the classical PI strategy. This means that the final converged imitation policy is similar to the PI strategy. Figures 6, 7 show the imitation learning curves through the root mean-squared error (RMSE) and mean absolute error (MAE) indices. The black line represents the loss of the test dataset, and the red line represents the loss of the training dataset. It is observed that the training of the imitation network converges efficiently after approximately 200 epochs. Although the test loss is larger than the training loss initially, the converged strategy achieves training and test losses that are both less than 0.01. The training convergence is achieved at 10.04 s. These results indicate that the imitation learning process is both fast and stable, making it feasible for practical implementation.

To demonstrate the convergence of the proposed ISAC method, Figure 8 shows the cumulative reward of each episode through the blue lines and moving average reward via orange lines. It is seen that the episode reward decreases from −55 initially to approximately −15 after over 2,000 episodes, which implies successful convergence. Note that the cumulative reward of each episode still shows oscillations even after later training. This is because the wind power fluctuations are uncertain and vary widely over different periods, causing the optimal cost to be dynamic and inconstant. However, from the moving average curve, the episode reward is seen to decrease with continuous training until convergence.

To investigate the effects of the imitation network, we applied a random initialization policy in the SAC approach for comparison. Figure 9 shows the CPS1 training results based on two RL methods (i.e., SAC and ISAC). The blue line shows the CPS1 result for the ISAC method, and the orange line denotes the CPS1 result for the conventional SAC method. First, we see that the CPS1 results in both methods converge to 2 after 1,500 episodes. At the beginning of the training process, the value of the CPS1 index in the conventional SAC method is much lower than −2, with a minimum value of almost −12. This means that the system frequency stability is unacceptable for real-world grids. However, in the proposed ISAC method, the CPS1 result is always maintained within a small fluctuation range of [-2, 2] even at the beginning. This is because the proposed imitation network provides a satisfactory initialization policy by collecting effective samples for better policy optimization. It is noted that the sampling efficiency is sacrificed in the ISAC method to ensure safe exploration, where the CPS1 value is stable after approximately 1,500 episodes. In the conventional SAC method, the CPS1 value reaches 2 at nearly 1,000 episodes; this is because a large number of unsafe samples accelerate the learning process. Therefore, these results indicate that the proposed imitation network effectively prevents the unsafe random explorations observed in conventional RL methods, which can help cope with safety issues in real-world grids, especially for safety-critical AGC problems.

We applied two well-trained agents and a PI controller for comparison in the AGC environment. Figure 10 presents the system frequency fluctuations over 15 min with the same wind inputs based on the three control strategies, where the purple line denotes PI control, blue line denotes the SAC agent, and red line represents the proposed ISAC agent. At time step t = 0, the initial frequency deviations are all 0. It is seen that the maximum frequency deviation of the PI controller is 0.076, which violates the upper limit of the frequency. However, the maximum values are only 0.045 and 0.028 in the SAC and proposed ISAC methods, respectively. This means that the two RL methods can adhere to the limits by effectively predicting wind fluctuations and taking actions in advance, with the proposed ISAC slightly outperforming the conventional SAC method with smaller frequency deviations.

Table 2 presents a detailed comparison of the system control results, including the ancillary service costs, CPS1 index values, maximum frequency deviations, and average frequency deviations. We can see that the ancillary costs in the SAC and ISAC methods are significantly higher than that with the PI method because accurate responses to wind fluctuations require more power regulation in the AGC units. Correspondingly, the CPS1 values are effectively improved from 192.58% to 196.46% in the SAC and 199.27% in the proposed ISAC approaches. This means that there exists a tradeoff between the frequency stability and economic benefit depending on the system preference. In addition, the average frequency deviation with PI is 0.029, which decreases to 0.009 in the proposed ISAC to achieve a more stable system frequency. Hence, the converged SAC agent can also achieve satisfactory control results, with the main disadvantage being the unsafe training process. With the proposed imitation learning approach, the training process is safer and the final converged policy is improved through consideration of expert experiences. Figure 11 shows the AGC power regulation curves of the three AGC units. Although the regulation capacities and ramp limits are different, we can see that the curve trends for the three AGC units are quite similar. Specifically, for AGC units 1 and 2, the power deviations start at approximately −3 MW and decrease, reaching approximately −13 MW. At the end time, slight improvements are observed, with the deviations moving toward −5 MW. For AGC unit 3, the power deviation starts at approximately −5 MW and decreases to approximately −20 MW, showing more significant regulation than the other two units. This is because AGC unit 3 has the largest upper output power of 1,100 MW. In summary, all three AGC units experience significant negative power deviations initially, possibly owing to load changes or outdoor environment conditions, before recovering toward the end.