Tertiary regulation of cascaded run‐of‐the‐river hydropower in the islanded renewable power system considering multi‐timescale dynamics

To enable power supply in rural areas and to exploit clean energy, fully renewable power systems consisting of cascaded run-of-the-river hydropower and volatile energies such as pv and wind are built around the world. In islanded operation mode, the primary and secondary frequency control, i.e., hydro governors and automatic generation control (AGC), ensure the frequency stability. However, due to limited water storage capacity of run-of-the-river hydropower and river dynamics constraints, without coordination between the cascaded plants, the traditional AGC with fixed participation factors cannot fully exploit the adjustability of cascaded hydropower. When imbalances between the volatile energy and load occur, load shedding can be inevitable. To address this issue, this paper proposes a coordinated tertiary control approach by jointly considering power system dynamics and the river dynamics that couples the cascaded hydropower plants. The timescales of the power system and river dynamics are very different. To unify the multi-timescale dynamics to establish a model predictive controller that coordinates the cascaded plants, the relation between AGC parameters and turbine discharge over a time interval is approximated by a data-based second-order polynomial surrogate model. The cascaded plants are coordinated by optimising AGC participation factors in a receding-horizon manner, and load shedding is minimised. Simulation of a real-life system shows a significant improvement in the proposed method in terms of reducing load shedding.


INTRODUCTION
In mountain areas such as in Tibet and Sichuan Province in China as well as the plateaus in South Asia and Africa, to enable electrical power supply to the local residents and to exploit clean energy, cascaded run-of-the-river hydropower plants have been built along river valleys [1][2][3][4][5]. Moreover, in recent years, photovoltaic (pv) and wind plants are built in these areas to further exploit the renewable energy and to compensate for power shortages in dry seasons [2,6,7]. The hydroelectric, solar, and wind power make up fully renewable local power systems. A In islanded operation mode, considering the volatility of solar, wind, and load power, energy balance and frequency stability are the most important issues [10,11]. In this situation, due to the lack of support from an external power grid and energy storage systems (ESSs) in rural areas, the primary and secondary control of cascaded hydropower, that is the hydro turbine governors and automatic generation control (AGC), are essential to stabilise the frequency against pv/wind and load volatilities.
But with a high penetration of renewables, frequency stability is challenging, as the pv and wind volatility may exceed the adjustability of hydropower units which may cause stability issues. Although by simply disconnecting pv or wind from the islanded grid instability can be completely avoided, it may cause power supply shortage and consequent loss of load. Considering socioeconomic benefits, applying controlled curtailment for smoothing renewable power [12][13][14] and using the hydropower to regulate the system seems to be a better option.
Until now, on the topic of using cascaded hydropower to mitigate solar and wind volatility, elaborations have been made by the community, including scheduling [15][16][17] and online control [18][19][20], but these works focus on grid-connected operations. On the other hand, studies on the frequency control of islanded power systems with hydropower focus on the timescale of the electromechanical dynamics, relating to the primary and secondary control [10,11].
In addition to the primary and secondary frequency control, in a longer time scale the mismatch between the stochastic pv power and load accumulates. If the accumulated energy is not properly allocated to the cascaded hydropower plants, due to river dynamics some of them may lose adjusting ability as explained later, endangering stability or causing loss of load. To ensure normal operation over a longer timescale, a study on the tertiary regulation that considers the hydraulic coupling between the cascaded plants is also needed, which is, to the best of the authors' knowledge, still lacking. This paper aims to fill this gap.
In stabilising the frequency against the volatilities of wind and solar energy, the primary and secondary frequency control adjusts the power output of the cascaded hydropower plants. Furthermore, in addition to the capacity and ramping limits, the adjustment of cascaded run-of-the-river hydropower is also limited by river dynamics. This is because in contrast to the conventional dam hydropower, in run-of-the-river hydropower, the water storage capacity is very limited, and the water energy is spatially distributed along the river. The utilisation of water is subject to the river dynamics, and the cascaded plants are therefore hydraulically coupled [18][19][20]. Moreover, river operation and ecological regulations often require that the water stage along the river is limited within an allowed interval [21], which further limits the adjustment of hydropower.
If we do not consider the river dynamics but only use a traditional AGC that allocates incremental power with fixed proportions to the cascaded plants, due to the impact of solar and wind volatility on hydropower generation and then the water distribution along the river, unacceptable violation of the river operation constraints or a large amount of load shedding may occur, as exemplified in Section 4.5. In contrast, if we consider the river dynamics and accordingly coordinate the cascaded plants in a tertiary regulation scheme by dynamically adjusting the AGC participation factors, as later illustrated in Figure 3, the adjustment ability of the hydropower can be maximised, and consequently load loss can be minimised. However, power system frequency dynamics and river dynamics, which is usually characterised by the shallow water equations [18][19][20][21][22], have very different timescales, as illustrated in Figure 2(a). If we directly combine them to establish a model predictive controller (MPC) for the tertiary regulation with a time resolution compatible with the frequency dynamics and a horizon that can accommodate the river dynamics, the curse of dimensionality will arise. To address this issue, in this work a data-based polynomial surrogate model describing the relation between the tertiary control variables (including AGC participation factors, system frequency reference and load shedding) and turbine discharges over a certain time interval is developed. Thus, the detailed model of power system dynamics can be replaced by simple algebraic functions and then easily incorporated into the MPC formulation with the dynamic river model without causing the curse of dimensionality, as illustrated in Figure 2(b).
Following the above presented ideas, this paper proposes a coordinated tertiary regulation approach for islanded power systems with cascaded run-of-the-river hydropower and volatile generations. The framework of the proposed approach is shown in Figure 3. Specifically, the following two contributions are made: 1. A data-based polynomial surrogate model is first established to describe the hydro turbine discharge as a function of AGC participation factors, power system frequency reference, and load shedding. This model can be easily incorporated into an MPC formulation for the tertiary regulation of cascaded run-of-the-river hydropower.  Framework of the proposed regulation approach. The orange, blue, and purple colours represent the electromechanical side, the hydraulic side, and the controller, respectively 2. Based on the polynomial surrogate models of power system dynamics and a dynamic river model based on the shallow water equations, a tertiary regulation approach based on MPC jointly considering the multi-timescale dynamics is proposed by updating the AGC participation factors to coordinate the cascaded plants and to minimise load shedding.
Simulation of a real-life system with actual real-time pv data verifies that compared to other regulation schemes the proposed approach significantly reduces load loss.
Specially note that our previous work [20] focuses on grid-connected operation of the cascaded hydro-solar system and does not consider any multi-timescale dynamics. Despite both works adopt the same river model, the scope, methodology, and contribution of [20] and this work are totally different. This paper is organised as follows. Section 2 develops the polynomial surrogate model for the electromechanical side of the islanded system and introduces the dynamic river model. Section 3 presents the optimal control formulation of the tertiary regulation. Finally, in Section 4, case studies are presented.

DYNAMIC MODELING OF ISLANDED POWER SYSTEMS WITH CASCADED RUN-OF-THE-RIVER HYDROPOWER
To coordinate cascaded hydropower plants in the tertiary regulation, the electromechanical side, that is the power system dynamics, and the hydraulic side, that is the river dynamics, need to be jointly considered. As noted in the Introduction, we establish a polynomial surrogate model in Section 2.1 to describe the relation between the water discharge of the hydropower plants and the tertiary control variables, which represents the impact of the primary and secondary frequency control on the hydraulic side. Section 2.2 uses the shallow water equations to depict the river dynamics. Then, in Section 2.3, these models are combined to depict the overall system.

2.1
Electromechanical side: Polynomial surrogate model of power system dynamics With proper smoothing curtailment control of renewable generations, the frequency of an islanded power system with small hydropower and volatile energy can be stabilised by turbine governors and AGC [8][9][10][11], known as the primary and secondary control. However, considering the river dynamic constraints in the cascaded run-of-the-river hydropower, the AGC arguments, for instance the participation factors, should be repeatedly optimised to coordinate the cascaded plants over longer time periods, which is referred to as the tertiary regulation or tertiary control in this work.
In establishing the optimal tertiary regulation model, the turbine discharges Q turb Hi at the hydropower plants connect the electromechanical and hydraulic sides. When pv or load fluctuates, the PI controller in the AGC calculates the incremental power reference for each hydropower plant based on system frequency deviation; see detailed AGC model in literature such as [23]. The AGC participation factors c Gi of the hydropower plants allocate the energy used in the secondary control and thus affect the turbine discharge Q turb Hi of the plants. Moreover, in islanded operation mode, load can be adjusted by deviating the system frequency reference ref Over a time interval t ∈ [0, T ], given time-domain trajectories of the volatile generation of solar or wind energy {P S (t )} t ∈[0,T ] and load {P D (t )} t ∈[0,T ] , the tertiary control variables p determine the electrical power P Gi of each generator, which is denoted by the following function, With a time horizon compatible with the river dynamics, the mean value of P Gi (t ) over the time interval t ∈ [0, T ] is considered, as Due to the non-linearity and complexity of the power system, the analytical expression of (3) does not exist. To deal with this problem to establish the system model as Figure 2(b), we instead employ an approximate surrogate model to facilitate modeling the optimal model predictive control problem for the tertiary regulation.
First, rewrite the trajectories of the volatile generation and the load as the products of preset normalised trajectories and the corresponding mean values, as T ] represent the typical trajectories of the volatile energy and load demand with a unit mean value, respectively. Substituting (4) into (2) and finally into (3), we can observe that the mean electrical power of a hydro plant is a function of the mean values of the volatile generationP S , the load demand P D over t ∈ [0, T ], and the tertiary control variables p, denoted asP Followingly, polynomial approximation technique is used to construct the surrogate model forP Gi (⋅). The obtained polynomial surrogate model has the form as where i (⋅) is the multivariate Legendre polynomial basis function; f i is the coefficient, obtained by the collocation method based on power system dynamic simulation results. See detailed procedure of constructing the polynomial surrogate model in the Appendix A.
To facilitate the formulation of the proposed tertiary regulation as a model predictive controller (MPC) that is solvable as a mixed-integer quadratic programming (MIQP), the approximation order of (7) is limited to 2; see Section 3.2 for details.
Further, generally load shedding is realised by tripping feeders. The total amount of load shedding is the sum of the products of the binary variable j and the capacity of feeders P sh D j , as and the tertiary control vector is re-denoted as Substituting (8) into (7) and rearranging, we can approximate the mean power of each hydropower generator as a function of P S ,P D , c G , ref , and , as The efficacy of the polynomial surrogate model (10) (or (7)) is validated numerically. In Section 4.6, simulations on PSS/E show that (10) accurately gives the mean values of the hydropower outputs over each 10-min interval, which is adequate for establishing a receding-horizon controller.
Finally, we establish the relation between the mean hydropower generationP Gi  Hi , known as the production function [24], which is expressed asP To facilitate the formulation of the controller, by linearising (11) the turbine discharge is approximated by a linear function aŝ Substituting (10) into (12), the turbine dischargeQ turb Hi of each hydropower plant during the primary and secondary frequency control process is finally represented as a polynomial function ofĤ up Hi ,Ĥ down Hi , c G , ref , and j . In a word, by adjusting the tertiary control variables (9), the turbine discharge of the hydropower plant over a time interval can be controlled.
Note that the linearisation of (11) has also been adopted in many existing works [18,20]. Because model (12) is only used in the tertiary regulation, such linearisation will not affect the primary and secondary control and therefore the frequency stability of the islanded system will not be affected. Further, in online operation the linearisation point is repeatedly updated based on the current operation point, through which the linearisation error can be alleviated in a receding-horizon regulation manner.

2.2
Hydraulic side: State-space model of river dynamics According to previous researches by the community [18,19,25] as well as our previous work [20], for modeling the river dynamics in the run-of-the-river hydropower the shallow water equations is a suitable model. This section only gives a brief introduction to this model, and details can be found in the above mentioned literature.
The shallow water equations are partial differential equations of the water volume and momentum conservation [22]: Qi-1

Hi+1
Hi-1 Spatial discretisation framework of the water flow. Detailed formulation can be seen in [20] where y denotes the position. The equivalent friction slope I f is empirically modeled by the Manning-Strickler formula, which is a function of Q and H (see [22]). The parameters such as the river width and slope can be measured in field or estimated via data assimilation [26].
In normal operation, water stage varies within only a small range compared to the depth. Thus, the river dynamics can also be linearised [22]. Then, discretising the water flow into non-overlapping cells of length h as Figure 4 and using the finite-difference format, the linear dynamic river model can be obtained [20,22] compactly aṡ where A R , B R , and C R are constant matrices. The river control variables u R (t ) including the turbine dis-chargeQ turb Hi are determined by the operation of the hydropower plants, which are essentially determined by the primary and secondary frequency control and the tertiary regulation. The hydraulic coupling of the cascaded plants is naturally modeled in (15). The river operation constraints can also be formed with x R . See detailed deduction of (15) in [18][19][20].
In online operation, a state estimator, for example the Kalman filter [27], can be employed to provide river state estimation using available measurements such as the turbine discharges and water stages at the plants. This provides full-state feedback for the proposed controller, as shown in Figure 3.

The overall state-space model
Substituting the second-order polynomial surrogate model (10) of the power system dynamics and the turbine discharge model (12) into the dynamic river model (15), the overall state-space system model is obtained, denoted compactly aṡ where A, B(⋅), C, and F (⋅) are coefficient matrices, where A and C are respectively identical to A R and C R in (15); x(t ) and z(t ) are the same as x R (t ) and z R (t ); is the vector of decision variables of the tertiary regulation; the quadratic term u(t ) T F (P D (t ))u(t ) represents the quadratic terms of control variables in the second-order polynomial surrogate model (10); and B(P S (t ),P D (t )) and F (P D (t )) indicate these matrices are functions ofP S (t ) andP D (t ).
The framework of this model is illustrated in Figure 2(b).

Discussion on the dimensionality of modeling
In the proposed model (16), because the frequency dynamics are modeled by algebraic functions (10), the state variables are related only to the river dynamics. In the case study in Section 4, the size of the dynamic river model (15) is 200. Considering the MPC step length T = 600 s and horizon N = 12 or 7,200 s, the size of the discrete system model (24) is 200 × 12 = 2400, which is appropriate for online rolling optimisation.
Otherwise, if additionally considering the detailed power system dynamics with total n (n varies from tens to thousands depending on system complexity) state and algebraic variables, with a step length of 0.02 s that is compatible with the electromechanical dynamics, the overall size of the state reaches (n + 200) × 50 × 7200. Even if the power system and river dynamics are discretised in different time resolutions, the size of the states reaches n × 50 × 7200 + 12 × 200, much larger than the proposed model. Considering using interior point method (IPM) with time complexity O(n 4 ) to solve the MPC, the efficiency of the proposed model is significantly better.

MATHEMATICAL FORMULATION OF THE COORDINATED TERTIARY REGULATION
As shown in Figure 3, for the tertiary regulation an MPC is employed to repeatedly optimise the tertiary control variables including the AGC participation factors to coordinate the cascaded plants and to give commands of the frequency reference and load shedding in a receding-horizon manner.
There are two main reasons for this work to adopt the MPC as the tertiary controller. First, the cascaded run-of-the-river hydropower system is a dynamic system modeled by the state space equation (16) with various operational inequality constraints; MPC is very suitable for the optimal control of such a constrained dynamic system. Second, due to the stochastic nature of pv, the operation trajectory of the system may deviate from the predicted one, which can be well alleviated by the well-known receding horizon scheme of the MPC.
Given the prediction horizon N and step length T , the objective and constraints of the MPC are formulated, as explained below.

Objective function
The overall control objective involves coordinating the cascaded plants and reducing load shedding. Meanwhile, the river state, such as the water stage, should not deviate too far from the nominal. The detailed objective includes the following components.

Deviation in the frequency reference
According to (10), the load power can be adjusted by the frequency reference, which assists power balance in the islanded system. However, the frequency bias should not deviate from zero if unnecessary. Therefore, the following quadratic function is included in the objective:

Load shedding
On the premise of power balance and frequency stability, load shedding should be minimised to avoid loss of load, which is expressed as:

River stage deviation
During the regulation process, the water stage and discharge along the river and channels of the plants should not deviate far from the nominal. This can be achieved by minimising the following quadratic function:

Water spillage
When the upstream inflow exceeds the power demand and the upper limit of the water stage is encountered, plant water spillage is needed to ensure operational security. Spillage occurs only when needed, represented as minimising 3.1.5 Quadratic terms of the AGC participation factors Finally, to avoid oscillation in the AGC participation factors, the following quadratic term is included in the objective: The overall control objective is defined as a weighted sum of the above terms with positive weight parameters: The weights in the objective function (23) are chosen based on a tradeoff between different components. For example, when power supply and reliability are more preferred, larger J 1 and J 2 can be chosen; and when river operational constraints are more critical, J 3 should be increased. We set 1 = 10, 2 = 10,

Equality constraints
The equality constraints in the proposed controller include two components, as listed below.

System dynamics model
The state-space model (16) is temporally discretised into the equality constraints as for k = 0, … , N − 1, whereĀ,B,C andF are coefficient matrices; the solar/wind powerP S (kT ), load powerP D (kT ), and river inflow z(kT ) in the future are given as forecasts by the dispatching system, and are constant parameters in the optimisation problem. When using second-order polynomials to approximate the mean power output of the hydropower plants, as introduced in Section 2.1 and exemplified in Table 1, bilinear terms in terms of the tertiary control variables, that is c Gi ref and j c Gi are included as the last term of (24). Directly incorporating them into an optimisation may cause non-convexity, making the optimal tertiary regulation problem difficult to solve.
To address the bilinear terms to make the optimisation problem convex, ref is discretised as where K > 0 is an integer; j is a binary variable. Thus, the bilinear term c Gi ref is replaced with the linear combination of c Gi j .
Then, the big-M method is used to address c Gi j as Similarly, the bilinear term c Gi j is addressed as Substituting (26) and (29) into (24)

Sum of the AGC participation factors
Obviously, the AGC participation factors sum up to 1, as

Inequality constraints
The inequality constraints for k = 0, … , N − 1 include the following components.

AGC participation factors
The AGC participation factors should be within the interval

Hydropower generation limits
The power output of each hydropower plant should not exceed the limits, as follows: The bilinear terms inP Gi (⋅) are replaced with (26) and (29), and the resulting constraints are linear.

Electrical network constraints
Network constraints such as the branch limits can be considered as dc power flow:

Water state limits
The operation regulations require that the water stage stay within an allowed interval. This is considered at monitoring points along the river and channels: If additional constraints such as the switching counting limit of load shedding need to be considered, they can also be easily included in the MPC. To save space, this is not discussed here.

MPC formulation
Summarising all the above, the overall MPC problem in the proposed coordinated tertiary regulation is established, denoted as min U (23), subject to (24) where U = [u(0) T , … , u((k − 1)T ) T ] T is the sequence of the control variables. The MPC problem (37) is a typical mixedinteger quadratic programming (MIQP), which can be solved using commercial solvers such as IBM ILOG Cplex.
Once the MPC problem is solved, the AGC participation factors, frequency reference, and load shedding are updated based on the first entry of U . When shifted to the next step, the MPC is solved again based on the updated state estimation and forecasts in a receding-horizon manner. The power references of the plants are updated in real time by the PI controller and the repeatedly updated participation factors of the AGC, shown in Figure 3.
Note that although modeling errors are introduced in approximating the power system dynamics in Section 2.1 and linearising the river dynamics in Section 2.2, they can be alleviated by the receding-horizon scheme.

3.5
Online smoothing control of PV to reduce volatility As noted in Introduction, with a high pv penetration, the stochastic volatility may exceed the adjustability of the primary and secondary frequency control of the hydropower units and therefore poses an adverse impact on stability. Hence, before applying the above proposed tertiary control, a smoothing control of pv is needed to reduce the impact of stochastic variation on frequency stability. Supposing if the pv power smoothly follows its prediction trajectory, its impact on system stability can be well counterbalanced. However, we cannot simply track the prediction because the available pv power may drop below it. Hence, a controlled curtailment that considers the real-time volatility of pv is adopted in this work.
Specifically, the controlled curtailment of pv at minute k is as, where P pred S (k) is the prediction and P cur S (k) is the curtailment at minute k; P S (k) is the smoothed output; P avail S (k) is the available pv power determined by irradiance, usually the maximum power point (MPP); (40) calculates the volatility of pv for the past minutes; (39) is an low-pass filter with parameter to make the curtailment varies smoothly; is the level of curtailment.
When irradiance varies randomly to cause the available pv power varies significantly, v(k) becomes large and the curtailment of pv is increased to counter the volatility. In contrast, when pv volatility is mild, v(k) becomes small and the pv output largely follows the prediction. Thus, the output of pv becomes relatively smooth under different circumstance. Examples are illustrated in Section 4.3.

Simulation platform
To verify the proposed regulation approach, a detailed simulation platform is established jointly based on PTI PSS/E 34 and Wolfram Mathematica 11.3, as shown in Figure 5. The electrical side is based on PSS/E, including electrical network, GENROU generator, SCRX exciter, PSS2A stabiliser, and HYGOVM governor models which include detailed penstock, turbine, and governor dynamics. The river model is based on Mathematica, including the shallow-water-equation-based river model (13)- (14), solved by the finite difference method. The electrical and hydraulic sides communicate with each other via the PSSPY interface. The proposed controller (37) is modeled on Mathematica and solved by IBM ILOG Cplex via the NETLink interface. The step length of dynamic river simulation is 10 s, and the control period of AGC is 4 s.

Case setting
The real-life system located in Xiaojin County, Sichuan Province, China, shown in Figure 1, is used in this case study. Modelling parameters for simulation are specified in Appendix B. Specifically, data for the ten river sections are listed in We choose to test the proposed regulation method in the dry season, as in the wet season, there is always abundant water to generate enough electrical power. In contrast, in the dry season, water resources are limited, and the total hydropower generation cannot satisfy the load demands without solar power generation and load shedding. In this situation, the adjustability of the hydropower should be fully exploited, where the proposed coordinated regulation shows its value.
Specifically, for this case study, the upstream water inflow and its forecast are as shown in Figure 6(a), being much smaller than the rated turbine discharge of 30 to 40 m 3 ∕s of the whole plants each with 3 turbines. Thus, of each plant only one turbine is in operation. The pv and load power and their forecasts are respectively shown in Figures 6(b) and 6(c). Specifically, the pv data in Figure 6(b) is the real-time pv power collected in the Xiaojin system on May 30, 2018.
By examining historical data, we find the maximal volatility of pv is 11.59 MW/min, and in 99.9% of the operation time pv volatility is within ±6.99 MW/min. To alleviate such random volatility, the smoothing control introduced in Section 3.5 is adopted with result is demonstrated in Section 4.3.
In the MPC of the tertiary regulation, the prediction step length is set as T = 600 s. Since the water wave travels through the cascaded plants for one more hour, we set the prediction horizon of the controller to be two-times longer, that is 2 h or N = 12. The objective function is set as (23) with 1 = 10, 2 = 10, 3 = 1, 4 = 1, 000, and 5 = 10. The river operational constraints include water stage limits at the monitoring points 800 m upstream of the dams and plants. On natural river reaches and channels, the limits are ±0.2 m and ±0.5 m around the nominal point, respectively. The frequency reference limit is set as ±0.1 Hz.

Result of PV smoothing control
Simulation result of the smoothing control introduced in Section 3.5 with = 15 and = 3 is illustrated on different days of available pv power collected in the Xiaojin system shown in Figure 7. As can be seen, when the volatility is mild as shown in Figure 7(a), the curtailment is very small, and the pv output largely follows the prediction. In contrast, when intensive volatility occurs as Figures 7(b) and 7(c), the curtailment of pv is controlled according to the intensity of volatility, and the output follows a similar trend to the prediction but in a lower profile, which significantly reduces the variations of pv and therefore alleviates the challenge of frequency stability. Further, the smoothing effect and energy loss due to curtailment control are investigated based on the pv data for the whole year of 2018. The maximum, and 99.9% and 99% largest minutely volatility of pv, and the average daily energy loss with different are shown in Figure 8. Here, = 0 means no smoothing control. Given larger , the pv output is smoother but energy loss becomes larger. In practice, can be flexibly selected based on the robustness of the power system. In the case study we choose = 3.

Polynomial surrogate model for the power system dynamics
Using the polynomial approximation method presented in the Appendix A with the error tolerance = 0.01 in the Smolyak approximation algorithm, polynomial surrogate models to characterise the mean power of the three cascade plants over a time interval of 10 min are obtained. The normalised trajectories of solar power and load demand, that is P norm S (t ) and   Table 1.
As can be seen, the quadratic terms with respect to the control variables in the polynomial functions include c Gi ref and c Gi P sh D (or c Gi j after discretising the continuous variable of load shedding by substituting (8) in). By using the transformation (26)-(29) introduced in Section 3.2, these quadratic terms can be replaced by linear terms, making the optimal tertiary regulation model (37) an MIQP.
The polynomial surrogate models are verified numerically. After the simulation later performed in Section 4.6, the actual power outputs of the three cascaded plants and the values obtained by the polynomial surrogate models are compared in Figure 9. We can see that the marks lie close to the line y = x, indicating that the employed surrogate models have good accuracy. Further, the decent controller performance shown in Section 4.6 again verifies these models.
Moreover, we calculate the surrogate models using different normalised trajectories of P norm S (t ) and P norm D (t ), and the accuracy of the obtained models are almost the same and good as well. Output of the polynomial surrogate models shown in Table 1 versus the actual generation of the three cascaded plants, where each mark represents the mean generation over a 10-min interval  Figure 10. The decision period of load shedding is 10 min.

Benchmarking regulation methods and simulation
BM3: Mostly the same as BM2, expect that the frequency reference is set as ref = −0.1 Hz to minimise load energy consumption and to minimise load shedding as the result.
BM4: A receding-horizon generation scheduling with N = 12 and T = 600 s is used to optimise the nominal power references of the hydropower plants, the frequency reference, and the load shedding every 10 min. The traditional PI-based AGC with fixed participation factors calculates the incremental power references to stabilise frequency in real time within each period of 10 min. The scheduling model is similar to (37) but without consideration of the adjustment of the AGC participation factors.
The benchmarking methods BM1 to BM3 are based on common idea of power system regulation. Specially note that this paper proposed optimisation-based BM4 first. As briefly mentioned in the Introduction, to the best of the authors' knowledge, there is no publication on the tertiary regulation of islanded power systems with cascaded run-of-the-river hydropower. Due to space limit, the detailed model of BM4 is not presented here.
To quantify the performances of the different controllers, the following indices are defined:

1) Total loss of load (in MWh)
: 2) Root-mean-square frequency deviation (in Hz): By simulation, the performance indices of the benchmarking methods and the feasibility of the river stage constraints are listed in Table 2. In detail, because load shedding is not considered in BM1 and the energy use in stabilising the frequency is not appropriately allocated to the three cascaded plants, water used for power generation significantly excesses the river upstream inflow at the MP plant, and as the result the water stage at the MP plant descended to an unacceptably low value (< −1.0 m), as shown in Figure 11(a). This violates the river operation constraints and may exhaust the water storage; therefore, it is strictly forbidden in operation. The load shedding control ( Figure 10) in BM2 ensures that the water stage limits are not violated. However, from the water stage curves in Figure 11(b), we can see that the cascaded plants are not coordinated at all. This causes the adjustability of the cascaded hydropower not fully exploited, leading to the largest load shedding, that is 78 MWh, as shown in Table 2. Setting the frequency reference to the lower limit, that is −0.1 Hz, to reduce the load demand in BM3, load shedding still reaches 58.83 MWh.
BM4 coordinates the cascaded plants by scheduling the base power references every 10 min. Thus, the adjustability is significantly improved, indicated by the decreased load shedding shown in Table 2. However, as shown in Figure 12(a), the power outputs of the plants deviate from the scheduling within each 10-min interval due to the solar power and load volatility and the consequent control actions of the AGC. This portion of the power adjustment is not coordinated in BM4, in contrast to the proposed method, again revealed by the differences between the three water stage curves of different plants, as plotted in Figure 11(c). In other words, the adjustability of the hydropower can be further exploited.

Simulation result of the proposed regulation method
By simulation of the proposed method, the AGC participation factors, as the outputs of the proposed tertiary regulation method, are plotted in Figure 13. The system frequency reference and the actual deviation are shown in Figure 14. The power generations of the plants are shown in Figure 12(b), and a comparison of the total power generation and load shedding to the load demand is given in Figure 15. The water stages of the plants are given in Figure 11(d).
The overall performance indices of the proposed method are listed in Table 2 compared to those of the benchmarking methods. As can be seen, the proposed method further reduces the total load loss and frequency deviation compared all the benchmarking methods. The maximal transient frequency deviation over the whole day is +0.006 Hz and −0.262 Hz, which satisfies the system's operational margin, that is ±0.5 Hz. If we do not apply the pv smoothing curtailment control but taking the risk of instability caused by pv volatility, the total load shedding is 17 MWh. The result with smoothing curtailment control, that is 32 MWh, is to some extent larger. Alternatively, if we simply disconnect the pv plant from the grid, due to the power supply shortage, load shedding jumps to 148 MWh, which is very socioeconomically costly. This result proves that tracking the pv prediction curve by smoothing curtailment is a sweet-point tradeoff between frequency stability and power supply reliability. The following phenomena exemplify how the proposed method coordinates the cascaded plants. From Figures 13 and  11(d), at 10:00 am, the stage at the MP plant descends to the lower limit, causing the plant to lose its power adjustment ability. This is caused by the drop in the upstream inflow shown in Figure 6(a). In this situation, the proposed controller lowers Then, comparing Figures 11(c) and 11(d), especially the magnified parts, the difference between the three water stage curves of the cascaded plants becomes smaller than that under BM4. This outcome reveals that the proposed method better coordinates the cascaded plants. As a result, the load loss and frequency deviation are reduced, as compared in Table 2.
Moreover, comparing Figures 12(a) and 12(b), we can see that the trajectories of hydropower generation under the proposed method are smoother than those of BM4. This result again shows that the proposed method offers a better coordination of the cascaded plants over short time periods. Obviously, this improvement is achieved by jointly considering the multitimescale dynamics.
Finally, to analyse the impact of changing AGC participating factors on system frequency, we magnify the factors and frequency deviation between 13:00 and 15:00 in Figure 14 with the events of changing AGC participating factors labelled. We can see the frequency fluctuates around the reference smoothly, and no observable disruption is observed at the events of changing the participating factors. This can be attributed to the fact that the frequency a power system is mainly stabilised by the primary control, and the AGC is designed to removes all stationary-state frequency deviation of primary control as long as the adjusting abilities of the hydro plants are preserved, which is guaranteed by the inequality constraints in the proposed tertiary controller.

4.7
Numerical discussions 4.7.1 Impact of objective weights on control performance Section 3.1 briefly discusses on how to choose the weights in the objective function (23) of the proposed tertiary controller. Here we numerically demonstrate how these

FIGURE 16
Events of changing AGC participating factors and system frequency deviation between 13:00 and 15:00 under the proposed tertiary control weights affect the control performance. With the weights perturbed, the corresponding components of the objective and performance indices are summarised in Table 3. As we can see, with one weight increased, the corresponding objective component decreases but others increase. The result shows that in engineering practice, the tradeoff between different objective components can be flexible adjusted by these weights.
Note that despite given different objective weight settings, the frequency stability of the system is not impacted. This is attributed to the fact that frequency stability is guaranteed by the primary and secondary control. With inequality constraints, the proposed tertiary control maintains the primary and secondary frequency control ability for in the longer time scale, which is not affected by choosing different objective weight factors.

4.7.2
Robustness of the proposed controller against perturbed modeling parameters To test the robustness of the proposed control method against inaccurate system model, we independently perturb the modeling parameters in the MPC (including river section geometry and hydropower production efficiency) by ±5% with uniform distributions. Overall 100 sets of perturbed parameters are created and used for simulation. Except for the parameters in MPC formulation, system parameters for simulation remain untouched. Empirical densities of the obtained performance indices of load shedding and RMS frequency deviation are shown in Figure 17. Despite significantly disturbed parameters, the performance of the proposed control method remains relatively consistent. This can also be attributed to the inherent robustness of the MPC.

Simulation results of the proposed regulation method under various scenarios
To validate the improvement in the proposed method more comprehensively, especially considering various operational conditions in a real-life system, we tested it under 100 different scenarios. Each scenario represents one day from 6:00 to 22:00. The solar power curves used are empirical data recorded in the Xiaojin system from April to July, 2018. For visualisation, 20 out of the 100 scenarios used for simulation are plotted in Figure 18. Specifically, the data of available pv power shown in Figure 18(b) are recorded on 100 consecutive days from January to April, 2018. Smoothing control in Section 3.5 is applied in simulation.
Simulation results of the proposed method are compared to those of BM3 and BM4. The differences between the total load loss and frequency deviation of the proposed method and those of BM3 and BM4 are respectively plotted as one mark per scenario in Figures 19 and 20. Seen from Figure 19, compared to BM3 the proposed controller performers dramatically better in terms of reducing total load loss. Comparing the proposed method and BM4, as observed in Figure 20, most points appear on the left halfplane, meaning that in most scenarios the proposed method still outperforms BM4 in terms of load loss while keeping frequency deviation almost constant.
Based on all the results above, we can conclude that the proposed regulation method that considers multi-timescale dynamics better coordinates the cascaded plants and exploits the overall power adjustability compared to the benchmarking methods, and thus improves the power supply ability of islanded power systems with cascaded run-of-the-river hydropower.

CONCLUSIONS
A coordinated tertiary regulation approach for islanded power systems with cascaded run-of-the-river hydropower and volatile generations is proposed. An MPC is established to dynamically adjust the AGC participation factors to coordinate the cascaded plants. A simulation of a real-life system shows that the proposed regulation method significantly reduces load loss compared to that yielded by the other regulation methods. Currently, the proposed method has not taken the statistical characteristics of solar and wind power into modeling and optimisation, but recent work has shown that considering stochastic characteristics improves the performance [20]. In future studies, taking the stochastic characteristics of the volatile generations and load into the tertiary regulation could be a promising work direction.   (7) are used to approximate the mean hydropower outputs over the time interval of t ∈ [0, T ], which are exemplified in Table 1. This appendix elaborates on how to obtain the polynomial surrogate models. Because Legendre polynomials are orthogonal over a finite input interval [28,29], leading to the optimal approximation in terms of L 2 -norm, this work chooses Legendre polynomials as the basis in approximation.
The j th-order Legendre polynomial P j (x) for x is defined as which admit the orthogonality over interval = [−1, 1], as where ⟨⋅, ⋅⟩ L 2 ( ) represents the inner product over ; jk is the Kronecker delta function.
Specifically, the Legendre polynomials up to order 2 are In finding the surrogate models for the power system dynamics in Section 2.1, the function to be approximated, that isP Gi (⋅), has multidimensional inputs, denoted as vector q, We first rescale all inputs into interval = [−1, 1], and denote the dimension of q as m. Given the approximation order N i for each scaler input q i , the multivariate Legendre polynomial basis { j (q)} for the vector input q is constructed as the following tensor product, { j (q) where j ≜ [ j 1 , j 2 , … , j M ] T is the multi-dimensional index. Then, the functionP Gi (q) is approximated as the following polynomial expansion, denotedP * Gi (q), as P Gi (q) ≈P * Gi (q) ∶= where f j are coefficients to be calculated.
Generally, the Galerkin method (GM) or the collocation method (CM) can be used to find the coefficients of a polynomial approximation [28][29][30]. Because the CM can treat the function to be approximated, that isP Gi (q), as a black box in a simulator such as the PSS/E whereas the GM relies on a white-box model with detailed mathematical formulation, we adopt the CM to find the coefficients in the polynomial surrogate model in this work.
First, define the collocation point set, as know as the Gaussian quadrature point set [28,30], to contain all the zeros of the product of the (N i + 1)th Legendre polynomials of the inputs in q, as where N b is the size of collocation point set, which equals the size of the basis (A.3).
For each collocation pointq i in (A.5), using (4) to create the corresponding trajectories of solar power and load demand and set the tertiary control variables according toq i . Then, using power system simulation software to find the values of P Gi (q i ). After the simulations of all collocation points are finished, the coefficients in the polynomial approximation model (A.5) are computed by where A is a constant matrix determined by the basis (reindexed ) and the collocation points (A.6), as Alternatively, due to that the high-order cross terms in the polynomial expression generally have less impact on the precision of the approximation [28], we can use the Smolyak adaptive sparse algorithm [31] to incrementally add higher-order cross terms in an adaptive way instead of adding all high-order cross terms at once. The precision of the approximation in terms of L 2 -norm is monitored by the iterative changes in the coefficients, and only the terms with a noticeable impact are added. This process terminates when the preset error tolerance is reached. Interested readers are refereed to the detailed theory of adaptive polynomial approximation in [31] as well as its applications in [32].