Economic Management Based on Hybrid MPC for Microgrids: A Brazilian Energy Market Solution

This paper proposes a microgrid central controller (MGCC) solution to the energy management problem of a renewable energy-based microgrid (MG). This MG is a case study from the Brazilian energy market context and, thus, has some operational particularities and rules to be obeyed. The MGCC development was based on a hybrid model predictive control (HMPC) strategy using the mixed logical dynamic (MLD) approach to deal with logical constraints within the HMPC structure, which results in a mixed integer programming (MIP) problem. The development of the solution is done through economic and dynamic modeling of the MG components; furthermore, it also takes into account the energy compensation rules of the Brazilian energy market and the white energy tariff. These conditions are specified through a set of MLD constraints. The effectiveness and performance of the proposed solution are evaluated through high-fidelity numerical simulation.


Introduction
One of the great advantages of some renewable energies is the possibility of generating energy directly in the region where it is consumed, with great emphasis on solar photovoltaic and small wind generation. Thus, this context enables the use of renewable energy in a distributed way, where each consuming unit is able to produce energy for self-sustenance. Considering this new scenario, the concept of the microgrid (MG) appears as a key solution, enabling the guarantee of operational reliability of electrical systems and, at the same time, presenting cost economies to the consumer units [1]. Microgrids are defined as "a cluster of loads and microsources operating as a single controllable system that provides energy and heat to their local area" [2].
The energy management problem of microgrids is carried by units which have been progressively called microgrid central controllers (MGCCs), which are responsible for implementing different control strategies that ensure the adequate energy generation efficiency and also enhances economic concerns. Many different approaches of algorithms for MGCCs, especially the ones based on model predictive control (MPC) [3] and its variants, are available in the literature [4][5][6][7][8][9][10][11][12].
Recently, [4] depicted the main state-of-the-art techniques of MPC applied to energy management in microgrids. In [5], the development of an optimal control for renewable energy microgrids with hybrid energy storage system (ESS) is presented using a hybrid MPC [6] aiming to maximize the economic benefit of the microgrid and to minimize the degradation causes of the storage systems. In other hand, the optimal load sharing of a renewable MG with hybrid ESS through an advanced optimization-based control technique is the research object in [7]. A hierarchical MPC structure acting in different time scales aiming to optimize the economic profit and the electric vehicles charging

Spot Market Rules and Proposed Solution
This section presents the Brazilian spot market rules, as well as the main ideas of the proposed solution and a description of the µGridLab microgrid where the simulation study is developed.

Spot Market Rules
Recently, the Brazilian market rules changed and the so-called "white tariff", which varies according to the time of day and penalizes the cost at peak times, became an interesting option for the users that produce and consume energy, the prosumers. Another important mechanism that has been created is the compensation rule, which allows the units with production capacity to inject energy to the main grid and use it with no purchase costs in other time periods when the load demand is bigger than the energy generation. These conditions allow small consumers, including those operating in low voltage networks, to generate their own energy and dispatch into the main grid autonomously, which creates an energy credit (compensation) that can be used (compensated) in up to 36 months [16].
At this point, it is important for the reader to understand the concept of energy injection and compensation. On one hand, the first concept is, in essence, the excess of energy generated-in other words, the amount of energy that is not instantaneously consumed, but is dispatched/supplied/injected to the main grid. According to the market rules, this energy will result in a credit that can be consumed in the future. On the other hand, the energy compensation concept resides in the usage of this energy credit-for example, in periods where the consumption is greater than the energy production capacity.
We also stress that for one to take benefit out of this energy market scheme, the price of the produced energy should be equal or lower than the price of the credit in each tariff spot. Another characteristic of these rules is that in peak periods, the energy price is bigger than in the other spots, which adds a degree of freedom to produce more energy than needed and inject the excess into the grid in peak periods. Thus, this allows one to generate credits with high values that can be translated in a bigger amount of energy in off-peaks and intermediate spots because of the credit conversion factors.
However, the energy credits are computed by the electricity provider, and to take them into account in the MGCC, it is necessary to model the market rules in the form that they can be dealt in the optimization problem. To solve this issue, the concept of virtual stocks (VS) proposed in this work arises as a powerful tool that can represent with accuracy the market operation, adding better possibilities of decision to the MGCC.
It is important to note here that, as the Brazilian market is different from other markets (e.g., European markets), the Brazilian prosumer cannot have economic benefits by selling energy to the DSO (there is no possibility to sell the energy to the DSO); the only possibility is to convert the injected energy into credits to be used in the future. Another important difference is that energy prices do not oscillate daily and remain constant for weeks or months, thanks to most of the electricity generated coming from large-scale hydroelectric power plants that, in most of the cases, have storage capacities that guarantee the fixed spot prices.
The main contribution of this work is the proposal of a new MGCC that takes the compensation rules in its optimization problem into account through the virtual stocks modeling framework. This concept allows the MGCC take the advantage of deciding the optimal solution for the energy management over the prediction horizon using the rules of energy injection/compensation and the different conversion factors among the three tariff spots. The solution presented hereafter uses the Brazilian energy market as a use case, but it is important to bear in mind that the methodology discussed here can be applied to every energy market with similar rules.

White Tariff
According to [17], in order to use the white tariff, a consumer unit must have the sum of the nominal powers of the transformers equal to or less than 112.5 kVA, which is the case of the µGridLab MG used as a case study in this paper. The white tariff is composed by 3 different spots, off-peak, Energies 2020, 13, 3508 4 of 20 intermediate and peak, and its costs are presented in Table 1. These prices will be used as weights in the main grid cost function, as well as converter factors in the virtual stocks modeling, as explained in following sections.

Compensation Rules
Hereafter, the modeling of the compensation rules and its integration in the MGCC, which is the novel contribution of this paper, will be discussed. For such, the concept of the virtual stocks (VS) is introduced. The control objectives, the dynamic models and the operational constraints are proposed under the following assumptions: 1.
Since there are three spots within the white tariff, three different VSs were proposed, one for each spot.

2.
The energy stored in the VSs has a minimum value of zero and no maximum value constraint.

3.
The energy injected to the grid at each tariff spot is stored in its respective VS. 4.
The energy compensation occurs primarily in the respective tariff spot. In cases where the VS is depleted in its respective tariff spot, it is possible to use the energy of other VS. In this case, a conversion factor γ given by the relation among the tariff spot prices is applied. This modeling approach is compatible with the energy compensation rules [19].

5.
The cost function to be minimized ponders the energy exchanged with the grid-this way, the main objective is to meet the produced energy with the demand as much as possible using the MG renewable sources. 6.
The compensated power is the sum of all powers extracted from the VSs. 7.
When the MG is using the compensated power, in other words, the energy from the VSs, it is not allowed to inject energy in the main grid. 8.
The MG is only allowed to buy energy from the grid if there is no more energy available in the VSs.
The following sections aims to detail how the previous assumptions were implemented in the optimization problem.

Proposed Control Solution
In order to solve the presented problem, a MGCC composed of two hierarchical MPC levels is proposed. The high-level MGCC is responsible for the economic management of the micro-network and it is the focus of this work. It operates on a time scale of hours and determines the operation points for the lower level. The low level is responsible to ensure the energy balance of the MG following the targets defined by the high-level MGCC and it is out of the scope of this work. It operates on a time scale of minutes and guarantees the operation of the MG as close as possible to the economic optimum but allows eventual deviations from this point when necessary. For the design and tuning of the high-level MGCC, correct operation of the low-level MGCC is assumed. A schematic of how this strategy works is shown in Figure 1. The main objective of the high-level MGCC is to operate the MG attending the demands and obtaining the minimum operation cost. Note that the efficient operation of the microgrid has to consider not only direct costs associated to the energy consumption, but also other indirect costs related to the several equipment operations, as will be detailed in the formulation of the cost functions of each one of the MG components in Section 3.
The manipulated (decision) variables of the high-level MGCC are the desired values (setpoints/references) at each sample time of 15 min of the generated power, startup and shutdown of the dispatchable sources, the charge/discharge power of the storage systems, the power extracted/injected from/to the grid, and a set of variables related to the market rules and virtual stocks modeling, as depicted in Section 3. Moreover, some of these variables can be continuous such as the generated power of a dispatchable source or the grid power, as well as binary-for example, the startup or shutdown operations of a dispatchable energy source.
The high-level MGCC is based on a hybrid model predictive control (HMPC) framework, which is an optimization-based control technique where the main goal is to minimize an objective function though a prediction horizon subjected to dynamic models and operational constraints. The term hybrid is used because the optimization problem has real and binary decision variables; in other words, it is a mixed integer optimization problem (MIOP). As in other MPC approaches, a receding horizon strategy is used, and only the first control action of the horizon is applied to the plant, recalculating the optimal solution at the next sample time with the new information available from the process [4]. In a general form, the optimization problem proposed in this work is presented in Equation (1), where the objective function is the sum of dispatchable sources, external grid and storage systems cost functions expressed over the prediction horizon ( ) and subjected to the dynamic models of the storage systems and virtual stocks ( ( , )), to the dispatchable sources and storage systems operational constraints, and to the energy market rules and to the energy balance constraints (respectively, ( , ) and ( , )). In this problem, is the vector of all decision continuous and binary variables related to each equipment operation, is the vector of state variables, in Equations (1) and (2) represents the time and represents the instant in the prediction horizon . The detailed descriptions of all parts of the optimization problem are provided in Section 3.
The general energy balance is presented in the following equation and aims to ensure that all produced energy at each time instant of the prediction horizon will be used by the loads, stored in The main objective of the high-level MGCC is to operate the MG attending the demands and obtaining the minimum operation cost. Note that the efficient operation of the microgrid has to consider not only direct costs associated to the energy consumption, but also other indirect costs related to the several equipment operations, as will be detailed in the formulation of the cost functions of each one of the MG components in Section 3.
The manipulated (decision) variables of the high-level MGCC are the desired values (set-points/references) at each sample time of 15 min of the generated power, startup and shutdown of the dispatchable sources, the charge/discharge power of the storage systems, the power extracted/injected from/to the grid, and a set of variables related to the market rules and virtual stocks modeling, as depicted in Section 3. Moreover, some of these variables can be continuous such as the generated power of a dispatchable source or the grid power, as well as binary-for example, the startup or shutdown operations of a dispatchable energy source.
The high-level MGCC is based on a hybrid model predictive control (HMPC) framework, which is an optimization-based control technique where the main goal is to minimize an objective function though a prediction horizon subjected to dynamic models and operational constraints. The term hybrid is used because the optimization problem has real and binary decision variables; in other words, it is a mixed integer optimization problem (MIOP). As in other MPC approaches, a receding horizon strategy is used, and only the first control action of the horizon is applied to the plant, recalculating the optimal solution at the next sample time with the new information available from the process [4]. In a general form, the optimization problem proposed in this work is presented in Equation (1), where the objective function is the sum of dispatchable sources, external grid and storage systems cost functions expressed over the prediction horizon (N) and subjected to the dynamic models of the storage systems and virtual stocks (F 3 (x, u)), to the dispatchable sources and storage systems operational constraints, and to the energy market rules and to the energy balance constraints (respectively, F 1 (x, u) and F 2 (x, u)). In this problem, u is the vector of all decision continuous and binary variables related to each equipment operation, x is the vector of state variables, k in Equations (1) and (2) represents the time and j represents the instant in the prediction horizon N. The detailed descriptions of all parts of the optimization problem are provided in Section 3.
The general energy balance is presented in the following equation and aims to ensure that all produced energy at each time instant of the prediction horizon will be used by the loads, stored in the Energies 2020, 13, 3508 6 of 20 storage systems or injected in the grid. This is an equality constraint where the sum of the power of all network components must equal to zero: where nl is the number of loads, nnd is the number of non-dispatchable energy sources, ns is the number of storage systems, and nd is the number of dispatchable energy sources. The nomenclature P sub means the power of each MG component expressed by the sub-indices sub.
Remark 1. The studied system has both AC and DC loads, but here, it was considered with a power factor near 1, and only the active power was used in the energy balance. Although this assumption simplifies the equations, it does not conceptually affect the results of the study.

GridLab Microgrid
The microgrid (MG) studied in this work is located in the µGridLab at the Federal University of Santa Catarina (Brazil) and is compound by a battery bank, a gas microturbine emulator, photovoltaic panels, a wind turbine emulator, DC loads and connection with the main grid, as can be seen in Figure 2.
Energies 2020, 13, x FOR PEER REVIEW 6 of 21 the storage systems or injected in the grid. This is an equality constraint where the sum of the power of all network components must equal to zero: where is the number of loads, is the number of non-dispatchable energy sources, is the number of storage systems, and is the number of dispatchable energy sources. The nomenclature means the power of each MG component expressed by the sub-indices .
Remark. The studied system has both AC and DC loads, but here, it was considered with a power factor near 1, and only the active power was used in the energy balance. Although this assumption simplifies the equations, it does not conceptually affect the results of the study.

GridLab Microgrid
The microgrid (MG) studied in this work is located in the μGridLab at the Federal University of Santa Catarina (Brazil) and is compound by a battery bank, a gas microturbine emulator, photovoltaic panels, a wind turbine emulator, DC loads and connection with the main grid, as can be seen in Figure  2. The battery bank consists of 10 lithium-ion battery modules (model Beckett 8224S) [20] with approximately 3000 life cycles. The total capacity of the bank is 10 kWh and the charge and discharge powers are 5 kW and 10 kW, respectively. In order to simulate the behavior of a gas microturbine, the model Capstone C3 [21] was developed, which is an emulator based on power electronics converters with an apparent power of 30 kVA.
The photovoltaic array consists of 10 panels (2 kW each), resulting in a maximum power of 20 kW. The wind turbine is emulated through the coupling of an electric motor and a generator with a capacity of 11 kW. The motor receives a profile to simulate the mechanical torque generated by the wind and the generator transforms the mechanical torque in an electric energy profile. The adjustable DC loads allow the reception of a power reference and can then vary between 0 and 35.5 kWh. The microgrid also has a connection to the main grid, which can be switched on and off, thus allowing operation in grid-connected and island modes. The battery bank consists of 10 lithium-ion battery modules (model Beckett 8224S) [20] with approximately 3000 life cycles. The total capacity of the bank is 10 kWh and the charge and discharge powers are 5 kW and 10 kW, respectively. In order to simulate the behavior of a gas microturbine, the model Capstone C3 [21] was developed, which is an emulator based on power electronics converters with an apparent power of 30 kVA.

Microgrid Economic Modeling and Control
The photovoltaic array consists of 10 panels (2 kW each), resulting in a maximum power of 20 kW. The wind turbine is emulated through the coupling of an electric motor and a generator with a capacity of 11 kW. The motor receives a profile to simulate the mechanical torque generated by the wind and the generator transforms the mechanical torque in an electric energy profile. The adjustable DC loads allow the reception of a power reference and can then vary between 0 and 35.5 kWh. The microgrid also has a connection to the main grid, which can be switched on and off, thus allowing operation in grid-connected and island modes.

Microgrid Economic Modeling and Control
This section presents the economic modeling of each MG component, as well as the inclusion of compensation rules in the optimization problem of the MPC. The modeling is based on the operational costs and it represents (roughly) the real cost of using each component, thus, minimizing the total cost function allows us to obtain the optimal operation point of the microgrid components. Hereafter, each of the component modeling and operational aspects will be depicted.

Energy Balance
The energy balance ensures that all energy produced is used by the loads, stored in the batteries or injected in the grid. This is an equality constraint where the sum of the power of all network components must equal to zero. It is important to note that the solar panels and wind turbines are non-dispatchable energy sources, and, from the control theory point of view, they act as disturbances in the system (MG).
In Equation (3), P i loads (k) is the power consumed by the loads, P solar (k) is the power produced in the photovoltaic panels, P wind (k) is the power generated by the wind turbine, P grid (k) is the power consumed/injected from/to the grid, P bat (k) is the charge/discharge battery power and P turb (k) is the microturbine generated power.

Costs of Operation Points
In order to model the costs of microturbine operation points, the datasheet of the capstone turbine model C30 was used [21]. The data presented in the datasheet can be approximated, with small error, by the following linear model, which represents the cost per hour (Cost turb ) of the microturbine as a function of the generated power (P turb ).
The linear coefficient LC = 1527.1428 l h represents the microturbine fuel hourly consumption without energy generation, while the angular coefficient AC = 311.4286 l kWh represents the fuel hourly consumption rate during the generation process. Finally, NTG is the natural gas tariff NGT = 0.0015719 R$ l , defined according to Resolution No. 098 of the Public Services Regulation Agency of Santa Catarina [22,23], used in Equation (4) in order to express the economic turbine operation cost in Brazilian currency. It is important to take the fact that the coefficients in Equation (4) were obtained for environmental conditions of 15 • C and 1 atm into account [24].

Microturbine Conversion Efficiency
It is well-known that the pressure varies with altitude; however, the MG is located in a coast city with an altitude close to zero and the variation of the microturbine conversion efficiency related to the pressure will be neglected. On the other hand, considerable ambient temperature oscillations directly influence the microturbine efficiency. In order to consider the ambient temperature variations, the weight factor Eff(Temp) is introduced and it varies according to the efficiency vs. temperature curve presented in Figure 3. The curve has been normalized so that Eff(Temp) is 1 at 15 • C. This normalization is given by multiplying the original efficiency curve by (1/0.26), where 0.26 represents the efficiency at 15 • C [21].

Startup and Shutdown Costs
In this work, the microturbine startup � = , ( $)� and shutdown � = , ( $)� total operation costs were considered. To obtain these values in local currency, the maximum startup time (2 min) and shutdown time (10 min) are multiplied by NTG and LC, as during these periods, the microturbine consumes gas without generating energy.

Objective Function
The objective function of the microturbine is a composition of the operating point, startup and shutdown costs. In this work, the sample time of 15 min for the high-level MGCC was used, which is bigger than the startup and shutdown times. Thus, when the high-level MGCC decides to change the operational state of the microturbine, the low-level MGCC is responsible to implement it, and this state is only possible to be changed in the next execution of the high-level MGCC.
Finally, the microturbine objective function is defined as: Equation (5) considers three terms: one for the normal operation, one for the shutdown process and one for the startup process. In Equation (5), ( ) is a binary variable that represents the microturbine on/off state at time instant , ( ) is the power of the turbine, is the sampling period and the variables ( ) and ( ) represent startup and shutdown of the turbine, respectively, being defined by the following constraints: To understand the behavior of Equation (6) Remark. Note that for the startup process, the cost will include the SU cost and the normal operation consumption for one sample. This is not exact, but the introduced error was not significant. This occurs because the sampling time is bigger than the startup period, and both terms depend on the variable ( ). The way

Startup and Shutdown Costs
In this work, the microturbine startup (SUC = 0.037023 (R$)) and shutdown (SDC = 0.185115 (R$)) total operation costs were considered. To obtain these values in local currency, the maximum startup time (2 min) and shutdown time (10 min) are multiplied by NTG and LC, as during these periods, the microturbine consumes gas without generating energy.

Objective Function
The objective function of the microturbine is a composition of the operating point, startup and shutdown costs. In this work, the sample time of 15 min for the high-level MGCC was used, which is bigger than the startup and shutdown times. Thus, when the high-level MGCC decides to change the operational state of the microturbine, the low-level MGCC is responsible to implement it, and this state is only possible to be changed in the next execution of the high-level MGCC.
Finally, the microturbine objective function is defined as: Equation (5) considers three terms: one for the normal operation, one for the shutdown process and one for the startup process. In Equation (5), δ turb (k) is a binary variable that represents the microturbine on/off state at time instant k, P turb (k) is the power of the turbine, T s is the sampling period and the variables SU(k) and SD(k) represent startup and shutdown of the turbine, respectively, being defined by the following constraints: To understand the behavior of Equation (6), let us first assume the startup scenario with δ turb (k − 1) = 0 and δ turb (k) = 1. In this case, we have the first and fourth constraints active, and the startup of the turbine is performed at a cost of SUC. In the opposite scenario, δ turb (k − 1) = 1 and δ turb (k) = 0, the second and third constraints are active and the turbine shutdown is carried out with cost SDC. In other words, minimization of the variables SU(k) and SD(k) implies, in minimization, changes in operational conditions in the turbine weighted by its economic costs. Note that if δ turb (k) = δ turb (k − 1), these costs are not considered in the sample time.

Remark 2.
Note that for the startup process, the cost will include the SU cost and the normal operation consumption for one sample. This is not exact, but the introduced error was not significant. This occurs because the sampling time is bigger than the startup period, and both terms depend on the variable δ turb (k) . The way to avoid this effect it to introduce one more binary decision variable with a set of MLD constraints, but since the error is infimum and startup procedures just happen during a few specific periods of the day, it does not justify the cost to augment the complexity of the optimization problem. Thus, the authors understand that the current solution is the one that provides the better compromise between the modeling and computational complexity and desired operation of the MG.
The microturbine power limits are given by: where M turb is the maximum power value (30kVA) and m turb is the minimum power value (0kVA). When δ turb (k) = 1, the power stays among the minimal and maximum limits, while in the case of δ turb (k) = 0, the turbine power should be zero.

Virtual Stocks Dynamic Model
In order to represent the behavior of the virtual stock, the following models were used: E I (k + 1) = E I (k) + γ 4 T s P OP I (k) + γ 5 T s P I I (k) + γ 6 T s P P I (k) + T s P injected grid I (k) E P (k + 1) = E P (k) + γ 7 T s P OP P (k) + γ 8 T s P I P (k) + γ 9 T s P P P (k) + T s P injected grid P (k) (8) where E sub1 (k) represents the energy of each stock related to the three white tariff spots. The subscript sub1 was used to generalize the problem since the same logic will apply to tariff spots periods OP-off-peak; I-intermediate, P-peak). P sub2 sub3 (k) represents the power consumed from the virtual stock with sub-index sub3 during the time of spot with sub-index sub2, taking into account that sub2 sub3. The power injected to the grid at each time spot is expressed by P injected grid sub1 (k). The constants γ n are the conversion factors given by the relation between the cost (C sub1 ) of the tariff spots: A simplified scheme of the virtual stocks operation is shown in Figure 4.

Minimum Value Constraints
To ensure that stocks will assume only positive values, the following constraints were added to the problem: Energies 2020, 13, 3508 10 of 20 Energies 2020, 13, x FOR PEER REVIEW 10 of 21

Minimum Value Constraints
To ensure that stocks will assume only positive values, the following constraints were added to the problem:

Usage of the VS Correspondent to the Actual Spot Period
The virtual stocks, in accordance with the energy market rules and to guarantee the compensation of the injected energy at the current tariff spot, must be subject to the following constraints according to the tariff spot: Off-peak: Intermediate: Peak: To understand the operation of these constraints, let us take a look at the scenario where the MGCC is operating in the off-peak period. In this case, if the MGCC decide to consume the energy from the off-peak VS, it will be done by the variable ( ); if it is needed to use energy from the intermediate and peak VSs, the MGCC will use the variables ( ) and ( ), respectively. During this spot period, the constraint imposed by Equation (11) should work in order to guarantee that the operation will be done as expected. It is important to ensure that the right decision variable will be used, since each of them is related to the conversion factor that represents the energy flow among the VSs. The same reasoning holds for the other two spots periods. This modeling framework is directly connected with assumption 4 of Section 2.1.2.

Priority Among the VSs
To represent the existing priority in energy compensation, given by the consumption of the virtual stock relative to the current spot period, the use of the following logical connectives is required:

Usage of the VS Correspondent to the Actual Spot Period
The virtual stocks, in accordance with the energy market rules and to guarantee the compensation of the injected energy at the current tariff spot, must be subject to the following constraints according to the tariff spot: Off-peak : P I OP = P I I = P I P = P P OP = P P I = P P P = 0 Intermediate : P OP FP = P OP I = P FO = P P OP = P P I = P P P = 0 Peak : P OP OP = P OP I = P OP P = P I OP = P I I = P I P = 0 To understand the operation of these constraints, let us take a look at the scenario where the MGCC is operating in the off-peak period. In this case, if the MGCC decide to consume the energy from the off-peak VS, it will be done by the variable P OP OP (k); if it is needed to use energy from the intermediate and peak VSs, the MGCC will use the variables P OP I (k) and P OP P (k), respectively. During this spot period, the constraint imposed by Equation (11) should work in order to guarantee that the operation will be done as expected. It is important to ensure that the right decision variable will be used, since each of them is related to the conversion factor γ n that represents the energy flow among the VSs. The same reasoning holds for the other two spots periods. This modeling framework is directly connected with assumption 4 of Section 2.1.2.

Priority among the VSs
To represent the existing priority in energy compensation, given by the consumption of the virtual stock relative to the current spot period, the use of the following logical connectives is required: Off-peak : E OP (k) > 0 → P OP I (k) = 0 and P OP P (k) = 0 Intermediate : E I (k) > 0 → P I OP (k) = 0 and P I P (k) = 0 Peak : E P (k) > 0 → P P OP (k) = 0 and P P I (k) = 0 This existing logical priority defines which storage will be used and is represented by the multiplexers showed in Figure 4. These constraints only allow using the VS referred to other spots if the VS of the current spot is out of energy. The presented conditions are related to assumption 4 of Section 2.1.2.
In order to enable the use of logical propositions by the MPC, the MLD framework is used [25]. In addition, since some constraints must have different behaviors according to the current tariff spot, the vector of binary variables defined in Equation (17) will be used, which will be decided out of the optimization problem and sent as a parameter. These vectors, which have a size equal to the prediction horizon, are conveniently used in some constraints that will be presented later. If at time instant k, the MGCC is operating in the spot period OP, for example, vector OP (k) = 1, otherwise vector OP (k) = 0.
An auxiliary binary variable δ sub1 (k) is used to force E sub1 (k) to zero when the condition is satisfied. Thus, it is desired that this variable is equal to zero if and only if E sub1 (k) is greater than zero: To model this logical condition in the form of a constraint, the constants M e and m e that represent the maximum and minimum values of E sub1 (k), respectively, are used. It should be emphasized that in this particular case, m e = 0 since the value of the stored energy credit may not be negative and M e ∞ because there is no ceiling for the value of the energy credit in the compensation system. In order to ensure the condition of equivalence, it must be verified that the logical implication is taken care of in both directions.
To guarantee that E sub1 (k) > 0 → δ sub1 (k) = 0 , the following inequality is considered: where ε is a constant of positive value close to zero. Note that this constraint meets this implication by forcing δ sub1 (k) to zero when E sub1 (k) assumes positive values. There is ambiguity when E sub1 (k) = ε, but it does not interfere with the desired implication, since E sub1 (k) − ε ≤ M e is still valid.
To ensure that δ sub1 (k) = 0 → E sub1 (k) > 0 , the following inequality is considered: (20) Note that this constraint meets this implication by forcing E sub1 (k) to assume positive values when δ sub1 (k) is equal to zero.

s P OP OP (k) + γ 2 T s z I OP (k) + γ 3 T s z P OP (k) + T s P injected grid OP
(k) E I (k + 1) = E I (k) + γ 4 T s z OP I (k) + γ 5 T s P I I (k) + γ 6 T s z P I (k) + T s P injected grid I (k) E P (k + 1) = E P (k) + γ 7 T s z OP P (k) + γ 8 T s z I P (k) + γ 9 T s P P P (k) + T s P injected grid P (k) (21) It is important to note that both the δ sub1 (k), P sub2 sub3 (k) and z sub2 sub3 (k), are decision variables, so that this multiplication between variables will generate a bilinearity in the objective function. To solve this problem, we adopt the approach introduced in [21], where the auxiliary variable z sub2 sub3 (k), subject to some constraints, is used to represent the bilinearity behavior, as follows: δ sub1 (k)).vector sub1 (k) (22) where M stock and m stock represent, respectively, the minimum and maximum values of P sub2 sub3 (k). When δ sub1 (k) = 0, we have z sub2 sub3 (k) = 0 and P sub2 sub3 (k) = 0, which represent that, at spot period sub2, there is no consumption from the VS sub3. The case of δ sub1 (k) = 1 implies that z sub2 sub3 (k) = P sub2 sub3 (k) and m stock ≤ P sub2 sub3 (k) ≤ M stock , which means that, at spot period sub3, consumption from the VS sub2 exists. At this point, it is important to point out that the possible combinations of the values of variables sub1, sub2 and sub3 are depicted in Table 2.

Cost Function
The portion of the cost function relative to the grid can be expressed by: where P grid (k) represents the grid power and C sub1 represents the cost of electricity at the time that depends on the current tariff spot presented in Table 1. The positive values of P grid (k) imply purchasing energy from the grid, while negative values imply injecting power to the grid. To relate the injected power P injected grid sub1 (k) with the grid power P grid (k), a set of MLD constraints is defined as follows: where δ aux (k) is a binary auxiliary variable; z aux (k) is a continuous auxiliary variable; and M grid and m grid represent, respectively, the minimum and maximum values of P grid (k). When δ aux (k) = 0, we have z aux (k) = 0 and P injected grid sub1 (k) = 0, which represent no power injection in the grid, but it allows consumption of the power from the grid. The case of δ aux (k) = 1 implies that z aux (k) = P grid (k) and P injected grid sub1 (k) = −z aux (k), which means that a power injection in the grid is happening.
Energies 2020, 13, 3508 13 of 20 In order for the energy injection not be accounted in all stocks simultaneously, the tariff spot vectors were used. By using the tariff spot vectors, only the VS referent to the current spot period will receive the injected power. Finally, Equation (21) must be rewritten as: (k)vector OP (k) E I (k + 1) = E I (k) + γ 4 T s z OP I (k) + γ 5 T s P I I (k) + γ 6 T s z P I (k) + T s P injected grid I (k)vector I (k) E P (k + 1) = E P (k) + γ 7 T s z OP P (k) + γ 8 T s z I P (k) + γ 9 T s P P P (k) + T s P injected grid P (k)vector P (k)

Compensated Power and Variables Interlocking
The variable related to total compensated power P comp (k), which represents the total power used from the VSs, does not appear in the cost function since its cost is zero. However, it appears in the constraints as a decision variable. P comp (k) = P OP OP (k) + z OP I (k) + z OP P (k) + P I I (k) + z I OP (k) +z I P (k) + P P P (k) + z P OP (k) + z P I (k) Furthermore, to ensure the interlocking between energy compensation and injection, that is, P comp (k) ≥ 0 → P injected grid sub1 (k) = 0 and P comp (k) = 0 → P injected grid sub1 (k) ≥ 0 , the following constraints were added: where δ comp (k) is a binary variable and M grid is the maximum value of P grid (k). When δ comp (k) = 0, there is no compensation from the VS, but power injection is allowed; however, in the case δ comp (k) = 1, the opposite statement is true.
In the same way, to ensure the interlocking between energy compensation and energy purchase, that is, E OP (k) + E I (k) + E P (k) > 0 → P grid (k) ≤ 0 , the following constraints are added: Here, when δ e (k) = 0, there is no compensation and the purchase from the grid is allowed; otherwise, the opposite situation rules.

Objective Function
To estimate the cost of using the batteries, two important aspects were considered: the loss of energy in the conversion and the degradation of the battery. The cost of using the battery considering purchase cost, efficiency and total number of cycles can be modeled as: where CB total is the purchasing cost of the battery, N cicles is the number of complete cycles, η C and η D are the battery charge and discharge efficiencies, P bat (k) is the battery power at instant k, and T s the sampling period. Moreover, in order to allow the use of different values for the charge/discharge efficiencies, the auxiliary variables z bat (k) (continuous) and δ bat (k) (binary) were introduced, so that z bat (k) = δ bat (k)P bat (k). The MLD constraints that guarantee the desired operation for the battery modes switching are given by: where M bat (positive value) and m bat (negative number) are, respectively, the maximal and minimal values allowed for P bat . Note that z bat (k) and δ bat (k) allow us to switch the cost function according the battery charge/discharge modes. The variable δ bat (k) is associated with the charge (δ bat (k) = 1) or discharge (δ bat (k) = 0). Note that when P bat (k) > 0, we have δ bat (k) = 1 and z bat (k) = P bat (k), which means that the battery is charging, therefore the charge efficiency is used. Otherwise, P bat (k) < 0 implies that δ bat (k) = 0 and z bat (k) = 0, and the discharge efficiency is used.
The efficiency of a battery depends on the charge or discharge speed. In this work, we chose to use near-normal efficiencies for lithium-ion batteries [26,27], then considered a charge efficiency of 92% and discharge efficiency of 95%.

Constraints
The optimization of the batteries cost must be subject to the dynamics of the behavior of the batteries, which is expressed by constraints. In order to model the dynamics of the battery, the MLD framework was used, since it allows representing logical variables in a suitable way to represent the restrictions in MPC. The modeling was based on the one presented by [4]. The battery bank state of charge is represented by: where SOC(k) is the state of charge at time k and Cap is the battery capacity.

Simulation Results
The simulations were developed in MATLAB software [28], where the controllers were described using the YALMIP modeling language [29] and the solver Gurobi [30] was used. Real data referring to the year 2019 were used to generate the vectors of solar incidence, temperature and wind speed. The sampling time used for the MPC was 15 min and the prediction horizon adopted was 24 samples (equivalent to 6 h). The control horizon is equal to the prediction horizon since the optimization problem does not have severe temporal and processing constraints, considering the computation time of the solution (less than 1 s).
In this section, three different simulation scenarios are evaluated. The first scenario seeks to represent the nominal operating conditions of the µGridLab microgrid, while the second scenario analyses the results if there is a 50% increase in expected demand. Finally, in the third scenario, an island operation of the MG in nominal conditions is studied and compared to scenario 1.

Scenario 1
In this scenario, a one-day period simulation is performed to observe the results graphically and a complete month simulation (period of 30 days) is simulated with the intention of performing the monthly economic analysis. Figures 5-8 show all the interesting variables during one day of operation. From Figure 5, where the generated and demanded power (kW) are shown, it can be observed that in many periods, the MGCC chooses to use the microturbine to generate power. This is due to the low cost of the microturbine use (compared to the grid power exposed in Figure 7b), thus, the controller prioritizes the energy generation through the microturbine instead of consuming power from the network. The microgrid also uses all generated power by the wind turbine and the PV (photovoltaic) panels. It is important to point out that in some periods of the day, the load demand is bigger than the generated power, but in these moments, as can be seen in Figure 7b, the compensated power from grid helps to feed the loads.

Scenario 1
In this scenario, a one-day period simulation is performed to observe the results graphically and a complete month simulation (period of 30 days) is simulated with the intention of performing the monthly economic analysis. Figures 5-8 show all the interesting variables during one day of operation. From Figure 5, where the generated and demanded power (kW) are shown, it can be observed that in many periods, the MGCC chooses to use the microturbine to generate power. This is due to the low cost of the microturbine use (compared to the grid power exposed in Figure 7b), thus, the controller prioritizes the energy generation through the microturbine instead of consuming power from the network. The microgrid also uses all generated power by the wind turbine and the PV (photovoltaic) panels. It is important to point out that in some periods of the day, the load demand is bigger than the generated power, but in these moments, as can be seen in Figure 7b, the compensated power from grid helps to feed the loads.  Figure 6 shows the SOC (state of charge) and power in the batteries. As can be seen, the battery power presents some oscillations near the operating point, as the battery plays the role of aiding the system in damping the effect of the renewable energy variations and load changes so that the energy balance is guaranteed and the variation of the energy exchanged with the external grid is minimized.
It is important to note that in the period between approximately 12 and 19 h, these oscillations are large, and are in accordance with the ones observed in the solar generation due to the strong presence of clouds, which causes abrupt variations of irradiation. On sunny days, it is not possible to observe this phenomenon. The constraints imposed in the controller ensure that the SOC is between 100% and 20% as recommended for the correct operation of the batteries.  Figure 7a shows the energy of the virtual stocks for the three tariff spots. Initially, the stocks have zero energy, so, it is assumed that in this scenario, there is no previous energy credit from previous months. This would be the case for the MG first month operation in the compensation system. It should be noted that in this scenario, there are several moments that happens the injection of energy for subsequent use. The energy injection at peak hours is the one that generates more offset  With respect to economic costs, the one-day period of the obtained operating data is presented in the second column of Table 3. It is noted that at the end of the day, there was an excess balance of compensated energy. The total operating cost of the microgrid during this day was R$ 221.84, but a credit of R$ 18.78 was generated for later use; this credit is calculated as the sum of the products between the power stored in each virtual store by its respective tariff (see Table 1).
In order to analyze the monthly economic operation of this scenario, a simulation is performed for a period of 30 days under the same conditions, and the data referring to the first 30 days of the month of January 2019 were used. In this case, the operation data are presented in the third column of Table 3.  Figure 6 shows the SOC (state of charge) and power in the batteries. As can be seen, the battery power presents some oscillations near the operating point, as the battery plays the role of aiding the system in damping the effect of the renewable energy variations and load changes so that the energy balance is guaranteed and the variation of the energy exchanged with the external grid is minimized. It is important to note that in the period between approximately 12 and 19 h, these oscillations are large, and are in accordance with the ones observed in the solar generation due to the strong presence of clouds, which causes abrupt variations of irradiation. On sunny days, it is not possible to observe this phenomenon. The constraints imposed in the controller ensure that the SOC is between 100% and 20% as recommended for the correct operation of the batteries. Energies 2020, 13, x FOR PEER REVIEW 18 of 21  Note that in this period, there was more energy injection than compensation; in addition, the MGCC has never decided to purchase energy once the compensated energy has zero cost. At the end of 30 days, there was R$ 19.57 of credit to be compensated in the upcoming months. The total cost of operation during the month considering the use of the battery and the turbine is R$ 8116.78.
For comparison purposes, the energy cost is calculated if all demand is supplied by the grid. In this case, the total operation cost is obtained by applying the cost of the conventional tariff to the energy demand (R$ 8755.27). This gives a 7.86% increasing in the energy bill with respect the distributed generators plus white tariff compensation system. Of course, the comparison with the use of the grid power only serves to give the reader a better perception of the values, since a rigorous economic analysis would have to take into account the purchase cost of the turbine, the wind generator and the PV panels. This analysis is disregarded in this work since the studied MG already has these components.

Scenario 2
As a second scenario to be analyzed, a 50% increase in energy demand is proposed. The curves obtained for this case have profiles similar to those obtained in the first scenario-this way, they will not be presented. In this case, simulating again the periods of 1 (one) day and 30 days give the operation data shown in Table 4.   Figure 7a shows the energy of the virtual stocks for the three tariff spots. Initially, the stocks have zero energy, so, it is assumed that in this scenario, there is no previous energy credit from previous months. This would be the case for the MG first month operation in the compensation system. It should be noted that in this scenario, there are several moments that happens the injection of energy for subsequent use. The energy injection at peak hours is the one that generates more offset credits.
The relationship between the compensated power, the grid power and the injected power in this scenario is shown in Figure 7b. It can be seen that the controller does not choose to purchase power from the grid, but only compensates for credits obtained, since grid power is always negative, which means that the power is being injected into the grid. This fact can be explained by the high price of energy in relation to the energy generated by microturbine and renewable sources. Figure 8 shows the value of the objective function during the first day of the scenario 1. As expected, the minimum values of the objective function occur in the period between 11 h and 19 h, where the microturbine is not operating most of the time and there is a peak in the photovoltaic generation as well as in the wind turbine.
With respect to economic costs, the one-day period of the obtained operating data is presented in the second column of Table 3. It is noted that at the end of the day, there was an excess balance of compensated energy. The total operating cost of the microgrid during this day was R$ 221.84, but a credit of R$ 18.78 was generated for later use; this credit is calculated as the sum of the products between the power stored in each virtual store by its respective tariff (see Table 1). In order to analyze the monthly economic operation of this scenario, a simulation is performed for a period of 30 days under the same conditions, and the data referring to the first 30 days of the month of January 2019 were used. In this case, the operation data are presented in the third column of Table 3.
Note that in this period, there was more energy injection than compensation; in addition, the MGCC has never decided to purchase energy once the compensated energy has zero cost. At the end of 30 days, there was R$ 19.57 of credit to be compensated in the upcoming months. The total cost of operation during the month considering the use of the battery and the turbine is R$ 8116.78.
For comparison purposes, the energy cost is calculated if all demand is supplied by the grid. In this case, the total operation cost is obtained by applying the cost of the conventional tariff to the energy demand (R$ 8755.27). This gives a 7.86% increasing in the energy bill with respect the distributed generators plus white tariff compensation system. Of course, the comparison with the use of the grid power only serves to give the reader a better perception of the values, since a rigorous economic analysis would have to take into account the purchase cost of the turbine, the wind generator and the PV panels. This analysis is disregarded in this work since the studied MG already has these components.

Scenario 2
As a second scenario to be analyzed, a 50% increase in energy demand is proposed. The curves obtained for this case have profiles similar to those obtained in the first scenario-this way, they will not be presented. In this case, simulating again the periods of 1 (one) day and 30 days give the operation data shown in Table 4. In fact, the MGCC behaves in an expected way as well in this more severe scenario, using all the power of the compensation system and virtual stocks in order to take advantage of the conversion factors among the spots and achieve a good economic efficiency. For comparative purposes, the equivalent cost of purchasing grid power for the 30 days' case with a conventional tariff is R$ 12845.85. Thus, in this case, the proposed control strategy gives a 3.72% reduction in the energy bill.

Scenario 3
Finally, in order to validate the controller and the sizing of the microgrid, a third scenario is simulated with the MG operating in islanded mode, which means without the connection to the power grid. The 1 (one) day and 30 days' results are given in Table 5 for the same loads as in scenario 1. Here, it is possible to compare scenarios 1 and 3, and as expected, the MG is capable of operating in both modes. Analyzing Tables 3 and 5, is it possible to conclude that the controller works more efficiently in the connected mode with a 4.76% reduction in the energy bill compared with the islanded mode. This happens because there is an extra economic element, the virtual stocks and compensation system, allowing more freedom in the optimization. The cost of battery use was lower in the island mode because the MGCC opted for the use of the turbine to assist in damping short-term oscillations.

Conclusions
This work proposed a MGCC strategy for simulation and control of microgrids considering the Brazilian energy compensation system. The Brazilian energy compensation rules are different from the ones used in other countries; thus, it is necessary to develop specific control strategies to take advantages of the particular processes and reduce the electricity bill. Therefore, the proposed MGCC considers in its formulation all the particularities of the local market and allows the user to take more profit of the installed EMS. The considered microgrids include renewable energy sources and conventional ones. The proposed strategy consisted of a MPC using mixed-integer programming, being the main contribution and innovation, the modeling and implementation of a set of virtual stocks that allow one to optimally manage the power consumption considering the Brazilian energy compensation rules. Through the analysis of the results, it is possible to conclude that the behavior of the proposed MGCC is in agreement with the expected one, taking management decisions that allow the operation of the microgrid in a more efficient way, reducing the energy bill in the studied scenarios. Although the obtained results were obtained in a particular microgrid, the methodology is general and can be applied to other microgrids and could also be simply adapted to other energy compensation systems. Regarding applicability, for companies selling EMS systems, the use of an advanced control system such as the one proposed here gives a competitive edge without increasing costs too much. On the other hand, users can reduce their energy bills and have faster payback. As for future works, the authors plan on implementation of a demand management layer and stochastic energy management techniques can be taken into account.

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