Model Predictive Control for Photovoltaic Plants with Non-Ideal Energy Storage Using Mixed Integer Linear Programming

Abstract

This paper proposes a model-based predictive control strategy based on mixed-integer linear programming for a photovoltaic power plant with battery energy storage. The control objective is to maximize the revenues from energy delivered from both photovoltaic panels and batteries to the grid in a deregulated electricity market. For each control interval, the proposed algorithm incorporates information on solar radiation, market prices, and the state of charge of the batteries to determine the intervals of energy injection into the grid to maximize the economic benefits. The proposed strategy considers the rate-based variable efficiency in the battery model and time-varying energies prices, thus providing a more general implementation than previous schemes proposed in the literature for the same purpose. Simulations considering the operational procedures of the Spanish market as a case study show that, by integrating the battery efficiency in the model, the proposed control strategy increments the economic benefits in 21% compared to previous results reported in the literature for the same operational conditions. Additionally, the proposed approach reduces the number of charge and discharge cycles, potentially extending the lifespan of batteries.

1. Introduction

Electricity generated from large-scale photovoltaic (PV) plants currently represents a significant part of the total generation in many countries. The proliferation of PV plants has been particularly prominent in regions with high solar irradiance, such as eastern Turkey [1] and northern Chile [2]; however, technological advances such as solar trackers for PV panels can improve the performance of PV systems and help on expanding PV generation to regions with less favorable irradiance conditions [3]. From an economic perspective, PV plants have become attractive not only in energy markets with feed-in tariffs such as Iran [4], but also in deregulated energy markets such as Spain [5,6].

In the deregulated Spanish market, electricity is traded in regular settlement periods, during which participating energy producers commit to delivering a certain amount of energy for future periods in hour blocks. The energy commitment price for each block is settled as pay-as-cleared before the actual time of delivery. At the time of delivery, i.e., in real-time, the power produced by a given generator may deviate from the plan considered at the moment of commitment due to unforeseen contingencies or variations with respect to the forecast of energy production considered during the settlement period. Producers failing to deliver the committed energy may generate a net imbalance in the system. System regulators must cover this imbalance by buying energy available from other generators participating in the market [7]. Energy bought to cover an imbalance is paid at an imbalance price, which is higher than the settled price agreed to in advance for the same block, and it is charged to the market participants that are generating the imbalance via a penalty. Therefore, system imbalances are economically detrimental for producers unable to fulfill their commitment and beneficial for producers with surplus energy available for selling at the time of the imbalance [8].

The variable nature of solar irradiation during a day prevents PV plants from providing a constant power supply [9]. Therefore, it is difficult to guarantee a constant energy delivery for a block of several hours without the risk of breaching this commitment. A typical solution for this issue is integrating battery energy storage systems (BESS), which enable energy accumulation during periods in which real-time production is higher than the commitment. The battery storage, can be of different technologies and be integrated into the PV system either in DC or AC side, with AC coupling being easier to use in retrofit projects and providing higher flexibility during operation [10]. On the other hand, lithium-ion batteries are considered more suitable for market-oriented applications [11,12]. Other applications include the suppression of power fluctuations according to the grid codes’ requirements [13], frequency regulation [14] or integration of different renewable energy sources [15,16].

The energy accumulated in the batteries can be, in turn, delivered to the system during periods in which the real-time energy production from the PV panels is lower than the committed energy. Moreover, producers can sell energy stored in the battery at a higher imbalance price. Considering that the price of committed energy, imbalance price, and the penalties vary for each settlement period, the energy stored in the batteries will have a dynamic relative value, and then the profits obtained by PV energy producers will depend on the adequate choice of intervals to charge or discharge the batteries. Therefore, to maximize the economic benefits obtained from PV plants, producers must implement optimized control strategies to manage the state of charge (SoC) in real-time, considering forecasts of both PV production and energy prices in the corresponding electricity market [14,17].

Multiple authors have tackled the problem of maximizing profits from photovoltaic power plant with battery energy storage system (PV+BESS) in deregulated electricity markets using Model Predictive Control (MPC) [6,16,18]. MPC strategies use prediction models of the system and operational variables to decide the times to charge or discharge the batteries considering estimations of future energy generation and market prices. To this end, MPC schemes must optimize a cost function subject to operational constraints at each sampling interval [19,20]. In our target application, battery efficiency and sometimes market prices are modeled as piece-wise functions, which naturally lead to a Mixed-Integer Linear Programming (MILP) formulation for the optimization problem.

In [6], the authors formulate an optimization problem to maximize profits from a PV+BESS plant in the Spanish market. In this case, the cost function depends on the imbalance price and penalty, which have a piece-wise representation that leads to an MILP formulation. Authors argue that solving an MILP problem in each control interval is computationally prohibitive for real-time operation, and thus they apply ad hoc heuristics to derive an equivalent LP formulation, showing the effectiveness of the control approach using simulations that consider a sampling time of four minutes. Although effective for incrementing profits, authors only consider a simplified model for the battery that omits the piece-wise function that represents the charging and discharging efficiency, which we argue is a relevant aspect for accurate decision-making to maximize profits from energy sales. The work in [18] proposes an MPC scheme for a PV+BESS system operating in the PJM Energy Market (www.pjm.com, accessed on 30 August 2022) with a price scheme that results in a linear objective function. This work does consider the BESS efficiency in the optimization model; however, the authors propose a configuration that uses two storage devices with independent charging signals generated through power spectral density (PSD) decomposition, which allows for simultaneous charging and discharging of the different storage devices. Hence, the proposed configuration avoids needing a binary variable and MILP formulation. The resulting optimization problem is solved in two steps, using a standard LP solver, considering the operational constraints, to generate references that are followed using MPC with a quadratic cost function and a standard Quadratic Programming (QP) solver. From the optimization perspective, the problem is simpler as no binary variables arise, but it requires a complex and unconventional BEES with multiple storage devices. In [16], the authors derive an MILP formulation to manage the BESS in a hybrid plant containing both photovoltaic and wind power generators. This work optimizes profits using a model that integrates the income from energy sales and the cost associated with under or oversupply, battery efficiency, and forecasting models for energy generation. This case considers fixed ratio between positive and negative imbalance costs; moreover, the complexity of the plant and the aggregation of forecasts in the optimization problem do not allow quantifying the benefit associated with modeling the battery’s efficiency.

In a related application, multiple works have used MILP solvers to optimize microgrids’ design, operational costs, and energy consumption; however, these problems consider simpler tariff schemes than our target scenario. The work in [21] proposes an MPC strategy that uses MILP to optimize the operational efficiency of a microgrid comprising PV panels, wind turbines, a diesel generator, and BESS. Results using data from an experimental plant in Greece show that the control strategy optimizes the use of diesel generators to reduce operational costs. The work in [22] tackles a similar problem, but using a more sophisticated optimization approach that considers MILP and particle swarm optimization. In [23,24], the objective is to design a control algorithm to satisfy the electrical demand and maximize self-consumption. In [25], the focus is on optimal asset (i.e., PV panels and BESS) planning to minimize the total associated energy-consumption cost for prosumers along the project’s horizon.

This paper presents a novel MPC scheme with an MILP formulation to maximize revenues from a PV+BESS generation plant participating in a deregulated electricity market. The proposed approach enables integrating piece-wise functions in the cost functional and restrictions of the optimization problem, facilitating the direct incorporation of a rate-based variable efficiency model of the battery efficiency and tariff schemes represented with integer variables. Using simulations that consider the Spanish market rules described in [6] as a case study, we show that the proposed MILP formulation that accounts for the efficiency of the BESS allows the MPC scheme to make better decisions about the time intervals to charge or discharge the batteries. In the simulated scenario, the proposed method increments the revenues by 21% with respect to the formulation in [6] that assumes ideal batteries. Evaluations of the execution time for the proposed MILP-based MPC scheme using the CPLEX solver report a worst-case execution time of around 45 milliseconds for a control loop, which is orders of magnitude below the typical required periods in the target application [6,18].

The results provide evidence of the effectiveness of MILP-based MPC schemes in optimizing the real-time operation of PV plants participating in electricity markets without requiring ad hoc simplifications [6] or unconventional BESS infrastructure [18]. Moreover, the reported execution times indicate that the approach can scale to incorporate more integer variables to accommodate more complex objective functions and different market rules.

The rest of the paper is organized as follows: Section 2 presents the MPC strategy for a PV+BESS plant operating in the Spanish market and its solution using an LP formulation with an ideal battery as described in [6]; Section 3 presents the proposed MILP formulation for the same problem that allows incorporating the rate-based efficiency of the battery in the functional cost; Section 4 compares the projected economic benefits obtained with the existing and the proposed approach using a reference scenario; finally, Section 5 and Section 6 discuss the main results and conclude the paper.

2. Materials and Methods

This section provides an overview of the main concepts related to MPC and its use in the optimization of economic gains from PV+BESS plants participating in electricity markets. In specific, the following paragraphs formalize the problem of maximizing revenues from PV+BESS plants following the assumptions presented in [6], which uses the regulations of the Spanish electricity market as a case study. The main theoretical contribution of [6] that allows presenting the associated optimization problem as an LP problem is described in detail. This formulation provides the baseline for validating and comparing the advantages of using the novel approach proposed in this paper.

2.1. Overview of Model Predictive Control

MPC is an advanced control strategy that, using a model of the target system and the (measured or estimated) value of the system’s states at a given sampling time t, is able to determine the optimal value of the control inputs to make the system’s output follow a reference along a number of N future samples while respecting operational constraints in the states and outputs.

The basic operation of MPC for a linear system can be summarized as follows, see, e.g., [19,20]:

• At each sampling time t, the current output is determined, and a dynamic model of the target system is used to forecast the future outputs along a prediction horizon of N sampling intervals.

• The forecast of system behavior at current interval t is used to calculate the optimal control signals that optimize a cost function for the entire prediction horizon. The cost function also imposes constraints on system states and outputs.

• An optimization algorithm is used to determine the sequence of control input signals that optimize the cost function for the prediction horizon; however, only the first element of the optimal control sequence is applied to the process. This approach is called Receding horizon.

• At the next sampling interval $t+1$, a new forecast is calculated, and the process is repeated.

The cost function to be optimized at each iteration of the MPC algorithm is specific for each problem and must consider the model of the controlled system and constraints.

2.2. Problem Specification

Plants participating in deregulated electricity markets must regularly commit a constant-by-hours power production some hours before the real-time delivery instant. The commitment is performed when the actual PV production is still unknown and can only be approximated with a certain degree of accuracy. In this context, when the delivery time arrives, it might happen that it is not possible to deliver the committed power if the actual PV generation is insufficient since low solar irradiation or the BESS are completely discharged. Therefore, it is convenient to have adequate decisions to manage the energy in the BESS to supply the committed power.

The objective of our application is to define a control strategy to manage the generation of a grid-tied PV power plant that takes part in the electricity market. Using short-term predictions for a horizon of N time intervals starting from current time t, we expect the control algorithm will indicate the actions of injecting energy from the battery to the grid, storing energy in the BESS, or leaving the battery idle, considering the optimization of economic return.

The problem of interest is represented in the diagram of Figure 1. Here, and in the rest of this document, the subscripts $t+k$, with $k=0,\dots ,N$, will symbolize a set of values for a given variable that expands $N+1$ samples starting from the current time interval t. In the diagram, $P^{\mathrm{grid}}$ is the total power delivered to the grid, corresponding to the sum of the power generated from the PV plant ($P^{\mathrm{PV}}$) and the power from the battery system ($P^{\mathrm{ES}}$). In the system, a positive value of $P^{\mathrm{ES}}$ indicates a power flows from the BESS into the connection node and, therefore, into the grid. Conversely, a negative value of $P^{\mathrm{ES}}$ indicates a power flow into the BESS. The value and direction of $P^{\mathrm{ES}}$ depend on the decision to either store energy or inject power into the grid ($P^{{\mathrm{ES}}^{*}}$) reached by the MPC algorithm. For simplicity, we consider $P^{{\mathrm{ES}}^{*}}=P^{\mathrm{ES}}$.

The MPC block uses the information of the commitment of constant hourly power to be injected into the grid $P^{\mathrm{ref}}$, the prediction of the energy production from the PV plant ${\widehat{P}}^{\mathrm{PV}}$, the amount of energy stored in the battery $E^{\mathrm{ES}}$, and energy prices. The variable $P^{\mathrm{ref}}$ and the prices are known in advance and $E^{\mathrm{ES}}$ can be measured. The prediction ${\widehat{P}}^{\mathrm{PV}}$ can be obtained from statistical information about the behavior of the solar radiation and PV system performance. Particularly, in Europa and Asia, a proper tool is the PhotoVoltaic Geographic Information System (PVGIS) [26]. The energy prices incorporate information about the market structure and the scheduling of settlement periods, which depend on the specific regulations for the given target market.

In the case of the Spanish electricity market, there is a daily market (one day ahead) where an amount of delivered power committed $P^{\mathrm{ref}}$ paid at a market price $c^{\mathrm{m}}$ is agreed upon, which can be adjusted in an intraday market divided into six sessions of three, four, and five hours to consider variations in the demand forecast and supply contingencies. On the one hand, the intraday market defines the commitment of $P^{\mathrm{ref}}$ and $c^{\mathrm{m}}$ for the following session. On the other hand, imbalances due to excess or shortage of generation are assigned two different prices around $c^{\mathrm{m}}$ that are fixed by the system operator and apply to the difference between the scheduled and actually generated energy. When there is a deviation in the sense of an excess of injected power, the price assigned to this imbalance is $c^{\mathrm{p}}$, with $c^{\mathrm{p}}\le c^{\mathrm{m}}$. On the other hand, when there is a negative difference between the injected power and the committed power, i.e., the producer fails to fulfill his contract, a penalty price $c^{\mathrm{n}}$ is applied, with $c^{\mathrm{n}}\ge c^{\mathrm{m}}$. To adjust the imbalances, the system operator purchases from generators capable of supplying the missing energy, and the penalty price $c^{\mathrm{n}}$ is imposed on the producer who generated the imbalance to cover this cost. The programming of the production market and the price of energy can be obtained in real-time through information provided by the system operator.

As is shown in the zooming part of Figure 1, for every time instant t, the block MPC uses the energy prices, the commitment power, and the prediction of the PV production along the entire prediction horizon, i.e., $c_{t+k}^{\mathrm{p}}$, $c_{t+k}^{\mathrm{n}}$, $P_{t+k}^{\mathrm{ref}}$, and ${\widehat{P}}_{t+k}^{\mathrm{PV}}$. Combining this information with the amount of energy stored in the battery and a set of constraints over the decision variables, the MPC strategy obtains an optimal control vector by solving a constrained optimization problem. Then, only the first element $P_t^{\mathrm{ES}}$ of this vector is applied to the system.

2.3. Mathematical Formulation for the Cost Function for the MPC Scheme

To maximize the economic revenue from PV generation, the cost function should account for the following situations:

• If the PV plant generates more than the committed energy, then the excess energy can be either delivered to the grid or stored in the BESS for future delivery. The energy stored in the BESS can be delivered to the grid in excess during periods when the energy prices are higher or used to attenuate economic penalties when the PV production is lower than the committed.

• If the PV is anticipated to generate less than the committed energy during some time intervals, then the control algorithm should use the energy stored in the BESS to cover the deficit or shift the failure to moments when the economic penalties are lower.

• The decision to charge the battery with energy from the PV plant, inject energy from the battery to the grid, or leave the battery idle is taken at every time interval targeting maximization of the economic revenue considering the forecast of energy prices, expected energy production, and state of charge of the battery along the prediction horizon.

Considering the behavior above, authors in [6] specify the following cost function:

(1)$\begin{array}{cc}\hfill \underset{P_{t+k}^{\mathrm{ES}}}{max}& {\mathcal{J}}_1=c_{t+N+1}^{\mathrm{soc}}{E}_{t+N+1}^{\mathrm{ES}}+\sum _{k=0}^N{\alpha }^{-k}{c}_{t+k}T(P_{t+k}^{\mathrm{grid}}-P_{t+k}^{\mathrm{ref}}),\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill \end{array}$

(2)$\begin{array}{cc}\hfill & c_{t+k}=\left\{\begin{array}{ccccc}{c}_{t+k}^{\mathrm{p}}& & \mathrm{if}& & P_{t+k}^{\mathrm{grid}}-P_{t+k}^{\mathrm{ref}}\ge 0,\\ c_{t+k}^{\mathrm{n}}& & \mathrm{if}& & P_{t+k}^{\mathrm{grid}}-P_{t+k}^{\mathrm{ref}}<0,\end{array}\right.\hfill \end{array}$

(3)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{grid}}={\widehat{P}}_{t+k}^{\mathrm{PV}}+P_{t+k}^{\mathrm{ES}},\hfill \end{array}$

(4)$\begin{array}{cc}\hfill & E_{t+k+1}^{\mathrm{ES}}=E_{t+k}^{\mathrm{ES}}-T{P}_{t+k}^{\mathrm{ES}},\hfill \end{array}$

(5)$\begin{array}{cc}\hfill & P_{\min }^{\mathrm{ES}}\le P_{t+k}^{\mathrm{ES}}\le P_{\max }^{\mathrm{ES}},\hfill \end{array}$

(6)$\begin{array}{cc}\hfill & E_{\min }^{\mathrm{ES}}\le E_{t+k}^{\mathrm{ES}}\le E_{\max }^{\mathrm{ES}},\hfill \end{array}$

where $c_{t+k}$ is the imbalance settlement price that at any given instant adopts a different value depending on the imbalance sign (see Equation (2)); T is the sampling period of the system in which the MPC’s control signal should be available; $\alpha \in \left[0,1\right]$ is a weighting parameter, $E_{t+N+1}^{\mathrm{ES}}$ represents the energy remaining in the BESS after the prediction horizon and $c_{t+N+1}^{\mathrm{soc}}$ is the price assigned to it. The first term of ${\mathcal{J}}_1$ sums the expected economic revenues obtained from power delivered to the grid for each interval in the prediction horizon, giving more relevance to the projected economic benefits obtained from time intervals further away from the current time. The second term in the cost function that gives an economic value to the remaining energy stored in the BESS after the prediction horizon. If the parameter $c_{t+N+1}^{\mathrm{soc}}$ is close to $c^{\mathrm{n}}$, then the control action will favor the accumulation of energy for future sale; conversely, if the value is closer to $c^{\mathrm{p}}$, then the control action will favor the delivery of energy from the BESS to the grid in the current time interval.

The restriction in Equation (3) indicates that the power delivered to the grid corresponds to the sum of the estimated power of photovoltaic panels and the power exchanged with the BESS. Equation (4) models the energy stored in the BESS on the whole prediction horizon. Equation (5) indicates that the power taken from the panels or delivered to the grid from the BESS is limited by the rated power capacity of the converter, and Equation (6) restricts the state of charge of the battery bank to avoid its accelerated degradation or damage. The selection of the sampling time (T) varies depending on the application of interest. The main factors that intervene in the appropriate selection of T depend on the size of the problem and the time constraints. In previous works that apply MPC strategies for controlling similar power generation systems, authors have proposed values of $T=15$ min [21,23] and $T=4$ min [6] to satisfy the control objectives.

The piece-wise expression for the cost in Equation (2) precludes using standard Linear Programming or Quadratic Programming solvers to find the solution. The authors in [6] identify this as a big limitation to achieving the required solving time of 4 min. Then they propose ad hoc simplifications to represent the cost function in an LP form that can be efficiently solved with standard solvers. Notice that the piece-wise constraints in Equation (2) are common in many applications such as [27,28], where there are variables that can take different values depending on whether a condition is met.

2.4. Representation as an LP Problem

In [6], the authors introduce two different variables to represent the power delivered to the grid depending on its value compared to the committed power:

(7)$P_{t+k}^{\mathrm{grid}}=\left\{\begin{array}{ccccc}{P}_{t+k}^{\mathrm{p}}& & \mathrm{if}& & P_{t+k}^{\mathrm{grid}}\ge P_{t+k}^{\mathrm{ref}},\\ P_{t+k}^{\mathrm{n}}& & \mathrm{if}& & P_{t+k}^{\mathrm{grid}}<P_{t+k}^{\mathrm{ref}},\end{array}\right.$

Using the change of variables and considering the recursive computation of Equation (4), it is possible to rewrite the optimization problem in the following equivalent form [6]:

(8)$\begin{array}{cc}\hfill \underset{P_{t+k}^{\mathrm{p}},P_{t+k}^{\mathrm{n}}}{max}& {\mathcal{J}}_1^{*}=c_{t+N+1}^{\mathrm{soc}}{E}_{t+N+1}^{\mathrm{ES}}+\sum _{k=0}^N{\alpha }^{-k}T\left[c_{t+k}^{\mathrm{p}}(P_{t+k}^{\mathrm{p}}-P_{t+k}^{\mathrm{ref}})+c_{t+k}^{\mathrm{n}}(P_{t+k}^{\mathrm{n}}-P_{t+k}^{\mathrm{ref}})\right],\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill \end{array}$

(9)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{p}}+P_{t+k}^{\mathrm{n}}-P_{t+k}^{\mathrm{ref}}={\widehat{P}}_{t+k}^{\mathrm{PV}}+P_{t+k}^{\mathrm{ES}},\hfill \end{array}$

(10)$\begin{array}{cc}\hfill & E_{t+k+1}^{\mathrm{ES}}=E_{t+k}^{\mathrm{ES}}-T{P}_{t+k}^{\mathrm{ES}},\hfill \end{array}$

(11)$\begin{array}{cc}\hfill & P_{\min }^{\mathrm{ES}}\le P_{t+k}^{\mathrm{ES}}\le P_{\max }^{\mathrm{ES}},\hfill \end{array}$

(12)$\begin{array}{cc}\hfill & E_{\min }^{\mathrm{ES}}\le E_{t+k}^{\mathrm{ES}}\le E_{\max }^{\mathrm{ES}},\hfill \end{array}$

(13)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{p}}\ge P_{t+k}^{\mathrm{ref}},\hfill \end{array}$

(14)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{n}}\le P_{t+k}^{\mathrm{ref}},\hfill \end{array}$

(15)$\begin{array}{cc}\hfill & (P_{t+k}^{\mathrm{p}}-P_{t+k}^{\mathrm{ref}})(P_{t+k}^{\mathrm{n}}-P_{t+k}^{\mathrm{ref}})=0,\hfill \end{array}$

Although the non-linear constraint in Equation (15) still precludes us from directly using traditional LP solvers, the authors in [6] argue that this product constraint is redundant when the other constraints are satisfied and $c_{t+k}^{\mathrm{n}}\ge c_{t+k}^{\mathrm{p}}$ for all time intervals within the prediction horizon, which is a natural condition in the target scenario. Therefore, it is possible to just omit this constraint, and the resulting optimization problem involves only linear cost function and linear constraints, i.e., an LP problem that can be efficiently solved using traditional tools.

Nevertheless, this solution involves an ad hoc transformation and is not directly applicable for a related problem that includes some other non-linearities such rate-based variable battery efficiency. Additionally, the work in [6] does not show the development with sufficient detail to directly reproduce the reported results. For completeness and, in order to compare our proposal with [6], we provide details in Appendix B to solve the optimization problem in Equation (8). It is noted that this problem cannot be solved directly with an available linear optimization solver, so it is necessary to rewrite it in the standard form of linear programming problems as follows:

(16)$\begin{array}{cc}\hfill \underset{\mathrm{s}}{max}& {\mathbf{h}}^T\mathbf{s},\hfill \end{array}$

(17)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& \mathbf{Ds}\le \mathbf{d},\hfill \end{array}$

where $\mathbf{h}$ summarizes the imbalance costs and the stored energy price over the prediction horizon. The vector $\mathbf{s}$ contains the decision variables $P_{t+k}^{\mathrm{p}}$ and $P_{t+k}^{\mathrm{n}}$ with $k=0,\cdots ,N$. The vector $\mathbf{d}$ and the matrix $\mathbf{D}$ define the constraints in Equations (9)–(14). The definitions of $\mathbf{h}$, $\mathbf{s}$, $\mathbf{d}$, and $\mathbf{D}$ are given in Appendix B. As can be seen, the problem posed by Equations (16) and (17) has the form of an LP optimization problem, where the upper and lower bound constraints of the decision variables have been rewritten as inequality constraints.

2.5. Implementation of the MPC Schemes and Measurement of Execution Time

To implement and compare the reference LP formulation for the optimization problem described previously and the new MILP formulation introduced in the next section, we use Matlab-Simulink to simulate the system described in Figure 1, including the CPLEX optimization toolbox [29] to solve the optimization problem for each MPC iteration. To estimate the execution time, we use the tic-toc functionality in Matlab to measure the elapsed time between two points of interest in the code. We measure the required time of the solver to get the optimal input vector for 165 MPC iterations. The worst-case value of the elapsed time is reported to show that our proposal is competitive to solve each MPC optimization problem in the required time. All tests were executed on a desktop computer with an Intel Core i7-8700 CPU @ 3.20 GHz processor, 16 GB of RAM, running Windows 11 and MATLAB 2019b.

3. Formulation of the New MPC-MILP Approach

We propose rewriting the problem of the optimal economic management of the storage system in a PV generation plant using an MILP approach. One of the advantages of the MILP approach is that the proposed control strategy can be adapted to the current conditions of the electricity market, as will be discussed later in Section 5.

3.1. Mathematical Formulation Considering an Ideal Battery

In order to solve the problem in (1) using MILP approach, the equality constraints in Equations (3) and (4) were introduced in the cost function ${\mathcal{J}}_1$, yielding the following optimization problem:

(18)$\begin{array}{cc}\hfill \underset{P_{t+k}^{\mathrm{ES}}}{max}& {\mathcal{J}}_1=\sum _{k=0}^N{c}_{t+k}T\left[P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]-c_{t+N+1}^{\mathrm{soc}}T{P}_{t+k}^{\mathrm{ES}},\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill \end{array}$

(19)$\begin{array}{cc}\hfill & c_{t+k}=\left\{\begin{array}{ccccc}{\alpha }^{-k}{c}_{t+k}^{\mathrm{p}}& & \mathrm{if}& & P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\ge 0,\\ {\alpha }^{-k}{c}_{t+k}^{\mathrm{n}}& & \mathrm{if}& & P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})<0,\end{array}\right.\hfill \end{array}$

(20)$\begin{array}{cc}\hfill & P_{\min }^{\mathrm{ES}}\le P_{t+k}^{\mathrm{ES}}\le P_{\max }^{\mathrm{ES}},\hfill \end{array}$

(21)$\begin{array}{cc}\hfill & E_{\min }^{\mathrm{ES}}\le E_t^{\mathrm{ES}}-T\sum _{s=0}^{k-1}{P}_{t+s}^{\mathrm{ES}}\le E_{\max }^{\mathrm{ES}}.\hfill \end{array}$

In Appendix C, we explain how to transform an optimization problem involving logical rules into an MILP problem. We consider a binary variable to deal with the condition $c_{t+k}$ that takes different values as $P_{t+k}^{\mathrm{ES}}$ is higher or lower than $(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})$. The condition $P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\ge 0$ can be associated to a binary variable ${\delta }_{t+k}\in \left\{0,1\right\}$ such that

(22)$\begin{array}{cc}\hfill P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\ge 0& \leftrightarrow {\delta }_{t+k}=1.\hfill \end{array}$

Since $P_{t+k}^{\mathrm{ES}}$ is bounded by $P_{\min }$ and $P_{\max }$, the logical statement in Equation (22) can be rewritten as linear constraints. Additionally, introducing a new continuous variable $z_{t+k}$ as follows:

(23)$z_{t+k}=\left[P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k},$

which also can be rewritten as linear constraints, then the optimization problem in Equation (18) is equivalent to:

(24)$\begin{array}{cc}\hfill \underset{P_{t+k}^{\mathrm{ES}}}{max}& {\mathcal{J}}_1^{+}=\sum _{k=0}^N{\alpha }^{-k}T\left(c_{t+k}^{\mathrm{p}}-c_{t+k}^{\mathrm{n}}\right)z_{t+k}+({\alpha }^{-k}{c}_{t+k}^{\mathrm{n}}-c_{t+N+1}^{\mathrm{soc}})T{P}_{t+k}^{\mathrm{ES}},\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill \end{array}$

(25)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{ES}}\ge \left[P_{\min }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right](1-{\delta }_{t+k})+(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}}),\hfill \end{array}$

(26)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{ES}}\le \left[P_{\max }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k}+(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}}),\hfill \end{array}$

(27)$\begin{array}{cc}\hfill & z_{t+k}\ge \left[P_{\min }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k},\hfill \end{array}$

(28)$\begin{array}{cc}\hfill & z_{t+k}\le \left[P_{\max }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k},\hfill \end{array}$

(29)$\begin{array}{cc}\hfill & z_{t+k}\le P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})-\left[P_{\min }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right](1-{\delta }_{t+k}),\hfill \end{array}$

(30)$\begin{array}{cc}\hfill & z_{t+k}\ge P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})-\left[P_{\max }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right](1-{\delta }_{t+k}),\hfill \end{array}$

(31)$\begin{array}{cc}\hfill & {\delta }_{t+k}\in \left\{0,1\right\},\hfill \end{array}$

(32)$\begin{array}{cc}\hfill & E_{\min }^{\mathrm{ES}}\le E_t^{\mathrm{ES}}-T\sum _{s=0}^{k-1}{P}_{t+s}^{\mathrm{ES}}\le E_{\max }^{\mathrm{ES}},k=1,\cdots ,N\hfill \end{array}$

where ${\mathcal{J}}_1^{+}$ stands for an equivalent problem of ${\mathcal{J}}_1$. Then, using the definitions in Appendix A, this problem can be rewritten in short form as follows:

(33)$\begin{array}{cc}\hfill \underset{\mathrm{x}\in \mathcal{X}}{max}& {\mathbf{f}}^T\mathbf{x},\hfill \end{array}$

(34)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& \mathbf{Ax}\le \mathbf{a},\hfill \end{array}$

where $\mathbf{f}$ summarizes the imbalance costs and the value of the stored energy over the prediction horizon. The vector $\mathbf{x}$ contains the continuous decision variables $P_{t+k}^{\mathrm{ES}}$ and $z_{t+k}$, and the logical decision variable ${\delta }_{t+k}$ over the prediction horizon N. The vector $\mathbf{a}$ and the matrix $\mathbf{A}$, contains the constraints in Equations (25)–(32). The set $\mathcal{X}$ is the space that contains all possible values of $\mathbf{x}$. The definitions of $\mathbf{f}$, $\mathbf{x}$, $\mathbf{a}$, and $\mathbf{A}$ are given in the Appendix D. Note that the change of variable proposed in Equation (22) requires an algorithm that allows the use of binary variables in the formulation, which is not possible with the LP solution proposed in [6]. Therefore, the optimization problem in Equation (33) is solved using standard solvers for MILP.

3.2. Validation Considering an Ideal Battery

To corroborate the functional equivalence between the solution with MPC-LP reported in Section 2 from previous literature and the proposed solution using MPC-MILP, we simulate the same system and conditions reported for the experiments in [6]. Table 1 shows the parameters of the PV+BESS plant and the MPC control parameters.

In this scenario, an ideal battery model with $100\%$ efficiency is considered. In other words, it is assumed that there are no losses in the charging and discharging processes in the battery, so the energy stored in the battery is equal to the energy available for subsequent discharge.

The data used for the MPC control algorithm explained in Section 2.2 are shown in Figure 2a–c, which correspond to the prices of one simulated day of the Spanish market, the committed power, and the estimated power, respectively. The moment of interest to sell energy in this context is when the committed power has been settled, where the MPC control algorithm acts optimally to generate economic profit or, in the worst case, reduce possible penalties due to the stochastic character of the power coming from the solar panels. As shown in Figure 2c, from 8:00 to 19:00 hours is the period that the estimated power generation from the panels is available. Therefore, it is possible to sell the total power of the plant or accumulate energy in the batteries.

The test scenario considers that the BESS has 25% of its total energy capacity at 8:00 hours. Figure 3a shows how both algorithms decide to completely discharge the batteries, equivalent to supplying all the available energy to the grid. This makes sense since both control algorithms have the same prior knowledge to act optimally to increase profits. In this case, the power generated by the PV-BESS system is greater than the committed power, as shown in Figure 3b, which increases the economic revenue according to the imbalance price. It can also be seen that throughout the operation of the system, the total power supplied to the grid is bigger than the committed power, which verifies that both algorithms avoid penalties. Starting at 11:00 hours, the BESS stores energy from the solar panels as indicated by negative values of $P_t^{\mathrm{ES}}$ in Figure 3c. The simulation results also show in Figure 3d that the heuristic solution employing MPC-LP proposed in [6] is valid and equivalent to the formulation MPC-MILP proposed in this work when considering an ideal battery. As can be seen, the results obtained with the proposed MPC-MILP formulation are validated by simulation results, offering the same performance as the strategy available in the literature.

Additionally, the measurements of the execution time for both MPC-LP and MPC-MILP, in the worst case, were about 18 and 46 milliseconds, respectively. As we explained above, the time required for the control algorithm to make the optimal decision in the PV-BESS application is $T=4$ min. It is noted that although the time required by the algorithms is not deterministic, the control signal is found quickly enough to actuate over the system. This fact shows that the MILP approach is not computationally prohibitive for real-time operation as the initial assumption in [6].

So far, it has been shown how to perform the MPC control considering an idealized battery model. The described problem has not considered the losses of the charging and discharging of the batteries. Neglecting these parameters can result in poor decision-making by the control algorithm. Consequently, in Figure 3a, a high amount of charging and discharging can be seen. Next, an alternative solution considering a model for the efficiency of the charging and discharging of the battery is proposed.

3.3. Formulation Considering BESS with Rate-Based Variable Efficiency

In this section, the same problem of interest as [6] that has been already discussed in Section 2.3 is considered, but introducing the charge and discharge efficiency of the battery. In this work, it is considered that the BESS is charged with constant rate efficiency ${\eta }^{\mathrm{c}}$ and discharged with $1/{\eta }^{\mathrm{d}}$, and it is not directly solvable using [6]. We note that both of these efficiencies are positive quantities, and smaller than one, as expected $(0<{\eta }^{\mathrm{c}}\le 1$$\mathrm{and}0<{\eta }^{\mathrm{d}}\le 1)$. However, as the objective function values the stored energy in the battery, the definition of efficiency must correspond to the natural definition for the charging efficiency, regardless of the actual operational condition. This leads to an efficiency definition as Equation (37) where ${\eta }^{\mathrm{c}}\le 1$ for charging $(P_t^{\mathrm{ES}}<0)$ and $1/{\eta }^{\mathrm{d}}\ge 1$ for discharging $(P_t^{\mathrm{ES}}\ge 0)$. Then, the problem to solve is the following:

(35)$\begin{array}{cc}\hfill \underset{P_{t+k}^{\mathrm{ES}}}{max}& {\mathcal{J}}_2=\sum _{k=0}^N{c}_{t+k}T\left[P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]-c_{t+N+1}^{\mathrm{soc}}T{\eta }_{t+k}{P}_{t+k}^{\mathrm{ES}},\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill \end{array}$

(36)$\begin{array}{cc}\hfill & c_{t+k}=\left\{\begin{array}{ccccc}{\alpha }^{-k}{c}_{t+k}^{\mathrm{p}}& & \mathrm{if}& & P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\ge 0,\\ {\alpha }^{-k}{c}_{t+k}^{\mathrm{n}}& & \mathrm{if}& & P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})<0,\end{array}\right.\hfill \end{array}$

(37)$\begin{array}{cc}\hfill & {\eta }_{t+k}=\left\{\begin{array}{ccccc}1/{\eta }^{\mathrm{d}}& & \mathrm{if}& & P_{t+k}^{\mathrm{ES}}\ge 0,\\ {\eta }^{\mathrm{c}}& & \mathrm{if}& & P_{t+k}^{\mathrm{ES}}<0,\end{array}\right.\hfill \end{array}$

(38)$\begin{array}{cc}\hfill & P_{\min }^{\mathrm{ES}}\le P_{t+k}^{\mathrm{ES}}\le P_{\max }^{\mathrm{ES}},\hfill \end{array}$

(39)$\begin{array}{cc}\hfill & E_{\min }^{\mathrm{ES}}\le E_t^{\mathrm{ES}}-T\sum _{s=0}^{k-1}{\eta }_{t+k}{P}_{t+s}^{\mathrm{ES}}\le E_{\max }^{\mathrm{ES}}.\hfill \end{array}$

To solve this problem, an approach similar to Section 3 is used. Consider a binary value to deal with the logical condition on ${\eta }_{t+k}$:

(40)$\begin{array}{c}\hfill P_{t+k}^{\mathrm{ES}}\ge 0\leftrightarrow {\beta }_{t+k}=1.\end{array}$

Additionally, a continuous auxiliary variable $w_{t+k}$ is introduced so that $w_{t+k}=P_{t+k}^{\mathrm{ES}}{\beta }_{t+k}$. Then, following Appendix C the optimization problem in Equation (35) can be rewritten as follows:

(41)$\begin{array}{cc}\hfill \underset{P_{t+k}^{\mathrm{ES}}}{max}{\mathcal{J}}_2^{+}=\sum _{k=0}^N{\alpha }^{-k}T& \left(c_{t+k}^{\mathrm{p}}-c_{t+k}^{\mathrm{n}}\right)z_{t+k}+T({\alpha }^{-k}{c}_{t+k}^{\mathrm{n}}-c_{t+N+1}^{\mathrm{soc}}{\eta }^{\mathrm{c}})P_{t+k}^{\mathrm{ES}}-c_{t+N+1}^{\mathrm{soc}}T({\eta }^{\mathrm{d}}-{\eta }^{\mathrm{c}})w_{t+k},\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill \end{array}$

(42)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{ES}}\ge \left[P_{\min }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right](1-{\delta }_{t+k})+(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}}),\hfill \end{array}$

(43)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{ES}}\le \left[P_{\max }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k}+(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}}),\hfill \end{array}$

(44)$\begin{array}{cc}\hfill & z_{t+k}\ge \left[P_{\min }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k},\hfill \end{array}$

(45)$\begin{array}{cc}\hfill & z_{t+k}\le \left[P_{\max }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right]{\delta }_{t+k},\hfill \end{array}$

(46)$\begin{array}{cc}\hfill & z_{t+k}\le P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})-\left[P_{\min }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right](1-{\delta }_{t+k}),\hfill \end{array}$

(47)$\begin{array}{cc}\hfill & z_{t+k}\ge P_{t+k}^{\mathrm{ES}}-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})-\left[P_{\max }-(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\right](1-{\delta }_{t+k}),\hfill \end{array}$

(48)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{ES}}\ge P_{\min }(1-{\beta }_{t+k}),\hfill \end{array}$

(49)$\begin{array}{cc}\hfill & P_{t+k}^{\mathrm{ES}}\le P_{\max }{\beta }_{t+k},\hfill \end{array}$

(50)$\begin{array}{cc}\hfill & w_{t+k}\ge P_{\min }{\beta }_{t+k},\hfill \end{array}$

(51)$\begin{array}{cc}\hfill & w_{t+k}\le P_{\max }{\beta }_{t+k},\hfill \end{array}$

(52)$\begin{array}{cc}\hfill & w_{t+k}\le P_{t+k}^{\mathrm{ES}}-P_{\min }(1-{\beta }_{t+k}),\hfill \end{array}$

(53)$\begin{array}{cc}\hfill & w_{t+k}\ge P_{t+k}^{\mathrm{ES}}-P_{\max }(1-{\beta }_{t+k}),\hfill \end{array}$

(54)$\begin{array}{cc}\hfill & {\delta }_{t+k}\in \left\{0,1\right\},\hfill \end{array}$

(55)$\begin{array}{cc}\hfill & {\beta }_{t+k}\in \left\{0,1\right\},\hfill \end{array}$

(56)$\begin{array}{cc}\hfill & ({\delta }_{t+k}-{\beta }_{t+k})(P_{t+k}^{\mathrm{ref}}-{\widehat{P}}_{t+k}^{\mathrm{PV}})\le 0,\hfill \end{array}$

(57)$\begin{array}{cc}\hfill & E_{\min }^{\mathrm{ES}}\le E_t^{\mathrm{ES}}-T\sum _{s=0}^{k-1}{\eta }_{t+k}{P}_{t+s}^{\mathrm{ES}}\le E_{\max }^{\mathrm{ES}},k=1,\cdots ,N,\hfill \end{array}$

where ${\mathcal{J}}_2^{+}$ stands for an equivalent problem of ${\mathcal{J}}_2$. Then, using the definitions in Appendix A, this MILP problem can be rewritten in short form as follows:

(58)$\begin{array}{cc}\hfill \underset{\mathbf{y}\in \mathcal{Y}}{max}& {\mathbf{g}}^T\mathbf{y},\hfill \end{array}$

(59)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& \mathbf{B}\mathbf{y}\le \mathbf{b},\hfill \end{array}$

where $\mathbf{g}$ summarizes the imbalance costs and the value of the stored energy over the prediction horizon. The vector $\mathbf{y}$ contains the continuous decision variables $P_{t+k}^{\mathrm{ES}}$, $z_{t+k}$ and $w_{t+k}$, and the logical decision variables ${\delta }_{t+k}$ and ${\beta }_{t+k}$ over the prediction horizon N. The vector $\mathbf{b}$ and the matrix $\mathbf{B}$, contains the constraints in Equations (42)–(57). The set $\mathcal{Y}$ is the space that contains all possible values of $\mathbf{y}$. The definitions of $\mathbf{g}$, $\mathbf{y}$, $\mathbf{b}$, and $\mathbf{B}$ are given in the Appendix E.

Next, the results of the simulation of the proposed formulation applied to a BESS with non-ideal charge and discharge efficiency are shown and compared with the results achieved when solved with the technique presented in Section 2.4, i.e., when assuming ideal efficiency in the control. In other words, the objective is to identify the economic benefit of using MPC-MILP in PV-BESS, considering a rate-based variable efficiency model for the battery, and compare it to the solution obtained with the MPC-LP strategy reported in the literature.

4. Results

The response of the MPC-LP, as described in Section 2.4, and the proposed MPC-MILP control algorithms, when applied to a system model with non-ideal charge and discharge efficiency, are shown. The initial conditions are the same as in the previous section. Figure 4 shows the comparison of the results obtained using the MPC-LP and MPC-MILP control algorithms. As can be seen, the system’s dynamic performance with the MPC-MILP control algorithm and a BESS model with a rate-based variable efficiency is notably better. Figure 4a shows how the performance of the MPC-MILP algorithm does not completely discharge the battery, unlike MPC-LP. Although both algorithms satisfy the output power needed to avoid a penalty, i.e., $P^{\mathrm{grid}}\ge P^{\mathrm{ref}}$, the power transients are less aggressive using the strategy MPC-MILP, as shown in Figure 4b. This is a direct consequence of making decisions about the charging and discharging cycles in the batteries, as shown in Figure 4c. It is also noted that reducing charge and discharge cycles contributes to extending the BESS’s health. The improvements obtained in the simulation results are understandable from a physical point of view. When considering a simplified battery model, the algorithm described in Section 2.4 estimates the charging and discharging processes with less precision because it does not consider the losses associated with the power exchanged with the BESS. This fact translates into a poor decision to store energy in the BESS or sell it to the power grid. While the MPC-LP algorithm acts aggressively, the MPC-MILP algorithm performs better in making decisions to supply power, charge the batteries or disconnect them from the system. Figure 4d shows how the proposed control algorithm obtains a bigger ${\mathcal{J}}_2$ over time compared with the MILP-LP, which is the main control objective, i.e., it makes the PV plant’s market participation more profitable. The computation of the cost function ${\mathcal{J}}_2$ using the optimal control signal of each time instant shows that the increment in the economic revenues obtained using our proposal is close to 21% with respect to the MPC-LP approach.

Furthermore, the execution time of the model that considers the rate-based variable efficiency increases with respect to the highly efficient MPC-LP method proposed in [6]. However, the solution of the MPC-MILP algorithm in the worst case was found in 55 milliseconds, which is still small enough for $T=4$ min. Additionally, the measurements for both MPC-LP and MPC-MILP, in the worst case, were about 18 and 46 milliseconds, respectively.

Influence of the Solar Radiation Estimation in the Control Scheme

In this section, a study of the influence of solar estimation on the performance of the control algorithm is carried out. As mentioned above, the database provided by PVGIS is used in this work. However, the estimated data could differ from the actual radiation due to the presence of meteorological phenomena, causing the stochastic behavior of this energy source. Estimated solar radiation data will become more accurate as parameter estimation techniques evolve. It should be considered that an inaccurate ${\widehat{P}}_t^{\mathrm{PV}}$ estimation model could lead to penalties in the market under the aforementioned conditions. An experiment is then carried out to assess how much the difference between the estimated ${\widehat{P}}_t^{\mathrm{PV}}$ power used in the control algorithm and the actual ${\widehat{P}}_t^{\mathrm{PV}}$ power can affect. Figure 5a shows an example of the difference in solar radiation on a slightly cloudy day.

Under these conditions, the actual ${\widehat{P}}_t^{\mathrm{PV}}$ is less than the model used in the algorithm. In this test, the approach described in the Section 2.4 and the one proposed in Section 3.3 are compared. The Figure 5b shows the objective function for both control algorithms. As can be seen, the MPC-MILP approach proposed in this work continues to generate higher economic returns compared to the MPC-LP approach.

5. Discussion

The use of the MPC strategy has increased in recent years for electric power generation due to its advantages. In particular, it allows optimal decision-making considering the physical and operational constraints of the system. The proposal incorporated can implement more realistic models of the storage system, including its efficiency, leading to better performance of the control loop than a previous method found in the literature, translating into higher earnings for the market participant.

This article provides a detailed methodology for the mathematical formulation of the optimization problem using the paradigm of Mixed Integer Linear Programming. Additionally, some measurements of the computational burden of this method are provided, showing that although higher than the alternative method found in [6], it is sufficiently low for the problem of interest, given the computational capability of modern desktop computers. Future work may consider integrating other energy sources, whether renewable or non-renewable, and applying the control algorithm to plants integrating different generation technologies to minimize energy production costs. Further, the proposed methodology can be adapted to different problems with similar characteristics.

This is very important, as rules in regulated markets are different from country to country and may change with new regulations. For example, the operating procedure through the provisions of P.O.-14.1 in the Spanish market has undergone recent modifications, where now there are four situations in which participating agents can incur in the market. However, we emphasize that the proposed control scheme can be directly applied to the new circumstances because it can adequately consider the value of $c_{t+k}$ in the constraints of the optimization problem for each situation.

6. Conclusions

This study investigates how to obtain economic benefits by participating in the Spanish electricity market. This work presents a new control strategy that incorporates the charging and discharging efficiency of the BESS model in PV plants with storage. The proposed strategy maximizes the economic benefit of a given power scheduling in electricity market participation by optimizing the use of the battery storage system. The proposal is validated by simulation results and compared with important previous works that do not consider the battery efficiency model. Three simulation examples show that the MPC-MILP algorithm can maintain the power above or equal to the committed power thanks to the integration of the BESS. The simulation results also show how the rate-based variable efficiency battery model can further improve the economic benefit using an MPC-MILP approach. It is also shown that with modern desktop computers, the proposed control is feasible within a sampling period of $T=4$ min sufficient for the intended application. Consequently, we believe that the proposed control methodology will allow an increased number of agents to be more competitive in the market, eventually potentially reducing the final cost to consumers and to society.

Author Contributions

Conceptualization, A.L.C., R.L.A., C.S. and J.C.A.; methodology, A.L.C., R.L.A., C.S. and J.C.A.; software, A.L.C.; validation, J.R. and G.C.; formal analysis, A.L.C., R.L.A., J.R., G.C., C.S. and J.C.A.; investigation, A.L.C., R.L.A., G.C., C.S. and J.C.A.; resources, G.C., C.S. and J.C.A.; writing—original draft preparation, A.L.C., R.L.A., J.R., G.C., C.S. and J.C.A.; writing—review and editing, A.L.C., R.L.A., J.R., G.C., C.S. and J.C.A.; visualization, G.C.; supervision, G.C., C.S. and J.C.A. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

MPC	Model Predictive Control
MILP	Mixed-Integer Linear Programming
PV	Photovoltaic
AC	Alternating Current
DC	Direct Current
BESS	Battery Energy Storage System
SoC	State of Charge
LP	LinearProgramming
QP	Quadratic Programming
PSD	Power Spectral Density
PVGIS	PhotoVoltaic Geographic Information System
CPLEX	IBM ILOG CPLEX Optimization Studio

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

Figures and Table

Figure 1. Model predictive control scheme.

Figure 2. (a) Synthetic energy prices resembling the ones reported in [6]; (b) committed power per hour; and (c) power production of photovoltaic panels considering ideal solar radiation.

Figure 3. (a) Exchange power with the BESS; (b) cost of the optimization algorithm; (c) energy stored in batteries; and (d) power injected into the grid.

Figure 4. (a) Exchange power with the BESS; (b) cost of the optimization algorithm; (c) energy stored in batteries; and (d) power injected into the grid.

Figure 5. (a) Difference between the estimated power and the real power for a slightly cloudy day (b) Cost of MPC-LP and MPC-MILP algorithms.

Table 1

PV+BESS plant and MPC control parameters.

N	T	${\mathit{c}}^{\mathrm{p}}$	${\mathit{c}}^{\mathrm{n}}$	${\mathit{\lambda }}_{\mathit{k}}$
60	4 min	as in Figure 2	as in Figure 2	$0.{999}^k$
$c_{t+N+1}^{\mathrm{soc}}$	$P_{\min }^{\mathrm{ES}}$	$P_{\max }^{\mathrm{ES}}$	$E_{\min }^{\mathrm{ES}}$	$E_{\max }^{\mathrm{ES}}$
$c^{\mathrm{n}}$	$-500$ kW	500 kW	0	$0.8$ MWh

Appendix A

Vector and matrices for the problem in Equation (1) considering a prediction horizon N:(A1)$$\mathbf{\alpha }=\left[\begin{array}{c}1\\ {\alpha }^{-1}\\ ⋮\\ {\alpha }^{-N}\end{array}\right],{\mathbf{c}}_{\mathrm{n}}=\left[\begin{array}{c}{c}_t^{\mathrm{n}}\\ c_{t+1}^{\mathrm{n}}\\ ⋮\\ c_{t+N}^{\mathrm{n}}\end{array}\right],{\mathbf{c}}_{\mathrm{p}}=\left[\begin{array}{c}{c}_t^{\mathrm{p}}\\ c_{t+1}^{\mathrm{p}}\\ ⋮\\ c_{t+N}^{\mathrm{p}}\end{array}\right],{\mathbf{c}}_{\mathrm{soc}}=\left[\begin{array}{c}{c}_{t+N+1}^{\mathrm{soc}}\\ c_{t+N+1}^{\mathrm{soc}}\\ ⋮\\ c_{t+N+1}^{\mathrm{soc}}\end{array}\right],$$ (A2)$${\mathbf{P}}_{\min }=\left[\begin{array}{c}{P}_{\min }^{\mathrm{ES}}\\ P_{\min }^{\mathrm{ES}}\\ ⋮\\ P_{\min }^{\mathrm{ES}}\end{array}\right],{\mathbf{P}}_{\max }=\left[\begin{array}{c}{P}_{\max }^{\mathrm{ES}}\\ P_{\max }^{\mathrm{ES}}\\ ⋮\\ P_{\max }^{\mathrm{ES}}\end{array}\right],{\mathbf{P}}_{\mathrm{ref}}=\left[\begin{array}{c}{P}_t^{\mathrm{ref}}\\ P_{t+1}^{\mathrm{ref}}\\ ⋮\\ P_{t+N}^{\mathrm{ref}}\end{array}\right],{\widehat{\mathbf{P}}}_{\mathrm{PV}}=\left[\begin{array}{c}{\widehat{P}}_t^{\mathrm{PV}}\\ {\widehat{P}}_{t+1}^{\mathrm{PV}}\\ ⋮\\ {\widehat{P}}_{t+N}^{\mathrm{PV}}\end{array}\right],$$ (A3)$${\mathbf{E}}_{\min }={\displaystyle \frac{1}{T}}\left[\begin{array}{c}{E}_t^{\mathrm{ES}}-E_{\max }^{\mathrm{ES}}\\ E_t^{\mathrm{ES}}-E_{\max }^{\mathrm{ES}}\\ ⋮\\ E_t^{\mathrm{ES}}-E_{\max }^{\mathrm{ES}}\end{array}\right],{\mathbf{E}}_{\max }={\displaystyle \frac{1}{T}}\left[\begin{array}{c}{E}_t^{\mathrm{ES}}-E_{\min }^{\mathrm{ES}}\\ E_t^{\mathrm{ES}}-E_{\min }^{\mathrm{ES}}\\ ⋮\\ E_t^{\mathrm{ES}}-E_{\min }^{\mathrm{ES}}\end{array}\right],\mathbf{F}=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 1& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \cdots & 1\end{array}\right],$$ where $\mathbf{\alpha }$, ${\mathbf{c}}_{\mathrm{n}}$, ${\mathbf{c}}_{\mathrm{p}}$, ${\mathbf{c}}_{\mathrm{soc}}$, ${\mathbf{P}}_{\min }$, ${\mathbf{P}}_{\max }$, ${\mathbf{P}}_{\mathrm{ref}}$, ${\mathbf{E}}_{\min }$, ${\widehat{\mathbf{P}}}_{\mathrm{PV}}$, ${\mathbf{E}}_{\max }\in {\mathbb{R}}^{N+1}$ and $\mathbf{F}\in {\mathbb{R}}^{N+1\times N+1}$. $\mathbf{I}\in {\mathbb{R}}^{N+1\times N+1}$ is the identity matrix. Additionally, $\mathbf{0}$ represents the zero vector or matrix, and its dimensions should be inferred from the context, i.e., $\mathbf{0}$ is a vector $\in {\mathbb{R}}^{N+1}$ or is a matrix $\in {\mathbb{R}}^{N+1\times N+1}$.

Appendix B. Vectors and Matrices for the Solution in [6]

In this section, the following vectors and matrices are introduced to solve the LP problem given in Equations (16)–(17). Considering Equation (16), the vector $\mathbf{h}$ is given by (A4)$$\mathbf{h}=\left[\begin{array}{c}{\mathbf{h}}_1\\ {\mathbf{h}}_2\end{array}\right],\begin{array}{c}{\mathbf{h}}_1=\mathbf{\alpha }\odot {\mathbf{c}}_{\mathrm{p}}-{\mathbf{c}}_{\mathrm{soc}},\\ {\mathbf{h}}_2=\mathbf{\alpha }\odot {\mathbf{c}}_{\mathrm{n}}-{\mathbf{c}}_{\mathrm{soc}},\end{array}$$ where ⊙ is the element-wise product operator, $\mathbf{\alpha }$, ${\mathbf{c}}_{\mathrm{p}}$, ${\mathbf{c}}_{\mathrm{n}}$ and ${\mathbf{c}}_{\mathrm{soc}}$ are defined in the Appendix A. The vector $\mathbf{s}$ in Equation (16) contains the decision vectors ${\mathbf{P}}_{\mathrm{p}}$ and ${\mathbf{P}}_{\mathrm{n}}$ that are given by (A5)$$\mathbf{s}=\left[\begin{array}{c}{\mathbf{P}}_{\mathrm{p}}\\ {\mathbf{P}}_{\mathrm{n}}\end{array}\right],{\mathbf{P}}_{\mathrm{p}}=\left[\begin{array}{c}{P}_t^{\mathrm{p}}\\ P_{t+1}^{\mathrm{p}}\\ ⋮\\ P_{t+N}^{\mathrm{p}}\end{array}\right],{\mathbf{P}}_{\mathrm{n}}=\left[\begin{array}{c}{P}_t^{\mathrm{n}}\\ P_{t+1}^{\mathrm{n}}\\ ⋮\\ P_{t+N}^{\mathrm{n}}\end{array}\right],$$ where ${\mathbf{P}}_{\mathrm{p}},{\mathbf{P}}_{\mathrm{n}}\in {\mathbb{R}}^{N+1}$. On the other hand, considering Equation (17), the matrix inequality $\mathbf{D}$ and the vector $\mathbf{d}$ are given by (A6)$$\mathbf{D}=\left[\begin{array}{cc}\hfill -\mathbf{I}& -\mathbf{I}\\ \hfill \mathbf{I}& \mathbf{I}\\ \hfill -\mathbf{F}& -\mathbf{F}\\ \hfill \mathbf{F}& \mathbf{F}\\ \hfill -\mathbf{I}& \mathbf{0}\\ \hfill \mathbf{0}& \mathbf{I}\end{array}\right],\mathbf{d}=\left[\begin{array}{c}\hfill -({\mathbf{P}}_{\min }+{\mathbf{P}}_{\mathrm{ref}}+{\widehat{\mathbf{P}}}_{\mathrm{PV}})\\ \hfill {\mathbf{P}}_{\max }+{\mathbf{P}}_{\mathrm{ref}}+{\widehat{\mathbf{P}}}_{\mathrm{PV}}\\ \hfill -({\mathbf{E}}_{\min }+\mathbf{F}({\mathbf{P}}_{\mathrm{ref}}+{\widehat{\mathbf{P}}}_{\mathrm{PV}}))\\ \hfill {\mathbf{E}}_{\max }+\mathbf{F}({\mathbf{P}}_{\mathrm{ref}}+{\widehat{\mathbf{P}}}_{\mathrm{PV}})\\ \hfill -{\mathbf{P}}_{\mathrm{ref}}\\ \hfill {\mathbf{P}}_{\mathrm{ref}}\end{array}\right],$$ where all vectors and matrices in $\mathbf{d}$ and $\mathbf{D}$ are defined in the Appendix A.

Appendix C. Two Equivalent Optimization Problems

Consider the problem to maximize a function upper and lower bounded by M and m, respectively, and subject to logical rules as follows: (A7)$$\begin{array}{cc}\hfill max& g\left(x\right)=\rho f\left(x\right),\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\hfill \end{array}$$(A8)$$\begin{array}{cc}\hfill & m\le f\left(x\right)\le M,\hfill \end{array}$$(A9)$$\begin{array}{cc}\hfill & \rho =\left\{\begin{array}{ccccc}a& & \mathrm{if}& & f\left(x\right)\ge 0,\\ b& & \mathrm{if}& & f\left(x\right)<0,\end{array}\right.\hfill \end{array}$$ where $f\left(x\right)=\mathcal{A}x+c$ and $m\le 0$. In this appendix, we show that the optimization in Equations (A7)–(A9) can be rewritten as an MILP. Following [30,31], the condition $f\left(x\right)\ge 0$ can be associated to a binary variable $\delta \in \{0,1\}$ such that:(A10)$$\begin{array}{cc}\hfill & m\le f\left(x\right)\le M,f\left(x\right)\ge 0,\iff \delta =1,\hfill \end{array}$$ then, Equation (A10) can be expressed by the following inequalities:(A11)$$\begin{array}{cc}\hfill & f\left(x\right)\ge m(1-\delta ),\hfill \end{array}$$(A12)$$\begin{array}{cc}\hfill & f\left(x\right)+ϵ\le (M+ϵ)\delta ,\hfill \end{array}$$ where $ϵ$ is a small positive scalar that represents a tolerance typically chosen as the machine precision. Then, the function $g\left(x\right)$ can be rewritten as:(A13)$$\begin{array}{c}\hfill g\left(x\right)=(a-b)\delta f\left(x\right)+bf\left(x\right).\end{array}$$

In fact, $\delta =1\Rightarrow g\left(x\right)=af\left(x\right)$ and $\delta =0\Rightarrow g\left(x\right)=bf\left(x\right)$, which is equivalent to Equations (A7)–(A9). In order to deal with the product of a continuous and a logical variable $\delta f\left(x\right)$ in Equation (A13), a new variable $z=\delta f\left(x\right)$ is introduced (see, e.g., [32]). This new variable can be rewritten as follows:(A14)$$\begin{array}{cc}\hfill & z\le M\delta ,\hfill \end{array}$$(A15)$$\begin{array}{cc}\hfill & z\ge m\delta ,\hfill \end{array}$$(A16)$$\begin{array}{cc}\hfill & z\le f\left(x\right)-m(1-\delta ),\hfill \end{array}$$(A17)$$\begin{array}{cc}\hfill & z\ge f\left(x\right)-M(1-\delta ).\hfill \end{array}$$

In fact, $\delta =1\Rightarrow (m\le z\le M,z=f(x\left)\right)$ and $\delta =0\Rightarrow (z=0,m\le f(x)\le M)$. Finally, the problem in Equations (A7)–(A9) can be written as (A18)$$\begin{array}{cc}\hfill max& g\left(x\right)=(a-b)z+bf\left(x\right),\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\hfill \end{array}$$ (A19)$$\begin{array}{cc}\hfill & z\le M\delta ,\hfill \end{array}$$ (A20)$$\begin{array}{cc}\hfill & z\ge m\delta ,\hfill \end{array}$$ (A21)$$\begin{array}{cc}\hfill & z\le f\left(x\right)-m(1-\delta ),\hfill \end{array}$$ (A22)$$\begin{array}{cc}\hfill & z\ge f\left(x\right)-M(1-\delta ),\hfill \end{array}$$ (A23)$$\begin{array}{cc}\hfill & f\left(x\right)\ge m(1-\delta ),\hfill \end{array}$$ (A24)$$\begin{array}{cc}\hfill & f\left(x\right)+ϵ\le (M+ϵ)\delta .\hfill \end{array}$$

This optimization problem is an MILP.

Appendix D. Vectors and Matrices for the MPC-MILP Solution

In this section, the following vectors and matrices are introduced to solve the MPC-MILP problem given in Equations (33) and (34). Considering Equation (33), the vector $\mathbf{f}$ is given by (A25)$$\mathbf{f}=\left[\begin{array}{c}{\mathbf{f}}_1\\ {\mathbf{f}}_2\\ \mathbf{0}\end{array}\right],\begin{array}{c}{\mathbf{f}}_1=\mathbf{\alpha }\odot {\mathbf{c}}_{\mathrm{n}}-{\mathbf{c}}_{\mathrm{soc}},\\ {\mathbf{f}}_2=\mathbf{\alpha }\odot ({\mathbf{c}}_{\mathrm{p}}-{\mathbf{c}}_{\mathrm{n}}),\end{array}$$ where ⊙ is the element-wise product operator, $\mathbf{\alpha }$, ${\mathbf{c}}_{\mathrm{p}}$, ${\mathbf{c}}_{\mathrm{n}}$ and ${\mathbf{c}}_{\mathrm{soc}}$ are defined in the Appendix A. The vector $\mathbf{x}$ in Equation (33) contains the continuous vectors ${\mathbf{P}}_{\mathrm{ES}}$ and $\mathbf{z}$, and the logical vector $\mathbf{\delta }$ that are given by (A26)$$\mathbf{x}=\left[\begin{array}{c}{\mathbf{P}}_{\mathrm{ES}}\\ \mathbf{z}\\ \mathbf{\delta }\end{array}\right],{\mathbf{P}}_{\mathrm{ES}}=\left[\begin{array}{c}{P}_t^{\mathrm{ES}}\\ P_{t+1}^{\mathrm{ES}}\\ ⋮\\ P_{t+N}^{\mathrm{ES}}\end{array}\right],\mathbf{z}=\left[\begin{array}{c}{z}_t\\ z_{t+1}\\ ⋮\\ z_{t+N}\end{array}\right],\mathbf{\delta }=\left[\begin{array}{c}{\delta }_t\\ {\delta }_{t+1}\\ ⋮\\ {\delta }_{t+N}\end{array}\right].$$ where ${\mathbf{P}}_{\mathrm{ES}},\mathbf{z},\mathbf{\delta }\in {\mathbb{R}}^{N+1}$. Notice that some components in the variable $\mathbf{x}$ are reals and other are integers, which defines the set $\mathcal{X}$ in Equation (33). Considering Equation (34), the matrix inequality $\mathbf{A}$ and the vector $\mathbf{a}$ are given by (A27)$$\mathbf{A}=\left[\begin{array}{ccc}\hfill -\mathbf{I}& \mathbf{0}& -\mathrm{diag}\left\{{\mathbf{P}}_{\min }-\mathbf{R}\right\}\hfill \\ \hfill \mathbf{I}& \mathbf{0}& -\mathrm{diag}\left\{{\mathbf{P}}_{\max }-\mathbf{R}\right\}\hfill \\ \hfill \mathbf{0}& -\mathbf{I}& +\mathrm{diag}\left\{{\mathbf{P}}_{\min }-\mathbf{R}\right\}\hfill \\ \hfill \mathbf{0}& \mathbf{I}& -\mathrm{diag}\left\{{\mathbf{P}}_{\max }-\mathbf{R}\right\}\hfill \\ \hfill -\mathbf{I}& \mathbf{I}& -\mathrm{diag}\left\{{\mathbf{P}}_{\min }-\mathbf{R}\right\}\hfill \\ \hfill \mathbf{I}& -\mathbf{I}& +\mathrm{diag}\left\{{\mathbf{P}}_{\max }-\mathbf{R}\right\}\hfill \\ \hfill \mathbf{F}& \mathbf{0}& \mathbf{0}\hfill \\ \hfill -\mathbf{F}& \mathbf{0}& \mathbf{0}\hfill \end{array}\right],\mathbf{a}=\left[\begin{array}{c}\hfill -{\mathbf{P}}_{\min }\\ \hfill \mathbf{R}\\ \hfill \mathbf{0}\\ \hfill \mathbf{0}\\ \hfill -{\mathbf{P}}_{\min }\\ \hfill {\mathbf{P}}_{\max }\\ \hfill {\mathbf{E}}_{\max }\\ \hfill -{\mathbf{E}}_{\min }\end{array}\right].$$ All vectors and matrices in $\mathbf{a}$ and $\mathbf{A}$ are defined in the Appendix A, except $\mathbf{R}={\mathbf{P}}_{\mathrm{ref}}-{\widehat{\mathbf{P}}}_{\mathrm{PV}}$.

Appendix E. Vectors and Matrices for the MPC-MILP Solution with Rate-Based

In this section, the following vectors and matrices are introduced to solve the MPC-MILP problem given in Equations (58) and (59). Considering Equation (58), the vector $\mathbf{g}$ is given by (A28)$$\mathbf{g}=\left[\begin{array}{c}{\mathbf{g}}_1\\ {\mathbf{g}}_2\\ \mathbf{0}\\ {\mathbf{g}}_3\\ \mathbf{0}\end{array}\right],\begin{array}{c}{\mathbf{g}}_1=\mathbf{\alpha }\odot {\mathbf{c}}_{\mathrm{n}}-{\mathbf{c}}_{\mathrm{soc}}\odot {\mathbf{\eta }}_{\mathrm{c}},\\ {\mathbf{g}}_2=\mathbf{\alpha }\odot ({\mathbf{c}}_{\mathrm{p}}-{\mathbf{c}}_{\mathrm{n}}),\\ {\mathbf{g}}_3={\mathbf{c}}_{\mathrm{soc}}\odot ({\mathbf{\eta }}_{\mathrm{c}}-{\mathbf{\eta }}_{\mathrm{d}}),\end{array}{\mathbf{\eta }}_{\mathrm{c}}=\left[\begin{array}{c}{\eta }^{\mathrm{c}}\\ {\eta }^{\mathrm{c}}\\ ⋮\\ {\eta }^{\mathrm{c}}\end{array}\right],{\mathbf{\eta }}_{\mathrm{d}}=\left[\begin{array}{c}1/{\eta }^{\mathrm{d}}\\ 1/{\eta }^{\mathrm{d}}\\ ⋮\\ 1/{\eta }^{\mathrm{d}}\end{array}\right],$$ where ⊙ is the element-wise product operator, $\mathbf{\alpha }$, ${\mathbf{c}}_{\mathrm{p}}$, ${\mathbf{c}}_{\mathrm{n}}$ and ${\mathbf{c}}_{\mathrm{soc}}$ are defined in the Appendix A, and ${\mathbf{\eta }}_{\mathrm{c}},{\mathbf{\eta }}_{\mathrm{d}}\in {\mathbb{R}}^{N+1}$. The vector $\mathbf{y}$ in Equation (58) contains the continuous vectors ${\mathbf{P}}_{\mathrm{ES}}$, $\mathbf{z}$, $\mathbf{w}$ and the logical vectors $\mathbf{\delta }$ and $\mathbf{\beta }$ that are given by (A29)$$\mathbf{y}=\left[\begin{array}{c}{\mathbf{P}}_{\mathrm{ES}}\\ \mathbf{z}\\ \mathbf{\delta }\\ \mathbf{w}\\ \mathbf{\beta }\end{array}\right],\mathbf{w}=\left[\begin{array}{c}{w}_t\\ w_{t+1}\\ ⋮\\ w_{t+N}\end{array}\right],\mathbf{\beta }=\left[\begin{array}{c}{\beta }_t\\ {\beta }_{t+1}\\ ⋮\\ {\beta }_{t+N}\end{array}\right],$$ with ${\mathbf{P}}_{\mathrm{ES}}$, $\mathbf{z}$, and $\mathbf{\delta }$ given in Equation (A26) and $\mathbf{w},\mathbf{\beta }\in {\mathbb{R}}^{N+1}$. Notice that some components in the variable $\mathbf{y}$ are reals and other are integers, which defines the set $\mathcal{Y}$ in Equation (58). On the other hand, consider Equation (59), then the matrix inequality $\mathbf{B}$ and the vector $\mathbf{b}$ are given by (A30)$$\mathbf{B}=\left[\begin{array}{ccccc}\hfill -\mathbf{I}& \mathbf{0}& -\mathrm{diag}\left\{{\mathbf{P}}_{\min }-\mathbf{R}\right\}\hfill & \mathbf{0}& \hfill \mathbf{0}\\ \hfill \mathbf{I}& \mathbf{0}& -\mathrm{diag}\left\{{\mathbf{P}}_{\max }-\mathbf{R}\right\}\hfill & \mathbf{0}& \hfill \mathbf{0}\\ \hfill \mathbf{0}& -\mathbf{I}& +\mathrm{diag}\left\{{\mathbf{P}}_{\min }-\mathbf{R}\right\}\hfill & \mathbf{0}& \hfill \mathbf{0}\\ \hfill \mathbf{0}& \mathbf{I}& -\mathrm{diag}\left\{{\mathbf{P}}_{\max }-\mathbf{R}\right\}\hfill & \mathbf{0}& \hfill \mathbf{0}\\ \hfill -\mathbf{I}& \mathbf{I}& -\mathrm{diag}\left\{{\mathbf{P}}_{\min }-\mathbf{R}\right\}\hfill & \mathbf{0}& \hfill \mathbf{0}\\ \hfill \mathbf{I}& -\mathbf{I}& +\mathrm{diag}\left\{{\mathbf{P}}_{\max }-\mathbf{R}\right\}\hfill & \mathbf{0}& \hfill \mathbf{0}\\ \hfill -\mathbf{I}& \mathbf{0}& \hfill \mathbf{0}& \mathbf{0}& -\mathrm{diag}\left\{{\mathbf{P}}_{\min }\right\}\\ \hfill \mathbf{I}& \mathbf{0}& \hfill \mathbf{0}& \mathbf{0}& -\mathrm{diag}\left\{{\mathbf{P}}_{\max }\right\}\\ \hfill \mathbf{0}& \mathbf{0}& \hfill \mathbf{0}& -\mathbf{I}& +\mathrm{diag}\left\{{\mathbf{P}}_{\min }\right\}\\ \hfill \mathbf{0}& \mathbf{0}& \hfill \mathbf{0}& \mathbf{I}& -\mathrm{diag}\left\{{\mathbf{P}}_{\max }\right\}\\ \hfill -\mathbf{I}& \mathbf{0}& \hfill \mathbf{0}& \mathbf{I}& -\mathrm{diag}\left\{{\mathbf{P}}_{\min }\right\}\\ \hfill \mathbf{I}& \mathbf{0}& \hfill \mathbf{0}& -\mathbf{I}& +\mathrm{diag}\left\{{\mathbf{P}}_{\max }\right\}\\ \hfill {\mathbf{V}}_1& \mathbf{0}& \hfill \mathbf{0}& {\mathbf{V}}_2& \hfill \mathbf{0}\\ \hfill -{\mathbf{V}}_1& \mathbf{0}& \hfill \mathbf{0}& -{\mathbf{V}}_2& \hfill \mathbf{0}\\ \hfill \mathbf{0}& \mathbf{0}& \hfill \mathrm{diag}\left\{\mathbf{R}\right\}& \mathbf{0}& \hfill -\mathrm{diag}\left\{\mathbf{R}\right\}\end{array}\right],\mathbf{b}=\left[\begin{array}{c}\hfill -{\mathbf{P}}_{\min }\\ \hfill \mathbf{R}\\ \hfill \mathbf{0}\\ \hfill \mathbf{0}\\ \hfill -{\mathbf{P}}_{\min }\\ \hfill {\mathbf{P}}_{\max }\\ \hfill -{\mathbf{P}}_{\min }\\ \hfill \mathbf{0}\\ \hfill \mathbf{0}\\ \hfill \mathbf{0}\\ \hfill -{\mathbf{P}}_{\min }\\ \hfill {\mathbf{P}}_{\max }\\ \hfill {\mathbf{E}}_{\max }\\ \hfill -{\mathbf{E}}_{\min }\\ \hfill \mathbf{0}\end{array}\right],$$ where the matrices ${\mathbf{V}}_1$ and ${\mathbf{V}}_2$ can be computed following: ${\mathbf{V}}_1=\mathbf{F}{\mathbf{\eta }}_{\mathrm{c}}$ and ${\mathbf{V}}_2=\mathbf{F}({\mathbf{\eta }}_{\mathrm{d}}-{\mathbf{\eta }}_{\mathrm{c}})$ and the remainder terms in $\mathbf{b}$ and $\mathbf{B}$ are defined in the Appendix A, except $\mathbf{R}={\mathbf{P}}_{\mathrm{ref}}-{\widehat{\mathbf{P}}}_{\mathrm{PV}}$. Notice that, if ${\eta }^{\mathbf{d}}={\eta }^{\mathbf{c}}=1$, the optimization problem becomes the original problem in Section 2.4.

References

1. Cubukcu, M.; Gumus, H. Performance analysis of a grid-connected photovoltaic plant in eastern Turkey. Sustain. Energy Technol. Assess.; 2020; 39, 100724. [DOI: https://dx.doi.org/10.1016/j.seta.2020.100724]

2. Fuentealba, E.; Ferrada, P.; Araya, F.; Marzo, A.; Parrado, C.; Portillo, C. Photovoltaic performance and LCoE comparison at the coastal zone of the Atacama Desert, Chile. Energy Convers. Manag.; 2015; 95, pp. 181-186. [DOI: https://dx.doi.org/10.1016/j.enconman.2015.02.036]

3. Rebolledo, F.; Mendoza-Araya, P.; Carvajal, G.; Ramírez, R. Performance evaluation of different solar modules and mounting structures on an on-grid photovoltaic system in south-central Chile. Energy Sustain. Dev.; 2022; 68, pp. 65-75. [DOI: https://dx.doi.org/10.1016/j.esd.2022.02.003]

4. Mohammadi, K.; Naderi, M.; Saghafifar, M. Economic feasibility of developing grid-connected photovoltaic plants in the southern coast of Iran. Energy; 2018; 156, pp. 17-31. [DOI: https://dx.doi.org/10.1016/j.energy.2018.05.065]

5. Khan, B.; Mahela, O.P.; Alhelou, H.H.; Padmanaban, S. Deregulated Electricity Market: The Smart Grid Perspective; Apple Academic Press: Palm Bay, FL, USA, 2022.

6. Perez, E.; Beltran, H.; Aparicio, N.; Rodríguez, P. Predictive Power Control for PV Plants With Energy Storage. IEEE Trans. Sustain. Energy; 2012; 4, pp. 482-490. [DOI: https://dx.doi.org/10.1109/TSTE.2012.2210255]

7. Ando, R.; Ishii, H.; Hayashi, Y.; Zhu, G. A Planned Power Generation for Battery-Assisted Photovoltaic System Using Short-Term Forecast. IEEE Access; 2021; 9, pp. 125238-125246. [DOI: https://dx.doi.org/10.1109/ACCESS.2021.3110950]

8. ESIOS. Red eléctrica de España. Sistema de Información del Operador del Sistema. 2022; Available online: https://www.esios.ree.es (accessed on 30 August 2022).

9. Ullah, Z.; Arshad,; Hassanin, H. Modeling, Optimization, and Analysis of a Virtual Power Plant Demand Response Mechanism for the Internal Electricity Market Considering the Uncertainty of Renewable Energy Sources. Energies; 2022; 15, 5296. [DOI: https://dx.doi.org/10.3390/en15145296]

10. DiOrio, N.; Denholm, P.; Hobbs, W.B. A model for evaluating the configuration and dispatch of PV plus battery power plants. Appl. Energy; 2020; 262, 114465. [DOI: https://dx.doi.org/10.1016/j.apenergy.2019.114465]

11. Bullich-Massagué, E.; Cifuentes-García, F.; Glenny-Crende, I.; Cheah-Mañé, M.; Aragüés-Peñalba, M.; Díaz-González, F.; Gomis-Bellmunt, O. A review of energy storage technologies for large scale photovoltaic power plants. Appl. Energy; 2020; 274, 115213. [DOI: https://dx.doi.org/10.1016/j.apenergy.2020.115213]

12. Arsalis, A.; Papanastasiou, P.; Georghiou, G.E. A comparative review of lithium-ion battery and regenerative hydrogen fuel cell technologies for integration with photovoltaic applications. Renew. Energy; 2022; 191, pp. 943-960. [DOI: https://dx.doi.org/10.1016/j.renene.2022.04.075]

13. Ai, X.; Li, J.; Fang, J.; Yao, W.; Xie, H.; Cai, R.; Wen, J. Multi-time-scale coordinated ramp-rate control for photovoltaic plants and battery energy storage. IET Renew. Power Gener.; 2018; 12, pp. 1390-1397. [DOI: https://dx.doi.org/10.1049/iet-rpg.2018.5190]

14. Hu, Y.; Armada, M.; Jesús Sánchez, M. Potential utilization of battery energy storage systems (BESS) in the major European electricity markets. Appl. Energy; 2022; 322, 119512. [DOI: https://dx.doi.org/10.1016/j.apenergy.2022.119512]

15. Souza, L.; Rotella Junior, P.; Aquila, G.; Maheri, A. Multiobjective optimization of hybrid wind-photovoltaic plants with battery energy storage system: Current situation and possible regulatory changes. J. Energy Storage; 2022; 51, 104467.

16. Yang, Y.; Bremner, S.; Menictas, C.; Kay, M. A Mixed Receding Horizon Control Strategy for Battery Energy Storage System Scheduling in a Hybrid PV and Wind Power Plant with Different Forecast Techniques. Energies; 2019; 12, 2326. [DOI: https://dx.doi.org/10.3390/en12122326]

17. Sánchez de la Nieta, A.A.; Paterakis, N.G.; Gibescu, M. Participation of photovoltaic power producers in short-term electricity markets based on rescheduling and risk-hedging mapping. Appl. Energy; 2020; 266, 114741. [DOI: https://dx.doi.org/10.1016/j.apenergy.2020.114741]

18. Kim, H.; Kwon, W.; Kim, S.; Lee, D. MPC-Based Optimal Operation for a PV Farm With Dual ESSs Using Spectral Density Analysis of Market Signals. IEEE Access; 2020; 8, pp. 215457-215466. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.3041593]

19. Borrelli, F.; Bemporad, A.; Morari, M. Predictive Control for Linear and Hybrid Systems; Cambridge University Press: Cambridge, UK, 2017.

20. Rawlings, J.B.; Mayne, D.Q.; Diehl, M. Model Predictive Control: Theory, Computation, and Design; Nob Hill Publishing: Madison, WI, USA, 2017; Volume 2.

21. Parisio, A.; Rikos, E.; Glielmo, L. A Model Predictive Control Approach to Microgrid Operation Optimization. IEEE Trans. Control Syst. Technol.; 2014; 22, pp. 1813-1827. [DOI: https://dx.doi.org/10.1109/TCST.2013.2295737]

22. Kim, R.; Glick, M.B.; Olson, K.R.; Kim, Y. MILP-PSO combined optimization algorithm for an islanded microgrid scheduling with detailed battery ESS efficiency model and policy considerations. Energies; 2020; 13, 1898. [DOI: https://dx.doi.org/10.3390/en13081898]

23. Delfino, F.; Ferro, G.; Minciardi, R.; Robba, M.; Rossi, M.; Rossi, M. Identification and optimal control of an electrical storage system for microgrids with renewables. Sustain. Energy Grids Netw.; 2019; 17, 100183. [DOI: https://dx.doi.org/10.1016/j.segan.2018.100183]

24. Sigalo, M.B.; Pillai, A.C.; Das, S.; Abusara, M. An Energy Management System for the Control of Battery Storage in a Grid-Connected Microgrid Using Mixed Integer Linear Programming. Energies; 2021; 14, 6212. [DOI: https://dx.doi.org/10.3390/en14196212]

25. Achiluzzi, E.; Kobikrishna, K.; Sivabalan, A.; Sabillon, C.; Venkatesh, B. Optimal Asset Planning for Prosumers Considering Energy Storage and Photovoltaic (PV) Units: A Stochastic Approach. Energies; 2020; 13, 1813. [DOI: https://dx.doi.org/10.3390/en13071813]

26. PVGIS. Photovoltaic Geographic Information System. 2022; Available online: https://re.jrc.ec.europa.eu (accessed on 30 August 2022).

27. Foroozandeh, Z.; Ramos, S.; Soares, J.; Lezama, F.; Vale, Z.; Gomes, A.; L. Joench, R. A Mixed Binary Linear Programming Model for Optimal Energy Management of Smart Buildings. Energies; 2020; 13, 1719. [DOI: https://dx.doi.org/10.3390/en13071719]

28. Tobajas, J.; Garcia-Torres, F.; Roncero-Sánchez, P.; Vázquez, J.; Bellatreche, L.; Nieto, E. Resilience-oriented schedule of microgrids with hybrid energy storage system using model predictive control. Appl. Energy; 2022; 306, 118092. [DOI: https://dx.doi.org/10.1016/j.apenergy.2021.118092]

29. CPLEX. IBM ILOG CPLEX Optimization Studio. 2022; Available online: https://www.ibm.com/cl-es/products/ilog-cplex-optimization-studio (accessed on 30 August 2022).

30. Williams, H.P. Model Building in Mathematical Programming; John Wiley & Sons: Hoboken, NJ, USA, 2013.

31. Bemporad, A.; Morari, M. Control of systems integrating logic, dynamics, and constraints. Automatica; 1999; 35, pp. 407-427. [DOI: https://dx.doi.org/10.1016/S0005-1098(98)00178-2]

32. Glover, F. Improved Linear Integer Programming Formulations of Nonlinear Integer Problems. Manag. Sci.; 1975; 22, pp. 455-460. [DOI: https://dx.doi.org/10.1287/mnsc.22.4.455]