Sustainable Energy, Grids and Networks

.


Introduction
Due to constantly stricter regulations towards increasing renewable penetration, maintaining a stable and efficient operation for isolated power systems has become challenging [1]. Based on recent advances in offshore technologies [2], more and more such isolated grids with high renewable penetration are about to be realized [3], and this highlights the need for efficient controls.
The stochastic nature of renewable energy sources (RES) and the low system inertia, induced by replacing traditional synchronous machines with converter interfaced generators, makes the operation of such grids vulnerable to large frequency variations [4][5][6]. A common and accepted way to mitigate such issues is the integration of energy storage systems (ESS) such as battery energy storage systems (BESS) [4,[7][8][9][10][11][12]. However, increasing the number of controllable components in the grid complicates the multi-input multi-output (MIMO) nature of these systems, exacerbating the chances of encountering multiple contradictory objectives, additional constraints due to safety, and optimal operation requirements. On top of that, uncertainty poses an additional level of difficulty for controlling such grids [13].
In [14,15] the fundamental methods to provide frequency control in isolated power systems by controlling the generator units, are explained in detail. Two main categories are identified, one being the well-known traditional droop control and the other the isochronous control. With the first, generators modify the  active power provision in correspondence to the system's frequency deviation and based on the set droop curve, whereas in isochronous mode the rotational speed and thus the frequency of a single or master unit are tightly regulated to their reference values. It is noteworthy that even though droop control facilitates the load sharing among several units, the frequency will not be restored to its nominal value. To achieve that in the absence of interconnection, frequency restoration will depend only on the local generator capabilities and requires that one of them be in isochronous mode. In other words, from a control perspective, droop control corresponds to a proportional controller whereas isochronous operation corresponds to controllers with integral action. In further detail and from a system dynamics perspective, droop control is the main source of damping, necessary for the stability of the system and isochronous control the term encapsulating error accumulation effects and responsible for achieving zero steady state deviation. For the isochronous case, a frequency variation will be immediately compensated by the generating unit, providing a service equivalent to the secondary frequency restoration in large interconnected systems. However, in contrast to the isochronous operation of a generator for an isolated system, in large interconnected systems such a service is achieved by the system operator, assigning different activation signals and participation factors to all available generators participating in the FRR (Frequency Restoration Reserve) market. Such services are typically activated after the provision of primary control to stabilize the system after a disturbance and on the scale of several seconds to minutes. We also highlight that the actual speed of frequency restoration will depend on the amount of power imbalance and the operational limits of the involved units (ramp constraints and technical maximum).
In [16] gas turbines dynamics are studied along with power systems frequency dynamics, demonstrating the fast-acting behavior of such generators and the isochronous operation capability (response time in the range of seconds). Such a practical case study for an islanded plant (industrial isolated power system) is also provided in [17].
Various types of advanced control methods for load frequency control of isolated power systems have been proposed in literature. In [18] optimal PID controllers were designed to provide load frequency control in microgrid clusters, in [19] a swarmoptimized fuzzy logic was used for robust secondary frequency control of islanded systems, in [20][21][22] the H ∞ synthesis concept was used to provide robust secondary frequency control in islanded microgrids, in [23] an MPC controller was investigated to provide load frequency control, in [1] a Grey Wolf Optimization was used to provide frequency support by a BESS for an island power system, in [24] a LQR stochastic based control was proposed to provide secondary frequency regulation in an independent microgrid, in [25,26] fractional order MPC and PID were designed correspondingly to provide frequency control for an islanded microgrids or for single area power systems, in [27] an IMC-PID design was investigated for load frequency control.
In the aforementioned studies, the (isolated) power system is modeled as a continuous time linear time-invariant dynamical system and the control performance is typically assessed through load step variations considering either average or worst cases disturbances. However, the actual empirical distributions are not derived, possibly leading to biased/over-conservative or rather optimistic conclusions. Nevertheless, all the proposed load frequency controllers include an integral part that brings the steady state error to zero in just a few seconds. Such action resembles the isochronous operation provided by single generators when speed and thus frequency are tightly regulated in stand-alone isolated systems. The system's response time highly depends on the various time constants involved (primarily the inertia) and the operational limits (power provision/absorption capability) of the related units. Therefore, for low inertia isolated power systems, response can be faster than large interconnected systems. We note that when secondary level control is designed and assessed (frequency restoration), primary control loops (i.e., droop) can be included in the dynamics of the power system, by modifying the system's damping.
Different approaches can be found in literature for dealing with uncertainty in frequency control of isolated grids. In [7] a standard model predictive control (MPC) algorithm based on a convex QP problem is employed to control an isolated power system containing critical and non-critical loads, diesel generators, BESS, and RES in the form of photovoltaics and wind power generation. The results demonstrates how MPC can effectively manage several objectives, like preserving power balance in the grid and reducing the fuel consumption of the diesel generators. However, this work integrates uncertainty in a simplistic way, by using scenarios with different level of accuracy for the load profile predictions ranging from a perfect forecast and up to approximately 10% deviation. In [10] authors develop a scheduling algorithm for an isolated power system with high penetration of RES which controls the energy production from fossil sources and the power transactions with the main grid in order to maintain power balance and maximize the RES penetration. This yields an alternative strategy to the MPC algorithm and shows good simulation results with a time frame of a day with uncertain forecasts of wind speeds and load profiles. In [28] a stochastic MPC (SMPC) approach that uses scenarios-based description of uncertainties is employed to optimize the fossil energy production and power transactions in the real time market. Such scenarios are generated from a scenario tree designed explicitly to capture the additive feature of uncertainty and avoid infeasibility. The authors of this work compare their algorithm to a so-called prescient optimal control strategy that assumes perfect knowledge of future realizations of the uncertainty, and a certainty equivalent MPC where the uncertain parameters are substituted by average values identified from historical data. Results demonstrate the superiority of the scenario based approach in decreasing costs. [29] performs a comparative study of what they call SMPC and scenario MPC (SCMPC). The difference between the two methods is that the SMPC method converts the probabilistic constraints into deterministic ones using knowledge of the co-variance of the random variables and their propagation along the prediction horizon; the SCMPC approach computes scenarios and forms a scenario tree from the probability distribution of the uncertain variables, as in [28]. The results show that SCMPC generates more realistic scenarios than SMPC because it uses information gathered online to adjust scenario predictions rather than exploiting knowledge on second moments. This accuracy comes though at the cost of increased computational requirements. In [30], the authors develop a scenario SMPC for hybrid vehicles with the goal of improving fuel efficiency while obeying constraints on the state of charge (SoC) of the battery and the power exchanged with it. However, the SCMPC is modified to only generate scenarios that are feasible and likely. The disturbance, i.e., the power requested by the driver, is instead estimated via a Markov chain that predicts the future driver inputs by learning the previous request pattern in real time. Results from [30] show that their SCMPC with learning may outperform classical MPC formulations, and that in many simulations the SMPC performs almost as an MPC with perfect knowledge of future realizations of the disturbances. Other methods consider then a robust control design approach to mitigate uncertain frequency variations [20]. For example, in [31] the authors apply both an H ∞ and µ synthesis approaches for robust frequency control in islanded microgrids. Results show that the two control algorithms may outperform an ''optimally tuned'' PID controller in the presence of structured uncertainty in the form of wind power generation, solar power generation and uncertain load conditions. The review above highlights the presence of a plethora of different MPC based methods for either optimal scheduling in isolated power systems or for frequency regulation under uncertainty, either based on robust or stochastic control algorithms. To the best of our knowledge there is though a lack of publicly available studies that analyze how to integrate both tasks into the same control strategy, and thus analyze how control strategies may satisfy optimal schedule tracking while simultaneously ensuring continuous, tight, and fast frequency regulation. The goals of this paper are thus two: 1. propose a series of control formulations and parametrizations (stochastic, robust and deterministic in primis) that are all capable of coordinating a master gas turbine and an energy storage to achieve fast frequency regulation for the local grid, while the gas turbine can still operate with minimum deviation from its optimal loading and the storage system can follow a pre-scheduled reference state of charge trajectory; 2. understand which one best handles uncertainties while guaranteeing control performance, and their tradeoffs.
For this purpose, the paper is organized as follows: Section 2 introduces the dynamic MIMO model of the isolated power system, Section 3 describes the design of the controllers we compare, and Section 4 performs the in silico analyses that lead to the main messages given by the paper, i.e., the comparative analysis of the proposed controllers. Lastly, Section 5 draws some final conclusions.

Power system modeling
We consider an isolated power system that includes a gas turbine generator in isochronous mode, a wind turbine generator, a battery energy storage system (BESS), and an aggregated load. The gas turbine is mainly responsible to stabilize the system by arresting the frequency deviation after a disturbance and simultaneously cooperate with the BESS to restore it to its nominal value. The following subsections describe then each of these elements in details, assuming the state of the system to be a sixdimensional vector where x 1 and x 2 are states related to the gas turbine (governor and turbine subsystems respectively), x 3 is the power deviation coming form the battery storage system, x 4 is the grid's frequency variation, and x 5 , x 6 are the states describing the internal dynamics of the batteries (relaxation and rate capacity effects).

Modeling of the grid dynamics
Given the scope of this paper, we choose to model the dynamics of the isolated power system as a first order transfer function from power balance to frequency deviation [14,15,32]. The transfer function, derived from the swing equation, is thus where M is the inertia constant (due to the generator's rotating mass), D is the damping constant (which encapsulates the combined effects of primary control layer (droop) and load damping, and ∆P g and ∆P ℓ the generated and consumed power deviations with respect to the operating point, respectively.

Gas turbine dynamics
To capture the dynamics of both the turbine and its governor, we model the dispatchable generation system as two first order low pass filters connected in series, giving the equivalent system where T g is the governor time constant, T t the turbine time constant, u 1 = P gt the power command to the gas turbine, relatively to its steady state operating point u 1|t=0 = ∆P * gt = 0. The input to Eq. (2) is defined as u 1|t = u 1|t−1 +∆u 1|t , from which we see that the part associated to the small signal dynamics of the gas turbine is ∆u 1 .

Wind turbine dynamics
The electromechanical conversion on the wind turbine generator is modeled as a first order filter [14,15] from the input (mechanical wind power) to the output (electric power injected into the grid). More precisely, we assume where T wt is the wind turbine generator time constant, ∆P w (s) is the uncertain wind power variation and w 1 (s) = ∆P wt (s) is the corresponding uncertain electric power injections from the wind turbine generator. As for modeling the stochastic process corresponding to the input, or in other words to implement the disturbance model used in our case, we consider wind speed scenarios from a wind speed generator that creates, starting from an average wind speed given as parameter, a realistic set of wind speed samples by means of a physics-driven model of the hydrodynamic effects occurring locally around wind turbines and rotor blades. Such a wind speed samples generator makes use of normal random variables and Kaimal distributions that better capture the small time scale wind intermittency. More details can be found in [33,34].

Battery energy storage system dynamics
To model the dynamics of the power converter interfacing the battery energy storage system (BESS) with the grid and capture the internal dynamics of the battery cells during charging and discharging processes, we consider the system where SoC = x 5 + x 6 is the state of charge, u 2 = P b = u 2|t=0 + ∆P b the reference power to the BESS (since u 2|t=0 = ∆P * b = 0 at the steady state), T B the time constant related to the power conversion and c r , c w the coefficients related to the linear modified KiBaM model, illustrated in Fig. 1 and which we assume as sufficiently detailed for our purposes [35][36][37]. Here it follows a detailed derivation of the proposed simplified battery model based on the modified KiBAM model. The battery charge dynamics (rate capacity and charge relaxation effects) can be approximated by an equivalent dynamic system of two interconnected water tanks with different volumes. If Q b is the total charge of the battery at full capacity, Q 1 and Q 2 the total charge of tank 1 and 2 and their widths c w and 1 − c w correspondingly we have that h 1 = q 1 cw and h 2 = q 2 1−cw , where q 1 and q 2 are normalized variables defined as and h 1 , h 2 represent the normalized water column heights (head) of each tank. Then, considering that the flow across the valve is proportional to the head difference between the two tanks, we can writė where c ′ r is the valve's coefficient and Remembering that the output current is defined as I b =Q b and based on the above definitions we can write the tanks system equations as which can be compactly expressed as where Then, under the assumption of constant (average) open circuit voltageV oc we can express the battery power as P b =V oc I b and that the battery energy capacity can be written as C b = Q bVoc [36,37] we have Finally, neglecting the charge/discharge efficiencies and considering a series-connected first-order delay for the power conversion stage (see Fig. 1) we end up in Eq. (4) where battery power −P b ≤ P b ≤P b is the system input.

Dynamics of the whole interconnected system
Combining subsystems Eqs. (1)-(4) as depicted in Fig. 2, letting x 4 = ∆f , and renaming the uncertain power demand as w 2 = ∆P ℓ , Eq. (1) can be rewritten aṡ This means that the system's model can be written in the state space forṁ From a block-scheme perspective, the overall system is thus modeled as in Fig. 2.

Deterministic model predictive control
To achieve optimal operation and reduced fuel usage in isolated power systems, supervisory power management is typically used to deliver the optimal scheduling set points to the local controllers of the different subsystems, such as the BESS and the Gas Turbines [8,38]. In this study we consider integrating the local control objectives of these subsystems along with the frequency regulation of the isolated power system. Local objectives mean that the gas turbine can still operate with minimum deviation from its optimal loading and provide primary grid stabilization while the BESS can still follow a SoC reference trajectory coming from a tertiary dispatch level, and at the same time work together to restore and tightly regulate the system's frequency. For this purpose we propose using a Model Predictive multi-objective multi-input multi-output (MIMO) control that, on top of the basic requirements, aims at minimizing the control effort so to promote reduced actuator wear, and cycling energy storage so to promote longer battery lifetime. To integrate the several control objectives considering the state space model developed Section 2 into a MPC formulation we first use a recursive elimination approach to express the states for the selected prediction horizon N p as where x 0 is the initial condition, is the initial operating point of the subsystems and X , ∆U , W , A e , B 0 , B e , and H e are defined in Appendix A. Hence, we define the following qualitative control objectives: • minimize the frequency deviation (to achieve fast regulation purposes); • minimize the amplitude of the control signal (to minimize the fuel consumption associated to gas turbine usage); • perform an optimal system operation (to follow the reference schedule); • reduce BESS degradation (to minimize replacement costs).
We then translate the above qualitative targets into quantitative ones as follows: first, penalize the states with the Q -norm x T Q x, with Q diagonal and positive definite, so to promote small frequency deviations. Then penalize deviations from operating the gas turbine at its maximum efficiency, and thus penalizing ∆u 1 through the norm ∆u T R∆u to minimize the gas turbine fuel consumption. To follow the reference schedule, penalize the deviations of the state of charge from the reference value using the affine plus quadratic cost ℓ ( Penalizing the BESS degradation can be then promoted in several ways.
Typically, the degradation of a battery is modeled as caused by two distinct effects [39], namely the calendar aging and the cycling effect. Since the cycling degradation is dependent on the number of cycles and the depth of the cycles, a common approach is to discourage cycling by minimizing the standard deviation of the state of charge for the prediction horizon as Then, by considering the augmented states X and control variables U for the prediction and control horizons (N p and N c respectively), and expressing the evolution of state of charge as X soc = SX , we can formulate the objective function for the finite horizon optimal control problem as J in Eq. (15) is convex in X and ∆U by construction (since a sum of basic convex functions. The constraints associated to the problem of optimizing J shall then include the physical limits of the components (like the allowable BESS SoC range, minimum and maximum governor opening, BESS power limits, and ramp rates for the changes of the manipulated variables), and the maximum allowable deviations from the nominal frequency and SoC. Summarizing, the constraints can be expressed mathematically as (17) ∆U min ⪯ ∆U ⪯ ∆U max (18) where P and c contain information of the hard limits on the states and are presented in Appendix A. Here we note that in practice, different battery cell types and technologies would have different C-rate limitations, resulting in different and tighter bounds on the allowed charge/discharge power. However, such limits can be directly integrated in Eqs. (17) and (18) without affecting the formulation of the problem proposed in this paper. To simplify our analysis while preserving generality, and since the focus of the paper is not on the comparison of different battery types and technologies, we chose to restrict charge/discharge power based on the rated power, as commonly done. This is a simplifying assumption which eventually does not alter the methodology proposed in this paper or the comparative analysis presented later, since all comparisons are performed under the same generic battery model and same constraints. In this way we intend to provide upper theoretical bounds on the performance, while further case specific studies are needed depending on the selected battery technology. As a complement, an additional sensitivity analysis focusing on the different factors that affect the resulting charge/discharge rate is provided in Section 4.5. Eqs. (14)- (18) define then the deterministic model predictive controller (DMPC ) that will then be compared against the stochastic one defined in the next subsection.

Stochastic model predictive control
One of the objectives of this manuscript is to investigate which strategy best handles uncertainties while guaranteeing control performance. To account for uncertainties in the control design we then adopt the commonly used scenario approach (SMPC ) [28][29][30]40,41], where the cost Eq. (15) from the DMPC scheme is replaced by an expectation over the possible outcomes defined from the scenarios. The implementation of such cost function is done through the sample average approximation (SAA) [42] as Where N is the number of scenarios considered. As explained in [41], different constraints are induced by Eqs. (16)-(18) for each random realization of the uncertainty from the randomly sampled uncertainty set. Those can be described as We highlight, that there is not specific requirement on the distribution of the random variables, or the disturbance model used for the scenario approach, as long as number of samples used is chosen based on the main theorem [43] and the same risk criteria. Another interesting feature related to the SMPC is the different choices for the parametrization of the control input [41], an aspect which is also investigated in this study. In this paper we thus employ different commonly reported controller parametrizations for the SMPC formulation Eqs. (19)- (22), and compare their performance in terms of handling uncertainties and system performance. More specifically, we consider the following parametrizations:

Open loop policy (no parametrization)
The direct open loop policy is also considered as a possible choice for the controller, meaning that the direct output of the solution of the numerical optimization problem Eqs. (19)-(22) is applied to the plant, i.e., The non-parameterized control action (SMPC-NP) is also included as a benchmark alternative to better understand and illustrate the effects of the control action parametrization, relatively to a less constrained open loop solution. Here we specify, that still under this parametrization, the controller is of the standard rolling horizon MPC type, which means that feedback of the system's response is given back to the controller in the form of the constantly updated initial condition for each step's optimization problem [44,45]. In that way, no substantial error is accumulated over time. This concept is illustrated in Fig. 3

Affine disturbance feedback
The disturbance feedback full parametrization (SMPC-FP) adjusts the control action directly based on past values of disturbance realizations. By using this strategy, the disturbance history is fed into the controller while preserving the convexity of Eq. (19). More in details, such full parametrization is defined as With such a parametrization, we include constant terms γ i and terms θ i,j proportional to past disturbance realizations w j , for the specified control horizon. In this way the resulting optimization problem is convex with respect to the control parameters which are iteratively updated on each MPC iteration. Note that, by reorganizing this into a recursive eliminated vector form, and by integrating the number of states, control actions, and disturbances, one obtains U = Γ + ΘW (25) where U , W , Γ , Θ are defined as in Appendix A.

State feedback
Another commonly used controller parametrization for state space formulations is the state feedback parametrization (SMPC-SF ), defined as where K lqr is a fixed constant. In contrast with SMPC-FP, this particular parametrization suffers from the fact that the feedback gain parameter cannot be treated as an optimization variable, since this would result in a non-convex problem [41]. To avoid this issue, K lqr may be chosen to be the infinite horizon optimal gain (see [41] for more details).

Robust control
To compare the stochastic control approach described in Section 3.2 against a robust control design approach, we also design a mixed synthesis H ∞ controller for the system Eqs. (10) and (11). More precisely we design the H ∞ controller by applying the principles of loop shaping, designing appropriate weights W p and W d for acceptable nominal performance and robust stability for the lumped unstructured output multiplicative uncertainty, and finally solving a linear matrix inequality (LMI) optimization problem that minimizes the cost where W p , Wu, and W d are the disturbance-rejection performance, controller effort and robustness weights respectively. Selecting proper weights for W p (s), W d (s) allows then to upper bound the sensitivity S(s) and complementary sensitivity T (s) functions, finding an overall satisfactory controller K (s). This paper follows the basic guidelines from [46] on how to select these weights for robust disturbance rejection; after an iterative procedure we set   ] .
From Fig. 4, we notice how the robustness weight selected as W d (s) = 0.32s + 0.018 s + 1.77 (31) bounds all realizations of the relative model error of the uncertain plant, validating the required specification for a robust controller design. Note that this type of robust control only considers parametric uncertainty being integrated as output multiplicative.

Simulation results and analysis
To analyze the effect of the proposed multi-objective MIMO control strategies, we compare the various control designs described in Sections 3.1 to 3.3 in terms of their dynamic response and capabilities of handling uncertainty. We note that our proposed formulations adopt the simplified but generic battery model described in Eq. (4) where no particular technologydependent maximum C-rates are enforced besides the power rating limits (Eqs. (17) and (18)). Since the focus of the study is on the comparative analysis and not on the different technologies' comparison, for the purpose of the analysis presented in Sections 4.1 to 4.4 we therefore assume no further technologydependent battery power constraints, but we ensure the same model and limits for all controllers, for a fair comparison. However, a detailed discussion on the various factors affecting the battery charge/discharge rate and the resulting C-rate is given in Section 4.5. In addition, we note that for the nominal case, the controller assumes no mismatch to the plant model, as to have a benchmark upper bound performance, for comparison with the proposed controllers that assume no knowledge of the nominal plant by considering both the uncertain plant parameters and the unmodeled dynamics. The following key performance indices are defined on the output signals to serve the relative comparison: • Max frequency deviation: • Average (expected) frequency deviation: • Total fuel usage: • Max SoC deviation: • Average (expected) SoC deviation: • SoC standard deviation (storage cycling indicator): where N sim is the number of simulated time steps.

Nominal dynamic performance
Initially, we compare the performance of designs DMPC, SMPC-FP, SMPC-SF, SMPC-NP, H ∞ at the nominal conditions, that is for the nominal set of plant parameters (see Table 1) and for a deterministic knowledge of the disturbance signal (i.e., assuming the uncertain wind power generation to be constant). This means that in this specific case the different versions of the SMPC controllers consider only one scenario implying that the DMPC, SMPC-FP, SMPC-NP formulations are in this case equivalent. However, the design of the SMPC-SF is different since it depends on the state feedback gain K lqr which is calculated independently of the uncertainty realization, thus resulting in a different performance.
The comparison of the dynamic performance of the various designs is illustrated through Fig. 5. In particular, from Fig. 5(a) we can observe the frequency response of the isolated power system, from Fig. 5(b) the SoC regulation to the scheduled value, and from Fig. 5(c) the control effort required from the two subsystems (BESS and gas turbine) to achieve the multiple objectives. From Fig. 5(a) it is clear that the DMPC/SMPC-FP/SMPC-NP design achieves the fastest regulation with the smallest peak deviation, the SMPC-SF design has the smallest damping and the largest peak deviation while the H ∞ design results in a slower but less oscillatory response. Then, from Fig. 5(b) we observe that H ∞ achieves the best tracking of the reference value with negligible dynamics, while DMPC/SMPC-FP/SMPC-NP once more performs better than SMPC-SF. However, the performance of H ∞ can be justified by Fig. 5(c) where we can observe the unbalanced use of gas turbine and battery power to achieve the multiple goals. As it is observed, DMPC/SMPC-FP/SMPC-NP and SMPC-SF split the power usage between the two units whereas H ∞ achieves the same goals by using almost exclusively the gas turbine in a way that avoids the power overshoot in contrast with the others. This sub-optimal response is associated with the tuning of the performance weights, a difficulty which inherently associated with the mixed sensitivity design procedure, since the calculation of the most appropriate weights is a notoriously challenging and tedious task.

Mixed model and disturbance uncertainty
To further compare the performance of the controllers introduced in Section 4.1 for the nominal conditions, we include a 10% uncertainty in the model parameters, so to capture imperfect knowledge one may have when identifying the system plant parameters. Dynamic simulations were performed for random values of the plant parameters and the corresponding responses as in Fig. 5 where recorded for each design: DMPC, SMPC-FP, SMPC-SF, SMPC-NP. The dynamic responses of the system under the various controllers are derived for random realizations of the uncertain parameters and disturbances and compared. By observing Figs. 6 to 10 we see the following patterns:

Monte-Carlo simulations and constraint violation
To assess the ability of the stochastic controllers SMPC-FP, SMPC-SF, SMPC-NP in integrating all kinds of uncertainties under the scenario approach, the various designs from Section 4.1 were compared under a Monte-Carlo simulation framework [23,47]. By performing numerous simulations and considering both parametric and disturbance uncertainty, the empirical distributions of the constrain violations and of the indices Eqs. (32)-(37) were calculated.
We note that from the scenario approach, the constraint satisfaction is meant on a probabilistic sense, that is by using specific number of scenarios we can have theoretical guarantees (bound) on the maximum violation probability when the original distribution is used. To demonstrate this effect, the SMPC-FP design was considered with different values of scenarios (small number: 250, big number: 750) and corresponding guaranteed violation probabilities. Increasing the number of scenarios generally results in a smaller violation probability, a fact which is numerically validated through Fig. 11 where we observe that for larger number of scenarios the constraint violation probability decreases. From the same figure it is also remarkable to notice the superiority of stochastic control over the deterministic version (DMPC ) which is always associated with a higher probability of constraint violation. In addition, we can observe that the combined effect of disturbance uncertainty and 10% parametric uncertainty has a greater impact on the max frequency deviation constraint Fig. 11(a) compared to the max SoC deviation constraint Fig. 11(b). From the latter we can see that even though enough cases would imply the maximum allowable deviation, none of them would cause an actual violation. Then, for the large number of scenarios (750) we compared the max constraint violation probabilities for the different controllers (deterministic, stochastic with different parametrizations, and robust). From Fig. 12 we observe that SMPC-FP and SMPC-NP achieve lower constraint violations for the max frequency and SoC deviations compared to DMPC whereas DMPC performs better for the max frequency constraint than SMPC-SF. This is an interesting result demonstrating (i) the incapability of state feedback to adapt its control law based on disturbance information and (ii) the unexpectedly good performance of the deterministic MPC, considering an expected behavior as the deterministic equivalent. However, the trends seem to be different for the max SoC deviation, where SMPC-FP and SMPC-SF seem to perform better than SMPC-NP and DMPC, meaning that the state feedback parametrization has a better capability of handling the parametric uncertainty related to the SoC variations relatively to DMPC. Finally, from Figs. 12(a) and 12(b) we observe that the robust design is associated with the highest violation probability in max frequency deviations and the smallest in SoC deviation respectively. This fact is in agreement with the fact that H ∞ considers the parametric uncertainty, minimizing the impact on the SoC deviation but does not considers disturbances that could cause higher frequency variations.

Expected performance and operation
Finally, the different control designs were compared in terms of their average constraint violation performance, the fuel usage of the gas turbine and the battery cycling. In particular from Fig. 13 we can observe the empirical cumulative distributions of the average frequency deviation defined in Eq. (33) and the average SoC deviation as defined in Eq. (36). From Fig. 13(a) we can see that SMPC-FP, SMPC-NP, DMPC significantly outperform the SMPC-SF, H ∞ , with DMPC having a bit more violations compared to the best stochastic designs. A similar trend is observed from Fig. 13(b), from which we can see that the average performance of the state feedback control is much worse compared to the other stochastic controllers. From these figures we can also see once more the pattern associated to robust control, where the average frequency violations are even worse than the worst stochastic design while the average SoC deviations are much better compared to the rest.
From Fig. 14 we can observe the effect of each controller on the battery degradation levels. As expected, since the robust control makes minimum use of the storage system, its degradation is the smallest with significant difference from the rest, while for the other controllers we can see again that SMPC-FP, SMPC-NP designs have a better degradation performance than SMPC-SF, DMPC with SMPC-SF again being the worst. However, interestingly the trends reverse for the total fuel usage defined in Eq. (34), where from Fig. 15 we can see that the SMPC-SF design is associated with the smallest total fuel usage probability compared to the other designs and the H ∞ is associated with the largest one with a big difference. The other designs (stochastic and deterministic) all have very similar behavior in terms of the total units of fuel distributions.

Effect of BESS related parameters on the dynamic performance
From the analysis above it is clear that different control methods and parametrizations lead to different dynamic performance of the system. However, it is not only these two factors that affect the system's behavior. Different weighting of the multiple objectives of the MPC optimization problem or different battery characteristics, such as its size, affect the resulting behaviors. Under this perspective, it is worth mentioning that the resulting behaviors shown in the previous sections where all derived for the same weightings and BESS size. The latter was set equal to unity (see Appendix A, Table 1) in order to consider a general case and in accordance with a per unit representation of the power deviations associated with the system operation (see Appendix A). Therefore, the proposed control methodology can be built upon a generic idealized system to remain as general as possible, since the target is not to provide case-specific results but a general control design framework where the control design procedure should be the same irrespectively of the battery capacity. . Even though the BESS capacity and the power rating are design parameters, whose selection is beyond the scope of this study but requires detailed analyses considering techno-economical perspectives, it is useful to study the effect of those parameters under the proposed MPC framework, for the nominal conditions. For this reason, we performed three individual sensitivity analyses which are explained in detail below. For all the following analyses we follow the same color convention where red represents the SoC deviation, blue the BESS power deviation and green the system's frequency deviation.

Sensitivity analysis for the SoC weighting
First, the effect of the SoC deviation weight ℓ was studied. We performed various simulations (load step disturbance) with increasing values ℓ, giving more significance to large state of charge deviations and thus better constraining the discharge/charge Crates of the battery. The results are illustrated in Fig. 16. From   Fig. 16(a) we observe that increasing ℓ not only leads to less SoC variation but also to decreasing discharging/charging rates, which can be of great importance when specific limits for the Crates of the battery have to be respected. For example, selecting the maximum value of ℓ = 80 and considering the battery's discharging at the load step (approximately 2% in 5 s, red curve in Fig. 16(a)) we can estimate a C-rate of around 0.1440, which is much lower than the one for ℓ = 2. However, putting a lot of weight to the SoC deviation comes at the cost of a higher overshoot at the systems frequency response, as can be observed    from Fig. 16(b) (bold green line). From this analysis, we can conclude that a value of ℓ = 20 would give a fair trade-off between the SoC deviation penalization and a realistic C-rate for the battery charging/discharging process. For this reason, we kept this value for the next analyses.

Sensitivity analysis for the BESS capacity
Another very important factor affecting the resulting charge/ discharge rates is the total available energy capacity of the battery system. It is expected that for larger capacities and same discharge levels, the C-rates will be decreasing. Again we performed various simulations (load step disturbance) where the BESS per unit scaled capacity C b is changing in multiples of the nominal one (C b = 1). In Figs. 17(a) and 17(b) we observe the resulting trajectories of the SoC, the BESS power and the system's frequency response. We see that increasing the battery size, the minimum value of SoC is increasing, leading to lower peak Crate and higher minimum frequency deviation, even though the frequency response is a less dumped. In other words, increasing the battery size, the system can tolerate higher discharge power rates, resulting in lower C-rates for the battery and improved frequency nadir.

Sensitivity analysis for the BESS power to energy ratio
Finally, the effect of the maximum power to energy ratio was studied by varying the λ parameter from 1 to 0.1. This means that a BESS with λ = 1 has 10 times higher power provision capability than the one for λ = 0.1, for the given (same) energy capacity. In other words, λ reflects the effect of further restricting the battery charge/discharge powerP b , resembling different technologies with different C-rate limits which are effectively power constraints. Then, from Fig. 18(a) we observe that by increasing λ we can effectively reduce the C-rate of the BESS not only by reducing its minimum SoC deviation but also by slowing down the BESS response, meaning that to reach its minimum SoC deviation value, more time is required. It is remarkable in that case, since the goalsẃeights were kept to the same values, that the controller requested more power from the BESS for lower λ values, resulting in almost identical system frequency response ( Fig. 18(b)). However, this is only available when the power rating of the BESS is large enough (as in our case) not to saturate the requested power.

Discussion and conclusions
This study proposed, analyzed and compared a set of multiobjective MIMO supervisory controllers for isolated power systems composed of a gas turbine, a wind turbine and a battery energy storage system. The proposed controllers are all based on the concept of integrating optimal frequency regulation for the isolated grid (isochronous operation) while at the same time tracking scheduling commands for the energy storage system. This was achieved though different model predictive controllers accounting for different types of uncertainty and different controller parametrizations. Overall, the goal has been of comparing deterministic MPC approaches against stochastic MPC ones and robust control. The results showed that stochastic control is not always associated with smaller constraint violation probability compared to the deterministic case. The occurrence of this event was, though, strongly correlated on the controller parametrization. In particular, we found that an affine disturbance parametrization is almost always associated with lower violations probabilities and very close to the open loop optimal policy (i.e., no parametrization of the controller). In addition, state feedback was found to be worse not only in terms of nominal performance but also in terms of constraint violation even when compared to the deterministic MPC. However a state feedback approach resulted in the smallest fuel consumption -at the cost though of slower regulation. Last but not least, the analyzed H ∞ robust controller was found resulting in higher constraint violations and sub-optimal coordination of the different subsystems, revealing the need for specialized tuning.
As for possible future work, higher parametric uncertainty could be considered for the comparison of the various proposed controllers and check to what extend the comparison results match the one presented in this study. In addition, higher fidelity simulation models could be used to further validate the results of the comparative analysis, while different specific battery technologies could also be examined for application suitability, given actual C-rate restrictions. Finally, a controller with risk aversion capabilities (i.e. by reducing the number of scenarios in the scenario approach) could be designed and be included in the comparison study.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. This paper presents a design methodology (proposed multiobjective controller) and a preliminary comparative analysis of various designs either with respect to the way uncertainty is considered (deterministic vs. robust vs. stochastic) or for different control parametrizations for the scenario-based MPC (SMPC). The scenario approach [40,41,43,48] gives the designer the unique capability to leverage risk of violating unseen constraints (from the randomly sampled scenarios). That means that we can reduce the number of scenarios considered by accepting higher risk for the sake of decreasing the controller's computational complexity and therefore the required calculation time. Under this perspective, we performed an analysis investigating the required time for the optimization problem to be solved, for different number of scenarios. The results are given in Table 2. Such times are recorded when using our available computational resources, a 2010 Mac-Book Pro with 2.4 GHz Intel Core 2 Duo processor, where just a single core was used (not parallel computing). We also not that the recorded execution times include several internal processes of the Matlab-based modeling system for convex optimization (CVX) and interactions with the selected simulation environment Matlab-Simulink. CVX is loaded with numerous additional tasks  other than solving the problem itself, such as interpreting the high-level modeling language into numerical matrices, loading searching and accessing those and many other staff that are not open to the user. In practice, increasing the number of scenarios we include more constraints to the optimization problem, making its solution harder. As we see from the table above, even for relatively low number of scenarios, the required solution time can be significant, giving rise to potential real-time implementation issues. The computational time to solve the optimization problem highly depends on the availability of computational resources which, from a system operator point of view, will not probably be a deadlock, since advanced computing may already be there as it is necessary for other tasks, too. While in the simulation environment there is always sufficient time for the optimization problem to be solved before the next simulation iteration, in a real time implementation non-sufficient time for the problem solution would result in a significant delay which is ignored in the simulations. Therefore, the results presented in this study can be thought of as an upper bound of the controller's performance and further studies are required to investigate what the required computational resources are and what parallelization techniques should be considered before the real-time implementation of the controllers proposed in this theoretical/numerical study.