This article presents a receding-horizon SMPC approach for energy management in a microgrid with local non-dispatchable renewable energy generation, a battery storage system, and an uncertain non-controllable electrical load. We use a modified linear forecasting model based on the fictitious input approach of Lauricella et al. (2020) to predict the future PV energy generation and the electrical load using only 2 weeks of data for model training. Then, we apply the modified DNCC approach of Ciftci et al. (2019a) to the residuals of the prediction models evaluated over 1 week of data. This allows us to obtain an algebraic constraint formulation that satisfies operational constraints in face of uncertainties to a pre-defined confidence level. The obtained algebraic constraints are then used to formulate a stochastic model predictive control for optimal energy management in a deterministic way.

The main contributions of this work can be listed as follows:

• We applied, for the first time, the forecasting method of Lauricella et al. (2020) based on the fictitious input in a receding-horizon fashion to the SMPC framework, showing both its accuracy on small datasets and its advantages for optimal control applications.

This paper is organized as follows. Section 2 presents the motivation and approach outline, and Section 3 describes the model learning procedure for time series forecasting. Section 4 presents our formulation of the data-driven non-parametric chance-constrained approach for energy optimization problems, while Section 5 introduces its practical implementation within the SMPC framework, and Section 6 discusses the numerical results. Finally, Section 7 concludes the paper and provides directions for future work.

We consider a multiple-input single-output ARX prediction model, where we have a linear dependency between future predictions of a signal with its past occurrences and with appropriate input signals. Let us denote with y ̂ θ ( k | t ) the predicted value of a time series at time step k, given the data available at time t, the model parameter vector θ, and the length of the prediction horizon K. Then, we obtain the following prediction model: y ̂ θ k | t = φ k − 1 T ⋅ θ , k ∈ t + 1 , … , t + K (1) The vector θ ∈ R n y + n u ⋅ d u + n f denotes the model weights, and φ ( k ) ∈ R n y + n u ⋅ d u + n f is the regressor, defined as φ k = Y r k T , U e k T , U f k T T (2) The regressor Eq. 2 consists of an auto-regressive part denoted by Y r ( k ) ∈ R n y and two exogenous components, namely, the external input U e ( k ) ∈ R n u ⋅ d u and a fictitious input U f ( k ) ∈ R n f . The auto-regressive component Y _r and the external input vector U _e are defined as Y r k = y k , … , y k − n y + 1 T U e k = u e k T , … , u e k − n u + 1 T T (3) Since we use the one-step-ahead predictor Eq. 1 recursively, to obtain a multi-step forecast of the load and PV power generation, Y _r is initialized with the most recent n _y past measurements for the first prediction step. Then, Y _r is filled with past predictions as the model is iterated over time until the end of the prediction horizon. For example, to obtain y ̂ ( t + 2 | t ) , the auto-regressive component would become Y r ( t + 1 ) = [ y ̂ θ ( t + 1 | t ) , y ( t ) , … , y ( t + 2 − n y ) ] . The external input vector U _e in (3) contains the most recent n _u past input signals of dimension d _u , with relevant information for the prediction of the considered time series. These can be external temperature and humidity measurements and/or forecast of weather and solar irradiation provided by an external service. To model the (nonlinear) periodic behavior of the considered time series while using a linear predictor, an additional fictitious signal is provided as input to the model. Such fictitious cyclic time encoding consists of sine and cosine functions having a given number of pre-selected time periods, and it is defined as follows: U f k = sin 2 π ⋅ t s k 60 ⋅ 60 ⋅ p 1 , … , sin 2 π ⋅ t s k 60 ⋅ 60 ⋅ p 6 , cos 2 π ⋅ t s k 60 ⋅ 60 ⋅ p 1 , … , cos 2 π ⋅ t s k 60 ⋅ 60 ⋅ p 6 T (4) The decision on the pre-selected time periods is based on the selected non-overlapping windows for the optimization objective in (7) and should be chosen by cross-validation. The variable p represents the period in hours, and for practical implementation on a computational unit, the time step k needs to be converted to t _s (k) having a UNIX timestamp representation in seconds to be correctly inserted in (4). This guarantees that the described model is actually time-invariant.

The prediction model Eq. 1 is trained on historical data by minimizing the simulation error obtained by recursive iteration of the one-step day-ahead predictor given the initial conditions Y r 0 = [ y ( t ) , … , y ( t − n y + 1 ) ] T over the prediction horizon K for each time step t ∈ N. This is done to train the model for forecasting in a receding-horizon fashion, i.e., generating new forecasts at every control step t of the SMPC. Assuming we have N available data points, the model training procedure can be described with the following optimization problem: θ * = arg min θ ∑ t = 0 N − K − 1 ∑ k = t + 1 t + K ‖ y t + k − y ̂ θ k | t ‖ 2 2 (5) Note that (5) is a nonlinear optimization problem since, by iterating the prediction model recursively, the resulting cost has a polynomial dependency on θ. Even if N is small, the dimension of the objective in (5) can grow very quickly depending on the prediction horizon K. To avoid long computation times during the model identification, we relax the objective to only contain D = ⌊ N K ⌋ non-overlapping prediction windows of the horizon K. The floor operation can lead to some information loss at the tail of the dataset; thus, ideally, the dataset should be chosen such that ⌊ N K ⌋ = N K . To simplify the notation, we define the prediction and the measurement vectors as Y ̂ θ t = y ̂ θ t + 1 | t , … , y ̂ θ t + K | t Y t = y t + 1 , … , y t + K (6) Due to the small size of the dataset, we introduce a regularizing term which is added to the objective function to reduce overfitting. The relaxation of (5) together with (6) and the regularization term result in θ * = arg min θ ∑ i = 0 D ‖ Y i ⋅ K − Y ̂ θ i ⋅ K ‖ 2 2 + λ ⋅ ‖ θ ‖ 2 2 (7) A larger value of the hyperparameter λ corresponds to a sparser structure of θ, and a smaller value corresponds to more flexibility in the magnitudes of θ. Larger values of λ avoid overfitting, but have a negative influence on training accuracy since they force a simpler model structure. On the other hand, smaller values of λ allow for closer resemblance of the training data, which could reduce model generality. The tuning of λ is a trade-off between good fitting performance and model generality; thus, it is recommended to tune its value via cross-validation. Since the optimization problem has become a ridge regression, we need to normalize our data to have a zero mean and unitary variance; otherwise, the regularization term would not function properly due to the different magnitudes of the signals. Thus, u ̃ e k = u e k − E u e Var u e , y ̃ k = y k − E y Var y (8) where E [ ⋅ ] and Var [⋅] represent, respectively, the expected value and the variance operators.

Let us denote an uncertain quantity as the random variable X and its distribution function as P X . In the classic chance constraint formulation, the expression that a randomly drawn value X ∼ P X will be less or equal to a value x ̃ with a probability 1 − α can be written as P X X ≤ x ̃ ≥ 1 − α (9) Here, we designate 1 − α the confidence level, and α is the risk level. We assume no prior knowledge of the true distribution function P X . Instead, we want to find an estimate f ̂ X of the true probability density function (PDF) using kernel density estimation (KDE). Such an estimate will then be used to formulate data-driven non-parametric chance constraints which can be used for stochastic optimization.

Given a dataset V = [ x 1 , x 2 , … , x n ] , x i ∈ R of n samples from a probability distribution function P X , an estimate f ̂ X of the true probability density function f _X can be found with KDE as f ̂ X x = 1 n h ∑ i = 1 n K x − x i h (10) Here, h ∈ R denotes the bandwidth, and K(x) is the kernel function. We adopt the Gaussian kernel function K ( x ) = 1 2 π e − 0.5 x 2 , but other kernels such as the Epanechnikov or uniform kernel can be used as well. We choose the Gaussian kernel due to its smoothness characteristic. It has shown smoother and more general distributions on small datasets than other kernels. The choice of the bandwidth has a significantly higher influence on PDF estimation than the choice of the kernel function, and it should be chosen according to the available samples. For simplicity, we chose the bandwidth as h = n − 1 4 according to Scott’s rule (Scott, 1992).

By law of large numbers, the estimated PDF f ̂ X will converge to the true density function f _X as n → ∞ (Tsybakov, 2009). Since we have only limited data samples from the unknown distribution available, the KDE will not converge to the true distribution, which results in an uncertainty in the estimation.

Considering that we can only estimate the true PDF using a finite dataset, we know that our estimation will have errors. Jiang and Guan (2016) propose a quantification of these errors and a method to include them in the chance constraint formulation. Here, the estimated distribution function P ̂ X is assumed to lie within a confidence set D of probability distribution functions P , whose ϕ-divergence to the estimated density function is less than a certain confidence set size d. D is defined as D = P ∈ M + : D ϕ f ‖ f ̂ X ≤ d , f = d P d x (11) In this work, we only consider 1-dimensional random variables; therefore, X ∼ P X , X ∈ R . We denote with M + the set of all probability distributions. Then, the divergence D ϕ ( f ‖ f ̂ ) is defined as D ϕ f ‖ f ̂ = ∫ R ϕ f x f ̂ X x f ̂ X x d x (12) To account for the uncertainty in the probability distribution estimation, the classic chance constraint Eq. 9 is reformulated to hold for every distribution within the confidence set: inf P ∈ D P X ≤ x ̂ ≥ 1 − α (13) To reformulate the chance constraint that holds for all P ∈ D , the risk level α is reduced to α′ according to the confidence set size d. Here, we choose the ϕ-divergence to be the χ-divergence of order 2, defined by ϕ(x)≔(x − 1)² for risk levels α ≤ 1/2. This choice allows for a closed-form derivation of the reduced risk level α′ (Jiang and Guan, 2016). Then, the reduced risk level α′ is given by α ′ = α − d 2 + 4 d α − α 2 − 1 − 2 α d 2 d + 2 (14) To avoid negative values for the risk level, the final risk level α + ′ is chosen as α + ′ = max ( α ′ , 0 ) . Note that in reality, the risk level rarely becomes negative, and in this work, it is not negative as well. This formulation is mainly made for consistency. With the reduced risk level, the chance constraint can be reformulated as P ̂ X X ≤ x ̃ ≥ 1 − α + ′ (15) Now, the chance constraint depends only on our estimated distribution and the confidence set of the estimated distribution. This chance constraint can now be easily reformulated to an algebraic representation by calculating the cumulative distribution function F ̂ X and then using its inverse F ̂ X − 1 (the quantile function) to determine x ̃ given a confidence 1 − α as F ̂ X x ̃ = P ̂ X X ≤ x ̃ ≥ 1 − α + ′ ⇔ x ̃ ≥ F ̂ X − 1 1 − α + ′ (16) Evaluating F ̂ X − 1 ( 1 − α + ′ ) yields a value x ̃ that will be greater/equal than any random variable X drawn from any distribution in D , with a probability of at least 1 − α. This makes sure that the user-defined confidence 1 − α is realized with respect to the uncertainty in the estimation P ̂ X .

As the last step of the DNCC approach, we need to determine a confidence set size d to construct the reduced risk level. Ciftci et al. (2019a) proposes using point-wise errors between the estimated density function f ̂ X and a histogram of the available data. Since the performance of point-wise errors with histograms can change significantly depending on the chosen histogram bins, we use a different approach to construct confidence intervals of the estimated PDF, as proposed by Fiorio (2004).

The general confidence interval for the estimated density function f ̂ X of the probability distribution of X can be written as E f ̂ X x ∈ f ̂ X x − u 1 − α / 2 ⋅ s , f ̂ X x − u α / 2 ⋅ s (17) Since the variance of the estimation is unknown, the sample variance s ² is used to construct the confidence interval. The values u _1−(α/2) and u _α/2 are the 1 − (α/2) and α/2 quantiles of the t-statistic defined in (18), respectively. t X x = f ̂ X x − E f ̂ X x s x (18) With the t-statistic, the quantiles are defined as P ( t X ≤ u α / 2 ) = α / 2 and P ( t X ≤ u 1 − ( α / 2 ) ) = 1 − ( α / 2 ) , respectively. It is important to make sure that the t-statistic in (18) is unbiased, such that u _α/2 ≤ 0 and u _1−(α/2) ≥ 0. In Fiorio (2004), the unbiased t-statistic is obtained by applying bootstrapping to the estimate f ̂ X as follows.

Let x _i ^* denote the bootstrapped data samples obtained by n times random sampling with replacement of data points from the original dataset V . With the total number of samples denoted as n, the bootstrap KDE can be formulated, similarly to Eq. 10, as f ̂ X * x = 1 n h ∑ i = 1 n K x − x i ′ h (19) The sample variance of the estimated density function f ̂ X is given by (20), and its bootstrap analog s X 2 * ( x ) by replacing the data points x _i with the bootstrapped data points x i ^*, and the estimated f ̂ X with the bootstrapped estimate f ̂ X * . s X 2 x = 1 n h 2 ∑ i = 1 n K x − x i h 2 − f ̂ X x 2 n (20) With the bootstrap estimate f ̂ X * and the bootstrap sample variance s X 2 * , the unbiased t-statistic t X * ( x ) is defined as t X * x = f ̂ X * x − f ̂ X x s X * x (21) Now, the quantiles can be obtained from the unbiased t-statistic with P ( t X * ≤ u α / 2 * ) = α / 2 and P ( t X * ≤ u 1 − ( α / 2 ) * ) = 1 − ( α / 2 ) . These quantiles together with the sample variance s X 2 can finally be used to construct the confidence interval for the estimate f ̂ X : E f ̂ X x U = f ̂ X x − s X x ⋅ u α / 2 * E f ̂ X x L = f ̂ X x − s X x ⋅ u 1 − α / 2 * (22) Finally, we use the obtained confidence intervals to construct the confidence set size of our distribution function estimate P ̂ X . The confidence set size d is defined as d = E f ̂ X x U − E f ̂ X x L 2 1 − α (23) Here, the subscript denotes the 1 − α quantile of the point-wise squared distances between the upper and the lower confidence bound.

Let t be the time instant for which a control decision has to be made, as described in Algorithm 1. Then, by implementing the receding-horizon SMPC with a fixed prediction horizon of length K, the resulting controller will plan the controllable inputs for all the next time steps k ∈ [t, t + K]. Then, only the first control decision corresponding to the time step k = t is implemented. The entire procedure is then iteratively repeated for all the following time steps within the current week of operation. All power variables are defined to be non-negative: P g k , P cur k , P b c k , P b d k ≥ 0 , ∀ k ∈ t , t + K (24)

Since the goal is to minimize the total grid energy import cost over the considered prediction horizon K, the objective function for the SMPC can be formulated as ∑ k = t t + K P g k ⋅ C g k (25)

The battery storage system here is described by a simple state-space model with one state representing the battery state of charge. It is assumed that the battery model is deterministic, i.e., control decisions regarding charging and discharging power will be executed exactly as mandated by the SMPC. We neglect the battery’s internal dynamics and all the external disturbances acting on the battery storage system. Let E _soc (k) denote the state of charge of the battery at time step k in kWh; η _c and η _d the charging and discharging efficiencies in percent, respectively; and Δt the duration of one time step in hours. Hence, the battery model is described by E soc k + 1 = E soc k + η c P b c k − 1 η d P b d k ⋅ Δ t , ∀ k ∈ t , t + K (26) The state of charge at the first time step k is initialized with the realized state of charge as available from the previous time step t − 1. To protect the battery and increase its lifespan, we add lower and upper bounds for the state of charge in the optimization problem: E soc min ≤ E soc k + 1 ≤ E soc max , ∀ k ∈ t , t + K (27) Furthermore, the battery charging and discharging capacities are limited as well, and to ensure that the battery is not charged and discharged at the same time, we introduce the binary decision variables δ _c and δ _d . Here, δ _c (k) = 1 means that the battery will be charged at time step k; therefore, δ _d (k) has to be 0 for the same time instant. These constraints are formulated as follows: 0 ≤ P b c k ≤ P b c , max ⋅ δ c k 0 ≤ P b d k ≤ P b d , max ⋅ δ d k 0 ≤ δ c k ≤ 1 , ∀ k ∈ t , t + K 0 ≤ δ d k ≤ 1 δ c k + δ d k ≤ 1 (28)

The MPC controller is required to maintain the balance between the overall power generation and consumption of the microgrid over the prediction horizon K. This operational requirement results in the following equality constraint: P g k + P b d k + P p v k = P l k + P b c k + P cur k , ∀ k ∈ t , t + K (29) In (29), P _g (k) and P _cur (k) are considered state variables, P b d ( k ) and P b c ( k ) are control variables, and P _pv (k) and P _l (k) denote the true photovoltaic power generation and the load consumption, respectively, which are both uncontrollable and unknown variables for the MPC. The equality constraint in (29) can be reformulated to include the forecasts of the uncertain variables instead of their true but unknown values. Here, we assume that P _pv (k) and P _l (k) can be expressed as their predictions plus a residual term as P p v k = P ̂ p v k + P ̃ p v k P l k = P ̂ l k + P ̃ l k (30) Here, P ̂ p v and P ̂ l denote, respectively, the PV and load forecast obtained with the prediction model described in Section 3, while P ̃ p v and P ̃ l denote the error terms (i.e., the residuals). To capture the stochastic nature of load and PV generation, we assume that the error terms are stochastic variables following unknown probability distributions. We adopt the procedure presented in Section 4 to estimate such probability distribution on historical error data obtained by comparing the predictions with their true realizations from the available validation dataset.

We assume that the prediction for time step k at the same time t will have similar errors over different days, i.e., P ̃ ( k | t ) ∼ P ̃ ( k | t + 24 h ) . Therefore, the error P ̃ ( k | t ) of the prediction for time k at time step t will follow the same probability distribution every day: [ P ̃ ( k | t ) , P ̃ ( k | t + 24 h ) , … ] ∼ P P ̃ ( k | t d ) . Considering the division of a day into T _d time steps per day and, therefore, t _d ∈ T _d , this results in a total of T _d ⋅ K probability distributions that need to be estimated from data.

Given this probabilistic definition of the prediction errors, we can now define a chance constraint to satisfy the power balance in (29) for the true load and PV realizations with a pre-defined confidence level. By substituting in (29) the true load and PV generation with their forecast plus residual counterpart from (30), we obtain the following equality constraint that depends on the random variables P ̃ l ( k ) and P ̃ p v ( k ) : P g k + P b d k + P ̂ p v k + P ̃ p v k = P ̂ l k + P ̃ l k + P b c k + P cur k (31) To ensure that this equality holds with a certain confidence 1 − α, the following chance constraint can be formulated using an inequality that will be later replaced in the optimization with an equality. P ( P g k + P b d k + P ̂ p v k + P ̃ p v k ≥ P ̂ l k + P ̃ l k + P b c k + P cur k ) ≥ 1 − α , (32) Notably, the time step notation t is omitted. Since P ̃ l ( k ) and P ̃ p v ( k ) have different underlying distributions, we cannot directly reformulate the upper chance constraint into an algebraic formulation. To mitigate this problem, we introduce a “dummy variable” d _pv , as suggested by Ciftci et al. (2019a). With d _pv , first, the probability distribution of the PV residuals is evaluated, and then the resulting value for d _pv can be inserted in the chance constraint in (32), which allows the evaluation of the probability distribution of the load residuals, resulting in a full evaluation of the chance constraint. P P ̃ p v k ≤ d p v k ≥ 1 − α P ( P ̃ l k − d p v k ≤ P g k + P b d k + P ̂ p v k − P ̂ l k − P b c k − P cur k ) ≥ 1 − α (33) Finally, the chance constraint can be reformulated as an equality algebraic constraint that can be easily included in the MPC optimization problem by using the DNCC method described in Section 4. F P ̃ l k − 1 1 − α + ′ − F P ̃ p v k − 1 1 − α + ′ = P g k + P b d k + P ̂ p v k − P ̂ l k − P b c k − P cur k (34) This chance constraint ensures that the SMPC will plan the battery operation such that there will be enough power generated or consumed in the system to match the true load and PV generation with a confidence of 1 − α while minimizing the overall energy cost.

By combining the objective function defined in Eq. 25, the battery dynamics and constraints of Eqs. 26-28, the power balance chance constraint Eq. 34 with p ′ = 1 − α + ′ , and the non-negativity constraints Eq. 24, the final SMPC formulation is min P g , P cur P b c , P b d , E soc δ c , δ d ∑ k = t t + K P g k ⋅ C g k s . t . P g k + P b d k + P ̂ p v k − P ̂ l k − P b c k − P cur k = F P ̃ l k − 1 p ′ − F P ̃ p v k − 1 p ′ . E soc k + 1 = E soc k + η c P b c k − 1 / η d P b d k ⋅ Δ t , E soc min ≤ E soc k + 1 ≤ E soc max , 0 ≤ P b c k ≤ P b c , max ⋅ δ c k 0 ≤ P b d k ≤ P b d , max ⋅ δ d k , ∀ k ∈ t , t + K 0 ≤ δ c k ≤ 1 0 ≤ δ d k ≤ 1 δ c k + δ d k ≤ 1 P g k , P cur k , P b c k , P b d k ≥ 0 (35)

The models for load and photovoltaic forecasting are trained on the publicly accessible dataset of a small microgrid by TronderEnergi (2021). This dataset includes measured values for weather, load, PV, and prices. Furthermore, it includes dimensioning hyperparameters for the battery. To evaluate the SMPC, we use historic weather forecast data instead of historic measured data to have a better representation of uncertainties. The historic weather forecast data are extracted from the publicly available data from the Norwegian Meteorological Institute METNorway (2023). Both datasets are sorted, cleaned, and merged together in one large dataset. The data are transformed to a time resolution of 15 min, and every data point is accompanied by a timestamp. A description of the final dataset can be found in Table 1.

The main assumption for the simulation is that the power balance equation in (29) must hold for the true load P _l and true PV generation P _pv at all times. Next, as mentioned previously in Section 5.3, the battery system is assumed to be deterministic, i.e., the control decision by the SMPC for time step t, denoted as P b c * ( t ) and P b d * ( t ) , will be exactly executed as such, and the state of charge of the battery will exactly follow the dynamics. As battery charging/discharging power and load and PV generation are fixed values, the simulation will decide the true grid import P _g and the true curtailed power P _cur at the time of simulation. The simulation results at time step t can be described by the following equations: P bal t = P l t − P p v t − P b d * t + P b c * t P g t = P bal t P bal t ≥ 0 0 P bal t < 0 P cur t = 0 P bal t ≥ 0 − P bal t P bal t < 0 (36)

All computations are carried out on a MacBook Pro Late 2013 with a 2 GHz Quad-Core Intel Core i7 and 8 Gb of RAM.