Accelerated Model Predictive Control for Electric Vehicle Integrated Microgrid Energy Management: A Hybrid Robust and Stochastic Approach

Abstract

A microgrid with an advanced energy management approach is a feasible solution for accommodating the development of distributed generators (DGs) and electric vehicles (EVs). At the primary stage of development, the total number of EVs in a microgrid is fairly small but increases promptly. Thus, it makes most prediction models for EV charging demand difficult to apply at present. To overcome the inadaptability, a novel robust approach is proposed to handle EV charging demand predictions along with demand-side management (DSM) on the condition of satisfying each EV user’s demand. Variables with stochastic forecast models join the objective function in the form of probability-constrained scenarios. This paper proposes a scenario-based model predictive control (MPC) approach combining both robust and stochastic models to minimize the total operational cost for energy management. To overcome the concern about the convergence time increasing from the combination of scenarios, the Benders decomposition (BD) technique is further adopted to improve computational efficiency. Simulation results on a combined heat and power microgrid indicate that the proposed scenario-based MPC approach achieves a better economic performance than a traditional deterministic MPC (DMPC) approach, while ensuring EV charging demands, as well as minimizing the trade-off between optimal solutions and computing times.

1. Introduction

Introduced to streamline the operation of increasing penetration of distributed generations (DGs) in distribution networks, microgrids are capable of providing a systematic scheme to integrate local dispatchable DGs, such as combined heat and power (CHP) units, renewable energy systems (RESs), energy storage systems (ESSs), electric vehicles (EVs), and require intelligent energy management systems subject to some criteria, such as economic performance, user satisfaction, reliability, etc. [1,2]. Energy management approaches are complex tasks which require dynamic assessment of current and future operations while considering time-varying energy consumption, production, and prices [3,4,5]. Meanwhile, by treating EV charging requests as postponable electrical loads under some appropriate demand-side management (DSM) strategies, the energy management of microgrids can be more flexible and cost-effective [6,7].

A promising domain for energy management is methods related to model predictive control (MPC). Being one of the most well-developed control techniques in the process industry, MPC deals with multiple-input-multiple-output systems quite well, and strong attention has been paid to its potential of reducing the operational costs of microgrids [8]. The MPC-based approaches employ a receding horizon strategy, which keeps the total operational cost low without violating the supply-demand balance and the pre-defined operation constraints. Prediction models for variables related to microgrids, such as RES outputs, real-time energy prices, load demands, etc., make it possible to implement MPC into energy management. While explicit prediction models have been used over a deterministic MPC (DMPC) approach, uncertainty models conveying detailed information about variable forecasts have drawn more attention lately. Two uncertainty models have been successfully integrated into MPC, namely robust models and stochastic models. Robust models treat uncertainties as one set-bounded scenario irrespective of scenario probabilities, providing rapid but conservative decisions while ensuring the robustness of the optimization when prediction models are not effective enough [9]. Stochastic models minimize the expected value of the objective function considering the probability-constrained scenarios of predictions, but are usually associated with large volumes of variables which raise challenges in the computing time for the realization of energy management [10]. To avoid intractable computation, continuous distributions in stochastic models are normally discretized into a finite number of scenarios, called scenario-based MPC, where each scenario represents a possible realization of the uncertain variables [11].

This paper presents a novel way of handling various uncertainty characteristics in MPC-based energy management combining both stochastic and robust models. For prediction models without sufficient information on probabilities, such as the charging load of diffusing EVs in a microgrid at the present stage, an original robust approach is created to decide an upper-bound scenario of each EV’s charging load, and to couple with the maximal postponement constraints. Variables with relatively mature prediction models, such as wind speeds and real-time electricity prices, are processed as stochastic models with scenario generation and combination. The objective function minimizing the total system’s operational cost in the prediction horizon is modified into a two-stage stochastic programming problem containing both robust scenarios and stochastic scenarios while the computational efficiency can be improved by implementing Benders decomposition (BD).

This paper is organized as follows: the remainder of Section 2 briefly reviews the present MPC approaches for energy management and introduces the scheme of the proposed approach; in Section 3, a general MPC-based energy management approach with uncertainties is marked as the inspiration of the hybrid robust and stochastic variables; Section 4 considers the implementation of the proposed approach in an exemplary combined heat and power microgrid model, as well as the incorporation of the BD technique; for demonstration purposes, Section 5 presents some comparative results of several approaches; and, finally, Section 6 provides a brief conclusion.

2. Motivation

2.1. Related Works

Different from the industrial control domain, where MPC is generally used to achieve the stability of a specified equilibrium point or track a reference signal, one of MPC’s modern philosophies is to bring down operational costs by calculating a sequence of optimal control inputs over a finite prediction horizon with future information while satisfying the input and state constraints [12]. Only the first control input of the sequence is applied to the system, and since MPC is implemented in a receding horizon manner, the same procedure is executed with updated measurements at the next sampling time. From a performance point of view, the uncertainties only influence the optimality of the objective cost function without affecting the stability and feasibility properties. Early DMPC energy management algorithms do not take into account prediction uncertainties directly; as research continues, further studies focusing on uncertainty-integrated MPC are further brought up for performance improvement.

Robust models have emerged to minimize the worst-case scenario cost with the inherent advantage that only limited information of the uncertainty, such as the lower or upper bound, is required, which is easier to obtain or estimate in practice [13]. Several robust models have been successfully used for energy management, and focus on representing forecast errors to enhance robustness in general. A DSM project is scheduled to incorporate robustness against the uncertainties in electricity prices and renewable resources to optimize power consumption in [14]. Robust parameters are introduced to control the extent of uncertainty resulting from the real-time electricity price in addition to the stochastic model, while minimizing the total operational cost in [15]. When discussing EV charging load profiles, although quantitative models are generated based on statistical information, the rather small sample size and the ongoing fast diffusion of EVs at the microgrid scale lead to significantly varying characteristics, making it impossible to implement a proper parametric forecast of EV fleets [16]. However, most DSM strategies considering EV demands assume that a parametric distribution or registration of the traveling schedule is needed in the prediction horizon [17,18]. Therefore, robust models are highly applicable where the explicit evaluation of the EV charging load is difficult to obtain.

On the other hand, forecast models, such as wind speeds, solar radiation, real-time electricity prices, etc., are recognized as relatively accurate, and often of a probabilistic nature [19,20,21]. Stochastic models minimizing the total expected cost have emerged as an alternative methodology when probabilistic descriptions of uncertainties are readily accessible. A stochastic energy management model for microgrids demonstrating a promising result in applying uncertainty scenarios associated with renewable energy resources is shown in [22]. A stochastic operational scheduling method for a distribution system considering a combination of wind power and demand uncertainties is demonstrated in [23]. Concerns of a reasonable computation time for these scenario-based MPC applications gradually increase since the computation size grows exponentially with the number of uncertainty scenario selections, as well as the length of a designated prediction horizon. Several scenario reduction techniques have been developed in the past decade to reduce the number of scenarios without sacrificing their accuracy to a large extent, such as the probability-based scenario reduction method in [24], and the forward selection in wait-and-see clusters method in [25]. From another perspective, the computational pressure can be released by reducing the convergence time by decomposition-and-parallelization–based solution techniques. In [26], a widely cited example of applying the primal BD from [27] for scheduling problems is proposed without enumerating the full set of states.

2.2. Novel Contributions

In this paper, we attempt to overcome those limitations mentioned in Section 2.1 by incorporating scenario-based MPC with forecasts on the various uncertainties. This paper presents an innovative energy management approach for optimizing the economic performance of a combined heat and power microgrid with EV integration. The key figure of merit is to include the predictable but uncertain variables in a receding horizon manner while retaining the computational complexity within the scope of acceptability. By building on a modified scenario-based MPC algorithm together with stochastic, robust, and deterministic models, different feasibilities of predicted variables are able to match with an MPC-based approach.

For variables where no accurate prediction models are appropriate, i.e., the EV charging requests with a small base and fast diffusion dynamics, a novel robust approach is created to ensure the charging demand satisfaction of every EV user under the demand-side management strategy, while achieving a minimized worst-case scenario cost for the microgrid. A state criterion is brought up to establish the upper-bound scenario of EV charging requests considering actual-connected EVs and virtual-connected EVs to all of the installed EV service equipment (EVSE) in the modeled microgrid. The maximal postponement constraint for the DSM of each postponable EV charging load is included as the boundary to ensure robustness. For variables with prediction models consisting of probability distributions, i.e., the real-time electricity price and RES outputs, stochastic models are used to represent these uncertainties with scenarios and their combination as well. The electricity demand is divided into two parts: the general electricity demand and the EV charging demand, respectively. A deterministic approach is also adopted for general electricity and heat demands of the microgrid to show the generality of MPC-based energy management. These forecasted quantities feed the scenario-based MPC algorithm with the solving process accelerated by BD.

The procedure of the proposed energy management scheme in discrete time with a sampling time step Δt is presented in Figure 1. After measuring the working parameters and environment variables in the microgrid at each time step t, these outputs are sent into the central controller, where a prediction horizon [t,t+N] is applied. After obtaining operation instructions of every time step in [t,t+N], perform these decisions to the controllable components only for the current time step. When information is updated in the next time step, the central controller adjusts the real-time energy dispatch decisions based on updated states and the associated future uncertainties.

3. A Hybrid Stochastic and Robust MPC Approach for Energy Management

3.1. Structure Model with Operation Constraints

An MPC-based energy management incorporated with uncertainties ξ is formulated to optimize the economic performance with the objective of minimizing the expected sum of the fuel purchasing cost C_f, the electricity purchasing cost from the external grid C_ex, and penalties P, such as the self-discharge and wear losses for ESSs or start-up costs for CHPs during the prediction horizon [t,t+N]. At every time step t, the following problem is solved:

$$\begin{array}{cc}\inf & \mathrm{E}\left[{\displaystyle \sum _{k=0}^N\left(C_{f,t+k|t}+C_{ex,t+k|t}+P_{t+k|t}\right)}|{\xi }_{t+k|t}\right]\end{array}$$

(1)

where the symbol U_t+k|t stands for the variable U at time step t + k predicted at time step k. The forecasted uncertainties can be correlated over time; thus, the variable U_t+k|t need not be the same as that of U_t+k+1|t+1. The production costs of RESs (i.e., wind and solar energy) are assumed to be negligible.

Although uncertainties are best modeled as continuous variables for presenting an ideal precise prediction as in Equation (1), integrals are always turned into probability sums based on a finite set of scenarios for the purpose of decreasing high computational complexity based on common sense.

3.2. Robust Approach for EV Charging with Demand-Side Management

3.2.1. Robust Approach for EV Charging Uncertainties

Current approaches specifying EV fleet characteristics neglect the feature that the EV user base is relatively small and of a fast dynamic process since we are still undergoing the early stage of EV diffusion progress. These commonly used parametric distribution characteristics of EV fleets may vary significantly in different microgrids, both at temporal and spatial scales. To overcome the statistical insignificance of EV charging prediction at the current stage of development, a novel robust approach is developed to address EV charging uncertainties for energy management.

Full awareness of prediction information is required for MPC-based energy management in the prediction horizon. For those unknown charging requests, a worst-case scenario is applied by setting a boundary of maximum possible charging loads, where all EVSEs with no EV connection at time step t + k will be plugged in a virtual EV in the following time step t + k + 1. The virtual EV will stay connected to this EVSE until t + N with a requested total charging time that is equal to N − (k + 1). The other worst-case scenario adopted is to deal with uncertainties beyond the prediction horizon for the actually-connected EVs. Despite the fact that the charging loads may be postponable for EVSEs in State 4, since the prediction horizon is less than the parking duration remaining, the charging requests are still treated as unpostponable in the whole [t,t+N]. These two types of worst-case scenario serve as extreme predictions to enhance the robustness.

Despite focusing on the changing numbers of EVs, we turn our attention to the settled EVSEs in a microgrid, with a total number of N_z which remains the same in the prediction horizon. Assume that all of the EVs connected in the microgrid accept DSM as postponable electrical loads. EV owners set up their desired state of energy SOE_z,set with an estimated departure time T_z,set when plugging the EV into the zth EVSE. Additionally, a charging behavior binary variable x_z,t+k|t is defined to represent whether there is a charging load of the zth EVSE at time step t + k predicted at current time step t.

At each time step t, every EVSE measures its EV connection state and parameters automatically. A state criterion with six possible states of the charging process can be generated in the prediction horizon [t,t+N] as shown in Figure 2. Notice that the charging behavior binary labeled in Figure 2 is arranged by the original order before DSM. To begin with, an EVSE detects whether there is an EV connected. If not, this EVSE is regarded to be in State 1. When there is a connected EV, measure the current state of energy SOE_z,t, the charging power P_z, and fetch the value of SOE_z,set and T_z,set. Since time is separated into discrete Δt, T_z,set sometimes need to get rounded to an integer by executing N_z,set = ceil(T_z,set/Δt). If the desired SOE_z,set is achieved, which means the connected EV is fully charged, compare whether the planned departure time step N_z,set precedes the prediction horizon N or not, where State 2 represents that the EV is going to leave the connected EVSE within the prediction horizon, and State 3 otherwise. If SOE_z,set has not been achieved, the difference of SOE_z,set and SOE_z,t denotes the needed sum of charging time step N_z,c before the EV departures as:

$$\begin{array}{l}{T}_{z,c}=\frac{SO{E}_{z,set}-SO{E}_{z,t}}{P_z}=\frac{\Delta SO{E}_z}{P_z}\\ N_{z,c}=ceil(T_{z,c}/\Delta t)\end{array}$$

(2)

If N_z,c is greater than or equal to N, label the EVSE as in State 4. If not, compare the departure time N_z,set and N again; define that an EVSE with N_z,set less than N is in State 5, and otherwise is in State 6.

At the end of this procedure, all EVSEs are able to send their respective states as well as the original charging behavior arrays containing actual demand and worst-case scenarios.

3.2.2. Demand-Side Management for EV Charging

To determine the role of each EVSE in the DSM, the six states are further divided into two groups:

1. Unpostponable group: including State 1, State 2, State 3, and State 4, with a total number of N_z1. The charging loads in the four states are unpostponable within the prediction horizon.

2. Postponable group: including State 5 and State 6, with a total number of N_z2. The charging loads in the two states are postponable without affecting the realization of original charging demand at the end of the prediction horizon. In other words, one or more x_z,t+k|t in these states can be scheduled from 1 to 0 or from 0 to 1 without violating the total incremental charging energy needed in [t,t+N]. The scheduled charging behavior E_z,t+k|t should satisfy the following constraints:

$$\begin{array}{l}{\displaystyle \sum _{t=t}^{t+N}{E}_{z,DSM,t}}=N_{c,z}{E}_z\forall z\in N_{z2}\\ E_{z,DSM,t}\in \left\{0,E_z\right\}\end{array}$$

(3)

Table 1. The assignment and circumstances of the charging behavior binary x_z,t+k|t.

Value	Assignment	Circumstances
0	No charging load	The desired SOEz,set has been achieved; Original charging behavior is arranged to be held by DSM.
1	Positive charging load	Unpostponable charging of actual-connected EVs; Unpostponable charging of virtual-connected EVs; Arranged charging behavior by DSM.

summarizes all possible circumstances of the charging behavior binary variable x_z,t+k|t due to DSM. When there is no charging behavior at time step t + k of the zth EVSE, such as when the desired SOE_z,set has been achieved of the actual-connected EV or the charging load is scheduled to be held by the schedule of the energy management with DSM, the value of x_z,t+k|t is attributed to be 0; when there is an unpostponable or arranged charging load, then x_z,t+k|t is attributed to be 1. In particular, for those EVSEs with virtual-connected EVs, i.e., x_z,t+1|t to x_z,t+N|t of EVSEs in State 1, and x_z,t+k+1|t to x_z,t+N|t of EVSEs in State 2 and 5, the charging demands are considered to be unpostponable; thus, the values of these variables above are also assigned to be 1.

Due to the receding horizon philosophy of MPC, all of the arranged EV charging postponements only apply to the current time step t, and the state criterion will be refreshed in the next time step. This allows the update of remotely setting changes of SOE_z,set and T_z,set from mobile phones by EV users and hypothetical worst-case scenarios by the system itself. In this way, EMS can achieve a minimized worst-case scenario cost while satisfying charging requirement of every EV user.

3.3. Stochastic Approach for Uncontrollable Inputs

For uncertainties with relatively accurate prediction models, stochastic approaches allow for the closer imitation of predictions. The so-called scenario tree method has been successfully applied to stochastic models, such as real-time electricity prices, outputs of WTGs, and other RESs. A scenario tree representing uncertainties with a tree of discrete scenarios is shown in Figure 3a with a general structure. In a scenario tree, each node represents a possible state at its particular time step and has only one predecessor, but multiple successors. The root node represents the initial state, and a forward branch from the root node to a leaf node describes a probability-weighted realization. The probability can be even or not, providing a causal but compact representation of the uncertainties acting on the system. A discrete-time formulation for a scenario tree can be written as:

$$\begin{array}{l}{x}_{s(t+k),t+k|t}=f\left(x_{p(t+k-1),t+k-1|t},u_{s(t+k-1),t+k-1|t},d_{r(t+k-1),t+k-1|t}\right)\\ {\displaystyle \sum {\pi }_{s(t+k),t+k|t}}=1\forall k=1,2,\mathrm{...},N-1\end{array}$$

(4)

where s(t + k) is the total scenario number at time step t + k when predicting at the current time step, and each state x_s(t+k),t+k|t is forecasted conditionally on the paths of the preceding time steps, depending on the parent state x_{p(t+k-1),t+k-1|t}, the corresponding control input u_{s(t+k-1),t+k-1|t}, and the corresponding disturbance d_{r(t+k-1),t+k-1|t}; and the scenario set is denoted by Ω_t+k|t with the sum of the probabilities of π_s(t+k),t+k|t at each state being 1. Note that this paper focuses on representing scenarios of forecast uncertainties without the intention to develop or compare their prediction models.

Multiple stochastic uncertainties are considered in this paper, and thus a scenario combination should be considered further. As shown in Figure 3b, a scenario set Ω_p,t+k|t of the real-time electricity price p_t+k|t, with probabilities of π_p,t+k|t and scenario sets Ω_wi,t+k|t of the ith RES output E_wi,t+k|t of π_wi,t+k|t (i = 1, 2), with the same prediction horizon N, are used to illustrate the scenario combination process at time step t + k. The combination happens at every time step, and each combined scenario consists of the value of a real-time electricity price and RES outputs with a multiplication of their respective probabilities. However, it is noticed that even with a properly reduced number of scenarios, the process of combination still brings about an exponential growth in the total scenario number compared to deterministic approaches, which is related not only to the scenario numbers of every stochastic uncertainty adopted, but also to the length of the prediction horizon. To deal with the computing space requirement and computing time beyond practicality, a parallel computing technique incorporated with the objective function will be introduced in Section 4.3.

4. Model Formulation for a Combined Heat and Power Microgrid

4.1. Combined Heat and Power Microgrid Model Description

A model of a combined heat and power microgrid is presented for case study purposes, with the basic structure with multi-energy flows depicted in Figure 4. For simplicity, the same type of units are represented as one module here.

The operation constraints of these DG units are introduced as follows:

Taking the wind turbine generator (WTG) as the RES modeled here, its output E_w,t is related with the cut-in speed v_in, rated speed v_r, and cut-off speed v_out. Suppose that it is fully used with no wind curtailment in the microgrid:

$$E_{w,t}=\{\begin{array}{cc}0\hfill & \hfill 0\le v_t\le v_{in},v_t>v_{out}\\ P_w\cdot \Delta t\cdot \left(v_t-v_{in}\right)/\left(v_r-v_{in}\right)\hfill & \hfill v_{in}<v_t\le v_r\\ P_w\cdot \Delta t\hfill & \hfill v_r<v_t\le v_{out}\end{array}$$

(5)

The gas input of a CHP unit G_g,t has a restricted minimum and maximum capacity, G_g,min and G_g,max, respectively. The outputs of electricity production part E_g,t and heat production part H_g,t are constrained by its electrical efficiency η_ge,t and thermal efficiency η_gh,t [28]. In a partial load operation, η_ge,t may drop approximately linearly by 5%–7% at a minimal part load of 50%, which is considered to be the lower bound G_g,min here, while the total efficiency η_g can be assumed to remain constant:

$$\begin{array}{l}{G}_{g,t}=H_{g,t}/{\eta }_{gh,t}=E_{g,t}/{\eta }_{ge,t}=H_{g,t}\cdot {\eta }_{ge,t}/\left({\eta }_g-{\eta }_{ge,t}\right)\\ G_{g,min}=50\%\cdot G_{g,max}\le G_{g,t}\le G_{g,max}\end{array}$$

(6)

The electrical heat pump (EHP) unit is modeled as the air-water type with minimum and maximum thermal capacity H_h,min and H_h,max, respectively. The relationship between heat production H_h,t and electricity consumption E_h,t of an EHP unit is restrained by the coefficient of performance COP_h, which can be seen as constant during the analysis horizon:

$$\begin{array}{l}{E}_{h,t}=H_{h,t}/CO{P}_h\\ H_{h,min}\le H_{h,t}\le H_{h,max}\end{array}$$

(7)

Boilers with a coefficient of performance COP_l are used as supplementary heat sources because of the wide selectable range of capacity H_l,max complementary to primary heat sources (CHPs and EHPs), with the flexibility of starting up/shutting down quickly. The heat production H_l,t follows:

$$H_{l,t}=G_{l,t}\cdot {\eta }_l, 0\le H_{l,t}\le H_{l,max}$$

(8)

The battery energy storage system (BESS) is used as the electricity storage here. There are recommended lower and upper limits for the stored energy SOE_b,t of the BESS of SOE_b,min and SOE_b,max, respectively, to extend the battery’s lifetime. Coefficient η_b features the internal self-discharge loss during the last time period [t−1,t]. Coefficients η_b,c and η_b,d feature the charging and discharging efficiency, respectively, with the charging energy donated by E_c,t with a peak value of E_c,max and the discharging energy by E_d,t with E_d,max. SOE_b,t is related to E_c,t and E_d,t by:

$$\begin{array}{l}SO{E}_{b,min}\le SO{E}_{b,t}=SO{E}_{b,t-1}{\eta }_b+E_{c,t}{\eta }_{b,c}-E_{d,t}/{\eta }_{b,d}\le SO{E}_{b,max}\\ 0\le E_{c,t}\le E_{c,max}, 0\le E_{d,t}\le E_{d,max}\end{array}$$

(9)

A water tank is adopted as the heat storage here with an energy capacity C_s. There is an upper limit H_b,max for the imported heat storage H_c,t and exported heat storage H_d,t. Coefficients η_s, η_s,c, and η_s,d are the converted storage, imported, and exported efficiencies, respectively. The heat storage—related behavior and limits are described as:

$$\begin{array}{l}0\le H_{s,t}=H_{s,t-1}{\eta }_s+H_{s,c,t}{\eta }_{s,c}-H_{s,d,t}/{\eta }_{s,d}\le C_s\\ 0\le H_{c,t}\le H_{c,max}, 0\le H_{d,t}\le H_{d,max}\end{array}$$

(10)

Electricity demand is divided into two parts: (a) general electricity demand E_e,t; and (b) EV charging demand E_ev,t, respectively. The supply side of the electricity subsystem covers the RES output, CHP-produced electricity, scheduled ESS discharging, and purchased electricity from the external grid E_ex,t. These supplements cover the total electricity demand, as well as the scheduled ESS charging energy and EHP-consumed electricity. Similarly, the supply side of the heat subsystem, including the heat production of the CHP, EHP, boiler, and scheduled thermal storage charging, should meet with the heat demand H_q,t and discharging of the thermal storage. The co-generation and re-allocation between electricity and heat flows, such as CHPs and EHPs, bring about an interaction in both power flows. The equilibrium of production and consumption in the two energy subsystems must be satisfied simultaneously at every time step t:

$$\begin{array}{l}{E}_{w,t}+E_{g,t}+E_{ex,t}+E_{d,t}=E_{e,t}+E_{ev,t}+E_{h,t}+E_{c,t}\\ H_{g,t}+H_{h,t}+H_{l,t}+H_{d,t}=H_{q,t}+H_{c,t}\end{array}$$

(11)

Note that time is separated into discrete periods by constant time step Δt, and all variables are considered constant within every Δt. The following assumptions are also applied during the formulation of the objective function:

• The microgrid is generally a price-taker in the market since the upper size of a microgrid is generally limited to a few mega volt-amperes (MVAs);
• The optimization is based on economic performance, and the general assumption is that voltage and frequency stability control are not considered on this level;
• The connection between the microgrid and external grid is unconstrained, and thus the output of the installed DGs can be fully absorbed to increase the revenue of the microgrid;
• Signal delays and transmission losses are not considered.

This model structure can also be expanded to include other DGs, such as photo-thermal (PT) units, combined cooling heating and power (CCHP) units, electrical heater units, etc.

4.2. Implementation of the Energy Management Objective Function

The general objective function in Section 3 minimizing the total operational cost is transformed into a two-stage stochastic programming framework here, and implemented to the microgrid model above. A symbol $V_{t+k|t}^{p,r}$ stands for the corresponding variable V at the pth scenario of the predicted real-time electricity price ${\lambda }_{ex,t+k|t}^p$ and the rth scenario of the predicted wind power $E_{w,t+k|t}^r$ at time step t + k predicted at time step k. The first stage denotes the optimized working state and the gas consumption of the CHP based on the expectation of all scenarios including purchased energies from the external grid, charging schedules for postponable EVs, as well as the scheduling control of EHPs, boilers, and ESSs within their constraints in the second stage. The energy management process is formulated as:

$$\begin{array}{l}\begin{array}{cc}\min & {\displaystyle \sum _{k=0}^N\left[{\lambda }_f{G}_{g,t+k|t}{X}_{g,t+k|t}+F\left(X_{g,t+k|t},X_{g,t+k+1|t}\right)+{\displaystyle \sum _p^{{\mathrm{\Omega }}_p}{\pi }_{p,t+k|t}{\displaystyle \sum _w^{{\mathrm{\Omega }}_w}{\pi }_{w,t+k|t}}}\left({\lambda }_{ex,t+k|t}^p{E}_{ex,t+k|t}^{p,r}+{\lambda }_f{G}_{l,t+k|t}^{p,r}\right)\right]}\end{array}\\ \mathrm{Subject}\mathrm{to}:\\ \left(\mathrm{a}\right)\mathrm{Energy}\mathrm{flow}\mathrm{balance}\\ E_{w,t+k|t}^r+{\eta }_{ge,t+k|t}{G}_{g,t+k|t}{X}_{g,t+k|t}+E_{ex,t+k|t}^{p,r}+E_{d,t+k|t}^{p,r}=E_{e,t+k|t}+{\displaystyle \sum _{z=1}^{N_{z1}}{E}_{z,t+k|t}}+{\displaystyle \sum _{z=1}^{N_{z2}}{E}_{z,DSM,t+k|t}^{p,r}}+E_{h,t+k|t}^{p,r}+E_{c,t+k|t}^{p,r}\\ {\eta }_{gh,t+k|t}{G}_{g,t+k|t}{X}_{g,t+k|t}+H_{h,t+k|t}^{p,r}+H_{l,t+k|t}^{p,r}+H_{d,t+k|t}^{p,r}=H_{q,t+k|t}+H_{c,t+k|t}^{p,r}\\ \left(\mathrm{b}\right)\mathrm{DSM}\mathrm{of}\mathrm{postponable}\mathrm{EV}\mathrm{charging}\mathrm{requests}\\ {\displaystyle \sum _{t=t}^{t+N}{E}_{z,DSM,t}^{p,r}}=N_{c,z}{E}_z\forall z\in N_{z2},E_{z,DSM,t}^{p,r}\in \left\{0,E_z\right\}\\ \left(\mathrm{c}\right)\mathrm{Operational}\mathrm{contraints}\mathrm{of}\mathrm{components}\\ SO{E}_{b,min}\le SO{E}_{b,t+k|t}^{p,r}=SO{E}_{b,t+k-1|t}{\eta }_b+E_{c,t+k|t}^{p,r}{\eta }_{b,c}-E_{d,t+k|t}^{p,r}/{\eta }_{b,d}\le SO{E}_{b,max}\\ 0\le E_{c,t+k|t}^{p,r}\le E_{c,max},0\le E_{d,t+k|t}^{p,r}\le E_{d,max}\\ 0\le H_{s,t+k|t}^{p,r}=H_{s,t+k-1|t}{\eta }_s+H_{c,t+k|t}^{p,r}{\eta }_{s,c}-H_{d,t+k|t}^{p,r}/{\eta }_{s,d}\le C_s\\ 0\le H_{c,t+k|t}^{p,r}\le H_{c,max},0\le H_{d,t+k|t}^{p,r}\le H_{d,max}\\ E_{h,t+k|t}^{p,r}=H_{h,t+k|t}^{p,r}/CO{P}_h,H_{h,min}\le H_{h,t+k|t}^{p,r}\le H_{h,max}\\ H_{l,t+k|t}^{p,r}=G_{l,t+k|t}^{p,r}\cdot {\eta }_l,0\le H_{l,t+k|t}^{p,r}\le H_{l,max}\end{array}$$

(12)

where the fuel price λ_f, as the natural gas price here, is seen as a constant in the predicted horizon; the binary variable X_g,t+k|t is given as 1 to denote that the CHP unit is on at time t + k predicted at time k, and 0 otherwise; and a penalty function F(X_g,t+k-1|t, X_g,t+k|t) is used to restrict frequent changes of the on/off working state, reflecting the nontrivial start-up and shut-down cost of the CHP.

Given the nature that CHPs should not change their on/off states too frequently, the variables to be scheduled in the first stage pertain to the CHP unit that interacts with the scenarios of all the stochastic processes instead of being explicitly linked to one specific scenario. Such a structure not only reduces the computing space and variables by separating the working binary state, but is also perfectly suitable for incorporating with the BD technique.

4.3. Benders Decomposition

As discussed in Section 2.1, the Benders decomposition is one of the most widely used decomposition methods with proven efficacy. It improves the convergence time by separating the original optimization problem into a two-fold structure: a master problem for optimizing the first-stage problem bounded by the Benders cut, and a sub-problem which is impacted by the trial solution from the master problem. When applied to stochastic problems, the BD performs parallel computing on solving a set of sub-problems as smaller deterministic-equivalent problems simultaneously, which reduces the computational burden and accelerates the convergence time [29].

The master problem determines the trial solution for the first-stage decisions subjected to the approximate solutions of the sub-problems as a set of constraints added on each iteration, namely the Benders cuts. In this particular case, the first-stage decisions are the binary variable of the CHP working state in [t,t+N]:

$$\begin{array}{l}{Z}_{\mu }=\begin{array}{cc}\min & {\displaystyle \sum _{k=0}^N\left[{\lambda }_f{G}_{g,t+k|t}{X}_{g,t+k|t}+F\left(X_{g,t+k-1|t},X_{g,t+k|t}\right)+{\theta }_k\right]}\end{array}\\ {\theta }_{\mu }\ge {\displaystyle \sum _p^{{\mathrm{\Omega }}_p}{\pi }_{p,t+k|t}{\displaystyle \sum _w^{{\mathrm{\Omega }}_w}{\pi }_{w,t+k|t}}}\left[Q_{\mu ,p,w,t+k|t}-{\epsilon }_{\mu ,p,w}{G}_{g,t+k|t}\left(X_{\mu ,g,t+k|t}-{\overline{X}}_{\mu ,g,t+k|t}\right)\right]\end{array}$$

(13)

where θ_μ is the Benders cut in iteration μ describing the second-stage value based on the previously-determined first-stage decisions, Q_μ,p,w,t+k|t is the value of sub-problems at iteration μ and the pth, ωth scenario, ε_μ,p,w is the extreme direction or sensitivity for the corresponding Q_μ,p,w,t+k|t, X_μ,g,t+k|t is the trial solution determined in iteration μ. In the initial step (μ = 1), optimality cuts are not considered.

The trial solution X_μ,g,t+k|t in iteration μ is passed as an input parameter to sub-problems subject to constraints in Equation (12); then a solution for the second-stage decisions is determined. Each scenario will generate a different sub-problem objective function which shares a similar formulation but different input parameters. The sub-problems minimize the electricity purchase and boiler consumption with the trial CHP working action by:

$$\begin{array}{l}{Q}_{\mu ,p,w,t+k|t}=\begin{array}{cc}\min & \left({\lambda }_{ex,t+k|t}^p{E}_{ex,t+k|t}^{p,r}+{\lambda }_f{G}_{l,t+k|t}^{p,r}\right)\end{array}\\ s.t.\mathrm{Equation}\left(12\right)\left(\mathrm{a}\right)–\left(\mathrm{c}\right)\end{array}$$

(14)

A convergence criterion is used to stop the algorithm when the gaps between the master problem and sub-problems are within a pre-defined tolerance. While the value of the master problem is set up as a lower bound (LB), the upper bound (UB) of the optimal solution is represented by the trial solutions of the sub-problems. The convergence criterion is formulated as:

$$\begin{array}{l}\mathrm{LB}-\mathrm{UB}\le \delta \cdot \mathrm{LB}\\ \mathrm{UB}={\lambda }_f{G}_{\mu ,g,t+k|t}{\overline{X}}_{\mu ,g,t+k|t}+F\left({\overline{X}}_{\mu ,g,t+k-1|t},{\overline{X}}_{\mu ,g,t+k|t}\right)+{\displaystyle \sum _p^{{\mathrm{\Omega }}_p}{\pi }_{p,t+k|t}{\displaystyle \sum _w^{{\mathrm{\Omega }}_w}{\pi }_{w,t+k|t}}}{Q}_{\mu ,p,w,t+k|t}\\ \mathrm{LB}=Z_{\mu }\end{array}$$

(15)

where it is assumed that the difference δ between the bounds will be lower than or equal to 0.1%LB.

The solution procedure of the stochastic BD with the receding horizon policy applied to MPC is represented below (see Algorithm 1).

Algorithm 1. Stochastic BD

At each time step t do Initialize μ = 1, UB = +∞, LB = −∞, θ = 0 For k = 1 to N Set μ = 0, do   Solve master problem in Equation (13) and determine a trial solution Xμ,g,t+k|t   Update the value for LB   Solve all sub-problems in Equation (14) with trial solution Xμ,g,t+k|t   Update the value for UB if LB/UB > δ then   Add the new optimality cut θ associated with iteration μ to the master problem   Set μ = μ + 1, and repeat the solution procedure else   an optimal solution is achieved, and execute the control signals at time step t end for

5. Case Study

5.1. Numerical Specifications

The proposed approach is applied to a hybrid heat and power microgrid model of campus size with the numerical parameters of the components specified in

Table 2. Parameters of the modeled microgrid components (half hour-based).

Type	Rated Power	Other Parameters
CHP	Gg,max = 2.45 MW	ηg = 0.75~0.85, ηge = 0.35~0.40
EHP	Hh,max = 2.5 MW	COPh = 3.0
Boiler	Hl,max = 2.65 MW	ηl,t = 0.85
BESS	Ec,max = 2.4 MW, Ed,max = 2.4 MW	ηb,c = 0.95, ηb,d = 0.95, ηb = 0.99, SOEb,max = 4 MWh, SOEb,min = 0.8 MWh
Water tank	Ec,max = 2 MW, Ed,min = 2 MW	ηs,c = 0.90, ηb,d = 0.90, ηs = 0.95, Cs,max = 26 MWh, Cs,min = 12 MWh

.

A total analysis horizon of 24 h is adopted here. The distributions of scenarios with probabilities for the WTG output and the real-time electricity price are shown in Figure 5a. The heat and general electricity demand predictions are treated as a forecast with a standard deviation range ε_d of 25%, as shown in Figure 5b. The gas price is treated as a constant value of $49.2 USD/MWh.

Ten EVSEs with a same charging power P_z of 7.2 kW are installed in the microgrid. For simplification, the total needed charging energy ΔSOE for each EV is calculated in advance as the difference between a desired SOE_z,set and the initial SOE when arriving. To display the EV charging demands intuitively, the actual EV admissions of each EVSE including total charged energy ΔSOE, arrival time T_z,arr, and departure time T_z,set are rounded up and are organized in the form of {ΔSOE, T_z,arr–T_z,set} where the unit (kWh) is omitted for simplicity. For instance, an array of {7.2, 0:00–6:00} signifies that this EV arrives at 12 a.m. and leaves at 6 a.m. with a total needed charging energy of 7.2 kWh. Every EVSE can be plugged with no more than one EV at the same time. EV charging demands for the simulation are listed in

Table 3. Actual EV charging energy, arrival, and departure settings.

EVSE No.	EV Charging Request Data
1	{7.2, 0:00–6:00}, {14.4, 9:00–15:30}, {10.8, 16:00–18:00}, {7.2, 21:30–23:30}
2	{10.8, 7:30–10:30}, {10.8, 11:30–15:00}, {14.4, 16:30–20:30}, {7.2, 22:00–24:00}
3	{10.8, 0:00–6:00}, {10.8, 7:30–10:00}, {14.4, 11:00–15:00}, {10.8, 17:00–19:00}
4	{10.8, 1:30–7:00}, {18.0, 8:00–11:30}, {10.8, 12:30–15:30}, {7.2, 17:00–19:30}
5	{3.6, 1:30–6:00}, {7.2, 6:30–8:00}, {10.8, 8:30–11:30}, {7.2, 12:30–15:30},{7.2, 21:00–24:00}
6	{10.8, 9:00–11:30}, {14.4, 13:00–16:30}, {7.2, 17:00–18:00}, {10.8, 20:30–23:00}
7	{7.2, 0:00–5:30}, {10.8, 8:00–11:00}, {10.8, 18:00–21:30}, {3.6, 22:30–24:00}
8	{10.8, 1:00–8:00}, {10.8, 9:30–13:00}, {3.6, 16:00–17:00}, {14.4, 19:00–22:00}
9	{3.6, 0:30–7:00}, {10.8, 8:00–10:00}, {10.8, 11:00–13:30}, {10.8, 16:30–20:30}
10	{10.8, 8:00–11:30}, {7.2, 14:30–16:30}, {10.8, 21:00–24:00}

. Although we regard these as non-perfectly-matched stochastic models for EV charging behavior in microgrids at this stage, the distribution sampled in Table 3 fits some accepted statistical information about travel patterns of vehicle users, such as in [30,31], other than randomized samples.

All numerical simulations are solved in MATLAB on an Intel i3 3.2 GHz Windows-based PC with 4 GB of RAM. The proposed approach is compared with some other approaches on the same model, although detailed outputs may differ significantly due to different settings of a microgrid.

5.2. Simulation and Discussion

5.2.1. Dealing with Robust and Stochastic Uncertainties

In order to exemplify the economic performance of the proposed scenario-based MPC approach, it is compared against DMPC, where these scenarios of stochastic variables are replaced by only one scenario in their predictions, and precise MPC (PMPC), where perfect deterministic foresight on the predicted inputs is assumed. Due to its ideality, PMPC is taken as the benchmark example for forecast accuracy with no uncertainties. For all three MPC-based approaches, the analysis horizon is conducted for each time period interval at 0.5 h, and the prediction horizon N = 16.

Figure 6 shows the results of these three approaches. PMPC with perfect future information attains a prominent behavior of overall expenditure, as expected. The accumulated operating cost for the proposed scenario-based MPC is $10,898.07 USD/day, which is the closest to PMPC, and saves the cost of 0.955% in the analysis horizon compared with DMPC. As to the computing time, no major difference is detected between DMPC and PMPC since both of them equivalently contain only one scenario. However, there is an observable deceleration of the convergence speeds for scenario-based MPC compared with the others due to the involvement of many more scenarios than the traditional one-scenario procedures. It is clear that when prediction uncertainty is inevitable in general, scenario-based MPC outperforms in reducing the total operating cost compared to DMPC at the expense of the computing time.

To further illustrate the DSM process of EV charging, Figure 7 shows the charging behavior schedule with accumulated SOE and the charging cost of each EV. Instead of executing a straightforward strategy as immediate and uncoordinated charging, postponements of EV charging are arranged by energy management to coordinate with the operational cost optimization of the microgrid as shown in Figure 7a. It can be seen in Figure 7b that each EV copes with a flat or declining charging cost after DSM, where the EV owners can also benefit from participating in the energy management. Since the objective function focuses on minimizing the total operational costs, the postponements often occur in the time steps with higher electricity prices or with lower local energy outputs in the prediction horizon, resulting in the reduced charging expenditure for each EV. The DSM cannot be achieved without the worst-case scenario for the unknown charging demand prediction while the robustness of the worst-case scenario ensures that the requested EV charging demand is fulfilled when an EV is leaving.

Figure 8 is supplemented to show the supply-demand balance of the electrical and heat flow by scenario-based MPC. Useful future information may affect the current decision more accurately at each time step, and the expected optimization can be scheduled more precisely by taking account of stochastic scenarios in advance, i.e., before peaks of the real-time electricity price and the demands.

5.2.2. Dealing with Increasing Scenario Sizes

As concluded, the proposed scenario-based MPC approach gains a lower operational cost but consumes more computing time. With the trend of progressive RES penetrations, this study case is simulated to investigate the importance of considering the computing performance while maintaining a good economic performance of scenario-based MPC. Assume that two WTGs with the same parameters are installed in the microgrid above, sharing dissimilar meteorological environments which lead to diverse wind speed scenarios. Figure 9 depicts the total operational cost and the peak calculating time for one time step in the whole simulated analysis horizon (N = 16). When uncertainties in the microgrid increase, the gap of the total operational cost between scenario-based MPC and DMPC also increases; in this particular case study, the cost savings rises from 0.955% to 1.249%. However, there is a significant drawback of the computing time observed, where the peak calculating time increase rises from 0.322 s to 2.077 s. A similar situation happens if more units of RES are integrated, since the problem dimension increases inevitably on account of more scenarios in the scenario combination.

The results of scenario-based MPC with BD are shown in

Table 4. Performances under different settings.

Setting Values	Total 1 Operational Cost (USD)			Peak Calculating Time for One Time Step 1 (s)
Installed WTG number	Scenario-based MPC		DMPC	Scenario-based MPC		DMPC
	without BD	with BD		without BD	with BD
1	10,898.07	10,902.76	11,003.12	0.404	0.095	0.082
2	8944.95	8957.33	9056.64	2.367	0.331	0.290

¹ In a simulated analysis horizon of 24 h.

to verify whether efficient solutions can achieved in a reasonable computational time. As expected, with the help of BD, the computing times of scenario-based MPC accelerate significantly, and the gap of the computing times between them is acceptable, which makes the proposed approach quite competitive. Meanwhile, the total operational cost is quite similar to that of the original scenario-based MPC approach without parallel computing. By using BD, energy management is capable of providing optimized solutions provided scenario-based MPC while ensuring a competitive computing performance with DMPC.

In addition to the scalability of integrated RES, different prediction models with various scenario generation and reduction methods can also result in increased combined scenario numbers, which lead to high volumes of variables. The adoption of BD makes the energy management systems able to choose applicable prediction models and suitable prediction horizons with more flexibility.

6. Conclusions

This paper presents a scenario-based MPC approach which enables the energy management system of a microgrid to improve the economic performance considering various characteristics of uncertainties. As for uncertainties of a robust nature, here, the EV charging loads at the present stage, a worst-case scenario, as well as the DSM are presented coping with the receding prediction horizon; as for the uncertainties of a stochastic nature, such as RES outputs and the real-time electricity price, the probabilistic scenarios are introduced while the computational efficiency of the objective function is improved by applying the BD technique. The numerical simulations indicate that the proposed approach offers cumulative cost savings while satisfying consumers' demands under a computationally tractable process. Moreover, the proposed approach can be tailored to microgrid applications with more DGs due to the scalability and flexibility to promote full-scale commercial microgrids. Future work could also be extended by incorporating the approach with available EV charging load prediction models when the penetration rate of EVs in a microgrid is adequate enough.