An Improved Model Predictive Fast Frequency Control for Power System Stability Against Unknown Time-Delay Switch Attack

A modern power grid is a cyber-physical system, which are vulnerable to cyber attacks. A recently found attack, the time-delay switch attack (TDSA), is made by inserting time delays into communication channels. A TDSA can be highly destructive to a power system as it can lead to instability. This paper presents a novel model predictive control (MPC) for fast frequency controller in a power system which can effectively mitigate the unknown TDSA. The MPC recently has received great attentions to be applied as FFC in a power system. Most of the MPC design are based on discrete-time model, whose future plant behaviour is calculated through iteration, rather than convolution. Nevertheless, one crucial step in the derivation of discrete-time MPC (DTMPC) is to capture the control trajectory over a finite prediction horizon. This imposes a challenge in designing a DTMPC to counteract the time-delay with unknown time length. Thus, a continuous-time MPC (CTMPC) is proposed to deal with TDSA. To overcome the unknown time delay, we synthesize an accurate time-delay estimator and sequential state predictor (SSP), are designed to accurately estimate and effectively counteract the unknown and random TDSA. All presented case studies are based on a real Taipower system and justification of the effectiveness of the proposed method was verified.


INDEX TERMS
change of mechanical power. P ESS power of energy storage system. P ESS rate of power of energy storage system. P e change of electrical power. 17 The associate editor coordinating the review of this manuscript and approving it for publication was Suman Maiti .f frequency deviation. K ob observer gain. K SPN state predictor gain. K mpc control gain. z (t) approximation of the predicted state. m(t) external signal. C 0 the set of all continuous functions. C 1 the set of all functions whose derivative are in C 0 . σ (t) sliding surface. l i (t) orthonormal basis functions i. r the number of Laguerre networks.

L(t)
Laguerre function. λ time scaling factor for the Laguerre functions. J expression of the cost function. Q weighting matrix of predicted state. T p horizon of predicted time.
ϕ an arbitrary future time. R weighting matrix of input. k step size. y max minimum frequency deviation. y min maximum frequency deviation. α 1 , α 2 parameters of super-twisting time-delay estimator.
works have been done by [5], [6], [7], [8], which employ PI 48 controllers for LFC applications. These works are more con-  Even though this method can enlarge the critical time-delay to 53 some extent, the robustness against the random TDSA has not 54 been ensured yet. A method to counteract random time-delay 55 attack based on perturbation observer has been devised by 56 [10]. In their work, the random time-delay is assumed as 57 an unknown disturbance, whose effect is necessary to be 58 minimized to some degree. Although this method does not 59 require to exactly estimate the time-delay values to tune the 60 control parameters, this method is not able to compensate 61 relatively long time-delay. It is to be noted that the methods 62 proposed by [4], [5], [6], [7], [8] merely evaluate the maxi-63 mum critical time-delay by selecting the appropriate control 64 parameters. Nevertheless, in practice, the TDSA can be 65 time-varying and random, which makes the work of [4], [5], 66 [6], [7], [8] not applicable to counteract the time-varying and 67 random TDSA. Two approaches for compensating the TDSA 68 are by redesigning the controller to satisfy the robustness 69 and by estimating the unknown TDSA without the need of 70 redesigning the controller. The first approach is prone to be 71 unstable if the TDSA is larger than the specified critical-time 72 delay. The second approach, on the other hand, is intended to 73 estimate the unknown TDSA and compensate the time-delay. 74 There are studies for estimating time delays, and they can be 75 generally classified into five approaches, which are 76 optimization-based, convolution-based algebraic, adaptive 77 backstepping-based, artificial intelligence-based, and slid-78 ing mode-based. In the optimization-based, the estimated 79 time-delay is an optimal solution obtained based on 80 a cost function minimization. Some studies related to 81 optimization-based approach can be found in [11], [12], [13] In the convolution-based algebraic approach, the convolution 92 method is applied to identify the unknown parameters and 93 constant delay of a linear system. This approach provides 94 high convergence, but it is not applicable to time-varying 95 or random time-delay estimation. Some studies related to 96 this approach are covered in [17], [18], and [19]. The third 97 approach, which is a class of nonlinear systems, applies 98 partial differential equation (PDE) transformation to estimate 99 the unknown time-delay. For instance, in [20], [21]

158
• It yields practicability as the proposed FFC is able to 159 compensate TDSA without the need of redesigning the 160 controller.

161
• It evaluates the robustness and effectiveness of the 162 improved model predictive fast frequency controller 163 under parameter changes and various TDSAs.

164
The paper is organized as follows. Section II describes 165 the structure of the proposed FFC. Section III describes how 166 various constraints are incorporated into the designed MPC to 167 be applicable in a power system. Section IV specifies control 168 parameters and the effect of system parameter changes on the 169 critical time-delay. Case studies are presented in Section V. 170 Conclusions are drawn in Section VI. Fig. 1 depicts the block diagram of the proposed system, 173 composed of the power system, acting as a plant, the unknown 174 time-delay, and the proposed control method. The unknown 175 time-varying time-delay constitutes lumps of input and output 176 delays and is represented as a RTT delay. Furthermore, the 177 proposed control method is composed of three subsystems, 178 which are the super-twisting time-delay estimator (ST-TDE), 179 sequential state predictor (SSP), and MPC, which consists of 180 an observer and optimizer. The unknown RTT delay, τ (t), 181 is estimated by ST-TDE to yield the estimated time-delay,τ . 182 Then, the estimated time-delay is fed to the observer and 183 SSP. The observer receives the measurement data, which 184 are specifically the measured frequencies, to yield the esti-185 mated state variablex (t). The estimated state variable is 186 subsequently processed by the SSP to predict the future state 187 variable under a given estimated time-delay, which results 188 in z N . Note that z N is an approximation ofx (t + τ ). In general, the aggregated PFR dynamics incorporated with 199 the ESS power is formulated as [35] 200ẋ

II. PROPOSED CONTROL SYSTEM DESIGN 172
Note that the PFR parameters, such as inertia constant (H ), To minimize the steady-state error, an integrator can be 221 embedded in (1), and an auxiliary variable is introduced as Note thatP e is basically a step function andṖ e = 0. Also, 226 the output state variable dynamics can be obtained as Augmenting (3) and (4) to the system yields where 232 Eq. (5) serves as an augmented system of (1) containing an 237 integrator. To consider the existence of time-delay in (5), the 238 dynamic representation of a time-delay system is expressed 239 as 240 where τ (t) is the unknown RTT time-delay. Note that since 242 τ (t) is an unknown value, the estimation is carried out to 243 obtain the estimated value of time-delay,τ (t). Section II.B. 244 discusses how the proposed ST-TDE can accurately estimate 245 the unknown TDSA. Furthermore, the observer is applied 246 to estimate the state-variables from the input and output 247 variables. The observer dynamics is represented as where K ob is the observer gain, whose value is determined via 251 pole-placement technique.

253
SSP is devoted as a countermeasure to deal with a relatively 254 long time-delay system. This method is also well-known as 255 a chain predictor. The term ''long'' is used to describe the 256 time-delay, which is larger than the critical time-delay. To cir-257 cumvent this problem, this approach extends a sequential 258 prediction containing copies of the system dynamics which 259 run at different time scale.

260
The concept of this approach resembles a basic predictor, 261 which essentially provides the estimation of the predicted 262 state, x (t + τ ). Without loss of generality, z (t) is defined as 263 the approximation of the predicted state, x (t + τ ). In [36], 264 the SSP is expressed as a set of successive estimation of 265 x t + i τ N , where i = 1, 2, . . . , N defines the number of set. 266 Accordingly, the SSP dynamics are represented as . . , N and K SP1 , K SP2 , 275 . . . , K SPN denote the predictor gain for i = 1, 2, . . . , N , 276 respectively. For simplicity, one may select K SP1 = K SP2 = 277 . . . = K SPN . As τ (t) is unknown, the estimated time-delay, 278 τ , yielded from ST-TDE is used. Moreover, as x (t) is not 279 readily available,x (t) obtained from the observer, i.e. (7) is 280 employed. Consequently, SSP in (8) can be modified into (9). 281 288 z N can be obtained by integrating (9) where m (t) is an external signal, which satisfies m (t) ∈ C 1 .

310
Note that C 1 is the set of all functions whose derivatives are 311 in C 0 , and C 0 is the set of all continuous functions. Thus, for 312 simplicity, we define m (t) = t.

313
ST-TDE dynamics can be represented as Recall that the input of the augmented system (5)  Laguerre function dynamics is expressed as ∈ R m×m and λ repre-361 sents the time scaling factor for the Laguerre functions whose 362 value is λ > 0. Thus, the exact solution of (16) can be 363 The cost function is formulated as a function of pre-368 dicted state variables and control, which are weighted by 369 Q = Q T ≥ 0 and R = R T > 0, respectively. Typically, these 370 matrices are diagonal matrices to simplify the representation.

371
The expression of the cost function is stated as The first term of (18) is the predicted state variables, whose 375 solution is obtained as Note that the orthonormality of Laguerre function allows Substituting (19) and (21) into (18) gives to η as Consequently, the optimal trajectory of derivative of control 404 at ϕ can be expressed as which after some algebraic manipulation, is also equivalent 413 to whereη denotes the estimated vector coefficient. Accord-416 ingly, the optimalη minimizing (28) is obtained as Note that the difference between η in (25) andη in (29) 419 lie on the use of all state variables measurement. Oncex (t) 420 converges to x (t), thenη resembles to η. In practical situations, the power of an energy storage sys-434 tem is restricted under specified limits. These limitations 435 are considered as actuator constraints and they need to be 436 incorporated to the control formulation to satisfy the control 437 objective. Furthermore, to maintain the system frequency 438 standard, the allowable range of the system frequency is 439 included in the control design. Accordingly, the constraints 440 are categorized into three terms, i.e. power limitation of 441 energy storage, rate of power limitation of energy storage, 442 and system frequency. These constraints are formulated into 443 linear inequalities, which are necessary to provide boundary 444 conditions for real-time optimization. it can be expressed as The allowable frequency deviation range is represented as  where M and β are the vectors whose elements are com-505 posed by the ESS power limitations and allowable system 506 frequency, defined in (32), (37), and (41) as In other words, the numerical solution of this method 510 is concerned with the problems of constrained minimiza-511 tion where the cost function is a positive definite quadratic 512 function and constraint functions are linear functions. A QP 513 approach can be applied to obtain such solutions.

516
According to (1), the relationship between the frequency 517 deviation to electrical power change can be described as From (45), one may observe that the lumped power system 520 and governor can be represented into an approximate SISO 521 system using generation and load aggregation. The unknown 522 parameters that need to be estimated are H , D, T g , and R g .
The estimated system parameters a 0 , a 1 , b 0 , and b 1 are 535 obtained via least squares method, which is describe as Since the element of the regressor is essentially known and is 554 based on the measurement data, the element of the regressor 555 containingf m (t) andḟ m (t) can be replaced byf (t) andḟ (t), 556 respectively. Consequently, it results in Defining the modeling error as the difference between the 559 actual data and model, stated as A cost function, which defines the function of modeling 562 error, is formulated as J LS , minimizing µ, can be obtained by taking the partial 568 derivative of J LS with respect to µ as Accordingly, the optimal µ can be obtained as The validation of the proposed method is conducted 573 for Taiwan power systems. Two events recorded in PMUs 574 occurred on September 14, 2018 (event 1) and December 1, 575 2018 (event 2) are taken as examples of generation-loss cases. 576 Event 1 dealt with three generator units trip whose total power 577 generation is 708 MW. On the other hand, one generator unit 578 is tripped with total power generation of 550 MW in event 2. 579 Figure 2 shows the frequency responses recorded from the 580 PMU, resulted from two actual generation-loss events. The 581 sampled time of PMU is 0.05 s. The figure shows that the 582 models obtained from system identification technique are in 583 great agreement with the real measurement data.

584
The total system loads of events 1 and 2 are 32607 and 585 24127 MW, respectively. The model for even 2 is stated 586 in (56). The corresponding estimated power system parameters are 589 obtained as H = 11.22, D = 0.48, R g = 0.22, and T g = 590 5.97 based on calculation in [33]. Table 1 lists parameters of the proposed method. The 592 weighting matrices Q = Q T ∈ R 3×3 and R ∈ R 1×1 are 593 associated with the predicted state variables and the time 594 derivative of the control input, respectively. The optimal Q 595 is often selected as Q = C T C [28]. With this choice, the 596 closed-loop eigenvalues are determined by R. The larger R 597 generates the smaller control input, which results in slower 598 settling time. On the contrary, the smaller R produces the 599 higher control input, which yields fast response, yet exhibits 600 overshoot. 601 Fig. 3 shows the frequency trajectory of the power sys-602 tem under various time-delay constants. As can be seen, the 603 maximum time-delay causing the system to be marginally 604 VOLUME 10, 2022  stable is τ = 0.568s. This indicates that MP-FFC is unable to 605 compensate the time-delay, which is greater than τ = 0.568s.

606
Changes of system parameters are also a key factor con-607 tributing to the shift of the critical time-delay. To validate this, 608 the MP-FFC designed for event 2 is tested to event 1, whose 609 system parameters are different from event 2. In event 1, two types of TDSA, i.e. constant and random TDSAs are 619 investigated. The first and second cases are to study the 620 robustness of the proposed controller against constant and 621 random TDSAs, respectively. On the other hand, the robust-622 ness of the proposed controller under parameter variations, 623 together with constant and random TSDAs, is investigated 624 in cases 3 and 4, respectively. It is to be noted that these 625 TDSAs occur at t = 0s, which are subsequently followed 626 by an under-frequency contingency (i.e. generator outage) 627 at t = 6s. Moreover, to justify the proposed method, three 628 controllers are tested under these cases. The first controller 629 is named as model predictive FFC (MP-FFC) without con-630 straints. It is based on the model developed in most of the 631 studies such as [38], [39] where the time-delay estimator 632 and corrector are not included in the control formulation. 633 Moreover, the system constraints such as the ESS limitation 634 and frequency constraints are also not taken into account. 635 The second control is called improved model predictive FFC 636 (IMP-FFC) without constraints. Different from MPC-FFC, 637 IMP-FFC employs the proposed time-delay estimator (i.e. 638 ST-TDE) and corrector (i.e. SSP). Finally, the third controller 639 is the complete proposed controller, which is called IMP-FFC 640 with constraints. It consists of the model predictive control 641 together with the proposed SSP and ST-TDE, and takes the 642 system constraints into account as well. The constant TDSA is set to τ (t) = 5s. The TDSA 646 value is selected greater than the critical time-delay (about 647 8.8 times) to prove the robustness of the proposed method. 648 Fig. 5a depicts the trajectory of frequency under constant 649 TDSA and generator outage. As can be seen, the MP-FFC 650 without constraints causes the frequency to diverge from its 651 nominal value and results in unstable response since the value 652 of TDSA is larger than the critical time-delay that can be 653 withstood by the corresponding controller. It is indicated 654 that at t = 6s, the frequency response of MP-FFC causes 655 the frequency uncontrollably oscillates. On the other hand, 656 the proposed method is able to regulate the frequency in 657 the presence of TDSA and under-frequency contingency by 658 accurately estimating the unknown TDSA and effectively 659    by the proposed method. As depicted in Fig. 8, the resulting 693 time-delay estimation produced by the proposed method can 694 quickly converge to the actual TDSA.

695
The proposed method is also able to allocate P ESS , satis-696 fying the range of the specified power limitation. As evident 697 in Figs. 7(b) and 9, the ESS power and ramp rate generated 698 by the proposed method are still within the specified limits. 699 This indicates that the proposed method is able to handle the 700 constraints determined by the ESS specifications. 701 VOLUME 10, 2022    Table 3. Figure 10(a) depicts 713 the trajectories of the system frequency. As can be seen, 714 the IMP-FFC with constraints is able to stabilize the power 715 systems. The frequency nadir governed by the IMP-FFC with 716   In this case, the system is tested under parameter variations, 731 whose setup is the same as that of Case 3. However, random 732 TDSA is injected at t = 0, whose value range is the same as 733 that of Case 2. The power limitations of ESS for this case are 734 listed in Table 3.  The system frequency trajectory is shown in Figure 12 without constraints are not able to stabilize the frequency. 744 The frequency trajectories yielded by these methods tend to 745 diverge and exhibit severe oscillations.

746
In terms of constraints handling, the IMP-FFC with con-747 straints is able to satisfy the constraints given by ESS 748 specifications. As can be seen from Figures 12(b) and 13, 749 the IMP-FFC with constraints is able to allocate P ESS and 750 P ESS within their specified ranges. On the other hand, both 751 MP-FFC and IMP-FFC without constraints are not able to 752 allocate P ESS andṖ ESS within their specified ranges since the 753 constraints are not considered in their control design.

754
Furthermore, the random TDSA can be accurately esti-755 mated via STA TDE as shown in Figure 14. By accurately 756 estimating the unknown TDSA, the controllers can provide 757 control actions corresponding to the estimated TDSA.

759
In this paper, a unified fast frequency control, based on 760 model predictive control has been proposed. The proposed 761 method is not only able to regulate the system frequency 762 to the nominal value during large disturbance, but also able 763 to counteract any types of unknown TDSA. By accurately 764 estimating the unknown TDSA, the proposed method can 765 compensate the attacks injected by the hacker. Thus, the 766 power system resiliency can be maintained. Furthermore, the 767 proposed method is also able to optimally allocate ESS power 768 based on ESS specifications, such as ESS power and ramp 769 rate, and frequency limits set by the grid code.