1 Introduction↩︎

The decarbonization of the economy as a whole is a significant challenge for modern societies. In particular, the sustainable transformation of both the energy and transport sectors poses significant technical and cultural challenges [1]. Both transitions meet/couple in the population’s homes where electricity, mobility, or heat are needed. Thus, possible synergies (if they exist) between the three systems can be exploited to achieve the desired decarbonization, freedom, resiliency, and cost savings at the local or aggregated level [2]. The successful exploitation of such coupling thus needs to be carefully tailored and built into the design of modern multicarrier energy systems [2]–[9]. This necessarily leads to advanced ems that schedule and control the der [6], [9]–[12]. Thus, the ems needs to handle uncertain and variable forecasts (electric vehicles [13], solar), aging assets (batteries [14]–[18]), and user preferences.

To dispatch and operate such systems, the literature suggests mpc [6], [12], stochastic optimization [19], rl [6], [7], [9] and many others. Nevertheless, the basis of such advanced systems is usually a day-ahead plan or dispatch that schedules the power of the assets along the day [4], [10], [11], [20]. This planner is usually an optimization-based system that uses approximated deterministic forecasts of certain inputs to schedule the different assets. The decisions taken are then implemented (and modified if necessary) in real time. For a residential mces such as the one presented in Fig. 1, the assets can generate, consume, or store energy from/in different carriers (electricity, heat, and transport). Controllable assets (loads or generation) combined with ess enable flexibility from the building, and thus, their operation and modeling are crucial in reducing safety margins and unlocking additional cost savings.

To coordinate hess, the different technologies must be modeled accordingly, properly representing their different power limits, dynamics, and particularities. In particular battery ess technologies are limited by their aging [5], [10], [15]–[18], [21], (lib and ev), as well as by their availability for mobility (ev) [5], [13]. Unfortunately, although aging mechanisms have been studied and modeled [5], [14], [22]–[25], they have not been entirely incorporated into ems design of residential mces [5], [10].

Energy management system formulations of residential mces tend to simplify the models used to linear or quadratic models, overlooking most technology particularities [2], [6], [9], [26]. White-box models are usually reserved for HVAC systems, focusing on thermal comfort [20], [26], and not to the ess. Ceusters et al. [6], [7] only use first-order linear models for both bess and tess, neglecting any differences between their dynamics. Similarly, Ye et al. [9] does not mention any difference between storage systems nor include evs in their system. Other works, only focus on the electrical carrier [4], [10], [16] using different types of battery models, used to describe the key variables such as state-of-charge \(SoC\), terminal voltage \(v_t\), state-of-health \(SoH\). However, a possible link to the thermal carrier is not explored.

Aging models fall within two categories: empirical or physics-based (PB) [5], [17]. The first are the most widely used in the literature due to their simplicity. They are obtained by performing long standardized calendar and cycling ageing tests [5], [21]. Unfortunately, empirical degradation models only have interpolation capabilities, usually use non-linear equations, represent a limited number of operating conditions (average C-rate, minimum \(SoC\), etc.), are prone to overfitting and are chemistry dependent [5], [25]. On the other hand, PB models are built through first-principles and specialized tests to identify individual degradation mechanisms[21], [22]. They have extrapolation features, can be expressed in the state-space form (or similar), account for several cathode chemistries and represent a wide range of operating conditions. However, they are also non-linear and, in general, non-convex [21]–[24], [27]. The integration of PB aging models into the operation of bess has been recently studied at the bms level for standalone and ev applications [14]–[16], [28], [29] usually through control-oriented pbrom. To the author’s best knowledge, their integration into wider schemes, such as transmission grids, microgrids, buildings, or industry, has not been extensively researched yet.

The contribution of this paper is the integration of pbrom aging models into ems day-ahead planning algorithms in the context of residential mces. pbrom aging models give the ems enough information to:

• handle different cathode chemistries.

• seamlessly operate new and used batteries

• identify dominant degradation mechanisms.

Our planning algorithm is an optimization-based secondary controller that minimizes energy cost and battery aging. A schematic of the mces is presented in Fig. 1. The system is composed of spv, bess, ev, pei, hp, st, tess, grid connection and loads.

This paper is organized as follows: section 2 presents the problem and modelling framework, section 3 presents the algorithm design and models used; section 4 explains the simulator; section 5 describes our case studies and validation; finally section 6 presents the conclusions and shortcomings of this approach.

2 Modeling and Optimal Planner↩︎

The following section describes the ems models, following the umf by Powell [30]–[32]. For a given system size the objective is to handle the operation cost, which is composed of three parts: the net cost of energy from the grid \(C_{\textrm{grid}}\), the degradation cost of losing storage capacity \(C_{\textrm{loss}}\), and a penalty for not charging the ev \(p_{\textrm{SoCDep}}\). The grid cost and the degradation cost are cumulative objectives because the goal is to optimize them through time, while the penalty for not charging the ev to the desired \(SoC\) is only a point reward at departure times \(t_{\textrm{dep}}\). The sdp is then:

\[\label{eq:sdp} \begin{align} \min_{P^*_{a,t}} \quad & \mathbb{E}_{W}[ C_{\textrm{grid}}+C_{\textrm{loss}}+p_{\textrm{SoCDep}}]\\ \textrm{s.t.} \quad & S_{a,t+1}=S^M_a (S_{a,t}, P^*_{a,t}, W_{t+1}|\theta_{a,t}) \\ & P_{a,t}^*=X^{\pi}_t(S_{a,t}) \in \mathcal{P} & \forall a \in \mathbb{A} \\ & S_{a,t} \in \mathcal{S} & \forall a \in \mathbb{A} \end{align}\tag{1}\]

\[\mathbb{A}=\{\text{SPV}, \text{grid}, \text{EV}, \text{BESS}, \text{HP}, \text{ST}, \text{TESS}\}\,.\] where the components of the objective are: \[\label{eq:cgrid} C_{\textrm{grid}}=w_{\text{grid}}\sum_{t=0}^{T}{(\lambda_{\textrm{buy},t} .P^{+}_{\textrm{g},t}-\lambda_{\textrm{sell},t} .P^{-}_{\textrm{g},t}).\Delta t}\tag{2}\]

\[C_{\text{loss}}=w_{\text{loss}}.c_{\text{loss}}.\sum_{t=0}^{T}\sum_{sa}{N_{s,sa} . N_{p,sa} . i_{\text{loss},as,t}.\Delta t},\;\forall\;sa \subset a,\]

\[p_{\textrm{SoCDep}}=w_{\textrm{SoC}}.||\varepsilon_{\textrm{SoC}, t_{\text{dep}}}||_2^2\]

The following definitions of the elements are considered:

• The state vector has 3 components, the physical state of the system \(R_t\), other (deterministic) information \(I_t\), and beliefs about uncertain quantities or parameters \(B_t\). All the observable physical quantities of our system, such as currents, voltages, and so on are included in \(R_t\). Other non-physical information, like prices are included in \(I_t\). Finally, our belief state \(B_{a,t}\) is composed of forecast medians. These are defined as: \[S_{a,t}=[R_{a}, I_{a}, B_{a}]_t^T \,,\]

  \[B_{a,t}=[\tilde{P}_{\text{PV}}, \tilde{P}_{\text{ST}}, \gamma_{n_{\text{EV}}}, \tilde{P}_{\text{load}}^{\textrm{e}}, \tilde{P}_{\text{load}}^{\textrm{th}}]_t^T \,.\]

• The actions or decision variables are \[P^*_{a,t}=[P_{\text{EV}}, P_{\text{BESS}}, P_{\text{HP}}^{\textrm{e}}]_t^T\,.\]

• The superscripts \(\textrm{e}\) and \(\textrm{th}\) refer to electricity or thermal carriers. They are used when the subscript is the same.

• Both the actions and state vectors have upper and lower limits denoted as \(\overline{P}_{a,t}^*\), \(\underline{P}_{a,t}^*\), \(\overline{S}_{a,t}\), and \(\underline{S}_{a,t}\).

• All bidirectional powers, either actions or states, are modeled with their converter efficiency \(\eta_a\): \[\eta_a S^+_{a,t} - \frac{1}{\eta_a} S^-_{a,t} = S_{a,t} \,,\]

• The order of the subscript is "name, device, time index".

• Capital \(C\) denotes total cost in , lowercase \(c\) denotes unit cost and \(w\) indicates tunning/scaling weight.

The operation’s inputs are the prices \(\lambda_{\text{buy/sell},t}\), ev availability \(\gamma_{n_{\textrm{EV}},t}\), the solar power \(P_{\text{PV},t}\), the electric demand \(P_{\text{load},t}^{\text{e}}\), and the thermal demand \(P_{\text{load},t}^{\text{th}}\). All of them can be considered to be exogenous processes that are known after making a decision.

In Eq. 1 the planner or policy \(\pi\) wants to minimize the likelihood of the operational cost \(\mathbb{E}[C]\) under the exogenous information process \(W\). The problem at hand is a state-dependent problem in which our decisions \(P_{a,t}^*\) are based on the current \(S_{a,t}\), and influence future states \(S_{a,t+1}\) (and thus, future decisions). Given the focus on future states and decisions, lookahead policies appear as attractive candidates for solving this sdp. Policy design and models are presented in the following Section 3.

3 Policy design↩︎

As mentioned before, the sdp in Eq. 1 is a state-dependent problem where current states influence future decisions. As such dla policies are commonly used in the literature to solve these problems. Two common examples of this policy family are optimal control strategies and stochastic dual dynamic programming. For this work, we focus on day-ahead planning, which is a subset of optimal control where deterministic inputs (forecast medians in this case) are used to decide the actions for the incoming day. The process is shown in Fig. 2.

More specifically a dla policy based on the devices’ approximated dynamic models. In this way, the ems plans future actions based on approximate predictive models of devices. Approximation is then denoted with \(\tilde{}\) . The policy is then solving an approximated economic ocp of the form:

\[\label{eq:appOCP} \begin{align} \min_{\tilde{P}^*_{a,tt'} \in \mathcal{P}} \quad & \tilde{C}_{\textrm{grid}}+\tilde{C}_{\textrm{loss}}+\tilde{p}_{\textrm{SoCDep}}+\tilde{p}_{\textrm{TESS}}\\ \textrm{s.t.} \quad & \tilde{S}_{a,t,t'+1}=\tilde{S}^M_a \left(\tilde{R}_{a,tt'}, \tilde{I}_{a,tt'}, \tilde{P}^*_{a,tt'}, \tilde{B}_{a,tt'}\right) \\ & \tilde{S}_{a,t} \in \mathcal{S} & \forall a \in \mathbb{A} \end{align}\tag{3}\]

\[\tilde{p}_{\textrm{TESS}}=w_{\textrm{TESS}}.\max\left(0, \tilde{SoC}_{\textrm{TESS}} - \overline{SoC}_{\textrm{TESS}} \right)\]

In our policy the upper limit constraint of the TESS, \(\overline{SoC}_{\textrm{TESS}}\), is implemented as a soft constraint to avoid infeasibilities during initialization or feedback. The penalty in the objective steers the \(SoC_{\textrm{TESS},t}\) towards the feasible region when the weight \(w_{\textrm{TESS}}\) is high enough.

The deterministic optimization problem in Eq. 3 approximates the real stochastic one by using forecasts, stored in \(B_{a,t}\), and approximated models for the transition function \(\tilde{S}_{a,t}^M\). In this model, the time \(t\) is the time at which the dla policy is created and \(t'\) is the time inside the policy itself. Note the subtle difference between the approximated dynamics \(\tilde{S}_{a,t}^M\) and the real ones \(S_{a,t}^M\). This is not to be overlooked because the assumption that the predictions done by the policy \(\pi\) hold true can lead to disappointing results in real-world applications. Making these distinctions early in design reveals important insights for future stages. In this work, the ema has an approximated model \(\tilde{S}_{a,t}^M\) to decide the setpoints \(\tilde{P}^*_{a,tt'}\) to be implemented in a simulator \(S_{a,t}^M\) containing detailed fidelity models. In the future, the simulator might as well grow enough to be considered a digital twin of the real building.

Thus the policy is: \[X^{\pi}_t(S_{a,t})= \arg \min_{P_{a,t}, ..., P_{a,t+H}} \tilde{C}_{\textrm{grid}}+\tilde{C}_{\textrm{loss}}+\tilde{p}_{\textrm{SoCDep}}+\tilde{p}_{\textrm{TESS}}\]

Notationally, we define a sequence matrix containing power setpoints and states from \(t\) to \(t+H\) with:

\[\mathcal{P}_{a,[t,t+H]}=[P_{a,t}, P_{a,t+1}, ..., P_{a,t+H}]\] \[\mathcal{S}_{a,[t,t+H]}=[S_{a,t}, S_{a,t+1}, ..., S_{a,t+H}]\]

The thermal assets are modelled in a linear way: \[P_{\text{HP}}^{\textrm{th}}=\eta_{\text{HP}}.P_{\text{HP}}^{\textrm{e}} \,,\] \[\tilde{P}_{\text{ST}}=\eta_{\text{ST}}.\tilde{P}_{\text{PV}} \,,\] \[\label{eq:socTESS} SoC_{\text{TESS},t+1}=SoC_{\text{TESS},t}-\frac{\Delta t}{Q_{\text{TESS}}.3600}.\eta_{\text{TESS}}.P_{\text{TESS},t} \,,\tag{4}\] \[\label{eq:thBalance} \tilde{P}_{\text{ST}}+P_{\text{HP}}^{\textrm{th}}+P_{\text{TESS}}=\tilde{P}_{\text{load}}^{\textrm{th}} \,.\tag{5}\] \[\label{eq:powerBalance} \tilde{P}_{\text{PV}}+P_{\text{BESS}}+ \gamma_{\text{EV}}.P_{\text{EV}}+P_{\text{grid}}=\tilde{P}_{\text{load}}^{\textrm{e}}+P_{\text{HP}}^{\textrm{e}} \,.\tag{6}\]

3.1 Batteries↩︎

The remaining devices in the mces are all battery-based ess. Batteries have complex nonlinear dynamics, and several modeling techniques are presented in the literature [17]. In this work, models coming from empirical and physics-based approaches are used. The modeling is divided into two different sub-models: performance and aging. Under the umf, this is represented in the transition function \(\tilde{S}^M_{sa,t}(\tilde{S}_{sa,t},x_{sa,t}|\theta_t)\), which contains both the perf. model \(p_{sa,t}^M(.)\) and the aging model \(d_{sa,t}^M(.)\). The performance model predicts stored energy \(SoC_{sa,t}\) and terminal voltage \(v_{t,sa, t}\). The aging model is used to update the parameters of \(p^M_{sa,t}(.)\), as shown in Fig. 3.

The perf. model is then: \[S_{sa,t+1}=p^M_{sa,t}\left(S_{sa,t}, P_{sa,t}, B_{a,t}|\theta_{sa,t}\right)\]

The aging model \(d_{sa,t}^M(.)\) is a set of equations that describe the dynamics of the performance parameters \(\theta_{sa,t}\). \[\theta_{sa,t+1}=d_{sa,t}^M\left(S_{sa,\;t}, P_{sa,t}, B_{a,t}, \theta_{sa,t}\right)\]

Finally, a terminal constraint is implemented to ease up feasibility and mitigate symmetries in the OCP of Eq. 3 , as in: \[SoC_{\textrm{BESS},t_1} = SoC_{\textrm{BESS},t_1+T}\]

3.1.1 Performance models \(p_{sa,t}^M\)↩︎

For the performance submodel, two alternatives have been implemented: a simple bm and a first order ecm. A basic bm of the operation of a battery assumes that its output voltage \(v_t\) is linear with the state of charge \(SoC\), assuming no voltage drop. Hence the only equations of this model are \[\label{eq:soc} SoC_{sa,t+1}=SoC_{sa,t}-\frac{\Delta t}{Q_{sa,t}.3600}.\eta_c.i_{sa,t} \,,\tag{7}\] \[\label{eq:ibranch} i_{sa,t}=\frac{P_{sa,t}}{v_{t,sa,t}.N_{s,sa}.N_{p,sa}} \,,\tag{8}\] \[\label{eq:ocv} OCV_{sa,t}=a_{OCV,sa}+b_{OCV,sa}.SoC_{sa,t} \,,\tag{9}\] \[\label{eq:vtBM} v_{t,sa,t}=OCV_{sa,t} \,,\tag{10}\] \[S_{sa,t}=[SoC_{sa}, v_{t,sa}, i_{sa}]^T_t\]

A first-order ecm has improved accuracy due to the incorporation of diffusion and series resistance, Fig. 4. The performance sub-model \(p_{sa,t}^M(.)\) is then modified by adding the equation: \[\label{eq:ecmRC} i_{R_1,sa,t+1}=e^{-\frac{\Delta t}{R_{1,sa}.C_{1,sa}}}.i_{R_1,sa,t}+ \left( 1-e^{-\frac{\Delta t}{R_{1,sa}.C_{1,sa}}} \right) .i_{sa,t}\tag{11}\] \[\label{eq:vtECM} v_{t,sa,t}=OCV_{sa,t}-i_{R_1,sa,t}.R_{1,sa}-i_{sa,t}.R_{0,sa} \,,\tag{12}\]

3.1.2 Degradation models \(d_{sa,t}^M(.)\)↩︎

For the aging models, the first alternative is an empirical sub-model presented by [25]. The empirical sub-model reduces all the degradation mechanisms into calendar and cyclic aging. \[\tag{13} \begin{align} \tag{14} & i_{\text{cycle},sa,t}=\frac{c_1.c_3}{c_4}.e^{c_2.|i_{sa,t}|}.(1-SoC_{sa,t}).|i_{sa,t}| \,, \\ \tag{15} & i_{\text{cal},sa,t}=c_5.e^{-\frac{\text{24 kJ}}{RT}}.\sqrt{t} \,, \\ & i_{\text{loss},sa,t}=i_{\text{cycle},sa,t}+i_{\text{cal},sa,t} \,, \end{align}\] \[\label{eq:Qdyn} Q_{sa,t+1}=Q_{sa,t} -\frac{\Delta t}{3600}.i_{\text{loss},sa,t} \,.\tag{16}\]

For the physics-based alternative, the reduced order model (pbrom) from [15] is used. It accounts for two degradation mechanisms: the sei and am. The author also presents a pbrom for Li-plating but its particular functional form prevents it from being implemented inside a MINLP solved by an interior-point solver.

The growth of the sei layer is modeled with a general reaction that aims to average all the different byproducts that compose the sei layer. This is synthesized in the reversible sei current \(i_{\textrm{SEI}}\): \[\label{eq:iSEI} i_{\textrm{SEI},sa,t} = \frac{k_{\textrm{SEI},sa}.e^{\frac{-E_{\textrm{SEI},sa}}{RT}}}{n_{\textrm{SEI}}.(1+\lambda_{sa}.\beta_{sa}).\sqrt{t}}\tag{17}\]

The system is completed with: \[\eta_{k,sa,t}=\frac{2.R.T}{F}.\text{sinh}^{-1}\left(\frac{i_{sa,t}}{n_{\textrm{SEI}}.a_s.A.L_n.i_0}\right)\]

\[z_{sa,t} = SoC_{sa,t}.(z_{100\%}-z_{0\%})+z_{0\%}\]

\[\beta_{sa} = e^{\frac{n_{\textrm{SEI}}.F}{R.T}.\left(\eta_{k,sa}+OCV_{n,sa,t}-OCV_s \right)}\]

The loss of active material due to the mechanical stress of the electrode is modeled with: \[i_{\textrm{AM},sa,t} = k_{\textrm{AM},sa}.e^{\frac{-E_{\textrm{AM},sa}}{R.T}}.SoC_{sa,t}.|i_{sa,t}|.Q_{sa,\;0}\] The total aging is the contribution of both mechanisms sei layer growth and am loss. The capacity fade current is: \[\label{eq:pbDeg} i_{loss,\;sa,t} = i_{\textrm{SEI},sa,t} + i_{\textrm{AM},sa,t}\tag{18}\]

Now, by carefully inspecting Eq. 12 , the reader will notice that if \(R_{0,sa,t}\) is incorporated as a variable in the OCP, Eq. 3 , this would add another non-convex constraint to it (since \(i_{sa,t}\) can be either positive or negative). Thus, its evolution is only included in the simulator \(S^M_{a,t}(.)\) updating the parameters without the policy \(X^{\pi}_t\) being directly aware of the process.

To model the power fade (i.e. the increase of \(R_0\)), the sei layer thickness \(\delta_{\textrm{SEI},\;sa,\;t}\) growth is described by: \[\delta_{\textrm{SEI},sa,t+1} = \delta_{\textrm{SEI},\;sa,\;t} + \frac{\Delta t}{M_{\textrm{SEI}}.n_{\textrm{SEI}}.F.\rho_{\textrm{SEI}}.A_n} i_{\textrm{SEI},sa,t}\] Hence the dynamics of the series resistance \(R_0\) are: \[R_{0,sa,t+1} = R_{0,\;sa,t} + \frac{\varepsilon_s}{\kappa_{eff}}.\frac{\Delta t}{M_{\textrm{SEI}}.n_{\textrm{SEI}}.F.\rho_{\textrm{SEI}}.A_n}i_{\textrm{SEI},sa,t}\] The solvent S leaves the electrolyte to form the SEI layer thus, the volume fraction of S evolves with: \[\varepsilon_{e,sa,t+1} = \varepsilon_{e,sa,t} - a_s. \frac{\Delta t}{M_{\textrm{SEI}}.n_{\textrm{SEI}}.F.\rho_{\textrm{SEI}}.A_n}i_{\textrm{SEI},\;sa,\;t}\]

3.2 Electric Vehicle↩︎

From the point of view of a residential building, the evs are a bess with availability constraints and certain requirements regarding their \(SoC\) at departure time \(t_{\text{dep}}\). For the availability \(\gamma\), the probability distributions of departure (\(t_{\text{dep}}\)) and arrival (\(t_{\text{arr}}\)) times can be described as random variables \(t_{\textrm{dep/arr}} \sim \mathcal{T_{\textrm{dep/arr}}}\), whose distributions \(\mathcal{T}_{\text{dep/arr}}\) are taken from Elaad [33]. The availability \(\gamma_t\) will then be: \[\gamma_t=\begin{cases} 0 & t \in [t_{\text{dep}};\;t_{\text{arr}}]\\ 1 & \text{otherwise} \end{cases} \,.\] The power balance of an ev is \[P_{\textrm{tot},\textrm{EV},t}=\gamma_{\textrm{EV},t}.P_{\textrm{EV},t}+(1-\gamma_{\textrm{EV},t})P_{\textrm{drive},\textrm{EV}}\]

At \(t_{\textrm{dep}}\) the ev is required to be delivered at \(SoC_{\textrm{dep}}^*\): \[SoC_{\textrm{EV}}(t_{\textrm{dep}})=SoC_{\textrm{dep}}^*\] This is implemented as a penalty in the objective function, Eq. 3 , as in any typical OCP. The deviation from the reference at the desired time is penalized with: \[\varepsilon_{SoC}=SoC_{\textrm{EV}}(t_{\textrm{dep}})-SoC_{\textrm{dep}}^*\]

4 MCES Simulator \(S_{a,t}^M(.)\)↩︎

The simulator is used to evaluate the policy \(\pi\) and close the loop with the state measurements. It is designed to:

• Provide high-accuracy simulation results that act as plant measurements.

• Adjust/reject setpoints that violate hard constraints.

• Re-balance power in case of rejections or infeasible optimizations.

The whole process is defined in Algorithm 5. First, the power setpoints must be adjusted for the grid and tess because the forecast used in \(X_t^{\pi}\) will never be the same as the actual exogenous inputs. Take a look at the balances, Eq. 5 and 6 , which contain the loads and solar generation. It is clear that \(\tilde{P} \neq P\) and a device must compensate for that difference. Thus, the simulator \(S_{a,t}^M\) recalculates: \[P_{\text{TESS},t}=P_{\text{load},t}^{\textrm{th}} - P_{\text{ST},t} -P_{\text{HP},t}^{\textrm{th}} \,.\] \[P_{\text{grid},t}=P_{\text{load},t}^{\textrm{e}}+P_{\text{HP},t}^{\textrm{e}} -P_{\text{PV},t} - P_{\text{BESS},t} - \gamma_{\text{EV},t}.P_{\text{EV},t} \,.\]

Second, once these powers have been adjusted, the simulator uses them to obtain the true/actual/fidelity state sequence \(\mathcal{S}_{a,[t,t+(H+\Delta t)/2]}\). For the tess it recalculates Eq. 4 . For the \(sa\), it uses LiiBRA.jl [34] to swiftly simulate pbroms of the performance of the battery [17], [35]. After that, the models from Jin [15] are used to calculate the true degradation outcome of the decisions \(P_{a,t}^*\). Again, the reader must remember that the capacity fade (decrease in \(Q_{sa,t}\)) is modeled in both the simulator \(S^M_{a,t}\) and the approximate model of the planner \(\tilde{S}_{a,t}^M\), whereas the power fade (increase in \(R_{0,as,t}\)) is only addressed in the simulator \(S^M_{a,t}\). Finally, if an action \(P_{sa,t}^*\) causes a future state to go out of bounds (\(S_{sa,t+1}\leq \underline{S}_{sa}\) or \(S_{sa,t+1}\geq \overline{S}_{sa}\)), the action is rejected and the \(sa\) remains in that state (either \(\overline{S}_{sa}\) or \(\underline{S}_{sa}\)) until the next day. Finally, the carriers are re-balanced if necessary.

The final state sequence is then fed back to the optimization-based planner. For practical implementation, in which the simulator is, in fact, an experimental setup, an online state observer is necessary to feed back the states to the ems. This is particularly important for the ess [17], [18], [36].

5 Case studies↩︎

The building has a grid connection with a smart meter with 15min resolution. The connection is also the physical link to the spot market in which the building participates. This is represented in the grid cost \(C_{\textrm{grid}}\) defined in Eq. 2 . The grid power \(P_{\textrm{grid}}\) is included in the state vector \(S_t\).

The system is composed of a 5kWp spv, a 20kWh bess with nmc or lfp cells, two 12.5kW ev charging points, a 4kWe heat pump, a 2.7kWth solar thermal collector, a 200kWh tess, a 6kWp electrical load, a 1kWp thermal demand, and 10kW LV grid connection. Power consumption profiles (\(P_{\textrm{load}}^{\textrm{e}}\)) were constructed for a year using data from 2021 to 2023 from the TU Delft’s Green Village smart meter data [37]. The output of the spv is taken from [38], [39], the market prices \(\lambda\) are taken from the EPEX day-ahead auction, and \(\lambda_{\textrm{buy}}=0.95\lambda_{\textrm{sell}}\) [40], and the heat demand \(P_{\textrm{load}}^{\textrm{th}}\) was modeled as [41].

The cells used are SANYO NCR18650 cells for nmc as in [15] and A123 cells for lfp. Their datasets where taken from PyBaMM [42] and LiiBRA [34]. To construct the ecm for both cell types sythetic/simulated cells were built and simulated in PyBaMM using standard 1D fom. The simulations followed the testing methodology from Plett, Chapter 2 Sections 2.9-.11 , [17]. Once the simulated profiles are ready the ecm parameters can be identified using subspace system identification as in [43]. The parameters can be found in 7, Table 5.

The OCP and simulations were modelled and run using Julia [44], JuMP [45], and InfiniteOpt [46]. The chosen solver was KNITRO from Artelys since it can handle MINLP formulations and NLP with complementarity constraints [47].

5.1 Case Study I: Day-ahead planning↩︎

To test and validate our ema the planner was simulated for two standard months (summer and winter) using 2023 data from the previously mentioned sources.

To quantify the impact of each performance and aging model 5 day-ahead planners were built:

• (BNoDeg) Including a bucket model and no degradation. \(\tilde{S}^{M1}_{a,t}\).

• (CEmpDeg) Including a first-order ecm and empirical aging for the \(sa\). \(\tilde{S}^{M2}_{a,t}\).

• (CPBDeg) Including first-order ecm with pbrom aging for the \(sa\). \(\tilde{S}^{M3}_{a,t}\).

  • (-w11), where weights are \(w_{\text{grid}}=w_{\text{loss}}=1\).

The simulation workflow is presented in Algorithm 7 and depicted in Fig. 2. First, the hyperparameters are initialized. This includes the time window to be optimized \(\mathcal{D}_t^{\pi}=[t;t+H]\), weights \(w\), number of days \(n_d\), user preferences, the initial state \(S_{a,0}\), etc. In our case for day-ahead planning, the time window is \(H=48\text{hr}-\Delta t\) and \(n_d=29\). At timestep \(t\) the OCP in Eq. 3 , is solved obtaining the actions sequence \(\mathcal{P}_{a,[t,t+(H+\Delta t)/2]}^*\). Together with the exogenous information sequence \(\mathcal{W}_{[t+1,t+(H+\Delta t)/2]}\) the actions are passed to the simulator \(S_{a,t}^M\) to get the feedback state sequence \(\mathcal{S}_{a,[t,t+(H+\Delta t)/2]}\). This feedback loop is repeated \(n_d\) times.

As a representative example Fig. 8 presents the results for the CPBDeg-w11 planner for a monthly period. It has the resulting power balances (electrical and thermal) and the use of the hess. The electric ess have daily cycles to minimize operating costs (energy arbitrage). This is particularly important for the ev since its mobility demand already establishes a daily periodicity. Thus, due to the ev’s battery pack size and its natural periodicity, it becomes the main electric storage of the system. This frees up the bess for energy arbitrage, which tries to capture prices variations when possible within the power balance. Thus price volatility incentivizes cycling. However, due to the sei model, batteries are also pushed downwards to the minimum aging state at \(\underline{SoC}_{sa,t}\). Thus the dispatch contemplates a trade-off between the 2 parts of the objective \(C_{\textrm{grid}}\) and \(Q_{\textrm{loss}}\). During the summer, price volatility is high with several hours with \(\lambda \leq 0\), the ems buys this energy to reduce costs. In winter, prices are less volatile and the load is higher leading to less opportunities for arbitrage and overall higher costs.

On the thermal side, the natural periodicity of the carrier is longer. Thus during the first week the initial \(SoC_{\textrm{TESS},0}\) influences costs greatly. After the first week the ems has already steered the buffer to its desired setpoint. This means a high \(SoC_{\textrm{TESS},t}\) during summer (high heat generation, low load) and a low setpoint during winter (low heat generation, high load). The high setpoint during summer entails 2 risks: overcharging the tess (i.e. activating soft-constraint) and not capturing negative prices/highly volatile prices due to past short-sight (tess starting a day with a high \(SoC\)).

The performance of the ess under the different planners is summarized in Figs. 9- 11. The first two present the 2D histogram between \(SoC\) and \(P\) of the electrical storages whereas the third presents the timeseries of the tess. In the histograms the left-hand side is for charge (\(P \leq 0\)) and the right-hand side is for discharge (\(P\geq 0\)). For the bess, Fig. 9, the BNoDeg and CEmpDeg cycle the battery pack more often, spreading towards the edges of \(P_{\textrm{BESS},t}\) and \(SoC\). The CPBDeg planner cycles less frequently, concentrating the operation near \(P=0\) and \(\underline{SoC}\) to reduce the ageing of the bess. This is true for both summer (top) and winter (bottom).

For the ev 2D histograms, Fig.10, the driving power \(P_{\textrm{drive},\textrm{EV}}\) is excluded and only the moments where the car is connected are analyzed (\(P_{\textrm{EV},t}.\gamma_{\textrm{EV},t}\)). The user’s mobility requirement leads to similar distributions for all the planners, with this denser distribution on the left half-plane for charging. Its shape can be explained by the fact that all planners try to delay the charging to extend the V2G as much as possible. V2G mode is less frequent in the BNoDeg/CEmpDeg than in the CPBDeg-w11. This occurs in both standard seasons. In summary, the operation pattern of the ev is dominated by its mobility demand. Finally, for the tess, Fig. 11 shows that in the summer the CPBDeg-w11 has the smallest st curtailment, because the chosen trajectory is the lowest of the three. BNoDeg and CEmpDeg chose almost the same trajectory for \(SoC_{\textrm{TESS},t}\).

The performance of each planner is summarized in Table 1. Starting with the grid cost \(C_{\textrm{grid}}\), the best performer is the CEmpDeg in the summer and BNoDeg in the winter. The differences between first and second are less than 1 . The worst performer is CPBDeg-w11. For the total capacity fade \(Q_{\textrm{loss}}\), CEmpDeg has the highest degradation in the summer and BNoDeg has the highest in the winter. The proposed aging aware CPBDeg achieves the lowest degradation in both seasons maintaining a reasonable \(C_{\textrm{grid}}\).

In general the planner that has the lowest \(C_{\textrm{grid}}\) is the one with the highest \(Q_{\textrm{loss}}\), but this is not always BNoDeg. The CEmpDeg fails to minimize the total capacity fade because of its model bias (\(S^M_{a,t} \neq \tilde{S}^M_{a,t}\)), in which it is rewarded to cycle at \(\overline{SoC}_{sa}\) and a calendar ageing independent of the \(SoC_{sa,t}\). Additionally the linear BNoDeg planner fails to fulfill its predictions because of the high number of rejected actions during summer. This is because when the BNoDeg plans for high-power discharges its predicted \(\tilde{SoC}_{sa,t}\) deviates from the true \(SoC_{sa,t}\) of the simulator \(S^M_{a,t}\), wrongly depleting the ess early (lower bound \(\underline{SoC}_{sa,t}\)). The linear BNoDeg is the planner with the highest percentage of rejects, roughly 45% of the total time in the summer.

However, when looking at the total storage usage, one must also analyze the number of cycles done by the \(sa\). Figure 12 presents the full equivalent cycles \(FEC\) over time \(t\), showing that the CPBDeg planners increase the ev usage \(FEC_{\textrm{EV}}\) at the expense of doing less cycles with the bess. When analyzing the capacity fade \(Q_{\textrm{loss}}\) against the full eq. cycles \(FEC\), Fig. 12, it is clear that the relative degradation per cycle (\(\frac{\partial Q_{\textrm{loss}}}{\partial FEC}\)) of the CPBDeg is smaller in the winter than the used benchmarks. In summer, the BNoDeg has the lowest degradation per cycle due to its high number of rejects. The use of the evs is the same for all planners, because the trajectory is mostly driven by user mobility. Finally, in both seasons CEmpDeg presents the highest \(\frac{\partial Q_{\textrm{loss}}}{\partial FEC}\). This appears to be as a risky strategy due to a lack of consistency across season and objectives (minimizing degradation or minimizing grid costs).

Finally, the computational time for the different strategies is presented in Fig. 13. Each sample is the total computational time it takes to solve Algorithm 7. Unexpectedly, BNoDeg has the lowest and most consistent \(t_{\textrm{comp}}\) distribution. However both CEmpDeg and CPBDeg planners have similar distributions, maintaining overall fast computational time, indicating that the modeling accuracy of CPBDeg is not prohibitively expensive.

5.2 Case Study II: Managing different cathodes↩︎

To demonstrate the PB models’ flexibility and extended capabilities, the CPBDeg scheduler is tested using two similar battery packs of the same rated capacity \(Q\) but using different cells. One is formed with nmc cells and the other with lfp. Since lfp cells have a lower rated capacity of \(Q_n=2.3\) Ah and a lower \(OCV\), the battery packs have more \(N_{s/p}\) to have roughly the same pack-rated capacity as its nmc counterparts. The power limits \(\overline{P}_{sa,t},\;\underline{P}_{sa,t}\) are also maintained to make an even comparison.

The PB aging models are suitable for both because they have graphite anodes [15], [23]. Nevertheless, they have different electrolytes. This is addressed by changing the electrolyte parameters in the model. Thus we use the same model equations but with different parameter values. This is a great advantage compared to the empirical fits presented in the literature. In the latter, the derived models are prone to overfitting to training conditions, delivering complex non-linear equations that can only be applied to specific chemistries, battery pack designs, and operating conditions. For our CPBDeg the physical parameters of the lfp cell used are taken from [24] and can be found in 7, Tables 4 and 5.

The simulation results are presented in Figs. 14 - 15 (a). Starting with the ev, Fig. 14 (a), the operation is similar except for a few days in summer and winter in which the lfp decides to have deeper discharges than its nmc counterpart. Moving forward to the bess, Fig. 14 (b), in summer the lfp scheduler is too ambitious generating early storage depletion roughly 40% of the time, 10% more than the nmc cells. During the winter total rejections are reduced to 25% and 10% for lfp and nmc respectively.

Continuing with the ageing analysis Fig. 15 (a) presents the \(FEC\) and \(Q_{\textrm{loss}}\) results. When looking at the ev, the the \(FEC_{n_{\textrm{EV}}}\) increases with the lfp cells. The opposite happens to the bess, which reduces its cycles in the lfp pack. In combination with its lower calendar ageing, represented in the parameter set, the lfp packs achieve lower degradation per eq. cycle \(\frac{\partial Q_{\textrm{loss}}}{\partial FEC}\) than their nmc counterpart. As such, the CPBDeg-LFP reduces bess total throughput (\(FEC\)) and capacity fade \(Q_{\textrm{loss}}\) while increasing ev throughput and reducing its degradation.

Table 2 presents the summary of performance for both cell types. Overall the CPBDeg-LFP achieve lower totals costs and ageing despite their high number of rejected actions. During winter \(C_{\textrm{grid}}\) stays the same with 37% less degradation. However, in the summer, just by changing the cell type from nmc to lfp \(C_{\textrm{grid}}\) is reduced by 25%, similar to BNoDeg in Section 5.1 1, and roughly 25% less capacity fade. This shows that CPBDeg rightly exploits its physical information of the system to achieve better performance.

5.3 Case Study III: Managing aged and fresh batteries↩︎

To demonstrate the flexibility and extended capabilities of the CPBDeg planner the scheduler is tested using two battery packs: the first is the fresh battery pack of nmc cells of Section 5.1 and the second is the same pack but with cells aged at \(SoH=90\%\). Only one benchmark is used: a BNoDeg with no \(SoH\) update. Thus, the BNoDeg ems sees a perfectly healthy cell with rated capacity when in reality the battery pack is aged 10%. The update is based on a 5% increase in series resistance \(R_{sa,0}\) and 10% decrease of the available Li content \(z_{100\%,sa}\).

Figure 15 (b) presents the degradation patterns for the 2 planners and different \(SoH\). In the aged battery, the share of calendar ageing (against the total) is much smaller and thus the percentual \(Q_{\textrm{loss}}\) is almost 35% smaller than in the new battery packs. In both seasons, the evs patterns are similar for both BNoDeg and CPBDeg. The change in ev trajectories between aged and fresh cells shows slight increase \(FEC\), due to smaller rated capacity, and a decrease in relative ageing, due to reduced calendar ageing. The impact on the bess is more pronounced. In winter, CPBDeg has smaller \(\frac{\partial Q_{\textrm{loss}}}{\partial FEC}\) (upper-right hand side of the graph) for both aged and new cells with BNoDeg doing more eq. cycles. In the winter CPBDeg still has a smaller \(\frac{\partial Q_{\textrm{loss}}}{\partial FEC}\) in the new cell but not in the used cell, where this is achieved by the benchmark controller.

However when adding the costs to the analysis, summarized in Table 3, it is noticed that:

• CPBDeg achieves lower capacity fade \(Q_{\textrm{loss}}\) than BNoDeg across all seasons and \(SoH\).

• CPBDeg improves its total grid cost \(C_{\textrm{grid}}\) as the cells degrade, in particular during summer.

• BNoDeg worsens its performance as the cells degrade. With higher \(C_{\textrm{grid}}\), \(Q_{\textrm{loss}}\) and model bias.

In summary, the performance of the proposed CPBDeg planner comparatively improves when using used cells, with respect to its linear BNoDeg benchmark.

6 Conclusions & Discussion↩︎

In summary, this paper presents an optimization-based day-ahead planner for residential multi-carrier energy systems that uses pbrom models to integrate battery ageing. The proposed CPBDeg planner can handle different cathode chemistries as well as batteries in different ageing states. The planner does this with lower degradation than the its benchmarks. This comes the expense of slightly higher computational times and grid cost.

In the first case study it is shown how advanced pbrom can be used to reduce battery ageing while maintaining a reasonable grid cost, in accordance to the literature for standalone utility-scale applications [21], [28], [29], [48], [49]. This is because the planner can control the degradation per full equivalent cycle better than the benchmarks.

In the second case study, the proposed planner is equipped with battery packs of different cathode chemistries and its performance is compared. This is out of the scope of most empirical degradation models. The lfp battery pack has a lower total \(Q_{\textrm{loss}}\) as per established knowledge and can achieve lower grid costs \(C_{\textrm{grid}}\) than its nmc counterpart. Even more so, when considering that the relative cost \(c_{\textrm{loss}}\) used for both packs was the same, when in reality lfp packs have a lower cost than nmc packs. Thus the 5.2 is a conservative approximation and lfp packs have the potential to enable even lower operation costs. This is achieved just by changing some simple parameter models taken from the literature. This is an important feature because it allows the ems to fully exploit any battery pack at hand. During the last case study it was shown how the proposed ema handles aged and new batteries seamlessly, even improving performance. Its linear counterpart however is unaware of the degradation, increasing the number of rejected actions. Summing up reducing the grid costs can be achieved by using LFP battery packs or modifying objective weights, the latter at the expense of higher degradation.

It is worth noting again that in this work only one optimization is done per day for all models (benchmarks and novel). If re-optimizations occur frequently as in mpc, model bias can be mitigated. Future works will seek to implement the proposed models in an economic MPC approach. Moreover, for physical setups the proposed physics-based approach requires a non-linear observer to identify internal states and a system identification algorithm to parametrize the models. Adaptive control techniques/ online learning techniques are crucial for scaling implementation. For the case of the bess initial tests and model identification can be done offline before start-up, and even offline during operation if historic data is continuously stored. However, for evs models, parametrization deems a challenge since a previously unknown car may appear or due to unknown driving conditions and profiles. Thus, the development of effective and accurate observers to identify and parametrize pbrom online automatically is crucial. This allows the ems predictions to be closer to reality, minimizing setpoint rejection by the bms [17], [18], [36].

Another limitation of the proposed ema is the modelling of the thermal carrier. In this paper, all thermal devices are assumed linear, and the horizon window \(H\) is only 48hs. Thus seasonal and monthly variations can not be handled properly [50] as in other studies, like long-term hydrothermal dispatch problems [51]. Future works will address this shortcoming and improve the thermal models for seasonal planning [20], [41].

7 Appendix: A↩︎

Graphite anode open-circuit voltage:

\[\begin{align} OCV_{n,sa,t}(z_{sa,t}) &= 0.6379 + 0.5416 \cdot e^{-305.5309 \cdot z_{sa,t}} \\ &\quad + 0.044 \cdot \tanh\left(-\frac{z_{sa,t} - 0.1958}{0.108}\right) \\ &\quad - 0.1978 \cdot \tanh\left(\frac{z_{sa,t} - 1.0571}{0.0854}\right) \\ &\quad - 0.6875 \cdot \tanh\left(\frac{z_{sa,t} + 0.0117}{0.0529}\right) \\ &\quad - 0.0175 \cdot \tanh\left(\frac{z_{sa,t} - 0.5692}{0.0875}\right) \end{align}\]

Table 4 presents the model parameters for the pbroms.

For the first-order ecms the model parameters are presented in Table 5.

References↩︎