*Article* **Resilient Predictive Control Coupled with a Worst-Case Scenario Approach for a Distributed-Generation-Rich Power Distribution Grid**

**Nouha Dkhili 1,2, Julien Eynard 1,2, Stéphane Thil 1,2 and Stéphane Grieu 1,2,\***


**Abstract:** In a context of accelerating deployment of distributed generation in power distribution grid, this work proposes an answer to an important and urgent need for better management tools in order to 'intelligently' operate these grids and maintain quality of service. To this aim, a modelbased predictive control (MPC) strategy is proposed, allowing efficient re-routing of power flows using flexible assets, while respecting operational constraints as well as the voltage constraints prescribed by ENEDIS, the French distribution grid operator. The flexible assets used in the case study—a low-voltage power distribution grid in southern France—are a biogas plant and a water tower. Non-parametric machine-learning-based models, i.e., Gaussian process regression (GPR) models, are developed for intraday forecasting of global horizontal irradiance (GHI), grid load, and water demand, to better anticipate emerging constraints. The forecasts' quality decreases as the forecast horizon grows longer, but quickly stabilizes around a constant error value. Then, the impact of forecasting errors on the performance of the control strategy is evaluated, revealing a resilient behaviour where little degradation is observed in terms of performance and computation cost. To enhance the strategy's resilience and minimise voltage overflow, a worst-case scenario approach is proposed for the next time step and its contribution is examined. This is the main contribution of the paper. The purpose of the min–max problem added upstream of the main optimisation problem is to both anticipate and minimise the voltage overshooting resulting from forecasting errors. In this min–max problem, the feasible space defined by the confidence intervals of the forecasts is searched, in order to determine the worst-case scenario in terms of constraint violation, over the next time step. Then, such information is incorporated into the decision-making process of the main optimisation problem. Results show that these incidents are indeed reduced thanks to the min–max problem, both in terms of frequency of their occurrence and the total surface area of overshooting.

**Keywords:** smart grid paradigm; distributed generation; model-based predictive control; robustness; worst-case scenario; min–max optimisation; intraday forecasting; Gaussian process regression; machine learning

### **1. Introduction**

Worldwide, the transition to renewable-energy-based power generation is in full swing. Because power grids were originally designed for centralised power generation with unidirectional power flow, large-scale deployment of renewable energy technologies, hereinafter referred to as distributed generation, ushers in numerous operational issues. The concept of "smart grid" was born out of the need to better monitor the behaviour of these evolving power grids, to more accurately anticipate the operational issues that could be caused by new components, and to more efficiently control them to ensure safety and service quality.

**Citation:** Dkhili, N.; Eynard, J.; Thil, S.; Grieu, S. Resilient Predictive Control Coupled with a Worst-Case Scenario Approach for a Distributed-Generation-Rich Power Distribution Grid. *Clean Technol.* **2021**, *3*, 629–655. https://doi.org/10.3390/ cleantechnol3030038

Academic Editor: Lieven Vandevelde

Received: 28 May 2021 Accepted: 27 August 2021 Published: 30 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Distributed power generators as well as storage devices are key players in modern power distribution grids. As a consequence, it goes without saying that these generators and devices can be used to balance power supply and demand. Based on this, the smart grid paradigm was mainly elaborated to tackle monitoring and control problems in power distribution grids [1]. Smart grids include advanced metering infrastructure and smart management schemes, which aim to boost grid observability and balance out power supply and demand, respectively. Regarding smart management schemes, optimal power flow (OPF) [2–7], demand-side management (DSM) [8–13], and multi-agent systems (MAS) [14–18] are some of the most numerous techniques used. One technique may be more appropriate than another, according to the target to be achieved, which can range from dimensioning and planning to real-time monitoring and control of power grids, and the technical and computational constraints to be taken into account. A survey of techniques for smart management of power distribution grids with prolific distributed generation is provided in [19]. The interested reader is referred to this paper.

This work falls within the scope of the "Smart Occitania" project (2017–2021), funded by the French agency for ecological transition (ADEME) and defined as a proof of concept for rural/suburban power distribution grids. An MPC-based strategy is proposed by PROMES-CNRS for smart management of a low-voltage power distribution grid with high levels of photovoltaics penetration using flexible assets, namely a biogas plant and a water tower. A simulated case study is carried out on a residential neighbourhood located in the Occitania region (southern France). The proposed strategy, which consists in an MPC scheme aiming to close the gap between power supply and demand, as stipulated by voltage constraints and the flexible assets' operational constraints [20], combines flexible asset management [21–23] and implicit model-based predictive control [24–27]. This strategy is fully explained in [20]; its principles are briefly outlined in Section 4.2. The suitability of model-based predictive control to the management of power distribution grids subject to disturbances—intermittent renewable-energy-based power generation and stochastic power demand—is plain. Because of their weakly meshed (often radial) structure, low-voltage power distribution grids are especially prone to cascading failures provoked by disturbances.

The main contribution of PROMES-CNRS is the a new problem formulation allowing the mixed-integer nonlinear programming (MINLP) [28–34] setting due to the ON/OFF controller the water tower is equipped with—this setting greatly increases the computational complexity of the optimisation problem—to be solved as a smooth continuous one without resorting to relaxation [20]. An in-depth analysis of the results has highlighted the MPC strategy's potential for upper-level power flow management and curtailment of voltage fluctuations in the considered low-voltage power distribution grid. Clearly, the strategy is a step towards low-voltage power distribution grids capable of integrating renewable-energy-based power generation while maintaining stability and quality of service.

The proposed MPC scheme incorporates intraday forecasts of grid load, water demand, and global horizontal irradiance (from which PV power generation is inferred; global horizontal irradiance is the total amount of shortwave radiation received from above by a surface horizontal to the ground), as well as their associated confidence intervals, in order to efficiently operate the aforementioned flexible assets towards balancing supply and demand within the grid. It does so using Gaussian process regression models. The stochasticity of the controller's inputs poses a threat to its efficacy, since it can incur unforeseen constraint violation (here, in the form of voltage undershooting/overshooting). The issue of stochastic data in MPC schemes is present in a wide array of problems, which makes it an active research field. There exists a significant body of literature handling stochastic MPC for linear problems [35–37].

As for nonlinear problems, which is the case addressed in this paper, more recent works have tackled them in different ways. A common rationale in the scientific literature is "multi-stage" control schemes. Recent examples include offline computation of an

incremental Lyapunov function, which is then used for an online construction of a "tube" to tighten the constraints [38], a decomposition into several deterministic sub-problems whose solutions are then aggregated using an operation-cost-based rule [39], and the modelling of the uncertainty through a tree of discrete scenarios, coupled with a tube-based MPC to balance the system's variability and its economic profitability [40,41].

The use of a multi-stage approach adds a layer of complexity into the control scheme. The advantages of such methods are their ability to combine several different techniques into a hierarchical scheme to tackle numerous difficulties in the problem. This often presents itself as an offline stage that feeds into an online one. Their obvious drawbacks is the added complexity and, in the case of scenario-based methods, significant computational burden, which make them ill-suited for real-time applications.

The method proposed in this paper is based on min–max MPC for uncertain nonlinear systems under constraints [42,43]. This technique's premise boils down to risk aversion by computation of a worst-case scenario upstream of a standard optimisation problem. The main optimisation problem upon which the predictive control strategy is based seeks optimal flexible assets' setpoints, in order to reduce the gap between power supply and demand in the power distribution grid. The solution proposed here adds a layer, to be called "the min-max problem", upstream of this main optimisation problem. The min– max problem will determine the values of possible PV power generation, grid load, and water demand for which constraint violation will be at a maximum, within the forecasts' respective confidence intervals. The scenario with these values is dubbed the worst-case scenario. Then, the predictive controller searches for optimal flexible assets' setpoints that would uphold constraints computed with values of PV power generation, grid load, and water demand corresponding to the worst-case scenario. The controller searches for worst-case scenario within a feasible space defined by the aforementioned confidence intervals associated to the GPR forecasts.

The paper is organised as follows: Section 2 presents the case study treated in this work. Section 3 provides a definition of Gaussian process regression, kernel compositions of the models developed to forecast all three stochastic quantities, and the obtained results. Then, Section 4 details the proposed model-based predictive control strategy, starting with a step-by-step explanation of the strategy's inner-workings, formulating the main optimisation problem upon which the strategy is based, and introducing the worst-case scenario approach which enhances the robustness of the control strategy to forecasting errors. In Section 5, the impact of forecasting errors on the performance of the model-based predictive control strategy is analysed and an evaluation is carried out of the contribution of the worst-case approach to enhancing the scheme's robustness to forecasting errors. Section 6 recapitulates the main findings of the paper and discusses possible improvements.

#### **2. Case Study**

#### *2.1. Low-Voltage Power Distribution Grid*

The case study [20] is carried out on a simulated low-voltage power distribution grid (a residential neighbourhood) in the Occitania region in southern France, composed of approximately 600 households, 200 household PV installations of 4 kW each, amounting to a total capacity of 800 kW, a small-scale biogas plant (100 kW), and a water tower (100 kW). Figure 1 displays some grid load and PV power generation data (over a week in April) used in this work. Grid load is measured at the MV/LV transformer level of the low-voltage power distribution grid whereas PV power generation is inferred from GHI measurements taken by a pyranometer installed at PROMES-CNRS laboratory, which is located a few kilometres from the residential neighbourhood.

In this section, models of the power distribution grid and the flexible assets used in this case study are presented, formulated over a forecast horizon *H*. Let *Hp* be the integer number of time slots within this forecast horizon. In the following, and for all time-dependant quantities, *t* ∈ {1, . . . , *Hp*}.

**Figure 1.** Grid load and PV power generation over a week in April (case study).

Measurements made at the MV/LV transformer level of the power distribution grid and used throughout this study correspond to means over each time slot (herein, a 10-min time step is considered). This study focuses on upholding contractual voltage bounds; voltage means must at all times remain within ± *δU* of the nominal value, i.e., ∀*q* ∈ {1, . . . , *N*} and ∀*t* ∈ {1, . . . , *Hp*}:

$$|\mathcal{U}\_{\emptyset}(t) - \mathcal{U}\_{n}| \lessapprox \delta \mathcal{U} \tag{1}$$

where *Un* is the nominal single-phase voltage value for all grid nodes, *N* is the number of nodes in the low-voltage power distribution grid, *Hp* is the integer number of time slots within the forecast horizon, and *δU* is the acceptable margin of voltage variations with respect to the nominal value. The acceptable margin fixed by the French distribution grid operator ENEDIS is 10% of the nominal value. However, for the purposes of this study, a lower margin of 3% is considered in order to better flesh out possible voltage overflow phenomena due to power supply/demand unbalance within the grid.

#### *2.2. Flexible Assets*

In this section, the biogas plant and water tower models are briefly described. Details about these models and the characteristics of the flexible assets are given in [20].

#### 2.2.1. Biogas Plant

The biogas volume in the storage unit (in m3) is described as:

$$\mathbf{V\_{b}}(t+1) = \mathbf{V\_{b}}(t) + \frac{T}{60} \left( \mathbf{Q\_{b,in}}(t) - \frac{\mathbf{P\_{b}}(t)}{\eta \cdot \text{LHV}} \right) \tag{2}$$

where *<sup>T</sup>* is the time step (*<sup>T</sup>* = 10 min), *Qb***,***in* is the flow rate of biogas production entering the storage unit (in m<sup>3</sup> h<sup>−</sup>1), *Pb* is the plant's active power output (in W), *η* is the generator's efficiency, and LHV is the lower heating value of the stored biogas (in kWh m<sup>−</sup>3). The output *Pb* is subject to the following constraint:

$$P\_{b,min} \lessapprox \mathbf{P}\_{\mathbf{b}}(t) \lessapprox P\_{b,max} \tag{3}$$

where *Pb*,*min* and *Pb*,*max* are the minimum and maximum power generation of the biogas plant, respectively.

The biogas volume in the storage unit is subject to the following constraint:

$$\mathbf{V}\_{b,min} \lessapprox \mathbf{V}\_{b}(t) \lessapprox \mathbf{V}\_{b,max} \tag{4}$$

where *Vb*,*min* and *Vb*,*max* are the minimum and maximum biogas storage capacities of the biogas plant, respectively.

#### 2.2.2. Water Tower

The water volume in the storage tank (in m3) is described as follows:

$$V\_{\rm av}(t+1) = V\_{\rm av}(t) + \frac{T}{60} \left( \frac{3600 \eta\_{\rm av}}{\rho \cdot g \cdot h} \mathbf{P}\_{\rm av}(t) - \mathbf{Q}\_{\rm av,out}(t) \right) \tag{5}$$

where *T* is the time step (*T* = 10 min), *Qw***,***out* is the flow rate of water demand (in m<sup>3</sup> h<sup>−</sup>1), *Pw* is the water pump's active power consumption (in W), *η<sup>w</sup>* is the water pump's efficiency, *ρ* is the water density (in kg m<sup>−</sup>3), *g* is the gravitational acceleration (in m s<sup>−</sup>2), and *h* is the water level in the storage tank (in m).

Let *Pw*,*min* and *Pw*,*max* be the minimum and maximum power consumption values of the water tower, respectively. The power consumption *Pw* can only be set following ON/OFF commands, i.e., it is subject to the following constraint:

$$P\_{\mathfrak{w}}(t) \in \{P\_{w,\min}; P\_{w,\max}\} \tag{6}$$

The water volume in the storage tank is subject to the following constraint:

$$\mathcal{V}\_{\text{uv},min} \lessapprox \mathcal{V}\_{\text{uv}}(t) \lessapprox \mathcal{V}\_{\text{uv},max} \tag{7}$$

where *Vw*,*min* and *Vw*,*max* are the minimum and maximum storage capacities of the water tank, respectively.

#### *2.3. PV Power Generation Inferred from Global Horizontal Irradiance*

PV power generation values are inferred from global horizontal irradiance values using the following Equation [44]:

$$
\hat{\mathbf{P}}\_{\rm PV}(t) = \eta\_{T\_{ref}} \cdot \mathbf{S} \cdot \bar{\mathbf{G}} \hat{\mathbf{H}} \mathbf{I}(t) \cdot \pi\_a \left[ 1 - \beta\_{ref} \left( T\_p - T\_{ref} \right) \right] \tag{8}
$$

where *ηTref* (herein, *ηTref* = 0.21) is the PV panel's efficiency, *S* is the total surface area of PV panels in the power distribution grid, **GHI** % are global horizontal irradiance forecasts, *τα* is the effective transmittance of the PV panels (herein, *τα* = 0.95), *βref* is the coefficient of power degradation due to high temperatures (*βref* = 0.004), *Tref* is the reference temperature (herein, *Tref* = 25 °C), and *Tp* is the PV panels' temperature, computed as follows [45]:

$$T\_p(t) = T\_a(t) + k \cdot \mathbf{GH} \mathbf{I}(t) \tag{9}$$

where *Ta* is the ambient temperature and *k* = 0.025.

#### *2.4. Initial Operation Strategy*

The initial operation of the flexible assets does not take into account power grid regulation purposes [20]. The biogas plant's power generation is a constant nominal value (herein, *Pb*,*<sup>n</sup>* = 100 kW). This is in line with the steady biogas flow generated by the bioreactor and coming into the storage unit. As for the water tower, its water pump is subject to an ON/OFF controller, which ensures that the water level remains between two threshold values at all times (herein, *Pw*,*<sup>n</sup>* = 100 kW). Both assets' priority is maintaining storage volume levels within pre-fixed extrema, which can sometimes be in conflict with the grid stability's best interest.

#### **3. Intraday Forecasting of Stochastic Quantities**

#### *3.1. Introduction*

In the context of the control strategy proposed in this paper, the three stochastic quantities that come into play are the following: power grid load, global horizontal irradiance (GHI) [46–48], and water demand. The power grid load represents agglomerated power demand of households in the studied suburban area. It is the grid operator's priority to

make sure this demand is met at all times, under adequate quality and security standards. PV power generation, inferred from GHI values as explained in Section 2.3, represents the agglomerated power generation of household PV panels in the studied area (the residential neighbourhood). One of the assets used by the control strategy is the water tower, whose operational priority is meeting the water demand. In this paper, the methodology used to forecast each of these quantities and its results are briefly presented.

As part of the research activities of the PROMES-CNRS laboratory, a comparative study of models developed for multi-step-ahead GHI forecasting is carried out for intrahour and intraday forecast horizons using a two-year database of GHI measurements sampled at a 10-min rate. Results demonstrate that, although all considered models outperform the persistence model, there is no clear frontrunner in terms of nRMSE values, with a slight advantage for LSTM (long short term memory) artificial neural network (ANN) and Gaussian process regression (GPR) models when taking into account the quality of the forecasts' associated confidence intervals.

Confidence intervals are a noteworthy perk of using GPR models in the forecast module of the smart management strategy since they are built-in in the models and do not require running a Monte Carlo simulation to be statistically inferred as is the case for other forecasting methods (like artificial neural networks). These confidence intervals give the predictive controller supplementary information it uses to achieve a more robust control strategy. The development of GPR models for intraday GHI forecasting are detailed in [46–48].

GPR models used herein are, in part due to the fact that the database used to update the model's parameters is a sliding one, sufficiently fast to be a valid candidate for the real-time control application at hand. On the downside, the predictive nature of the controller makes it dependant upon real-time access to measurements to be able to sustain the forecasting models' sliding databases and to update their parameters at each time step. This limitation is true for all machine-learning-based methods operating in real time, hence the growing interest in the development of state-of-the-art smart metering technologies that ensure quick and reliable data transfer.

For the reasons listed above, the choice has been made to develop GPR models to provide intraday forecasts for the smart management scheme developed herein. The associated confidence intervals are incorporated into the control scheme to improve its robustness to forecasting errors as explained in Section 4. This is the main contribution of the paper as it focuses on the developement of a resilient predictive control strategy for low-voltage power distribution grids with prolific distributed generation.

#### *3.2. Gaussian Process Regression*

A Gaussian process (GP) can be seen as a collection of random variables, any finite number of which have a joint Gaussian distribution [48,49]. A prior over functions is defined and can then be converted into a posterior over functions once some data has been observed (i.e., observations). *f*(*x*) ∼ GP(*μ*(*x*), *k*(*x*, *x* )), with *x* and *x* arbitrary input variables, *μ*(*x*) = E[ *f*(*x*)] the mean function (which is usually assumed to be null) and *k*(*x*, *x* ) = E (*f*(*x*) − *μ*(*x*))(*f*(*x* ) − *μ*(*x* )T) the covariance function (also known as *kernel*), indicates that this random function follows a GP. The interested reader is referred to [49] for a detailed explanation of Gaussian processes and Gaussian process regression (GPR).

The standard regression model with additive noise is formulated as follows:

$$y = f(\mathbf{x}) + \varepsilon \tag{10}$$

where *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>*D*×<sup>1</sup> is the input vector with a dimension of 1, *<sup>f</sup>* is the regression function, *<sup>y</sup>* is the observed value and *<sup>ε</sup>* ∼ N (0, *<sup>σ</sup>*<sup>2</sup> *<sup>ε</sup>* ) is an independent, identically distributed Gaussian noise.

The mean prediction *μ*<sup>∗</sup> can be written as a linear combination of kernel functions, each one centred on a data point:

$$
\mu\_\* = \mu(\mathbf{x}\_\*) + k\_\*^{\mathsf{T}} \cdot \boldsymbol{a} = \mu(\mathbf{x}\_\*) + \sum\_{i=1}^n a\_i \cdot k(\mathbf{x}\_i, \mathbf{x}\_\*) \tag{11}
$$

where *α* = *K* + *σ*<sup>2</sup> *<sup>ε</sup>* · *I* −<sup>1</sup> (*y* − *μ*(*x*∗)). The coefficients *α<sup>i</sup>* are referred to as parameters.

It is possible to combine several kernel functions to obtain a more complex one. The only requirement is that the resulting covariance matrix must be a positive semi-definite function. A detailed list of kernel functions is provided in [49]. Hereinafter, the kernel functions used in this work are defined.

The periodic kernel (Per) is given by:

$$k\_{\rm Per}(\mathbf{x}, \mathbf{x}') = \sigma\_1^2 \cdot \exp\left(-\frac{2\sin^2\left(\frac{\pi(\mathbf{x} - \mathbf{x}')}{P}\right)}{\ell\_1^2}\right) \tag{12}$$

where *σ*<sup>1</sup> is the amplitude, -<sup>1</sup> is the correlation length and *P* is the period.

The isotropic squared exponential (SE) kernel is given by:

$$k\_{\rm SE}(\mathbf{x}, \mathbf{x}') = \sigma\_2^2 \cdot \exp\left(-\frac{(\mathbf{x} - \mathbf{x}')^2}{2\ell\_2^2}\right) \tag{13}$$

where *σ*<sup>2</sup> is the amplitude and -<sup>2</sup> is the correlation length.

The isotropic rational quadratic (RQ) kernel is given by:

$$k\_{\rm RQ}(\mathbf{x}, \mathbf{x}') = \sigma\_3^2 \left( 1 + \frac{(\mathbf{x} - \mathbf{x}')^2}{2\ell\_3^2 \cdot \mathbf{a}} \right)^{-\alpha} \tag{14}$$

where *σ*<sup>3</sup> is the amplitude, -<sup>3</sup> is the correlation length and *α* defines the relative weighting of large-scale and small-scale variations.

In addition to forecasts, the regression model provides an associated confidence interval within which measurements have a probability of 95% of staying. This interval (CI) is computed as follows:

$$\text{CI} = \mu\_\* \pm 1.96 \sqrt{\sigma\_\*^2} \tag{15}$$

where CI represent the confidence interval bounds, *<sup>μ</sup>*<sup>∗</sup> is the predictive mean, and *<sup>σ</sup>*<sup>∗</sup> is the predictive variance.

During the training phase, the values of hyperparameters are optimised using a 2 week database [48]. During the test phase, at every time step, new measurements are incorporated in order to update values of the model's parameters using a sliding 24 h database.

#### *3.3. Data Description and Kernel Compositions*

The forecast horizons considered in this work are intraday: they range from 1 h to 24 h. For the training phase, two weeks of data are used. The testing phase is then performed over one week. Available data of power grid load, water demand, and global horizontal irradiance span a year for the three stochastic quantities. They are sampled at a 10-min rate. An examination of the signals' behaviours informs the following choices of kernel compositions for their intraday forecasting.

1. Power grid load

$$k^{\mathbb{g}l}(\mathbf{x}, \mathbf{x}') = k\_{\text{Per}}^{\mathbb{g}l}(\mathbf{x}, \mathbf{x}') + k\_{\text{SE}}^{\mathbb{g}l}(\mathbf{x}, \mathbf{x}') + k\_{\text{RQ}}^{\mathbb{g}l}(\mathbf{x}, \mathbf{x}') \tag{16}$$

Herein, the period of the data is 24 h and the structure of *k gl* Per translates the daily periodic pattern. *k gl* SE is used to capture the long-term trend of the data, namely the seasonal tendencies present in the power grid load. *k gl* RQ is used to capture the intraday variations observed in the data, namely intraday fluctuations due to behavioural patterns of end-users.

2. Water demand

$$k^{\rm vvt}(\mathbf{x}, \mathbf{x}') = k^{\rm vvt}\_{\rm Per}(\mathbf{x}, \mathbf{x}') + k^{\rm vvt}\_{\rm SE}(\mathbf{x}, \mathbf{x}') \tag{17}$$

*kwt* Per models the daily periodic shape of the data and *<sup>k</sup>wt* SE is used to fit the data's intraday fluctuations. A brief kernel study leads to the conclusion that the simple addition of a squared exponential kernel is sufficient to have satisfactory results without adding too much computational burden to the model.

3. Global horizontal irradiance

$$k^{\rm GHI}(\mathbf{x}, \mathbf{x}') = k\_{\rm Per}^{\rm GHI}(\mathbf{x}, \mathbf{x}') \cdot k\_{\rm RQ}^{\rm GHI}(\mathbf{x}, \mathbf{x}') \tag{18}$$

*k*GHI Per models the daily periodic shape of the data and *<sup>k</sup>*GHI RQ captures the intraday fluctuations present in these data. A thorough study has been conducted in [47] to determine the most suitable kernel composition for intraday GHI forecasting.

#### *3.4. Forecasting Results*

The evaluation metrics used in this work, evaluated over the considered week in April (see Figure 1), are the following:

1. The normalized root mean square error (nRMSE), expressed as follows:

$$\text{nRMSE} = 100 \frac{\sqrt{\frac{1}{n\_\*} \sum\_{i=1}^{n\_\*} \left( y\_{\text{best}}(i) - y\_{\text{forecast}}(i) \right)^2}}{\frac{1}{n\_\*} \sum\_{i=1}^{n\_\*} y\_{\text{best}}(i)} \tag{19}$$

where *<sup>y</sup>*test <sup>∈</sup> <sup>R</sup>*n*∗×<sup>1</sup> are the test data, *<sup>y</sup>*forecast <sup>∈</sup> <sup>R</sup>*n*∗×<sup>1</sup> are the forecasts given by the models, and *n*<sup>∗</sup> is the number of data points in the forecast horizon.

2. The coverage width-based criterion (CWC) [50], defined as a combination of two different criteria, i.e., the prediction interval normalized average width (PINAW) and the prediction interval coverage probability (PICP).

The PINAW criterion allows the surface area of the confidence intervals associated with the forecasts to be assessed:

$$\text{PINAW} = \frac{1}{n \cdot R} \sum\_{i=1}^{n} (\mathcal{U}\_{i} - L\_{i}) \tag{20}$$

where *R* is the difference between the maximum and minimum in test data and *Ui* and *Li* are the upper and lower bounds of the confidence interval, respectively. The PICP criterion informs on the probability that measurements would fall within

the confidence interval:

$$\text{PICP} = \frac{1}{n} \sum\_{i=1}^{n} \epsilon\_i \tag{21}$$

where *<sup>i</sup>* is used to detect whether the target value is in the confidence interval. The coverage width-based criterion (CWC) is then defined as follows:

$$\text{CWC} = \text{PINAW} \{ 1 + \gamma(\text{PICP}) \cdot \exp \left( -\eta \cdot (\text{PICP} - \mu) \right) \} \tag{22}$$

where *η* (*η* = 10) and *μ* (*μ* = 0.95) are parameters, and *γ* is defined such that:

$$\gamma(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} < \mu \\ 0 & \text{if } \mathbf{x} \gg \mu \end{cases} \tag{23}$$

Figure 2 displays nRMSE values of intraday forecasts of power grid load, water demand, and GHI. For intraday forecasts of grid load, results show that the GPR model performs well for short time horizons. However, the forecasting error rapidly increases as the forecast horizon does, which is an expected result. For a 1-h horizon, the GPR model achieves an error as low as 13%. As for the 6-h horizon and the 12-h one, nRMSE values are virtually constant (21% and 23%, respectively). This remains the case for a forecast horizon of 24 h, where the GPR model has an error of 26%.

**Figure 2.** Evaluation criteria (nRMSE (19) and CWC (22)) for intraday GPR forecasts of GHI, power grid load and water demand, with respect to the forecast horizon, over a one-week period.

For intraday forecasts of water demand, values of nRMSE start at 9% for a 10 min forecast and quickly stabilises around 13% for longer forecast horizons. Due to the data's regular nature, a simple combination of a periodic kernel and a rational quadratic kernel for intraday fluctuations proves enough to capture the data's behaviour quite faithfully. As for intraday forecasts of GHI (from which PV power generation is inferred), the selected GPR model performs well for short forecast horizons (15.74% for a 10-min forecast horizon), but its performance degrades as the forecast horizon grows and it stabilises around 32% for horizons beyond 6 h. This is reflected in the temporal evolution of the forecasts, where it can be seen that intraday fluctuations become increasingly difficult to fit. Further work is currently underway to enhance the model's performance for long intraday forecast horizons.

Figure 2 also displays values of CWC, which follow similar patterns to those of nRMSE values: starting with lower values for short forecast horizons, then quickly stabilizing around a constant value for longer horizons. It can be seen that the higher the nRMSE values, the higher the corresponding CWC values. They stabilize around 57% for water demand forecasts, around 65% for grid load forecasts, and around 68% for GHI forecasts.

#### **4. Model-Based Predictive Control Strategy**

*4.1. Inner-Workings of the Smart Management Scheme*

Figure 3 shows the synoptic scheme of the proposed MPC strategy. The following is a detailed explanation of the different steps taken by the proposed strategy at each time step.

**Figure 3.** Synoptic scheme of the amended MPC-based strategy for smart management of a low-voltage power distribution grid using flexible assets. Let **GHI**, *Pcons*, *Qw***,***out*, *Vw*, *Vb* be measurements of global horizontal irradiance, grid load, water demand, water volume, and biogas volume, respectively. Let *<sup>Y</sup>*\$ be forecasts of stochastic inputs for the following time steps. Let *δ*PV, *δc*, *δ<sup>w</sup>* be margins that define confidence intervals for the next time step of GPR forecasts of PV power generation, grid load, and water demand, respectively. Let *<sup>P</sup>*\$**PV**, *<sup>P</sup>*\$*cons*, *<sup>Q</sup>*\$*w***,***out* be forecasts of PV power generation, grid load, and water demand over the forecast horizon, respectively. Let *Yrisk <sup>s</sup>* and *Yrisk* be candidate values and optimal values of worst-case scenario stochastic inputs. Let *<sup>P</sup><sup>s</sup> <sup>b</sup>* and *<sup>P</sup><sup>s</sup> w* be candidate values of biogas plant setpoints and water tower setpoints, respectively. Let *<sup>V</sup><sup>s</sup> <sup>b</sup>* and *<sup>V</sup><sup>s</sup> w* be the biogas volume and the water volume, respectively, corresponding to the candidate optimisation variables of a given iteration.

1. Data acquisition: measurements are taken of stored biogas volume and stored water volume and are injected into the predictive controller in order to update the system's state. GHI and grid load are measured and water demand is inferred from measurements of water volume and incoming flow rate *Qw***,***in*. These values are injected into the forecast module.


An optimisation algorithm searches the feasible space defined in Figure 4 for worstcase scenario input values. It does so based on the optimisation problem formulated by Equations (46)–(49) and using the low-voltage grid model to evaluate volume bounds and potential voltage overshooting.

4. Reduction of power supply/demand unbalance: phase 2 of the controller's decisionmaking process consists in determining the flexible assets' optimal setpoints. This bloc receives GPR forecasts of PV power generation, grid load, and water demand over the entire forecast horizon, as well as worst-case scenario input values of the following time slot, produced by phase 1.

An optimisation algorithm searches for optimal setpoints of biogas plant power generation and water tower power consumption. It does so based on the optimisation problem defined in Section 4.2 and using the low-voltage grid model based on Kirchhoff laws to evaluate various constraints.

5. Implementation of flexible assets' setpoints: optimal setpoints of flexible assets are produced, the first step of which are implemented by the biogas plant and the water tower.

#### *4.2. Main Optimisation Problem*

The optimisation problem [20] is formulated such as the discrete setpoint values of the water tower are replaced with a real-values optimisation variable *<sup>t</sup>*¯ <sup>∈</sup> <sup>R</sup>*Hp* , with *Hp* the integer number of time slots within the forecast horizon, effectively avoiding the MINLP setting [28–34] without relaxing the problem. In fact, between sampling times *ti* and *ti*+1, the water tower can switch between its two states of operation (*Pw*,*min* and *Pw*,*max*) and the biogas plant can switch between two setpoints within the interval [*Pb*,*min*, *Pb*,*max*]. If the assets do make the switch within a given time slot, they do so at the same instant ¯*ti*. Let us denote *X* as follows:

$$X = \begin{bmatrix} P\_{b, \text{ON}} \ P\_{b, \text{OFF}} \ \text{if} \ \mathcal{U}\_{\text{ON},q} \ \mathcal{U}\_{\text{OFF},q} \end{bmatrix}^{\mathsf{T}} \tag{24}$$

where *Pb***,ON** <sup>∈</sup> <sup>R</sup>*Hp* and *Pb***,OFF** <sup>∈</sup> <sup>R</sup>*Hp* form the biogas plant setpoint like this:

$$P\_{\mathsf{b}}(\tau) = \begin{cases} P\_{\mathsf{b},\mathsf{ON}}(\tau), \tau \in \left[t\_i, t\_i + \overline{t}\_i\right] \\ P\_{\mathsf{b},\mathsf{OFF}}(\tau), \tau \in \left[t\_i + \overline{t}\_i, t\_{i+1}\right] \end{cases} \tag{25}$$

For every node *<sup>q</sup>* ∈ {1, ... , *<sup>N</sup>*}, *<sup>U</sup>***ON**,*<sup>q</sup>* <sup>∈</sup> <sup>R</sup>(*Hp*·*N*) and *<sup>U</sup>***OFF**,*<sup>q</sup>* <sup>∈</sup> <sup>R</sup>(*Hp*·*N*) form the voltages in the grid:

$$\mathcal{U}\_{\boldsymbol{q}}(\boldsymbol{\tau}) = \begin{cases} \mathcal{U}\_{\text{ON},\boldsymbol{q}}(\boldsymbol{\tau}), \boldsymbol{\tau} \in [t\_{i\prime}t\_{i} + \overline{t}\_{i}] \\ \mathcal{U}\_{\text{OFF},\boldsymbol{q}}(\boldsymbol{\tau}), \boldsymbol{\tau} \in [t\_{i} + \overline{t}\_{i\prime}t\_{i+1}] \end{cases} \tag{26}$$

The optimisation problem can be solved at each time step assuming that the first state of the water tower is always ON. In some extreme cases, this assumption may induce some issues of implementability given volume constraints, which are tackled in a posttreatment phase (the interested reader is referred to [20] for details about the post-treatment algorithm). However, this simplification reduces the complexity of the model without sacrificing much performance. The objective function *fobj* is formulated as follows:

$$f\_{obj}(\mathbf{X}) = \sum\_{i=0}^{H\_p - 1} \left[ \int\_{t\_i}^{t\_i + \overline{t}\_i} \mathbf{S\_{ON}}(\tau) \mathbf{d}t + \int\_{t\_i + \overline{t}\_i}^{t\_{i+1}} \mathbf{S\_{OFF}}(\tau) \mathbf{d}t \right] \tag{27}$$

with

$$\mathbf{S\_{ON}}(\tau) = \left| \mathbf{P\_{PV}}(\tau) + \mathbf{P\_b}(\tau) - \mathbf{P\_{counts}}(\tau) - \mathbf{P\_{w,max}} \right|^2 \tag{28}$$

$$S\_{\rm OFF}(\tau) = |P\_{\rm PV}(\tau) + P\_{\rm b}(\tau) - P\_{\rm cons}(\tau) - P\_{w, min}|^2 \tag{29}$$

The optimisation problem is then formulated as follows:

$$X^\* = \arg\min\_X f\_{obj}(X) \tag{30}$$

and is subject to the following bounds and constraints, ∀*i* ∈ {1, . . . , *Hp*}:

• Biogas plant power bounds

$$P\_{b,min} \lessapprox P\_{b, \mathbf{ON}}(\tau) \lessapprox P\_{b,max} \tag{31}$$

$$P\_{b,min} \lessapprox \mathbf{P}\_{b, \text{OFF}}(\tau) \lessapprox P\_{b,max} \tag{32}$$

• Switch time bounds

$$0 \ll\_{\mathbb{R}} \bar{t}\_i \ll\_{\mathbb{R}} T \tag{33}$$

• Biogas volume constraints

$$\mathbf{V}\_{b,min} \lessapprox \mathbf{V}\_{\mathbf{b}}(\mathbf{t}) \lessapprox \mathbf{V}\_{b,max} \tag{34}$$

• Water volume constraints

$$V\_{\text{uv,min}} \lessapprox V\_{\text{uv}}(t) \lessapprox V\_{\text{uv,max}} \tag{35}$$

• Voltage constraints

$$\ddot{\boldsymbol{\pi}}\_{i} \cdot \mathcal{K}\_{t}(\mathbf{P}\_{\text{b},\text{ON}}(\boldsymbol{\pi}), \boldsymbol{P}\_{\text{w},\text{max}}, \mathbf{U}\_{\text{ON},\text{q}}(\boldsymbol{\pi})) = \mathbf{0} \tag{36}$$

$$(T - \bar{t}\_i) \cdot K\_l(\mathbf{P\_{b,OFF}}(\tau), P\_{w,min}, \mathbf{U\_{OFF,q}}(\tau)) = 0 \tag{37}$$

$$|\mathcal{U}\_{\mathbf{ON},q}(\tau) - \mathcal{U}\_n| \leqslant \delta \mathcal{U} \tag{38}$$

$$|\mathcal{U}\_{\text{OFF},\mathcal{J}}(\tau) - \mathcal{U}\_{\text{n}}| \lessapprox \delta \mathcal{U} \tag{39}$$

*Kt* describes voltage variations across the low-voltage power distribution grid as a function of the power which is injected or absorbed at each node. Kirchhoff's law constraints are presented as two sets of constraints (Equations (36) and (37)), which guarantee that Kirchhoff's laws are upheld in both sub-intervals of each time slot. The equation set depicting voltage variations is multiplied by ¯*ti* (Equation (36)) or *<sup>T</sup>* − ¯*ti* (37), using appropriate values of biogas plant and water tower setpoints for each interval, in order to ensure that only one constraint is activated in case of extreme values of ¯*ti*. In fact, in case ¯*ti* = 0, Equation (36) is ignored, reflecting the fact that during time slot *i*, the water tower is turned off immediately at the beginning of the time slot. Likewise, in case ¯*ti* = *T*, Equation (37) is ignored since the water tower remains on for the duration of the time

slot. Voltage constraints are also written as two sets of constraints (Equations (38) and (39)) that account for voltage variations in both states of the low-voltage power distribution grid within each time slot. While depicting the same physical constraints, this problem formulation has a bigger feasible set than the MINLP one [28–34]. As a result, with such a formulation, the global optimum is guaranteed to be equal or better than the one of the MINLP formulation.

#### *4.3. The Worst-Case Scenario Approach*

In this section, amendments are made to the predictive control strategy in order to enhance the controller's performance by making it more robust to forecasting errors of the system's various stochastic quantities. The premise of the method is to use min–max optimisation in order to find and solve a "worst-case scenario" at each time step based on confidence intervals given by the inputs' respective GPR models. Eliminating, or at least minimising, the constraint violations corresponding to the worst-case scenario guarantee that these violations are also minimised for all other possible scenarios.

It should be noted that, at each time step, these amendments are only made to the decision-making of the subsequent time step, and not over the entire forecast horizon. This choice is motivated by two reasons. The first is that determining a worst-case scenario over the entire forecast horizon is a conservative and computationally expensive optimisation problem, which is incompatible with the real-time application at hand. As a matter of fact, a min–max problem over the entire forecast horizon has a feasible space of dimension (3 × *Hp*), as opposed to the problem concerned only with the following time step only having a 3-dimensional feasible space. The second is that, due to the closed-loop nature of MPC, computing robust setpoints for the entire forecast horizon is in high likelihood a waste of resources because, at each time step, only the first setpoint is applied and the whole procedure is reiterated at the next one. Ergo, the scheme only seeks to provide a setpoint robust to forecasting errors for the next time step.

First, let *<sup>P</sup>*\$**PV**, *<sup>P</sup>*\$*cons*, and *<sup>Q</sup>*\$*w***,***out* be forecasted values of PV power generation, grid load, and water demand, respectively, over the forecast horizon. Then, let *P***PV**, *Pcons*, and *Qw***,***out* be measured values of PV power generation, grid load, and water demand, respectively, over the forecast horizon. Finally, *δ*PV, *δc*, and *δ<sup>w</sup>* define confidence intervals, for the next time step, of GPR forecasts of PV power generation, grid load, and water demand, respectively, as follows:

$$P\_{\rm PV}(t+1) \in \left[\hat{P}\_{\rm PV}(t+1) - \delta\_{\rm PV}, \hat{P}\_{\rm PV}(t+1) + \delta\_{\rm PV}\right] \tag{40}$$

$$P\_{\rm cons}(t+1) \in \left[\widehat{P}\_{\rm cons}(t+1) - \delta\_c, \widehat{P}\_{\rm cons}(t+1) + \delta\_c\right] \tag{41}$$

$$Q\_{w,out}(t+1) \in \left[\hat{Q}\_{w,out}(t+1) - \delta\_{\mathcal{W}}, \hat{Q}\_{w,out}(t+1) + \delta\_{\mathcal{W}}\right] \tag{42}$$

Herein, the confidence intervals are computed such that there is a 95% probability of measurements remaining inside them (15). There exists a triplet (*Prisk* **PV** (*<sup>t</sup>* <sup>+</sup> <sup>1</sup>), *<sup>P</sup>risk cons*(*t* + 1), *Qrisk <sup>w</sup>***,***out*(*t* + 1)), contained in the feasible space illustrated in Figure 4, that induces a worst-case scenario vis-a-vis the optimisation problem constraints in the next time step. Finding this triplet and incorporating it into the predictive control strategy described in Section 4 allows the controller to adjust its decision-making in order to reduce, if possible eliminate, any constraint violation that could occur in the next time step as a result of the stochastic quantities' measured values being different from the forecasted ones, within the limits of the associated confidence intervals. Let *Y* be the vector containing measurements of PV power generation, grid load, and water demand, for the next time step:

$$\mathbf{Y} = \begin{bmatrix} P\_{\text{PV}}(t+1) \ P\_{\text{cons}}(t+1) \ Q\_{w,out}(t+1) \end{bmatrix}^T \tag{43}$$

Let *<sup>Y</sup>*\$ be the vector containing forecasts of PV power generation, grid load, and water demand, for the next time step:

$$\hat{Y} = \begin{bmatrix} \hat{P}\_{\text{PV}}(t+1) \ \hat{P}\_{\text{cons}}(t+1) \ \hat{Q}\_{w,out}(t+1) \end{bmatrix}^{T} \tag{44}$$

Let *Yrisk* be the vector comprised of critical values defining the worst-case scenario of the next time step:

$$\mathcal{Y}^{risk} = [P\_{\rm PV}^{risk}(t+1) \ P\_{cons}^{risk}(t+1) \ \mathcal{Q}\_{w,out}^{risk}(t+1)]^T \tag{45}$$

where *Prisk* **PV** , *<sup>P</sup>risk cons*, and *Qrisk <sup>w</sup>***,***out* are critical values of PV power generation, grid load, and water demand, respectively, corresponding to the worst-case scenario.

Based on Equations (40)–(42), *Yrisk* exists in the feasible space illustrated by Figure 4. At every time step, the following min–max problem is solved:

$$Y^{risk} = \arg\min\_{\mathbf{Y}} (-\Phi(\mathbf{Y})) \tag{46}$$

where Φ(*Y*) is the voltage undershooting/overshooting corresponding to input values of *Y*, as defined in Section 5.1, subject to:

$$
\hat{Y}(1) - \delta\_{\rm PV} \leqslant \mathcal{Y}(1) \leqslant \hat{Y}(1) + \delta\_{\rm PV} \tag{47}
$$

$$
\hat{Y}(2) - \delta\_{\mathfrak{c}} \leqslant \Upsilon(2) \leqslant \hat{Y}(2) + \delta\_{\mathfrak{c}} \tag{48}
$$

$$
\hat{Y}(\mathfrak{Z}) - \delta\_w \lesssim \Upsilon(\mathfrak{Z}) \lesssim \hat{Y}(\mathfrak{Z}) + \delta\_w \tag{49}
$$

**Figure 4.** Feasible space of the min–max problem, defined by the confidence intervals of one-stepahead forecasts of grid load, water demand, and PV power generation. The time indices are removed to avoid cluttering the illustration. All quantities in this figure correspond to values in the next time slot.

#### **5. Results and Analysis**

#### *5.1. Evaluation Metrics*

Various evaluation metrics used throughout this paper are defined hereinafter. Let *Ht* be the number of time slots in the simulation period.


$$
\Omega\_{\rm PV} = \frac{\sum\_{t=0}^{H\_l - 1} \mathsf{P} \mathbf{p} \,(t+1) - \mathsf{P} \mathbf{p} \,(t+1)}{H\_l} \tag{50}
$$

<sup>Ω</sup>*Pcons*, which is the mean deviation of grid load values (*Pcons*) from the ones forecasted at a one-step-ahead forecast horizon (in W), is given by:

$$
\Omega\_{P\_{\rm conv}} = \frac{\sum\_{t=0}^{H\_l - 1} P\_{\rm cons}(t+1) - P\_{\rm cons}(t+1)}{H\_t} \tag{51}
$$

<sup>Ω</sup>*Qw*,*out* , which is the mean deviation of water demand values (*Qw***,***out*) from the ones forecasted at a one-step-ahead forecast horizon (in m3 h<sup>−</sup>1), is given by:

$$
\Omega\_{Q\_{w,out}} = \frac{\sum\_{t=0}^{H\_l - 1} \mathbb{Q}\_{w,out}(t+1) - \mathbb{Q}\_{w,out}(t+1)}{H\_l} \tag{52}
$$


$$\Phi(\mathbf{Y}) = \sum\_{s=1}^{N} \left( \sum\_{t=1}^{H\_l} \max\left( \mathbf{U}\_{\mathbf{s}}(t) - \mathbf{U}\_{\min}, \mathbf{U}\_{\max} - \mathbf{U}\_{\mathbf{s}}(t), 0 \right) \right) \tag{53}$$

where *N* is the number of nodes in the power distribution grid, *Umin* is the lower voltage threshold, *Umax* is the upper voltage threshold, and *Us*, where *<sup>s</sup>* ∈ {1, ... , *<sup>N</sup>*}, are voltage values in various nodes of the power distribution grid. Let **<sup>Ψ</sup>** <sup>∈</sup> <sup>R</sup>*Hp*×<sup>3</sup> be the matrix containing the measured stochastic inputs of the control strategy, over the simulation period, defined as follows:

$$\mathbf{Y} = \begin{bmatrix} P\_{\text{PV}} \ P\_{\text{cons}} \ Q\_{w,out} \end{bmatrix} \tag{54}$$

Note that voltage values across the power distribution grid are linked to the values of **Ψ** through Kirchhoff laws.

6. Average voltage overshooting per time step *φ*: it corresponds to the mean of maximum voltage overshooting over the number of instances at which overshooting is observed during the simulation period. It is defined as follows:

$$\phi(\mathbf{Y}) = 100 \frac{\Phi\_{\text{max}}(\mathbf{Y})}{\nu \cdot H\_{\text{f}}} \tag{55}$$

with

$$\Phi\_{\max}(\mathbf{Y}) = \sum\_{t=1}^{H\_t} \max\_{s \in \{1, \dots, N\}} (\mathbf{U}\_{\mathbf{s}}(t)) \tag{56}$$

where values of *Y* and *Us* are related through Kirchhoff laws.

#### *5.2. Impact of Forecasting Errors on MPC Performance*

In this section, the impact of forecasting errors on the MPC strategy's performance is studied. The control scheme is tested over the week in April presented in the previous section and for sliding window sizes ranging from 1 to 24 h. Intraday GPR forecasts of grid load, water demand, and PV power generation (inferred from GHI forecasts), acquired as explained in Section 3, are used to run these simulations.

The performance of the predictive controller fed with GPR forecasts is evaluated in comparison with a controller fed with measurements, i.e., a case where no forecasting errors are made. This comparison will focus on three main aspects of the scheme's performance: the power supply/demand gap *fobj*,, *final* (Figure 5), computational cost *κ* (Figure 6), and surface area of voltage overshooting Φ (Figure 7).

**Figure 5.** The cumulative power supply/demand gap within the power distribution grid, per sliding window size. The gap before optimisation is 10,035 kW.

In Figure 5, the cumulative power supply/demand gap given by the MPC scheme over the considered week is displayed per sliding window size. Though these values degrade as the sliding window size increases in both cases, the ones given by the MPC scheme when it uses GPR forecasts are not significantly degraded with respect to those given by the MPC scheme that uses measurements. In fact, the maximum difference between the two curves is 198.39 kW, obtained for the 22-h sliding window, which constitutes only 1.98% of the initial value (10,035 kW).

**Figure 6.** Computational complexity, measured as the mean number of function evaluations per sliding window weighted by its size.

The computational cost, presented in Figure 6, steadily grows in both cases as the sliding window size does. It is different for the scheme that uses GPR forecasts, which is to be expected since the updated forecasts displace the optimisation problem's starting point with respect to the previous time step. That being said, the increase in computational cost due to the use of forecasts remains subdued. For window sizes between 1 and 10 h, its average value is a 12.3% increase from the scheme using measurements. For all window sizes up to 24 h, this average is evaluated at 16.4%.

Figure 7 displays the surface area of voltage overshooting, the initial value of which is 4371.4 kV. For sliding window sizes up to 3 h, the MPC scheme is unable to reduce voltage overshooting and actually makes matters considerably worse. Starting from a 4-h sliding window, voltage overshooting given by the MPC scheme decreases significantly from the initial value, in both the case where the scheme uses measurements and the one where it uses GPR forecasts. Then, for larger window sizes, it quickly stabilizes around the same level. As of the 4 h window size, the MPC strategy effectively eliminates more than 50% of voltage overshooting.

**Figure 7.** The total surface area of voltage overshooting per sliding window size.

When taking into account the fact that the forecasting errors are at their lowest for short forecast horizons (inferior to 3 h) and rapidly grow for longer horizons, it becomes apparent that the accuracy of forecasts for these short horizons is not enough to guarantee better management of voltage fluctuations on its own. In reality, the availability of forecasts over a longer forecast horizon is pivotal to better equip the MPC scheme to anticipate emerging voltage overshooting and work to prevent it. In light of this observation, it is recommended, for the purposes of this study, to prioritise reducing forecasts' error rates for forecast horizons up to several hours, rather than only focusing on short horizons.

Figure 8 displays the percentages of time steps during the simulated week where overshooting is observed, with respect to the size of the sliding window. As of the 2-h window, this percentage is significantly lower than the initial one (23.71%) for both MPC schemes, the one using measurements and the one using GPR forecasts. It quickly stabilizes at roughly 7% for both cases and reaches a minimum of 3.67% for the former and of 5.16% for the latter, both corresponding to the 24 h window.

Figure 9 displays the average voltage overshooting per time step (*φ*), with respect to the sliding window size. These values are lower than the initial value (18.29%) for both schemes, with the values corresponding to the MPC scheme using GPR forecasts being slightly lower than the ones corresponding to the MPC scheme using measurements. It is interesting to note that despite the percentage of overshooting occurrences (*ν*) being significantly higher for small windows than for larger ones, values of *φ* remain at roughly the same level regardless of sliding window size. Their average is 11.6 V per time step for the scheme using measurements and 10.7 V per time step for the one using GPR forecasts. This means that, for large window sizes, overshooting incidents are less frequent, but have a higher amplitude than for small ones.

**Figure 8.** Percentage of time steps where an overshooting is observed, with respect to the sliding window size used by the MPC scheme.

**Figure 9.** Average voltage overshooting per time step, with respect to the sliding window size used by the MPC scheme.

On another note, voltage overshooting is not remarkably impacted by the use of GPR forecasts as opposed to the case where measurements are used. The two sets of values, whether in terms of cumulative voltage overshooting (Figure 7) or percentage-wise (Figure 8), behave similarly and remain roughly at the same level. These results point to the MPC strategy proposed herein being inherently resilient to forecasting errors of PV power generation, grid load, and water demand. This is thanks to the closed-loop structure of model-based predictive control, which allows course-correction as new information comes in at every time step.

Although the MPC strategy succeeds in reducing voltage overshooting with respect to the initial case, it does not completely eliminate it. In order to further enhance the strategy's performance, a complementary module can be added upstream of the main optimisation problem upon which the predictive controller is based. This module attempts to limit overshooting due to erroneous GPR forecasts of grid load, water demand, and PV power generation.

#### *5.3. Contribution of the Worst-Case Approach*

It should be noted here that the main optimisation problem responsible for balancing power supply and demand in the power distribution grid, as constructed in Section 4.2, prioritises constraints so that the volume ones are always upheld. In other words, in cases where no feasible solution is found, voltage constraints are relaxed in order to guarantee that biogas volume constraints and water volume constraints are always upheld. For this reason, hereinafter, only voltage constraint violation is examined to evaluate the efficiency of the proposed min–max problem in enhancing the control strategy's robustness to forecasting errors.

In this section, the amended predictive control strategy as explained above is implemented, and its results are presented and discussed in comparison with three other cases:


For each case, a simulation is run over a week in April. This period is selected because it presents high PV power generation and therefore demonstrates significant voltage overshooting. Two sizes of sliding windows are used for the MPC scheme in the results that are presented hereinafter: a 4 h window and a 10 h window. These two sliding window sizes are chosen to examine the difference in effects of the min–max problem on the MPC strategy's performance for both short sliding windows and long ones. The evaluation metrics of the MPC scheme with both sliding window sizes are assembled in Table 1.

**Table 1.** Assessment of the min–max problem's contribution to the control strategy's robustness to forecasting errors, for a week in April. See below (Section 5.3) for details about the three cases. For Case 2 and Case 3, 4 h and 10 h sliding windows are considered.


When examining the inner-workings of the min–max problem, it can be deduced that there are noteworthy deviations between forecasted values of PV power generation and grid load and the ones corresponding to worst-case scenarios. For the considered week, for the tested MPC windows of 4 h and 10 h, the deviation of worst-case scenario PV power generation values from the forecasted values (Ω*P*PV ) represents, on average, 7% and 9% of the data's mean, respectively. When it comes to the grid load, the deviation of worst-case scenario values from forecasted ones Ω*Pcons* is less notable. For the tested MPC windows of 4 h and 10 h, it represents 3% an 2% of the data's mean, respectively.

Having said that, forecasted values of water demand and the ones corresponding to worst-case scenarios are virtually identical (Ω*Qw*,*out* is virtually null). This observation reaffirms the presumption that water demand, and by extension water levels in the water tower's storage tank, have no direct impact on voltage fluctuations in the power distribution grid. Their influence resides solely in determining the water tower's capacity in absorbing excess power off the grid at a given time, which can be properly foreseen through proper dimensioning of the storage tank.

The instances of voltage overshooting decrease steadily from Case 1 through Case 3. In fact, the amended MPC scheme with the min–max problem (Case 3) brings down their percentage (*ν*) to 4.96% and 6.35% for the 4 h window and the 10 h window, respectively, from an initial value of 23.71%. The total surface area of voltage overshooting Φ is also considerably reduced from the initial value. It brought down to 1464.5 kV and 1176.8 kV for the 4 h window and the 10 h window, respectively, from an initial value of 4371.4 kV.

The gain procured through the addition of the min–max problem to the MPC scheme is deduced by comparing the metrics of Case 2 and Case 3. As it happens, for the 4 h sliding window, voltage overshooting is further decreased from Case 2 to Case 3 by 168.2 kV, which amounts to 3.8% of the total surface area of voltage overshooting in the initial case. Percentage-wise, this decrease corresponds to 0.4% of the initial instances of overshooting.

For the 10 h sliding window, voltage overshooting is decreased from Case 2 to Case 3 by 287.8 kV, which amounts to 6.6% of the total surface area of voltage overshooting in the initial case. Percentage-wise, this decrease corresponds to 0.5% of the initial instances of overshooting. The small fraction of the eliminated instances of overshooting with respect to the corresponding percentage of reduced surface area suggests that the min–max problem is particularly apt in eliminating major overshooting incidents. Table 1 reveals that the drop in the total surface area of voltage overshooting observed between the scheme using a 4 h window and the one using a 10 h window is accompanied by an increase in the percentage of instances of overvoltage. Figure 10 illustrates the extrema of voltage fluctuations within the power distribution grid for the standard MPC strategy and the one using a min-max problem, for both sliding window sizes (4 h and 10 h). Voltage overshooting is considerably reduced in both cases with respect to the initial case. Voltage values mostly remain within the acceptable voltage bounds and veer closer to the nominal value (230 V). Unfortunately, this is achieved at the expense of the smoothness of the voltage curves. A possible solution to this issue could be the addition of a regularisation term to the objective function in order to penalise high-frequency voltage fluctuations. That being said, several voltage fluctuations are eliminated thanks to the addition of the min–max problem to the control strategy. This is especially noticeable for the MPC scheme using a 10 h window where notable overshooting incidents are removed during midday of 16 April and 18 April.

**Figure 10.** Extrema of voltage fluctuations within the power distribution grid for the standard MPC strategy compared to the one using a min–max problem and to the initial case, displayed for a 4 h sliding window (**top**) and a 10 h sliding window (**bottom**).

The addition of the min–max problem does not introduce additional computational burden to the control strategy. In fact, the computational complexity (measured by *κ*) decreases from Case 2 to Case 3 by 2.9% for the scheme using a 4 h window and by 6.8% for the scheme using a 10 h window. Figure 11 displays the evolution of the gap between power supply and demand for MPC schemes with and without the min–max problem. It is clear that neither scheme succeeds in reducing this gap significantly. However, they trim the peak occurring everyday around noon, due to the peak in PV power generation. This trimming effect is more visible for the scheme using a 10 h window than the one using a 4 h window. This is reflected in the the final values of the objective function. Though reduced from the initial case, they change very little between Case 2 (MPC using GPR forecasts) and Case 3 (MPC using GPR forecasts and the min–max problem). This suggests that the min–max problem does not provoke any degradation to the MPC strategy's ability to reduce the gap between supply and demand in the power distribution grid.

**Figure 11.** Gap between power supply and demand within the power distribution grid for the standard MPC strategy compared to the one using a min–max problem and to the initial case, displayed fora4h sliding window (**top**) and a 10 h sliding window (**bottom**).

In the case study considered here, it can be seen that the MPC strategy's reduction of the gap between supply and demand remains humble. This may be traced back to the dimensioning of the flexible assets, especially the dimensioning of these assets' storage units. One could argue that the power generation capacity of the chosen biogas plant and the power demand capacity of the chosen water tower are too small to have any meaningful impact on the reduction of the supply/demand gap within the power distribution grid. This observation further illustrates the importance of optimal dimensioning of flexible assets in order for the smart management scheme to yield efficient results.

#### *5.4. Discussion*

The introduction of the worst-case scenario approach, detailed in Section 4.3, is inspired by previous works on min–max optimisation for uncertain nonlinear systems under constraints, which is by definition a conservative risk aversion technique. It is chosen for its relative ease of implementation and is used in this paper to complement the previously developed predictive control strategy in order to anticipate the worst-case scenario, in terms of stochastic input values, and minimise resulting voltage overshooting. The conservativeness of the technique is lessened by its containment to the following time step instead of the entire forecast horizon. This is done primarily to reduce the computational burden added by the min–max problem, in light of the real-time aspect of the application at hand. Nevertheless, the merit of extending the min–max problem to the entire forecast horizon and a quantification of its added cost are valid questions worth investigating in a follow-up to this work. It is worth mentioning that, on top of reducing voltage overshooting, the min–max problem has virtually no effect on the reduction of the gap between power supply and demand.

It is clear that both the growing power demand and the deployment of distributed generation within power distribution grids are not slowing down in the near future. Therefore, configurations like the one studied herein will likely materialise in the upcoming years, accompanied by emerging constraints such as the ones observed in this study. This goes to show the pertinence of studies such as this one in order to prepare for this inevitability.

The damper on the control scheme's performance in the configuration studied in this paper can be traced back, at least partially, to the flexible assets' dimensioning and their suitability to the task. As grid load and PV power generation levels rapidly rise, the flexible assets' room for manoeuvring diminishes. Consequently, two critical questions transpire. The first question is that of optimal dimensioning of these flexible assets. In this study, the control strategy operates at the level of the MV/LV transformer. However, in an application with finer spatial granularity, this question also encompasses optimal placement of the assets within the grid. The second question is that of the assets nature and their suitability to the application.

The combination of a biogas plant and a water tower in this case study was selected in an attempt to utilise small-scale assets compatible with the rural setting of the "Smart Occitania" project and offer the possibility of deferment of the assets' operation. Having said that, the examination of this duo's potential reveals several flaws. To begin, the water tower's ON/OFF operation adds computational complexity to the optimisation problem and is therefore a handicap to the real-time aspect of the applications. Besides, it infers choppier setpoints, which not only worsen voltage fluctuations but also shorten the equipment's life expectancy.

This is especially problematic when the installation's latitude in terms of storage levels is limited. Furthermore, flexible assets need to be extensible in order to adapt their room for manoeuvring to the rapidly growing power demand and distributed generation levels within the grid with minimal cost. The assets considered herein are not easily extensible, particularly the biogas plant, whose energy source is based on a fairly delicate organic process.

In the case of model-based predictive control applications for power grids, choosing the time step is a pivotal task with no definitive answer. A compromise is always made between the high computational cost of this type of control scheme and the granularity of the model, which allows us to capture a maximum of phenomena occurring in the system. This type of control strategy is also dependent upon access to data, in real-time, which comes with its own set of technical issues. Solutions to these issues are starting to come together through the maturing of advanced metering infrastructures in recent years.

The 10 min time step considered in this work is very much an instance of the aforementioned compromise. It allows the necessary computations of both the forecast module and the optimisation problem to run their course. However, it limits the strategy's visibility into the high-dynamics of power grids and thus makes it impossible to intervene between two sampling times. This type of strategy can therefore be seen as an-upper level control scheme, to be coupled with longer-term planning strategies and lower-level operation methods that have the capacity to react to rapid electrical phenomena, namely methods that fall under the umbrella of electrotechnical engineering.

#### **6. Conclusions and Prospects**

The work presented in this manuscript falls within the scope of the "Smart Occitania" project, whose goal is to demonstrate the feasibility of the smart grid concept for rural and suburban low-voltage power distribution grids. To this end, a simulated case study is constructed, based on data collected during the project's run and made available by ENEDIS, in order to elaborate a predictive control strategy for more efficient management of power flows within a power distribution grid with prolific levels of distributed generation, namely PV power generation.

The premise of the proposed strategy is to use a model-based predictive controller to optimise setpoints of flexible assets present in the power distribution grid, in order to reduce the gap between supply and demand. This optimisation is constrained by pre-defined acceptable voltage margins, in addition to the assets' operational restrictions.

A forecast module is constructed using Gaussian process regression and provides intraday forecasts of grid load, water demand, and GHI from which PV power generation is inferred. The quality of the forecasts decreases as the forecast horizon grows longer, but quickly stabilizes around a constant error value. The GPR models also provide confidence intervals associated with the forecasts. Herein, the computed confidence intervals correspond to a probability of 95% of containing the measurements.

An evaluation is carried out of the MPC strategy's resilience to forecasting errors, and the induced errors are quantified. Results show that the predictive control strategy is inherently resilient to forecasting errors as the final objective function value varies little between the case where measurements are used and the one where GPR forecasts are used.

Finally, a min–max problem is added upstream of the main optimisation problem. Its purpose is to anticipate and minimise the voltage overshooting resulting from forecasting errors. In this min–max problem, the feasible space defined by the confidence intervals of the forecasts is searched, in order to determine the worst-case scenario in terms of constraint violation, over the next time step. Then, the main optimisation problem incorporates this information into its decision-making process. Results show that these incidents are indeed reduced thanks to the min–max problem, both in terms of frequency of their occurrence and the total surface area of overshooting.

There are several axes along which the present work can make headway. To begin, the optimisation problem could be modified to take into account the implementability of the flexible assets' setpoints. For example, this could take the form of a multi-objective optimisation that balances out sometimes-conflicting goals.

In addition, the min–max problem integrated into the predictive control strategy in order to improve the scheme's resilience to forecasting errors, can be extended from focusing on the next time step to span the entire forecast horizon. The length of the min–max problem's forecast horizon would be another parameter to be optimised. This will inevitably increase the computational burden of the control scheme, but should also enhance its performance through a better anticipation of issues that may arise along the forecast horizon, especially in light of the degradation of the forecasts' quality as the algorithm advances into the forecast horizon.

**Author Contributions:** Conceptualization, N.D. and J.E.; methodology, N.D. and J.E.; formal analysis, N.D. and J.E.; investigation, N.D., J.E., S.T. and S.G.; writing—original draft preparation, N.D.; writing—review and editing, J.E., S.T. and S.G.; supervision, J.E., S.T. and S.G.; project administration, S.G.; funding acquisition, S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the French agency for ecological transition (ADEME).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors would like to thank the French agency for ecological transition for its financial support. They also thank all the academic and industrial entities involved in the "Smart Occitania" project for their contribution to this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

ADEME French agency for ecological transition ANN Artificial neural network


#### **References**

