Optimal Day-Ahead Scheduling of Microgrids with Battery Energy Storage System

Abstract

Optimal scheduling is a requirement for microgrids to participate in current and future energy markets. Although the number of research articles on this subject is on the rise, there is a shortage of papers containing detailed mathematical modeling of the distributed energy resources available in a microgrid. To address this gap, this paper presents in detail how to mathematically model resources such as battery energy storage systems, solar generation systems, directly controllable loads, load shedding, scheduled intentional islanding, and generation curtailment in the microgrid optimal scheduling problem. The proposed modeling also includes a methodology to determine the availability cost of battery and solar systems assets. Simulations were carried out considering energy prices from an actual time-of-use tariff, costs based on real market data, and scenarios with scheduled islanding. Simulation results provide support to validate the proposed model. Data illustrate how energy arbitrage can reduce microgrid costs in a time-of-use tariff. Results also show how the microgrid’s self-sufficiency and the storage system’s capacity can impact the microgrid’s energy bill. The findings also bring out the need to consider the scheduled islanding event in the day-ahead optimization for microgrids.

1. Introduction

The introduction of distributed energy resources (DER) in the distribution system has brought several benefits to consumers, utilities, and society. These include increasing system reliability and resilience, reducing greenhouse gases, relieving distribution and transmission systems, reducing energy tariffs for consumers who also produce energy (referred to as prosumers in [1]), and the possibility of forming an energy market in which consumers and prosumers can participate actively. However, it also came with some challenges as the need to readapt distribution systems (infrastructure, automation, protection, control, operating systems, planning) to receive these new resources and to deploy modern energy tariffs to meet the equity criterion [1], finding a trade-off between encouraging renewable sources while avoiding cross-subsidization between consumers and prosumers. Microgrids (MG), as part of this new active distribution system, experience the same benefits and challenges; however, they have an additional purpose: being economically efficient, producing energy at the lowest possible cost while eliminating or minimizing waste.

A microgrid can contain loads and distributed energy resources such as distributed generators, storage devices, and controllable loads, and must be operated in a controlled way whether or not connected to the main grid [2]. With such resources, a microgrid can provide ancillary services to the distribution system operator (DSO), perform energy arbitrage (store energy when the price is low to sell it or use it when the price is high), actively participate in energy markets, in addition to participating in distribution systems with multiple microgrids. However, to make these activities feasible, a microgrid must optimize the use of its resources, seeking economic efficiency. Due to the renewable energy sources and load forecasting features, the time horizon for this optimization is the day-ahead. On the other hand, markets whose energy prices can vary hourly may require optimization with an equivalent time horizon. Thus, day-ahead and intra-day optimal schedulings are requirements for microgrids to participate in current and future energy markets.

The optimal scheduling of microgrids with battery energy storage system (BESS), solar and/or wind generation has been studied in [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. Although these works address the modeling of solar photovoltaic systems for microgrids, none of them discusses curtailment modeling in day-ahead scheduling. In addition, the costs related to the PV system can be modeled in the objective function by a fixed term, or by a term proportional to the power generated, or by a combination of them, depending on the methodology adopted to address such costs. The authors [7,14,17] model these costs through a term proportional to the power generated. In [18], these costs are modeled using a fixed term and another proportional to the power. Similar modeling is presented in [3] with a fixed term and two others which are proportional to the power and the square of the power. However, all of these works do not present a methodology for calculating costs related to the PV system. As a consequence, in the simulations presented in those works, the authors consider such costs as known values.

Regarding the modeling of battery energy storage systems, most authors present simplified modeling that considers only upper and lower limits for charging and discharging power and for the state of charge (SOC), as in [3,5,7,13,14]. Others present a recursive equation for the BESS state of charge but without clarifying how to transform such a recursive equation into a constraint of an optimization problem, as in [4,5,6,10,16,20]. The authors in [19] address the modeling of BESS costs considering the ratio between the cost of capital and the total number of cycles, however without considering the impacts of the depth of discharge and battery degradation. Finally, the authors in [18] address the modeling of BESS costs considering the cost of battery degradation, for a given depth of discharge, an exponential relationship that depends on a function fitting to the curve provided by the manufacturer. However, they do not address the method for obtaining the curve parameters.

Demand response programs in microgrids have been investigated through modeling of directly controlled loads (DCL) in [14,15,16,21]. Authors in [14,21] use a discrete DCL model with load blocks to analyze the problem of microgrids optimization with demand response, which is equivalent to a model of interruptible loads per block. Authors in [16] present a similar approach but considering the load to be interrupted as a continuous variable. However, these works do not address the modeling of shiftable loads as an option for DCLs. In [15], loads are classified as fixed and adjustable (DCL), including shiftable and curtailable (interruptible) loads. Both types of DCL are modeled only by an upper power limit. However, shiftable loads require more complete modeling due to their need for a continuous duty cycle, which was not addressed at that work.

Several studies address the problem of optimal day-ahead scheduling of microgrids by modeling uncertainties in load and generation forecasts. Authors in [6,10,16,17] divide the optimization problem into levels or stages to deal with uncertainties in load, generation, and price forecast. Thus, an intermediate stage of optimization has the role of finding, in a set of bounded uncertainties, the condition that maximizes the cost of the microgrid (the worst case). Once such a condition has been defined, the next stage seeks to minimize the cost of the microgrid. Authors refer to this technique as robust optimization. In [18], the authors model the uncertainties probabilistically, and in [19], the authors propose the interval mathematics technique to solve the same problem. However, considering uncertainties in the optimization problem comes at a price because, in addition to increasing the complexity of the problem, it leads to a sub-optimal solution compared to the deterministic solution. It is what the authors in [22,23] call the price of robustness in robust optimization. If on the one hand, there may be a mismatch between the forecasted and realized values as a deterministic solution is used, on the other hand, such a mismatch may even increase in solutions such as robust optimization. That is because in such a technique there is no correlation between the uncertainty considered in the day-ahead problem and that which should actually happen in the intra-day. These are the main reasons why the present study addresses this problem in a deterministic way, which has the advantage of simplicity and low implementation complexity.

Linear programming is the most used technique among the works consulted in the present literature review to solve the problem of day-ahead scheduling optimization for microgrids [6,10,13,15,16,21]. Furthermore, in an extensive literature review carried out in [24], it was found that mixed-integer linear programming (MILP) is the most appropriate approach for solving the scheduling problem for microgrids and virtual power plants. These works support the choice for linear programming in the present study, without prejudice to the mathematical modeling of the elements here presented.

The optimal dispatch of distributed energy resources in distribution systems and microgrids is of economic interest to the entities that operate these systems. Although there are some works published in the literature on the subject, there also are some gaps in the mathematical modeling that this work seeks to fill. For example, no previous work was found that formulated the day-ahead problem considering scheduled intentional islanding, as defined by the IEEE Std. 1547-2018 [25]. The same for the modeling of shiftable loads of continuous cycle and solar photovoltaic (PV) generation curtailment. Furthermore, there is no work in the present literature review that model in detail the BESS system and its cost of availability, maintaining the linearity of the problem.

Thus, this work presents detailed modeling of the optimal dispatch of microgrids with the distributed generation resources namely: a battery energy storage system; a PV system with the possibility of curtailment; demand response (DR) through directly controllable loads, such as the shiftable and interruptible ones; and the option of considering scheduled intentional islanding events with the possibility of performing load shedding (LS) when necessary.

The battery storage system is modeled according to the IEEE Std 2030.2.1-2019. Details of this modeling are presented, considering aspects such as the efficiency of the whole BESS, and state of charge as a non-recursive constraint. Furthermore, the ranges of some practical values for the model parameters are presented. In cost modeling, an exponential equation for the battery state of health (SOH) is presented in order to fit a possible BESS manufacturing data curve. The costs of BESS storage, charging, and discharging are then presented as a function of the SOH threshold, rated depth of discharge (DOD), cycle life, and the battery capacity.

In the cost modeling of the PV system, factors such as annual degradation of the photovoltaic panels, the lifespan of PV panels, and the solar generation capacity from the region of installation are considered. The methodology adopted in the modeling gave rise to a problem of mixed-integer linear programming, which was implemented in the MATLAB software. Simulation results seek to validate the mathematical model against all distributed energy resources of the microgrid. In addition, a case study based on the time-of-use tariff from the Brazilian energy market is presented considering actual tariff rates. Such a study examines the conditions for energy arbitrage performed by BESS in a microgrid.

The main contributions of this paper are:

• Presenting a practical methodology to compute the BESS availability cost in USD/kWh; modeling of BESS system considering battery, converter, and transformer useful values of efficiencies; modeling the BESS state of charge as a non-recursive constraint to the MILP problem;

• Modeling of scheduled intentional islanding in a microgrid with the possibility of curtailment and load shedding;

• Modeling of non-interruptible shiftable loads (continuous cycle) as washers and driers when participating in demand response programs;

• Presenting a practical methodology to compute the availability cost of PV systems.

The rest of this paper is organized as follows. Section 2 presents an overview of an energy management system (EMS) for microgrids considering the day-ahead and real-time modules. Section 3 presents the mathematical model of a microgrid regarding its dynamics of operation and costs. The simulation methodology is presented in Section 4, and simulation results in Section 5. A discussion of the main findings is presented in Section 6, and the main conclusions of the work in Section 7.

2. Microgrid Energy Management System

From a logical point of view, a microgrid EMS can be divided into day-ahead and real-time modules, as illustrated in Figure 1. Furthermore, the day-ahead EMS requires some input modules to update their data daily. The forecasting data input module is responsible for the daily forecast of the microgrid generation and load curves, which can be based on historical data and weather forecasting. The DSO inputs represent information passed on by the distribution system operator, as scheduled outages and distribution network operational limits, which can be established by a power flow analysis of the distribution network region in which the microgrid is located. The DR contracts data module has the role of updating technical information about consumers who participate in demand response programs. Finally, the energy price module is responsible for updating the energy market prices for the next day.

The day-ahead EMS module must contain the microgrid settings data, the mathematical model of all the microgrid devices and the desired operating conditions, and an optimizer algorithm. The decision-maker block calculates costs that depend on settings data and sets up the optimization problem with constraints and parameters for the optimizer. The output of this module is the optimal solution, which is formed by the states and powers that will be the setup of operating for the real-time EMS module throughout the next day. Therefore, the day-ahead EMS module, as illustrated in Figure 1, defines and limits the scope of this work.

In addition to daily receiving the optimal operating condition from the day-ahead module, the real-time EMS module is updated with data from local and central controllers, electronic meters, and the battery management system (BMS) in real-time. This module has the role of a real-time decision-maker as it must correct the deviations between the forecasted load and generation and the real-time acquisition data. Such a correction can be made, for example, considering a one-hour horizon. These data are then passed on to lower-level controllers. In addition, this module has the role of deciding the microgrid’s operating strategy if the schedule received from the day-ahead module is broken, as can happen in the event of intentional unscheduled islanding [25]. Finally, this module can daily pass data estimated by the BMS to the day-ahead module, such as the current BESS capacity and state of health, as well as the state of charge expected to the end of the day. It will be the initial state of charge in the day-ahead module. Real-time EMS modeling is outside the scope of this work.

3. Modeling Methodology

3.1. Microgrid Modeling

Figure 2 illustrates a schematic diagram of the microgrid reference model for the MILP problem addressed in the present study. The model comprises integer (binary) and real decision variables. The former ($s{t}^{\left(t\right)}$ type) has the function of virtual disconnect switches, and the latter ($P^{\left(t\right)}$ type) represents the active powers of each element of the model. Variables of type $\widehat{P}$ represent forecasted values, which are input variables for the model. Features of the model include: solar PV system, battery energy storage system, directly controllable loads, conventional loads, load shedding, energy transactions with the external grid, and the possibility of islanding. As this model requires a large number of decision variables, they are described in Table 1 and Table 2. Table 3 presents some input variables of the model in Figure 2.

Some elements in the model do not necessarily correspond to a real device present in the microgrid. The energy transacted with the external grid pass through virtual disconnectors (in practice, they are not required). Creating these branches connected to the exchange bus allows the model to identify the agents that are exchanging energy with the microgrid, in addition to quantifying that energy. In fact, such transaction can be a contractual obligation only, since all the energy that enters and leaves the microgrid through the point of common coupling (PCC) bus has a single source or destination: the external grid. Furthermore, in practical applications, the BESS is connected to the main bus through a single branch, and a battery management system is responsible for controlling the charge and discharge operations. On the other hand, without loss of generality, modeling charging and discharging by different branches (see Figure 2) may be more suitable for mathematical programming purposes.

In the MG model of Figure 2, DCL are loads in a demand response program that the MG operator can remotely disconnect or cycle upon contractual agreement and notification [26]. Furthermore called flexible loads, they can be interruptible or shiftable. The former are loads subject to curtailment or interruption under system contingencies, as air conditioners and water heaters [27]. The latter are continuous cycle loads that can be shifted in time as washers and dryers. On the other hand, when the microgrid operator needs to disconnect non-flexible loads, the conventional ones in Figure 2, due to a contingency or generation shortage, it must use the load shedding resource. As this type of load is not participating in a demand response program, there is a penalty associated with such a disconnection.

As a day-ahead scheduling is made up of a 24-hours time horizon, then the time slot for the model can be defined as,

(1)$\Delta t={\displaystyle \frac{24}{N}}\left(\mathrm{h}\right),\mathrm{for}N\in \mathbb{Z},N\ge 1$

where N represents the number of time slots. In addition, throughout this work, the independent variable $t\in {\mathbb{S}}_N$ represents the discrete-time of the model, where ${\mathbb{S}}_N=\{k\in \mathbb{Z}\mid 1\le k\le N\}$

The rate of change of the physical quantities handled by the model at each time slot should guide the choice of $\Delta t$. It should be small enough to capture significant changes in the forecasted load and solar generation curves, for example. On the other hand, it must be large enough to the meet the time requirements of microgrid’s tertiary control. Typical values for $\Delta t$ are from $1/4$ h [9,10,12] to 1 h [6,16,28]. However, 1/4 h may be a more appropriate time slot when the problem involves controllable loads, as it can better adjust to load cycles.

3.1.1. Objective Function

Consider a daily MG energy bill ($E_{bill}$) made up of the difference between daily costs and revenue, such that positive $E_{bill}$ values indicate an amount to be paid, and negative, an amount to be received or converted into energy credits, i.e.,

(2)$E_{bill}:=-\underset{revenue}{\underbrace{\Delta t\sum _{t=1}^N{P}_{sel}^{\left(t\right)}{p}_{sel}^{\left(t\right)}}}+\underset{costs}{\underbrace{\Delta t\sum _{t=1}^N{P}_{pur}^{\left(t\right)}{p}_{pur}^{\left(t\right)}+\left\{\mathrm{MG}\mathrm{costs}\right\}}}\left(\$\right)$

where $p_{pur}$ and $p_{sel}$ are the purchase and sale prices of energy, respectively. Although it is usual in other studies to name $E_{bill}$ a daily operating cost, Equation (2) shows its value does not always represent a cost. Therefore, optimizing the dispatch of a microgrid consists of minimizing its daily energy bill; and, the optimization problem can be formulated as

(3)$\begin{array}{cc}\hfill \underset{\begin{array}{c}\mathbf{x}\end{array}}{min}& E_{bill}\left(\mathbf{x}\right)=C_{pv}+\sum _{t=1}^N({st}_{{pur}_g}^{\left(t\right)}{c}_{{pur}_g}+{st}_{{sel}_g}^{\left(t\right)}{c}_{{sel}_g}+\sum _{j=1}^M\left({st}_{{pur}_{\mathsf{\mu }}}^{(j,t)}{c}_{{pur}_{\mathsf{\mu }}}+{st}_{{sel}_{\mathsf{\mu }}}^{(j,t)}{c}_{{sel}_{\mathsf{\mu }}}\right)+\hfill \\ \hfill & {st}_{shd}^{\left(t\right)}{c}_{shd}+{st}_{int}^{\left(t\right)}{c}_{int}+{st}_{shf}^{\left(t\right)}{c}_{shf}+{st}_{chr}^{\left(t\right)}{c}_{chr}+{st}_{dch}^{\left(t\right)}{c}_{dch}+\Delta t(P_{pvc}^{\left(t\right)}{p}_{pvc}^{\left(t\right)}+P_{chr}^{\left(t\right)}{p}_{chr}+P_{dch}^{\left(t\right)}{p}_{dch}+\hfill \\ \hfill & P_{shd}^{\left(t\right)}{p}_{shd}^{\left(t\right)}+P_{int}^{\left(t\right)}{p}_{int}^{\left(t\right)}+P_{shf}^{\left(t\right)}{p}_{shf}^{\left(t\right)}+P_{{pur}_g}^{\left(t\right)}{p}_{{pur}_g}^{\left(t\right)}-P_{{sel}_g}^{\left(t\right)}{p}_{{sel}_g}^{\left(t\right)}+\sum _{j=1}^M\left(P_{{pur}_{\mathsf{\mu }}}^{(j,t)}{p}_{{pur}_{\mathsf{\mu }}}^{(j,t)}-P_{{sel}_{\mathsf{\mu }}}^{(j,t)}{p}_{{sel}_{\mathsf{\mu }}}^{(j,t)}\right)\left)\right)\hfill \end{array}$

where $\mathbf{x}$, the vector of decision variables, is a feasible solution to the problem; and M is the number of external MGs with whom the microgrid carries out energy transactions.

The linear objective function in Equation (3) is a sum of products. Those of type $\{st\times c\}$, whose variables are shown in Table 1 and Table 4, represent fixed costs of energy transaction and operating DERs such as BESS, load shedding, and controllable loads. Products of type $\Delta t\{P\times p\}$, whose variables are described in Table 2 and Table 5, represent the revenues and expenses with the amount of energy transacted and the costs to operate each kWh of the DERs.

3.1.2. Energy and Power Exchanges

Regarding the energy transaction, there must be a logical constraint to deny that the microgrid can buy and sell energy simultaneously, that is,

(4)$0\le {st}_{{pur}_g}^{\left(t\right)}+{st}_{{sel}_g}^{\left(t\right)}+\sum _{j=1}^M{st}_{{pur}_{\mathsf{\mu }}}^{(j,t)}+{st}_{{sel}_{\mathsf{\mu }}}^{(j,t)}\le {st}_{con}^{\left(t\right)},\forall t\in {\mathbb{S}}_N$

However, if energy market rules allow the purchase of energy from two or more entities simultaneously, then this constraint can be decomposed into constraints that consider each term individually. The same goes for a market logic of simultaneous selling.

According to Equation (4), to insert a condition of islanding into the schedule, it is enough to make ${st}_{con}=0$ during the desired period of islanding, $T_{isl}\in {\mathbb{S}}_N$, and ${st}_{con}=1$, otherwise. In  fact, ${st}_{con}=0$ means that the microgrid can be either indeed operating in island mode (the real disconnect switch is open, e.g., due to a scheduled maintenance outage planned by the DSO) or connected but without exchanging power with the grid (see Equation (5)). Furthermore, if the microgrid cannot sell energy to the main grid (e.g., due to contractual restrictions with the DSO), then it is enough to make, ${st}_{{sel}_g}^{\left(t\right)}={st}_{{sel}_{\mathsf{\mu }}}^{(j,t)}=0,\forall t\in {\mathbb{S}}_N,\forall j\in {\mathbb{S}}_M$, where ${\mathbb{S}}_M=\{k\in \mathbb{Z}\mid 1\le k\le M\}$ is the set of external microgrids.

Regarding the power exchange, some technical features as line capacity, PCC sizing, or contracted limits can restrain the power exchange ($P_{pur}$ and $P_{sel}$) between the microgrid and the main grid or external microgrids [6]. Furthermore, if a state variable ${st}_{pur}$ or ${st}_{sel}$ is zero, then the respective power $P_{pur}$ or $P_{sel}$ must be zero, i.e.,

(5a)$\begin{array}{c}\hfill \begin{array}{c}0\le P_{{pur}_g}^{\left(t\right)}\le {st}_{{pur}_g}^{\left(t\right)}{\overline{P}}_{{pur}_g},\forall t\in {\mathbb{S}}_N\hfill \end{array}\end{array}$

(5b)$\begin{array}{c}\hfill \begin{array}{c}0\le P_{{sel}_g}^{\left(t\right)}\le {st}_{{sel}_g}^{\left(t\right)}{\overline{P}}_{{sel}_g},\forall t\in {\mathbb{S}}_N\hfill \end{array}\end{array}$

(5c)$\begin{array}{c}\hfill \begin{array}{c}0\le P_{{pur}_{\mathsf{\mu }}}^{(1,t)}\le {st}_{{pur}_{\mathsf{\mu }}}^{(1,t)}{\overline{P}}_{{pur}_{\mathsf{\mu }}}^{\left(1\right)},\forall t\in {\mathbb{S}}_N\hfill \\ ⋮⋮⋮\hfill \\ 0\le P_{{pur}_{\mathsf{\mu }}}^{(M,t)}\le {st}_{{pur}_{\mathsf{\mu }}}^{(M,t)}{\overline{P}}_{{pur}_{\mathsf{\mu }}}^{\left(M\right)},\forall t\in {\mathbb{S}}_N\hfill \end{array}\end{array}$

(5d)$\begin{array}{c}\hfill \begin{array}{c}0\le P_{{sel}_{\mathsf{\mu }}}^{(1,t)}\le {st}_{{sel}_{\mathsf{\mu }}}^{(1,t)}{\overline{P}}_{{sel}_{\mathsf{\mu }}}^{\left(1\right)},\forall t\in {\mathbb{S}}_N\hfill \\ ⋮⋮⋮\hfill \\ 0\le P_{{sel}_{\mathsf{\mu }}}^{(M,t)}\le {st}_{{sel}_{\mathsf{\mu }}}^{(M,t)}{\overline{P}}_{{sel}_{\mathsf{\mu }}}^{\left(M\right)},\forall t\in {\mathbb{S}}_N\hfill \end{array}\end{array}$

where ${\overline{P}}_{pur}$ and ${\overline{P}}_{sel}$ are individual upper limits imposed by external entities (index $\mathsf{\mu }$ for MGs, and g for the main grid). Besides, a power flow analysis considering the MG exchange bus as a node in the external grid can result in a constraint of physical upper limit ${\overline{P}}_{pcc}$ for the point of common coupling,

(6a)$\begin{array}{c}\hfill \begin{array}{c}\hfill P_{{pur}_g}^{\left(t\right)}+\sum _{j=1}^M{P}_{{pur}_{\mathsf{\mu }}}^{(j,t)}\le {\overline{P}}_{pcc},\forall t\in {\mathbb{S}}_N\end{array}\end{array}$

(6b)$\begin{array}{c}\hfill \begin{array}{c}\hfill P_{{sel}_g}^{\left(t\right)}+\sum _{j=1}^M{P}_{{sel}_{\mathsf{\mu }}}^{(j,t)}\le {\overline{P}}_{pcc},\forall t\in {\mathbb{S}}_N\end{array}\end{array}$

3.1.3. Power Balance

To satisfy the power balance at a given time t, the sum of the energy generated, imported and exported must equal the microgrid demand. Assuming the total load power is given by (see Figure 2)

(7)$P_{load}^{\left(t\right)}={\widehat{P}}_{load}^{\left(t\right)}-P_{shd}^{\left(t\right)}-P_{int}^{\left(t\right)}+P_{shf}^{\left(t\right)},\forall t\in {\mathbb{S}}_N$

then, the power balance equation is given by

(8)$\begin{array}{c}\hfill P_{int}^{\left(t\right)}+P_{shd}^{\left(t\right)}+P_{dch}^{\left(t\right)}+P_{{pur}_g}^{\left(t\right)}+\sum _{j=1}^M{P}_{{pur}_{\mathsf{\mu }}}^{(j,t)}-P_{{sel}_{\mathsf{\mu }}}^{(j,t)}-P_{chr}^{\left(t\right)}-P_{{sel}_g}^{\left(t\right)}-P_{pvc}^{\left(t\right)}-P_{shf}^{\left(t\right)}={\widehat{P}}_{load}^{\left(t\right)}-{\widehat{P}}_{pv}^{\left(t\right)}s{t}_{pv},\\ \hfill \forall t\in {\mathbb{S}}_N\end{array}$

where ${\widehat{P}}_{load}^{\left(t\right)}$ and ${\widehat{P}}_{pv}^{\left(t\right)}$ are the forecasted load and PV output curves in kW.

3.1.4. Battery Energy Storage System

A battery energy storage system for power system applications can comprise a coupling transformer, a converter system, and a battery pack. The efficiency of each one (${\eta }_{trf}$, ${\eta }_{conv}$, and ${\eta }_{bat}$, respectively) impacts the overall BESS efficiency as

(9)${\eta }_{bess}={\eta }_{trf}{\eta }_{conv}{\eta }_{bat}$

which is valid for both charging and discharging operations, once in practical applications, it is usual to consider the charge and discharge efficiency (${\eta }_{bat}$) are the same. In this case, ${\eta }_{bat}=\sqrt{{\eta }_{rte}}$, where ${\eta }_{rte}$ is the round trip efficiency. The storage efficiency (${\eta }_{sto}$) could also make up the value of ${\eta }_{bat}$; however, batteries present high ${\eta }_{sto}$. The lead-acid technology, for example, can achieve a self-discharge of 2–5% per month, and lithium-ion 1% per month [29]. Therefore, the storage efficiency can be neglected in models with a time horizon of a day.

According to the authors in [30], higher technologies lithium-ion batteries can achieve an ${\eta }_{rte}$ of 95%, which implies an ${\eta }_{bat}$ of $0.9747$. On the other hand, authors in [31] stand lead-acid and lithium-ion technologies can achieve a round trip efficiency of 85% and 90%, respectively, which leads to the values of $0.9220$ and $0.9487$ for ${\eta }_{bat}$. In addition, according to data in [32], a battery manufactured with the LiFePO${}_4$ technology can achieve an ${\eta }_{bat}$ of 0.9747. Regarding the efficiency of converters and transformers, manufacturing data [33] for systems from 500 kW to 5000 kW show that each can be greater than 97%. Table 6 presents a summary of these efficiency values for three battery technologies, and the respective ${\eta }_{bess}$ calculated from Equation (9).

The IEEE Std 2030.2.1™ [34] presents the main concepts related to BESS engineering. Some of them are illustrated in Figure 3. According to the standard, state of charge is “the degree to which a BESS is charged relative to the maximum possible amount of energy that can be stored by the system”, and depth of discharge is “the degree to which a BESS is discharged relative to the maximum possible amount of energy that can be discharged by the system”. Thus, let $E_{avl}$ be the BESS available capacity (the maximum possible amount of energy..., kWh), then

(10)$SOC={\displaystyle \frac{E_{sto}}{E_{avl}}}$

and $DOD=1-SOC$; where $E_{sto}$ is the stored energy.

Furthermore, the ratio of available capacity to rated capacity ($E_r$) defines the battery state of health,

(11)$SOH={\displaystyle \frac{E_{avl}}{E_r}}$

Therefore, the amount of stored energy during the BESS operation is $E_{sto}=SOCSOH{E}_r$. In real life systems, a battery management system can estimate the battery state of health at run time.

As the state of charge is a constraint in the optimization problem, Equation (10) needs to be rewritten to illustrate the dynamics of charging and discharging the battery. Adopting the main bus as the reference point, the next state of charge can be computed from the previous one as

(12)$\begin{array}{c}\hfill \begin{array}{cc}\hfill SO{C}^{\left(t\right)}{E}_{avl}& =SO{C}^{\left(t-1\right)}{E}_{avl}+{\eta }_{bess}{P}_{chr}^{\left(t\right)}\Delta t-{\displaystyle \frac{1}{{\eta }_{bess}}}{P}_{dch}^{\left(t\right)}\Delta t\hfill \\ \hfill \Rightarrow SO{C}^{\left(t\right)}& =SO{C}^{\left(t-1\right)}+{\displaystyle \frac{{\eta }_{bess}\Delta t}{E_{avl}}}{P}_{chr}^{\left(t\right)}-{\displaystyle \frac{\Delta t}{{\eta }_{bess}{E}_{avl}}}{P}_{dch}^{\left(t\right)},\forall t\in {\mathbb{S}}_N\hfill \end{array}\end{array}$

Although Equation (12) is recursive, it can also be written directly, bringing up the initial state of charge $SO{C}^{\left(0\right)}$, i.e.,

(13)$SO{C}^{\left(t\right)}=SO{C}^{\left(0\right)}+\sum _{\begin{array}{c}k=1\end{array}}^t\left({\displaystyle \frac{{\eta }_{bess}\Delta t}{E_{avl}}}{P}_{chr}^{\left(k\right)}-{\displaystyle \frac{\Delta t}{{\eta }_{bess}{E}_{avl}}}{P}_{dch}^{\left(k\right)}\right),\forall t\in {\mathbb{S}}_N$

Furthermore, the SOC may vary within a limited range with maximal ($\overline{\mathrm{SOC}}$) and minimal ($\underline{\mathrm{SOC}}$) values,

(14)$\begin{array}{c}\hfill \underline{\mathrm{SOC}}\le SO{C}^{\left(t\right)}\le \overline{\mathrm{SOC}},\forall t\in {\mathbb{S}}_N\end{array}$

As manufacturers usually define the battery cycle life ($L_{cyc}$) at a rated depth of discharge ($DO{D}_r$), then it can be useful to consider $\underline{\mathrm{SOC}}=1-DO{D}_r$, if the goal is to reach the amount of $L_{cyc}$ cycles.

Finally, the substitution of Equation (13) in (14) results in N MILP constraints for SOC, bringing out the decision variables $P_{chr}$ and $P_{dch}$,

(15)$\begin{array}{c}\hfill \underline{\mathrm{SOC}}-{SOC}^{\left(0\right)}\le \sum _{\begin{array}{c}k=1\end{array}}^t\left({\displaystyle \frac{{\eta }_{bess}\Delta t}{E_{avl}}}{P}_{chr}^{\left(k\right)}-{\displaystyle \frac{\Delta t}{{\eta }_{bess}{E}_{avl}}}{P}_{dch}^{\left(k\right)}\right)\le \overline{\mathrm{SOC}}-{SOC}^{\left(0\right)},\forall t\in {\mathbb{S}}_N\end{array}$

Furthermore, by setting $t=N$ in Equation (13), it is possible to add an equality constraint to set the desired final state of charge $SO{C}^{\left(N\right)}$, that is,

(16)$\sum _{\begin{array}{c}k=1\end{array}}^N\left({\displaystyle \frac{{\eta }_{bess}\Delta t}{E_{avl}}}{P}_{chr}^{\left(k\right)}-{\displaystyle \frac{\Delta t}{{\eta }_{bess}{E}_{avl}}}{P}_{dch}^{\left(k\right)}\right)=SO{C}^{\left(N\right)}-SO{C}^{\left(0\right)}$

However, to make a more balanced economic analysis, at the end of an N period, the state of charge should return to its initial value, as in [10], which can be reached by setting $SO{C}^{\left(N\right)}=SO{C}^{\left(0\right)}$ in Equation (16).

The state variables for charging and discharging are binary, i.e., $\{s{t}_{chr},s{t}_{dch}\}\in \{0,1\}$. As the system cannot perform charging and discharging operations at the same time, then

(17)$0\le {st}_{chr}^{\left(t\right)}+{st}_{dch}^{\left(t\right)}\le 1,\forall t\in {\mathbb{S}}_N$

When a state variable disables operation of charging or discharging, naturally, the respective power decision variable must be forced to zero. Otherwise, the power is limited by the BESS rated power ($P_r$). Then, for all $t\in {\mathbb{S}}_N$,

(18a)$\begin{array}{c}\hfill \begin{array}{cc}\hfill 0\le P_{chr}^{\left(t\right)}& \le {st}_{chr}^{\left(t\right)}{P}_r\hfill \end{array}\end{array}$

(18b)$\begin{array}{c}\hfill \begin{array}{cc}\hfill 0\le P_{dch}^{\left(t\right)}& \le {st}_{dch}^{\left(t\right)}{P}_r\hfill \end{array}\end{array}$

3.1.5. Directly Controllable Loads and Load Shedding

The state variable of load shedding is an integer of bounds $0\le {st}_{shd}^{\left(t\right)}\le 1$, $\forall t\in {\mathbb{S}}_N$, and the amount of load shedding is bounded by ${\overline{P}}_{shd}^{\left(t\right)}$,

(19)$0\le P_{shd}^{\left(t\right)}\le {st}_{shd}^{\left(t\right)}{\overline{P}}_{shd}^{\left(t\right)},\forall t\in {\mathbb{S}}_N$

Typically, ${\overline{P}}_{shd}^{\left(t\right)}={\mathsf{\gamma }}_{shd}{\widehat{P}}_{load}^{\left(t\right)}$, as in [10], where $\{{\mathsf{\gamma }}_{shd}\in \mathbb{R}\mid 0\le {\mathsf{\gamma }}_{shd}\le 1\}$ is a factor of shedding.

Since the optimizer algorithm could schedule load shedding to avoid buying energy when it is more expensive, it may be of interest to refuse this type of operation. For example, it may be useful to allow disconnecting the load only during island condition, which requires that

(20)$0\le {st}_{shd}^{\left(t\right)}\le 1-{st}_{con}^{\left(t\right)},\forall t\in {\mathbb{S}}_N$

In the model of directly controllable loads illustrated in Figure 2, the amount of interruptible loads $P_{int}^{\left(t\right)}$, included in the forecasted daily load ${\widehat{P}}_{load}^{\left(t\right)}$, is available for the optimizer algorithm to decide if it is necessary to disconnect it and when to do it. On the other hand, the amount of shiftable loads $P_{shf}^{\left(t\right)}$ has different requirements, and the optimizer must connect this load block optimally throughout the day and ensure that it has a continuous cycle.

The state variables for interruptible and shiftable loads are integers such that $0\le {st}_{int}^{\left(t\right)}\le 1$ and $0\le {st}_{shf}^{\left(t\right)}\le 1$, $\forall t\in {\mathbb{S}}_N$, respectively. Now assuming that, according to the demand response program, the period of time that the MG operator can keep interruptible loads disconnected from the system is limited to $T_{int}\in {\mathbb{S}}_N$, thus,

(21)$\sum _{t=1}^N{st}_{int}^{\left(t\right)}\le T_{int}$

Furthermore, the amount of interruptible loads is bounded by ${\overline{P}}_{int}^{\left(t\right)}$, i.e.,

(22)$0\le P_{int}^{\left(t\right)}\le {st}_{int}^{\left(t\right)}{\overline{P}}_{int}^{\left(t\right)},\forall t\in {\mathbb{S}}_N$

In fact, ${\overline{P}}_{int}^{\left(t\right)}$ can vary throughout the day according to the characteristics of the interruptible load group available for the demand response program. Hence, in this work, interruptible loads are modeled as a fraction of the forecasted load, such that ${\overline{P}}_{int}^{\left(t\right)}={\mathsf{\gamma }}_{int}{\widehat{P}}_{load}^{\left(t\right)}$, where $\{{\mathsf{\gamma }}_{int}\in \mathbb{R}\mid 0\le {\mathsf{\gamma }}_{int}\le 1\}$.

Assuming the block of shiftable loads has a cycle of period $T_{shf}\in {\mathbb{S}}_N$ and magnitude ${\overline{P}}_{shf}$, then Equations (23)–(25) are constraints required for power and time of connection,

(23)$0\le P_{shf}^{\left(t\right)},\forall t\in {\mathbb{S}}_N$

(24)$P_{shf}^{\left(t\right)}-{st}_{shf}^{\left(t\right)}{\overline{P}}_{shf}=0,\forall t\in {\mathbb{S}}_N$

(25)$\sum _{t=1}^N{st}_{shf}^{\left(t\right)}=T_{shf}$

However, to meet the continuous load cycle requirement, thus avoiding a sequence of turning the load on and off throughout the day, an additional restriction is required, which is

(26)$\begin{array}{cc}\hfill & \mathrm{for}k=1,2,\dots ,N-T_{shf}\hfill \\ \hfill & \mathrm{for}n=0,1,\dots ,N-k-T_{shf}\hfill \\ \hfill & {st}_{shf}^{\left(k\right)}+{st}_{shf}^{(k+T_{shf}+n)}\le 1\hfill \\ \hfill & \mathrm{end}\hfill \\ \hfill & \mathrm{end}\hfill \end{array}$

Algebraically, Equation (26) intents to ensure the time distance between any two enabled states is less than one period of load cycle $T_{shf}$. In this equation, a pseudocode technique was adopted as a mathematical resource to represent the inequality restrictions in a compact and readable way.

To add other blocks of shiftable loads with different magnitudes, periods, and states to the model, consider add to the problem the set of variables $\{P_{shf}^{(r,t)},T_{shf}^{\left(r\right)},{st}_{shf}^{(r,t)}\}$ for the r-th load block, which must meet the constraints of Equations (23)–(26) for each r. However, the computational cost can be high, since the number of restrictions per load block is $N+1+(N-T_{shf})(N-T_{shf}+1)/2$. For example, each block of two hours of shiftable loads requires 4013 constraints in the MILP problem ($\Delta t=1/4$ h).

3.1.6. PV System

In a day-ahead schedule problem, the PV output power variable is usually equal to the forecasted value for the next day. However, there can be conditions in which the microgrid has surplus generation but is unable to transfer energy to the grid (technical restrictions or islanding condition) or to the battery (SOC is at the upper limit, or the battery is out of service for maintenance). In such a situation there must be generation curtailment $P_{pvc}^{\left(t\right)}$ to satisfy the power balance restriction. In addition, $P_{pvc}^{\left(t\right)}$ must be limited by the amount of generation expected for the next day. Such requirements result in the restriction of Equation (27) for the PV system.

(27)$0\le P_{pvc}^{\left(t\right)}\le {\widehat{P}}_{pv}^{\left(t\right)}s{t}_{pv},\forall t\in {\mathbb{S}}_N$

The input variable $s{t}_{pv}$ must be set equal to zero if the PV system is out of service for maintenance; otherwise, it must be equal to 1. In fact, according to the modeling assumed in Equation (3), the daily cost of the PV System $C_{pv}$ remains for all days of the year, even when the condition $s{t}_{pv}=0$ is valid (for some days of the year). According to the authors in [35], solar panels suffer degradation mainly due to climatic and environmental conditions. Therefore, in this modeling, the degradation of the PV system is assumed to be a continuous process that occurs even if it is turned off.

3.1.7. Matrix Formulation

Formulating the optimization problem in matrix notation may be more useful for performing it in computer algorithms. To this end, the solution vector can be split into a state vector ${\mathbf{x}}_{st}$ and a power vector ${\mathbf{x}}_p$, i.e., $\mathbf{x}=[{\mathbf{x}}_{st};{\mathbf{x}}_p]$, where,

(28a)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{x}}_{st}=\left[1;\left\{{st}_{{pur}_g}^{\left(t\right)}\right\};\left\{{st}_{{sel}_g}^{\left(t\right)}\right\};\left\{{st}_{{pur}_{\mathsf{\mu }}}^{(j,t)}\right\};\left\{{st}_{{sel}_{\mathsf{\mu }}}^{(j,t)}\right\};\left\{{st}_{shd}^{\left(t\right)}\right\};\left\{{st}_{int}^{\left(t\right)}\right\};\left\{{st}_{shf}^{\left(t\right)}\right\};\left\{{st}_{chr}^{\left(t\right)}\right\};\left\{{st}_{dch}^{\left(t\right)}\right\}\right]\end{array}\end{array}$

(28b)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{x}}_p=\left[\left\{P_{pvc}^{\left(t\right)}\right\};\left\{P_{chr}^{\left(t\right)}\right\};\left\{P_{dch}^{\left(t\right)}\right\};\left\{P_{shd}^{\left(t\right)}\right\};\left\{P_{int}^{\left(t\right)}\right\};\left\{P_{shf}^{\left(t\right)}\right\};\left\{P_{{pur}_g}^{\left(t\right)}\right\};\left\{P_{{sel}_g}^{\left(t\right)}\right\};\left\{P_{{pur}_{\mathsf{\mu }}}^{(j,t)}\right\};\left\{P_{{sel}_{\mathsf{\mu }}}^{(j,t)}\right\}\right]\end{array}\end{array}$

Elements inside braces are row vectors; those with index $\left(t\right)$ are $N\times 1$ vectors, and with index $(j,t)$, are $NM\times 1$ vectors, i.e., a chain of vectors with index $\left(t\right)$; e.g., $\left\{{st}_{{sel}_{\mathsf{\mu }}}^{(j,t)}\right\}=\left[\left\{{st}_{{sel}_{\mathsf{\mu }}}^{(1,t)}\right\};\dots ;\left\{{st}_{{sel}_{\mathsf{\mu }}}^{(M,t)}\right\}\right]$. Elements in the state vector are integers (binaries) and in the power vector, real numbers. Furthermore, as $E_{bill}\left(\mathbf{x}\right)$ and the constraints are linear functions, the model fits a MILP problem. Furthermore, ${\mathbf{x}}_{st}$ and ${\mathbf{x}}_p$ are $\left(N(2M+7)+1\right)\times 1$ and $2N(M+4)\times 1$ vectors; therefore, the solution vector has $N(4M+15)+1$ variables. As an example, a scenario with one external MG and $\Delta t=1/4$ h ($M=1$ and $N=96$) requires a MILP problem with 1825 decision variables, of which 865 are integers and 960 real numbers.

Similarly, the parameters presented in Table 4 and Table 5 can be grouped into a vector, i.e.,

(29)$\mathbf{f}=[{\mathbf{f}}_{st};{\mathbf{f}}_p]$

where

(30a)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{f}}_{st}=\left[C_{pv};\mathbf{1}{c}_{{pur}_g};\mathbf{1}{c}_{{sel}_g};\left\{\mathbf{1}{c}_{{pur}_{\mathsf{\mu }}}^{\left(j\right)}\right\};\left\{\mathbf{1}{c}_{{sel}_{\mathsf{\mu }}}^{\left(j\right)}\right\};\mathbf{1}{c}_{shd};\mathbf{1}{c}_{int};\mathbf{1}{c}_{shf};\mathbf{1}{c}_{chr};\mathbf{1}{c}_{dch}\}\right]\end{array}\end{array}$

(30b)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{f}}_p=\left[\left\{p_{pvc}^{\left(t\right)}\right\};\mathbf{1}{p}_{chr};\mathbf{1}{p}_{dch};\left\{p_{shd}^{\left(t\right)}\right\};\left\{p_{int}^{\left(t\right)}\right\};\left\{p_{shf}^{\left(t\right)}\right\};\left\{p_{{pur}_g}^{\left(t\right)}\right\};-\left\{p_{{sel}_g}^{\left(t\right)}\right\};\left\{p_{{pur}_{\mathsf{\mu }}}^{(j,t)}\right\};-\left\{p_{{sel}_{\mathsf{\mu }}}^{(j,t)}\right\}\right]\Delta t\end{array}\end{array}$

where $\mathbf{1}$ is a $N\times 1$ vector of ones; also, by definition $\left\{\mathbf{1}{c}_{{pur}_{\mathsf{\mu }}}^{\left(j\right)}\right\}:=\left[\left\{\mathbf{1}{c}_{{pur}_{\mathsf{\mu }}}^{\left(1\right)}\right\};\dots ;\left\{\mathbf{1}{c}_{{pur}_{\mathsf{\mu }}}^{\left(M\right)}\right\}\right]$ and $\left\{\mathbf{1}{c}_{{sel}_{\mathsf{\mu }}}^{\left(j\right)}\right\}:=\left[\left\{\mathbf{1}{c}_{{sel}_{\mathsf{\mu }}}^{\left(1\right)}\right\};\dots ;\left\{\mathbf{1}{c}_{{sel}_{\mathsf{\mu }}}^{\left(M\right)}\right\}\right]$. The same notation rule for vectors used in ${\mathbf{x}}_{st}$ can be applied to the elements in Equations (30a) and (30b).

Finally, the MILP problem written in matrix notation is

(31)$\begin{array}{cccc}\hfill & \underset{\mathbf{x}}{\mathrm{minimize}}\hfill & \hfill & {\mathbf{f}}^{\mathrm{T}}\mathbf{x}\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \left\{\begin{array}{c}\mathbf{A}\mathbf{x}\le \mathbf{b}\hfill \\ {\mathbf{A}}_{eq}\mathbf{x}={\mathbf{b}}_{eq}\hfill \\ {\mathbf{l}}_b\le \mathbf{x}\le {\mathbf{u}}_b\hfill \\ {\mathbf{x}}_{st}\in \mathbb{Z}\hfill \end{array}\right.\hfill \end{array}$

where $\mathbf{A}$ and $\mathbf{b}$ are the matrix and vector of linear inequality constraints, respectively; ${\mathbf{A}}_{eq}$ and ${\mathbf{b}}_{eq}$ are the matrix and vector of linear equality constraints, respectively; and ${\mathbf{l}}_b$ and ${\mathbf{u}}_b$ are the vectors of lower and upper bounds, respectively. Table 7 identifies the equations required to define and create the matrices and vectors of constraints.

3.2. Cost Estimation

3.2.1. BESS Costs

The calculation of the availability cost of a BESS ($p_{bess}$) in $/kWh requires estimating the total charge and discharge energy of the battery over its lifetime ($L_{E,bess}$), and knowing its capital cost ($C_{cap,bess}$), i.e.,

(32)$p_{bess}=C_{cap,bess}/L_{E,bess}(\$/\mathrm{kWh})$

Theoretically, the battery life ends when it reaches the SOH threshold ($SO{H}_{tsh}$) for a given application. Thus, the curve SOH versus cycles $\left(SOH\left(n_{cyc}\right)\right)$ can bring out the cycle life $\left(L_{cyc}\right)$ of the battery, since $SOH\left(L_{cyc}\right)=SO{H}_{tsh}$. Although it is possible to know the true $SOH\left(n_{cyc}\right)$ only in time of operation from the BMS estimates, the curve provided by the BESS manufacturer can be used for cost estimation purposes.

The IEEE Std 2030.2.1™ [34] defines cycle as “the process of BESS discharging or charging from initial state of charge to the same state within a single discharge and charge”. A cycle can be full if it has $DOD=1$, or partial if $DOD<1$ [36]. In a characterization testing, a battery manufacturer must determine how many full cycles (or partial cycles at a rated depth of discharge) the battery will achieve [36], and then build the $SOH\left(n_{cyc}\right)$ curve. Therefore, the charge and discharge energy in a manufacturer’s test cycle is $2E_{avl}DO{D}_r$, or similarly, $2SOH\left(n_{cyc}\right)E_{r}DO{D}_r$. As a consequence, a lower bound on the battery charge and discharge energy over its lifetime can be

(33)$L_{E,bess}\ge 2E_{r}DO{D}_r{\int }_0^{L_{cyc}}SOH\left(n_{cyc}\right)d{n}_{cyc}$

as long as during the BESS operation, the condition $DOD\left(n_{cyc}\right)\le DO{D}_r$ is satisfied. Equation (33) is a lower limit because it considers that charging and discharging operations whose $DOD$ is less than $DO{D}_r$, as it happens under real-world conditions, prolong battery life as pointed out by authors in [37]. In addition, calculating the integral of that equation requires knowing an analytical $SOH\left(n_{cyc}\right)$ function that fits the empirical curve from the manufacturer. For lithium-ion batteries, an exponentially decreasing relationship between the state of health and the number of cycles can approximate a real $SOH\left(n_{cyc}\right)$. Thus, consider the base function $y\left(x\right)=e^{-x}$. Applying the expansion technique on the x-axis, and compression and displacement on the y-axis, this function can be adjusted to the desired exponential $SO{H}_{exp}\left(n_{cyc}\right)$, i.e.,

(34)$SO{H}_{exp}(n_{cyc},k_e)=k_1{e}^{-n_{cyc}/k_3}+k_2,n_{cyc}\ge 0$

where $k_1=(1-SO{H}_{tsh})/k_e$, $k_2=1-k_1$, $k_3=-L_{cyc}/ln(1-k_e)$, and $k_e\in \mathbb{R}$, $0<k_e\le 1$, is a factor of non-linearity. Note that $SO{H}_{exp}(0,k_e)=1$, and $SO{H}_{exp}(L_{cyc},k_e)=SO{H}_{tsh}$.

It is possible to show (through limits calculation) that Equation (34) is a family of curves between a linear approximation for SOH and the line $SO{H}_{tsh}$, for $k_e$ ranging from 0 to 1, respectively, as illustrated in Figure 4. In fact, the non-linearity factor allows a fine adjustment to the SOH curve of interest, as the one provided by the manufacturer. As an example, data in [38] show $SOH\left(n_{cyc}\right)$ curve of a Lithium-ion BESS of 6000 cycle life, and $SO{H}_{tsh}$ of 0.8, which fits $SO{H}_{exp}(n_{cyc},k_e)$ for $k_e\approxeq 0.55$.

Now, applying Equation (34) in (33),

(35)$\begin{array}{c}\hfill L_{E,bess}\ge 2E_{r}DO{D}_r{L}_{cyc}\left(\frac{SO{H}_{tsh}-1}{ln(1-k_e)}+\frac{SO{H}_{tsh}-1}{k_e}+1\right)\left(\mathrm{kWh}\right)\end{array}$

As an example, consider a ${\mathrm{LiFePO}}_4$ BESS of 140 kW/280 kWh, $SO{H}_{tsh}$ of 0.8, 6000 cycle life at $DO{D}_r$ of 0.90, ${\eta }_{bess}$ of 0.92 (according to Table 6), where the cost of capital $C_{cap,bess}$ is 91,000 USD (real market data). Thus, using Equations (32) and (35), and assuming $k_e=0.55$, the result is $L_{E,bess}\ge \mathrm{2,681,776}$ kWh, and $p_{bess}\le 0.0339$ USD/kWh, respectively.

As $p_{bess}$ represents the cost of charging and discharging each kWh at the battery cells node, then the BESS powers in the objective function, whose reference is the main (bus) node, must be weighted by the efficiency of BESS, which is similar to consider

(36a)$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathit{cost}\mathit{of}\mathit{charging}:\Delta t{P}_{sto}{p}_{bess}=\Delta t{P}_{chr}{\eta }_{bess}{p}_{bess}=\Delta t{P}_{chr}{p}_{chr},\Rightarrow p_{chr}={\eta }_{bess}{p}_{bess}\end{array}\end{array}$

(36b)$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathit{cost}\mathit{of}\mathit{discharging}:\Delta t{P}_{sto}{p}_{bess}=\Delta t\frac{P_{dch}}{{\eta }_{bess}}{p}_{bess}=\Delta t{P}_{dch}{p}_{dch},\Rightarrow p_{dch}=\frac{1}{{\eta }_{bess}}{p}_{bess}\end{array}\end{array}$

where $\Delta t{P}_{sto}$ is the portion of energy from $\Delta t{P}_{chr}$ that was stored.

Finally, in energy arbitrage applications, a microgrid can purchase energy during an off-peak period, at a price $p_{out}$ ($/kWh), and store it for use during an on-peak period, when the energy is more expensive ($p_{peak}$). Thus, the unit cost of energy after the storage process must be less than the on-peak price, i.e.,

(37)$\Delta t\left(P_{pur}{p}_{out}+P_{chr}{p}_{chr}+P_{dch}{p}_{dch}\right)<\Delta t{P}_{dch}{p}_{peak}$

Now, making $P_{pur}=P_{chr}=P_{dch}/{\eta }_{bess}^2$ in Equation (37), and using Equations (36a) and (36b),

(38)$\frac{1}{{\eta }_{bess}^2}{p}_{out}+\frac{2}{{\eta }_{bess}}{p}_{bess}<p_{peak},\Rightarrow p_{bess}<\frac{1}{2}{\eta }_{bess}\left(p_{peak}-\frac{1}{{\eta }_{bess}^2}{p}_{out}\right)$

If this relationship is not satisfied, the optimizer algorithm may decide not to use the storage system to perform energy arbitrage.

3.2.2. PV Costs

The capacity of solar photovoltaic generation can vary considerably according to the region where the solar panels are installed. For example, according to data in [39], the region of Curitiba (Brazil) has a generation capacity ${\lambda }_{reg}$ of about $1261.57$ kWh/year per kW, at standard test conditions (STC), of solar panel (tilt ${24}^{\circ }$, azimuth $0^{\circ }$), while the New York City has ${\lambda }_{reg}=1360.49$ (tilt ${37}^{\circ }$, azimuth ${180}^{\circ }$), considering system losses of 14% and inverter efficiency of 96% in both cases. As a consequence, given an annual amount of energy $E_{year}^0$ (kWh/year) to be generated by a PV system, its capital cost $C_{cap,pv}$ depends on the region of installation due to the local price and generation capacity, i.e.,

(39)$C_{cap,pv}=\frac{E_{year}^0}{{\lambda }_{reg}}{p}_{reg},\left(\$\right)$

where $p_{reg}$ ($/kW) is the local cost per kW of a installed PV system.

A PV system with an annual degradation rate of ${\kappa }_{pv}$ (%/year) should produce throughout its lifespan $L_{spv}$ (years) an annual amount of energy of,

(40)$E_{year}^{\left(n\right)}=E_{year}^0\left(1-\frac{{\kappa }_{pv}}{100}n\right),\left(\mathrm{kWh}\right),\mathrm{for}n=0,1,\dots ,L_{spv}-1$

Similarly, the annual cost of generation can be proportional to the yearly energy generated, that is,

(41)$C_{year}^{\left(n\right)}=C_{year}^0\left(1-\frac{{\kappa }_{pv}}{100}n\right),\left(\$\right),\mathrm{for}n=0,1,\dots ,L_{spv}-1$

required that

(42)$C_{cap,pv}=\sum _{n=0}^{L_{spv}-1}{C}_{year}^0\left(1-\frac{{\kappa }_{pv}}{100}n\right),\Rightarrow C_{cap,pv}=C_{year}^0{L}_{spv}\left(1-{\kappa }_{pv}\frac{L_{spv}-1}{200}\right)$

Now, applying Equation (42) in Equation (39) to obtain the value of $C_{year}^0$ and then using the result in Equation (41), results in

(43)$C_{year}^{\left(n\right)}=\frac{(E_{year}^0/{\lambda }_{reg})p_{reg}}{L_{spv}(1-({\kappa }_{pv}/200)(L_{spv}-1))}\left(1-\frac{{\kappa }_{pv}}{100}n\right),\left(\$\right),\mathrm{for}n=0,1,\dots ,L_{spv}-1$

For application in the day-ahead optimization model, the cost values must be daily. Thus, assuming an average daily cost of generation $C_{pv}^{\left(n\right)}=C_{year}^{\left(n\right)}/365$ and an average daily capacity of generation $E_{pv}=E_{year}^0/365$ in Equation (43) results in,

(44)$C_{pv}^{\left(n\right)}=\frac{(E_{pv}/{\lambda }_{reg})p_{reg}}{L_{spv}(1-({\kappa }_{pv}/200)(L_{spv}-1))}\left(1-\frac{{\kappa }_{pv}}{100}n\right),\left(\$\right),\mathrm{for}n=0,1,\dots ,L_{spv}-1$

Authors in [40] summarized photovoltaic degradation rates using field tests reported in the literature. Results show an average degradation rate of $0.8$%/year, which can be a practical value for ${\kappa }_{pv}$. Furthermore, manufacturers usually offer a 25-year standard solar panel warranty, which can be a practical value for $L_{spv}$.

3.2.3. Other Costs

Costs in Table 4 are related to the operating states of the virtual disconnect switches, whereas those in Table 5 are related to the amount of energy for each operation. When participating in a free energy market, a microgrid may incur third-party service costs due to energy purchase and sale operations, which can be modeled by $c_{pur}$ and $c_{sel}$. On the other hand, the market prices are covered in the model by $p_{pur}$ and $p_{sel}$. Furthermore, in the present model, $c_{shd}$ represents a penalty, imposed by the regulatory agency, for each operation of load disconnection performed by the microgrid, while it lasts. On the other hand, $p_{shd}$ may represent a penalty due to the amount of disconnected load.

Finally, according to the modeling in the present study, if a cost does not apply to a given approach, it can be set zero.

4. Simulation Methodology

All results of the present paper were obtained through the implementation of the MG model, described in Section 3.1, in the MATLAB^® software. The optimizer function intlinprog() was used to solve the MILP problem from Equation (31), subject to the constraints from Table 7. Furthermore, it was used a time resolution $\Delta t$ of 1/4 h, and $N=96$ time slots.

The BESS is a 140 kW/280 kWh LiFePO${}_4$ system with coupling transformer, where: ${\eta }_{bess}=0.92$ (see Table 6), $DO{D}_r=0.90$, $SOH=0.90$, $L_{cyc}=6000$, $SO{H}_{tsh}=0.80$, $SO{C}^{\left(0\right)}=SO{C}^{\left(N\right)}=0.40$, $\overline{SOC}=1$, $\underline{SOC}=1-DO{D}_r$, and $k_e=0.55$. Furthermore, according to the example from Section 3.2.1, $p_{bess}=0.0339$ USD/kWh, and from Equation (36), $p_{chr}$ and $p_{dch}$ are $0.0312$ and $0.0369$ USD/kWh respectively.

The solar PV system is the same from Section 3.2.2, for the city of Curitiba, with: ${\lambda }_{reg}=1261.57$ kWh/year per kW, $L_{spv}=25$ years, ${\kappa }_{pv}=0.8$%, $p_{reg}=2060$ USD (market price). The average daily PV generation $E_{pv}$ was defined according to the daily load $E_{load}$, i.e., $E_{pv}:={MG}_{\mathrm{SS}}{E}_{load}$, where ${MG}_{\mathrm{SS}}$ is the microgrid self-sufficiency, and $E_{load}$ is constant and equal to 2400 kWh. Thus, the daily cost of the PV system for the first year of use $(n=0)$, calculated from Equation (44) and used in the simulations of this work is a function of the $M{G}_{ss}$. $C_{pv}=173.40$ USD, for ${MG}_{\mathrm{SS}}=1$, was used in the first three simulation cases. Furthermore, the PV generation curve was retrieved from [41], whose data are from a region in Africa and were collected on a sunny day of summer.

The load curve comes from field measurements of a group of middle-income residential consumers, which were collected during a business day [42].

The Brazilian time-of-use tariff, named White Tariff, was used to set the prices used in the simulations of this work. In the Brazilian captive energy market, customers can only buy electricity from an agent authorized to do the distribution service, typically the DSO. A regulatory agency (ANEEL) sets the price of energy, and the customer cannot negotiate it. However, the White Tariff allows customers to manage their consumption voluntarily according to the energy price. The regulated ToU time blocks are on-peak, off-peak, and pre and post-peak. Each energy distribution company can set its time blocks according to the load demand in its operation area. Table 8 shows the time blocks with their respective average periods [43]. Furthermore, it presents the weighted average of rates charged in Brazil for the conventional and ToU tariff. They were calculated considering data from 11 major energy distributors in Brazil [43], weighted by the number of customers of each company. In this table, ToU prices and time blocks define $p_{{pur}_g}^{\left(t\right)}$ and $p_{{sel}_g}^{\left(t\right)}$ (USD/kWh) used in this model. The conventional tariff was used as the reference price $p_{ref}$ in Equation (45).

In load shedding simulations, it was adopted $p_{shd}=3p_{ref}$ ($/kWh) as a penalty for the amount of disconnected load. Furthermore, in directly controllable loads simulations, it was adopted $p_{int}=2p_{ref}$ and $p_{shf}=0$, with $T_{int}=4$ (one hour) and $T_{shf}=10$ (two hours and a half). All other costs whose values have not been mentioned have been set to zero.

Finally, although the model presented here is prepared to simulate external MGs supplying and acquiring energy from the microgrid, this feature was not explored in any simulation because it does not fit the adopted ToU tariff model.

4.1. Normalized Energy Bill

For results analysis purposes, it may be preferable to handle a normalized energy bill than to work with monetary values. To this end, consider a reference daily MG energy bill $R_{bill}$ in which all distributed energy resources are turned off, that is,

(45)$R_{bill}:=\Delta t\sum _{t=1}^N{P}_{load}^{\left(t\right)}{p}_{ref}^{\left(t\right)}$

where $p_{ref}$ stands for a reference energy price. Therefore, the normalized energy bill is

(46)$E_{\eta }:={\displaystyle \frac{E_{bill}}{R_{bill}}}$

$E_{\eta }$ can represent three distinct regions of operation for an MG: $E_{\eta }>1$, $0\le E_{\eta }\le 1$, and $E_{\eta }<0$, as Figure 5 shows. In the first, there was a relative increase in the energy bill; as an example, a microgrid could operate in this region when it is connected to the main grid with generation off, BESS in operation, and the purchase price of energy is constant throughout the day. Naturally, such a region should be avoided. In the second, there was a relative reduction in the energy bill; a microgrid can operate in this region when in the first mode of operation, however with energy prices favoring the practice of energy arbitrage. Finally, a microgrid can operate in the third region when in the previous operating modes, but with enough generation surplus to meet the load and still sell energy.

4.2. Simulation Cases

• Energy arbitrage as a result of the Brazilian White Tariff applied to a microgrid;

• MG scheduled intentional islanding using the resources of interruptible loads and PV output curtailment; two islanding periods: $T_{{isl}_1}=13$, from 02:00 to 05:15, and $T_{{isl}_2}=8$, from 10:30 to 12:30; $M{G}_{ss}=1$.

• MG scheduled intentional islanding using the resources of interruptible loads, PV output curtailment, load shedding, and shiftable loads; two islanding periods: $T_{{isl}_1}=16$, from 02:00 to 06:00, and $T_{{isl}_2}=8$, from 10:30 to 12:30; $M{G}_{ss}=1$.

• The impact of the microgrid self-sufficiency on the normalized energy bill; limit for energy import and export: ${\overline{P}}_{pcc}=800$ kW; curves for $E_r=12$, 16, 20, and 30% of $E_{pv}$.

5. Results

In power system buses with loads and generators connected, there is a convention in which the power generated is positive while the power consumed is negative [44]. Although in the modeling presented in this work all powers are positive, such a convention of signs was adopted in the figures throughout this section. This minimizes overlapping curves in the figures. Thus, the curves $P_{load}$, $P_{sel}$, $P_{chr}$, and $P_{shf}$ are shown with negative values. Furthermore, these quantities are included in parentheses in the legends of the figures. Figure 6 presents the main findings of simulation 1. The pre and pos-peak, and on-peak time blocks are represented by $T_{pre,pos}$ e $T_{peak}$, respectively. Figure 6a shows the curves that make up the balance of active powers in the microgrid. Figure 6b presents the BESS state of charge curve with its upper and lower bounds.

Results show that, before solar generation begins, the microgrid has to import the energy necessary to supply its load. On the other hand, during the period of sunlight, the microgrid exports the PV generation surplus in addition to storing energy to use it during the on-peak period. As expected for the energy arbitrage simulation, the storage system reaches its maximum state of charge during the off-peak period and then performs a complete discharge to the lower limit during the on-peak. At the end of the day, the BESS is recharged to reach its final state of charge. Furthermore, the BESS performs a complete cycle throughout the day and presents an energy loss ($E_{loss}$) of 38 kWh, which is the loss per full cycle of the storage system.

Figure 7 shows the results of simulation 2. Two islanding periods $T_{{isl}_1}$ and $T_{{isl}_2}$ were added to the simulation. The storage system has the role of avoiding load shedding in the first islanding, minimizing curtailment in the second, and performing energy arbitrage during the on-peak period. Although the results show that BESS reaches its maximum state of charge before the first islanding and performs a complete discharge during the same (Figure 7b), the stored energy is not sufficient to supply 100% of the microgrid load during that period. Thus, the interruptible load resource is required to complete the remaining energy and avoid load shedding, as illustrated in Figure 7a and expanded in Figure 7c; $P_{int}$ presents four pulses of 15 min long, at 02:45, 03:15, 04:00, and 04:15, during the first islanding, which contributes approximately 10 kWh to the energy balance for that period. In practice, each pulse can represent a temporary power reduction or even shutdown of interruptible load groups, such as air conditioners and water heaters.

According to results, during the second islanding period, the BESS performs a full charge in response to the surplus of solar generation that cannot be injected into the main grid. Even so, a curtailment of $223.85$ kWh in the solar generation is expected for the next day. Besides, according to the optimal solution, all surplus energy produced outside the islanding must be exported to the main grid. Finally, BESS discharges during the on-peak period all energy stored during $T_{{isl}_2}$, performing the energy arbitrage, and then it returns to its initial state of charge.

Figure 8 illustrates the results of simulation 3, in which the first islanding period was increased by 45 min, and the shiftable load resource was added to the simulation. Results illustrated in Figure 8a show that although MG uses all its storage and interruptible loads resources (respecting the $T_{int}$ limit), it will be necessary to perform $61.98$ kWh of load shedding the next day, during the first islanding period.

However, the use of shiftable loads allowed a drop in curtailment from 223.85 to 127.85 kWh during the second islanding. In the optimal solution, the optimizer reallocated 96 of the 120 kWh of shiftable loads taken from the on-peak period to the islanding period; the remaining 24 kWh were allocated in the half-hour before the second islanding. It is important to emphasize the optimizer performs the reallocation of shiftable loads in compliance with the continuous cycle requirement for this type of load. Finally, on this schedule (and the previous one), there should be a loss of 76 kWh due to the two full cycles to be performed by BESS the next day. In addition, the behavior of the BESS at the on-peak period is similar to simulations 1 and 2.

In the simulation 4, for a daily load of 2400 kWh (the same from previous simulations), the capacity of the PV system and the storage system were gradually increased. The storage system was configured as a percentage of the PV system capacity. Furthermore, the imported and exported energy at the microgrid PCC was limited by setting ${\overline{P}}_{pcc}=800$ kW in the model. A simulation was performed for each setting, and an optimal solution was reached in each case. Figure 9 shows the curves $E_{\eta }\times M{G}_{ss}$ obtained from simulations. No directly controllable load resources were used in those simulations.

All the curves start in the increased energy bill area where $E_{\eta }>1$, because the microgrid does not have sufficient distributed energy resources to make economically attractive the participation in the White Tariff. From $M{G}_{ss}=0.4$ onwards, the curves enter the reduced energy expenses area where $0<E_{\eta }<1$; but only the curves of 16%, 20%, and 30% reach the profit or energy credit area ($E_{\eta }<0$), for $M{G}_{ss}$ greater than or equal to 3.5, 3.3, and 3.1, respectively. In the light blue region ($P_{sel}<{\overline{P}}_{pcc}$), the optimal use of the distributed energy resources is resulting in an energy export to the main grid whose power is less than the limit of 800 kW. On the other hand, the light green region contains the points for which this limit was reached. In that region, the generation surplus that cannot be exported must be stored by BESS. However, when the storage limit is reached, the curtailment process begins, which can be identified by the turning point of curves of 12%, 16%, and 20%. Similar reasoning can be applied to the MG without BESS curve. From the turning point onwards, the higher the capacity of the PV system, the more energy will be wasted, and by consequence, the more will be the cost of the microgrid, as can be seen in the figure for those curves. The turning point of the 30% curve is not shown in the figure as it must happen for $M{G}_{ss}>4.0$.

Finally, although in this work, the MILP problem has thousands of variables and restrictions, the execution time in all simulations was less than 15 s on a personal computer of general use.

6. Discussion

Simulations presented in this work provide support to evaluate the proposed MG model. In general, results corroborate the consistency of the microgrid mathematical modeling proposed in this work. In practice, the MG used in the simulations can correspond to a medium to large industry, or to a residential area containing more than a hundred houses. Simulation 1 shows the Brazilian White Tariff can be economically attractive for a microgrid with solar generation and energy storage systems. In this example, a savings of 20% $(E_{\eta }=0.7992)$ in the energy bill was observed. It is necessary to emphasize the energy arbitrage is directly related to the BESS costs. According to Equation (38), $p_{bess}$ should be less than 0.0544 USD/kWh to enable energy arbitrage in the White Tariff; the value used in simulations, whose calculation is based on market data, was $p_{bess}=0.0339$ USD/kWh. Furthermore, the microgrid load profile used in this study presents an increase in consumption during peak hours, which helps to reduce the costs of a microgrid with BESS and to perform energy arbitrage.

The results of the second simulation highlight the role of the storage system when the microgrid is disconnected from the main system. In the first islanding, for more than three hours, the BESS can supply the microgrid load without requiring load shedding; in the second, the storage system can reduce the PV curtailment from 450.65 to 223.85 kWh (50.33%) in the next day.

Results also show directly controllable loads can play a complementary role to that of the storage system. In the second simulation, the optimizer algorithm makes use of interruptible loads during the first islanding to avoid load shedding. In the third, although it is not possible to prevent shedding, the results show the optimizer first exhausts the resource of interruptible loads and then uses the load shedding one; the cost to use the former is less than the penalty for using the latter. Regarding shiftable loads, the results of the third simulation show their use can also help to minimize the PV curtailment, which was reduced from 450.65 to 127.85 kWh (71.63%). Thus, even with the penalties of load shedding, the costs in the third simulation are lower than in the second.

Results of simulation 4 show that a BESS system can reduce the energy costs of a microgrid with operational limits (as in practice) at the PCC bus, in addition to the reduction achieved by the PV system. In this model, the energy balance of a full day of operation must be zero when the MG self-sufficiency is unitary, even when there is no BESS. This is because the microgrid uses the main grid as a battery of infinite capacity and power limited by ${\overline{P}}_{pcc}$ to exchange energy with the grid both when it generates surplus energy and when it only demands energy. Thus, from $M{G}_{ss}>1$ onwards, the region in which there is an energy credit or revenue for the microgrid begins. However, making a profit from the generation of energy requires the revenue must exceed the costs of operating the microgrid. Thus, the results show that in practice reaching the profit region for this microgrid may require self-sufficiency as high as $M{G}_{ss}>3$. Furthermore, the size of the storage system must be adequate for this purpose; otherwise, such a region may not be hit, as illustrated in Figure 9.

It is necessary to emphasize that in Brazil, the current captive energy market does not allow customers to sell the surplus energy, even in the White Tariff; however, it allows the accumulation of energy credits. In this market, a scenario with $M{G}_{ss}>1$ may not be economically attractive, since day after day there will be an accumulation of energy credits that cannot be traded. However, the model presented in this work has the purpose of being generic enough to be applied in other types of markets, including a free energy market and future markets that include microgrids and renewable energy resources.

Considering that load and solar generation forecasts can vary considerably over the year in a microgrid, then addressing the day-ahead scheduling problem for a single day may be a limitation of this work. A day-by-day analysis covering the four seasons as well as variations in load and solar generation forecasts could result in a richer set of simulation data. This type of analysis can be considered in future works. Furthermore, although in this work the decision variable for curtailment is continuous (a real number), in practice, it may be necessary to model this variable by ranges due to converter system technology limitations, which would make it an integer in the MILP problem. In addition, although the values used in this work are based on real market values, the results presented here must be interpreted with caution, as they are dependent on market prices for the PV system, the BESS, and the energy tariff.

7. Conclusions

This work aimed to model the optimal day-ahead scheduling problem in microgrids, filling some gaps found in mathematical modeling. Scheduled intentional islanding, shiftable loads of continuous cycle, the methodology for calculating the availability cost of PV and battery systems, and a comprehensive and detailed method for modeling battery energy storage systems are examples of gaps in modeling filled throughout this work. Although the microgrid resources have been modeled in detail, the linearity of the problem has been preserved, which contributes to simplifying its implementation in optimization software, and it helps to focus the work effort on modeling rather than optimization algorithms. Computational simulations validated the proposed model regarding all the modeled distributed energy resources and scheduled intentional islanding situations. A normalized energy bill approach was used as the variable to optimizing, and simulations were carried out to investigate how the microgrid self-sufficiency and the capacity of the storage system impact the normalized bill. Results showed a practical difficulty in reaching the region of profit operation.

The detailed mathematical modeling presented and validated in this work, covering the three types of distributed energy resources of a microgrid (renewables, storage, and controllable loads), can be relevant for microgrid operators and managers in the process of seeking economic efficiency through optimized use of their resources. Future works may examine the application of the model to an energy market with a framework of rules for free energy trading between entities, which can be multiple microgrids and as well as generation units distributed in an energy distribution network. Future work may also include an annual analysis of the optimal day-ahead scheduling problem. Finally, a study on the computational complexity and scalability of intra-day control of microgrids based on real-time data may also be a subject of future works.

Author Contributions

V.A.S. performed the literature review, conceptualization, methodology, software codification, simulations, validation, and writing. A.R.A. and G.L.-T. contributed to conceptualization, methodology, validation, and performed the review, corrections, supervision, advising, and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Companhia Paranaense de Energia - COPEL research and technological development (RTD) program, through the PD-02866-0511/2019 project, regulated by ANEEL.

Acknowledgments

Alexandre Rasi Aoki thanks the grant for Technology Development of CNPq—National Council for Scientific and Technological Development within the Ministry of Science, Technology, Innovations, and Communications. Authors also thank the support of CAPES—Brazilian Federal Agency for Support and Evaluation of Graduate Education within the Ministry of Education of Brazil.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Schematic diagram of the MG day-ahead EMS model addressed in the present work, showing its relationship with the input and output data and its interaction with the real-time EMS; CD, ED, FD, LD, PD, and TD stand for cost, estimated, forecasted, limit, price, and technical data, respectively.

Figure 2. Schematic diagram of the microgrid model used in the MILP problem of the present study. The model presents integer (binary) and real decision variables. The former (st(t) type) work as virtual disconnectors and the latter (P(t) type) represent the active powers of each element of the model. Type P^ variables represent forecasted values.

Figure 3. Clarifying the nomenclature of BESS used throughout the present study. The definitions of SOC and DOD are following the IEEE Std 2030.2.1™ [34].

Figure 4. Exponential approximation curves for state of health, where SOHtsh=0.8 and Lcyc=6000.

Figure 5. Regions of MG operation represented by the normalized energy bill.

Figure 6. Results of simulation 1. (a) Results of simulation 1 for the power balance curves; Ebill = 249.36 USD; Eloss represents energy losses throughout the day with BESS charging and discharging operations; Tpre,post stands for the pre and post-peak period; Tpeak is the on-peak period; (b) Results of simulation 1 for the BESS SOC.

Figure 7. Results of simulation 2. (a) Results of simulation 2 for the power balance curves; [Epvc = 223.85, Eint = 10.16] (kWh); Ebill = 294.82 USD; (b) Results of simulation 2 for the BESS SOC; (c) Expanding the curve of interruptible loads during the first period of islanding.

Figure 8. Results of simulation 3. (a) Results of simulation 3 for the power balance curves; [Epvc = 127.85, Eshd = 61.98, Eint = 14.92] (kWh); Ebill = 285.92 USD; (b) Results of simulation 3 for the BESS SOC.

Figure 9. Normalized energy bill for a sequence of simulations with increasing values of microgrid self-sufficiency; curves for BESS rated capacity (Er) from 12% to 30% of the average daily PV system capacity (Epv); P¯pcc=800 kW; Eload=2400 kWh; Epv:=2400MGss kWh.

Table 1

State variables for the model, in which the value 1 means state enabled and 0 disabled.

State Variables	Description
${st}_{{pur}_g}^{\left(t\right)}$, ${st}_{{pur}_{\mathsf{\mu }}}^{(j,t)}$	represent the state of purchasing transaction at time t; index g stands for trading with the main grid and $\mathsf{\mu }$ trading with the j-th microgrid;
${st}_{{sel}_g}^{\left(t\right)}$, ${st}_{{sel}_{\mathsf{\mu }}}^{(j,t)}$	represent the status of energy sale transaction at time t; index g stands for trading with the main grid and $\mathsf{\mu }$ trading with the j-th microgrid;
${st}_{chr}^{\left(t\right)}$, ${st}_{dch}^{\left(t\right)}$	represent the state of charge and discharge of the battery, at time t, respectively;
${st}_{int}^{\left(t\right)}$, ${st}_{shf}^{\left(t\right)}$	represent the connection state of interruptible and shiftable loads, respectively.
${st}_{shd}^{\left(t\right)}$	represents the activation state of the load shedding resource.

Table 2

Power variables for the model.

Power Variables (kW)	Description
$P_{{pur}_g}^{\left(t\right)}$, $P_{{pur}_{\mathsf{\mu }}}^{(j,t)}$	are the powers injected by the main grid into the microgrid ($pur$ index) through the $PCC$, at time t. The index g stands for trading with the main grid and $\mathsf{\mu }$ trading with the j-th microgrid;
$P_{{sel}_g}^{\left(t\right)}$, $P_{{sel}_{\mathsf{\mu }}}^{(j,t)}$	are the powers injected by the microgrid into the main grid ($sel$ index) through the $PCC$, at time t. The index g stands for trading with the main grid and $\mathsf{\mu }$ trading with the j-th microgrid;
$P_{chr}^{\left(t\right)}$, $P_{dch}^{\left(t\right)}$	are the battery charging and discharging power, respectively.
$P_{pvc}^{\left(t\right)}$, $P_{shd}^{\left(t\right)}$	are the amount of PV output power curtailment and load shedding, respectively.
$P_{int}^{\left(t\right)}$, $P_{shf}^{\left(t\right)}$	are the amount of interruptible and shiftable loads, respectively.

Table 3

Input variables for the model.

Non-Decision Variables	Description
${\widehat{P}}_{load}^{\left(t\right)}$, ${\widehat{P}}_{pv}^{\left(t\right)}$	are the day-ahead forecasted load and PV output power curves, respectively (kW).
${st}_{con}^{\left(t\right)},{st}_{pv}$	are the microgrid and PV system state of connection, respectively; state 0 means disconnected and 1 connected.

Table 4

State parameters.

State Parameter ($)	Description
$c_{{sel}_g}$, $c_{{sel}_{\mathsf{\mu }}}$	are fixed costs of service when selling energy to the grid or to the j-th microgrid.
$c_{{pur}_g}$, $c_{{pur}_{\mathsf{\mu }}}$	are fixed costs of service when purchasing energy from the grid or from the j-th microgrid.
$c_{chr}$, $c_{dch}$	represent a fixed cost for charging and discharging the battery, respectively;
$c_{int}$, $c_{shf}$	represent a fixed cost for disconnecting interruptible loads, and for connecting shiftable loads, respectively;
$c_{shd}$, $C_{pv}$	represent a fixed cost for disconnecting load (shedding), and the availability cost of the PV system, respectively;

Table 5

Power parameters of the MILP problem.

Power Parameter ($/kWh)	Description
$p_{{sel}_g}^{\left(t\right)}$, $p_{{sel}_{\mathsf{\mu }}}^{(j,t)}$	are the sale prices of energy, at time t, when trading with the grid (g index) or with the j-th microgrid ($\mathsf{\mu }$ index);
$p_{{pur}_g}^{\left(t\right)}$, $p_{{pur}_{\mathsf{\mu }}}^{(j,t)}$	are the purchase prices of energy, at time t, when trading with the grid (g index) or with the j-th microgrid ($\mathsf{\mu }$ index);
$p_{chr}$, $p_{dch}$	are the cost of each kWh for charging and discharging the battery, at time t, respectively.
$p_{pvc}^{\left(t\right)}$, $p_{shd}^{\left(t\right)}$	are the cost of each kWh of PV output power curtailment and load disconnected (shedding), at time t, respectively.
$p_{int}^{\left(t\right)}$, $p_{shf}^{\left(t\right)}$	are the cost of each kWh to disconnect interruptible loads and to connect shiftable loads, at time t, respectively.

Table 6

Values of energy storage system efficiency for three battery technologies.

Battery Technology	Coupling Transformer	Efficiency
		${\mathbf{\eta }}_{\mathit{t}\mathit{r}\mathit{f}}$	${\mathbf{\eta }}_{\mathit{c}\mathit{o}\mathit{n}\mathit{v}}$	${\mathbf{\eta }}_{\mathit{b}\mathit{a}\mathit{t}}$	${\mathbf{\eta }}_{\mathit{b}\mathit{e}\mathit{s}\mathit{s}}$
Lead-acid	no	−	>$0.97$	$0.9220$	$0.8943$
	yes	>$0.97$	>$0.97$	$0.9220$	$0.8675$
Lithium-ion	no	−	>$0.97$	$0.9487$	$0.9202$
	yes	>$0.97$	>$0.97$	$0.9487$	$0.8926$
LiFePO${}_4$, and Lithium-ion of higher technology	no	−	>$0.97$	$0.9747$	$0.9454$
	yes	>$0.97$	>$0.97$	$0.9747$	$0.9170$

Table 7

Summary of the equations, organized by features, required to build the matrices and vectors of constraints of the MILP problem; $\mathbf{0}$ is a $N\times 1$ vector of zeros.

Features	Equations for A and b	Equations for ${\mathit{A}}_{\mathbf{eq}}$ and ${\mathit{b}}_{\mathbf{eq}}$	Equations for ${\mathit{l}}_{\mathit{b}}$ and ${\mathit{u}}_{\mathit{b}}$
State variables in general			$\mathbf{0}\le {\mathbf{x}}_{st}\le \mathbf{1}$
Energy exchange and islanding	(4)		(4)
Power exchange	(5), (6)		(5)
Power balance		(8)
BESS, state variables	(17)		(17)
BESS, power variables	(15), (18)	(16)	(18)
Load shedding, state variable			(20)
Load shedding, power variable	(19)		(19)
Interruptible load, state variable	(22), (21)
Interruptible load, power variable	(22)		(22)
Shiftable load, state variable	(26)	(24), (25)
Shiftable load, power variable		(24)	(23), (24)
PV output power curtailment		(27)	(27)

Table 8

The Brazilian White Tariff; average time blocks were computed from 104 companies data; average prices were calculated from the 11 major energy distributors in Brazil, weighted by the number of customers of each company.

Tariff	Average Time Blocks			Average Prices		Percentage Prices (%)
	Start Time	End Time		BRL/kWh	USD/kWh
Conventional tariff	−	−		0.558	0.130	$100$
ToU, off-peak	00:00	17:00		0.466	0.109	$84$
ToU, pre-peak	17:00	18:00		0.679	0.158	$122$
ToU, on-peak	18:00	21:00		1.059	0.247	$190$
ToU, post-peak	21:00	22:00		0.679	0.158	$122$
ToU, off-peak	22:00	24:00		0.466	0.109	$84$