Hierarchical model predictive control-based electric vehicle fleet charging management

Hierarchical model predictive control-based electric vehicle fleet charging management ☆ Branimir ˇ Skugor * , Luka Grden , Joˇ sko Deur Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Zagreb, Croatia ARTICLE INFO Keywords: Electric vehicle fleet Optimal charging Renewable energy sources Model predictive control Hierarchical charging Mixed integer linear programming Dynamic programming ABSTRACT Due to the relatively long parking/grid-connection intervals, charging of electric vehicles is characterized by flexibility that could be exploited for different benefits such as charging cost minimization and better utilization of intermittent renewable energy sources. To this end, an optimal and predictive hierarchical charging management method of electric vehicle fleet, characterized by computational efficiency and good scalability, is proposed in the paper. The method relies on an aggregate electric vehicle fleet model, a model predictive control law that performs an on-line dynamic programming optimization of an aggregate charging power, and a heuristic algorithm that distributes the aggregate charging power to individual electric vehicles. The heuristic algorithm is set to prioritize electric vehicles with lower state-of-energy levels and sooner time-of-departure. The proposed charging management strategy is demonstrated for the case of virtually electrified delivery vehicle fleet of a local retail company and virtual electricity production from renewable energy sources, and it is verified against the offline globally optimal mixed integer linear programming benchmark and a baseline dumb-charging scheme in terms of charging cost, renewable energy utilization, and related optimization time execution. 1. Introduction Advanced information and communication technologies in combination with optimization tools facilitate coordinated charging of electric vehicle (EV) fleets in terms of grid load levelling [1], minimizing related charging costs [2] or/and power peaks [3], and maximizing renewable energy sources (RES) exploitation [4]. Offline charging optimizations are typically performed over longer time periods (months or years) to reveal techno-economic potential of conventional fleet electrification [5] or in the context of energy planning studies to assess impacts and opportunities related to EV charging on a larger scale (city or even country) [6]. For instance, different offline optimizations are conducted in [7] for large EV fleets to analyse potential trade-offs of reducing charging costs and excess of RES production. Offline optimizations are often employed for system configuration optimization, such as locations and sizing of photovoltaics plants and charging stations [8]. On the other hand, online optimization is employed within EV fleet charging management, performed in real-time under time-varying conditions and in the presence of different uncertainties. The charging optimization in this context is typically performed in a model predictive control (MPC) fashion [9], i.e., by repeatedly solving related optimization problem on receding time horizon in each sampling time step. To perform control trajectory optimization, a proper EV fleet dynamics model is needed. There are two distinctive modelling approaches: (i) distributed model, where each individual EV within the fleet is modelled separately [10], and (ii) aggregate model, where the whole EV fleet is modelled as a single aggregated battery characterized by time-varying input profiles [11]. The former yields better accuracy, but at the expense of increased computational complexity, i.e., due to large number of state and control variables depending on the number of EVs within the fleet. The latter is somewhat less accurate, as it represents a relatively crude approximation of EV fleet behaviour, but it is computationally very efficient and essentially insensitive to the number of EVs. This makes the aggregate modelling approach particularly appropriate in applications which involve numerous EVs. An alternative variant of aggregate modelling approach is a representation of EV fleet via aggregated lower and upper constraints of cumulative charging energy of individual EVs, where the lower ones are obtained by applying as-late-as-possible and the upper ones as-soon-as-possible charging while meeting target battery state of charge (SoC) values [12], ☆ This article is part of a special issue entitled: ‘SDEWES 2024_ECM’ published in Energy Conversion and Management. * Corresponding author. E-mail addresses: [email protected] (B. ˇ Skugor), [email protected] (L. Grden), [email protected] (J. Deur). Contents lists available at ScienceDirect Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman https://doi.org/10.1016/j.enconman.2025.120119 Energy Conversion and Management 342 (2025) 120119 Available online 4 July 2025 0196-8904/© 2025 Published by Elsevier Ltd. appropriate for machine learning-based prediction of EV fleet behaviour, charging optimization and demand response assessment [13]. For instance, such aggregate boundaries are predicted in [14] by using multiple linear regression model and used for optimal EV fleet charging. In the context of real-time charging management, authors in [10] propose a centralized event-based MPC law to handle real-time EV fleet charging, aiming at the charging cost minimization and tracking of predefined aggregate charging power profile (set, e.g., by a distribution system operator). To ensure IEC 61851 standard regarding compliant semi-continuous charging power (i.e., to be either zero or above some minimum charging power threshold), the resulting optimal control problem is formulated as a mixed integer linear program (MILP). A similar MPC-based charging method is applied in [15] as a basis of a smart home controller, aimed at load shifting of a household with smart appliances, storage unit, electric vehicles, and photovoltaic microgeneration. The MPC is also used in [16] to control power of electric energy storage system (ESS) based on forecasts of electric power demand of other users and RES production, with the aim to track a predefined grid power profile established in a day-ahead manner for grid efficiency. In [17], the ESS power is optimized to mitigate impact of EV charging at fast charging stations, i.e., to keep as low and smooth power flow as possible at the point of connection, because increased power fluctuations would lead to higher grid connection fees and charges. Stochasticity of EV arrivals and related charging power profiles is tackled therein by formulating a stochastic variant of the optimal control problem, being solved via Pontryagin Minimum Principle. Similarly, the MILP optimization is used in [18] to coordinate photovoltaic, battery ESS, and EV charging, and to optimally size ESS for maximal photovoltaics exploitation and profit. An appealing MPC-based charging approach is proposed in [19], where an adaptation of objective function weights is applied based on the current EVs status feedback. Namely, EVs having higher deviation of their current SoC values from the target values and lower remaining charging times are given higher priority via larger weights within the objective function. The main implication of this adaptation is in possibility of reducing MPC prediction horizon and thus improving related computational efficiency and scalability to larger EV fleets. Another promising aspect of EV charging is a bi-directional energy flow between the grid and EVs, known as vehicle-to-grid (V2G) functionality, which could further increase EV charging flexibility for various benefits. The main concern here arises due to faster battery degradation caused by elevated energy flow. To tackle this, authors in [20] included battery degradation cost within a charging optimization problem along with an energy cost concerning time-of-use tariffs, peak power prices, with a robust MILP algorithm employed to handle energy consumption uncertainties. A reference [21] similarly accounts for it within a two-stage charging, first determining charging profiles of individual EVs within a fleet and another performing charging scheduling over limited number of chargers, both tending to maximize remaining useful life of EV fleet batteries. Inherent uncertainties related to individual EVs, i.e., their connection and disconnection from the grid and respective energy levels, and RES production are often treated in literature by stochastic variant of MPC algorithm (SMPC). For instance, authors in [22] uses probabilistic one-day-ahead predictions of charging load, RES production, and electricity price and employs a scenario-based SMPC to tackle respective uncertainties. Authors in [23] similarly uses SMPC to manage uncertainties at fast charging stations accommodating RES production and ESS. As opposed to the above-outlined real-time charging strategies, optimizing the charging power directly at the level of individual EVs (single-level/distributed charging) and relying on distributed EV fleet model, an alternative hierarchical approach is proposed in [24]. The real-time MILP-based optimization of charging power is performed therein on the aggregate level jointly for EVs and static batteries with the aim to maximize local exploitation of energy produced by photovoltaics. The obtained optimal aggregate charging power is distributed over individual EVs by using a rule-based algorithm. Similar approach is addressed in [25], where charging powers of individual charging stations are optimized by distributed MPC to reduce energy cost while avoiding local grid overload. The hierarchical charging manner is employed on the level of individual charging stations, whose optimized aggregate charging powers are distributed in real-time over individual EVs connected therein, using a fuzzy logic accounting for charging urgency reflected in departure time and current SoC. In general, the hierarchical approach is characterized by the simplicity of implementation and invariance to the number of EVs within a fleet [24], thus leading to good scalability to relatively large EV fleets (e.g., for e-hubs applications involving thousands of EVs on the city level [26]). It represents an appealing solution for different subjects intermediating between EV fleets and grid operators (aggregators), e.g., offering grid ancillary services, participation in energy markets, and potentially overcoming privacy concerns of individual EV users [26]. However, the above studies do not establish a globally optimal benchmark, which would be obtained offline on the full-time horizon or online on receding or shrinking horizon. Thus, they do not include a strict and systematic assessment of the proposed charging management strategies. For instance, a suboptimal uncoordinated charging baseline is used in [24] as a reference assessment point. Furthermore, charging problem formulations are often too specific/narrow. For example, only the energy exchanged with grid is minimized and the problem is often restricted to linear fleet models and related constraints [24]. Finally, practical issues including computational efficiency, scalability to different fleet size, adaptability to different, generally nonlinear formulations, and robustness to prediction errors are usually not addressed in systematic and comprehensive manner. Apart from hierarchical approach, scalability issues encountered with large fleets can be alternatively mitigated by distributed charging approach, characterized by flexibility, modularity, data privacy preservation, and computational efficiency. For instance, authors in [27] apply an alternating direction method of multipliers for distributed EV fleet charging, performing optimization on three levels in parallel, minimizing the cost on the individual EV level, supressing power peaks at aggregator level, and providing voltage regulation on distribution system operator level. A deep reinforcement learning (DRL) is another powerful method widely applied nowadays, shown in [28] to successfully handle charge scheduling of large number of EVs for reduction of average queuing and chargers’ idling times [28]. It is applied in [29] as a part of federated learning framework, performing optimization EV charging profiles are calculated at the local level. A hierarchical EV fleet charging management method proposed in this paper is conceived to optimize charging power at two levels: (i) aggregate level, and (ii) distributed level of individual EVs, as illustrated in Fig. 1. It is motivated by the authors’ previous work on the offline EV fleet charging optimization [5], where it has been found that such computationally efficient approach can provide near globally optimal results. The charging power time profile on the aggregate level is optimized by MPC in a receding horizon manner, by using a simplified and numerically efficient aggregate battery-based EV fleet model and different optimization algorithms such as MILP and dynamic programming (DP). The obtained optimal aggregate charging power is then distributed over individual EVs in each time step by using a heuristic algorithm based on charging priorities, thus prioritizing EVs with a lower level of energy in battery and sooner time-of-departure. The proposed method is verified against the offline globally optimal MILP benchmark for the scenario of virtually electrified delivery vehicle fleet of a local retail company and two-tariff electricity price model; along with an alternative single-level MPC charging performing direct optimization of individual charging powers. The verification study includes demonstration of computational efficiency, scalability to different fleet sizes, and robustness to RES power production prediction errors. The main contributions of the paper are: (i) setting up a computationally efficient, optimal, robust, and scalable hierarchical EV fleet charging management framework, combining (a) MPC-based EV fleet B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 2 charging at the aggregate level, relying on DP (and LP/MILP) optimization and capable of handling nonlinear/nonconvex fleet models and constraints, and (b) carefully-tailored, near-optimal and computationally inexpensive heuristic algorithm for distribution of aggregate charging power over individual vehicles, and (ii) systematically assessing the proposed hierarchical charging management framework versus single-level MPC charging approach and a globally optimal single-level MILP-based charging benchmark obtained offline, including computational efficiency, sensitivity, and scalability analyses. The remaining part of the paper is organized as follows. Section 2 presents the distributed and aggregate EV fleet models employed. Section 3 deals with offline charging management optimization aimed at setting the globally optimal benchmark. Section 4 elaborates on online charging methods including high-level MPC and low-level heuristic charging, while Section 5 presents simulation results. Concluding remarks are given in Section 6. 2. Electric vehicle fleet models Two types of EV fleet models adopted from [11] are considered: (i) aggregate, and (ii) distributed ones; where the former considers all EVs within fleet as a single aggregated battery for the purpose of aggregatelevel charging optimization/control (see Fig. 1), while the latter models each EV battery separately for single-level optimization and simulation. The batteries are modelled as energy storages with the state-of-energy (SoE) and the charging power as their state and control variables, respectively. 2.1. Aggregate electric vehicle fleet model Dynamics of the aggregate EV fleet model is described by the following state equation [11]: SoEagg(k+1) = SoEagg(k) + SoEin,avg(k)nin(k) Nv − SoEout,avg(k)nout(k) Nv + η ch Pc,agg(k)ΔT NvEmax,ind ,(1) where k is the discrete time step (k=0,1,⋯,Nt−1), Nt is the total number of time steps, SoEin,avg and SoEout,avg are average SoE values of EVs connecting to the grid and disconnecting from the grid within k th step, respectively, with the corresponding number of EVs denoted by nin and nout, respectively, Nv is the total number of EVs within a fleet, Pc,agg is the aggregate charging power, Emax,ind is the energy capacity of the individual battery (expressed in Wh; NvEmax,ind is the energy capacity of all batteries within the fleet), and ΔT is the time step (expressed in hours; here set to ΔT =0.25 h which corresponds to 15 min). The aggregate SoE state variable, SoEagg, is defined as normalized average energy of connected EVs: SoEagg(k) = ∑Nv i=1Ec,i(k) NvEmax,ind ,(2) where Ec,i is the battery energy (in Wh) of i th EV, which equals the actual battery energy if EV is connected within the k th time step, while it is zero, otherwise. The lower limit on SoEagg is zero, while the upper limit is set to be dependent on the number of EVs connected to the grid (nc): 0≤SoEagg(k) ≤ nc(k)/Nv≤1.(3) The aggregate charging power is limited in the range from zero (only one-direction power flow is enabled, i.e., from a grid to EVs) to the charging power capacity of connected EVs: 0≤Pc,agg(k) ≤ nc(k)Pcmax,ind,(4) where Pcmax,ind is the maximum charging power of individual EV. Additionally, the aggregate charging power is limited by the fixed upper constraint: Pc,agg(k) ≤ Pc,agg,max,(5) to account for the grid power limit. Fig. 1. Concept of hierarchical EV fleet charging management framework. B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 3 2.2. Distributed electric vehicle fleet model for offline charging power optimization The model structure given by the state equation (1) may also be used to model the individual (i th ) EV battery within a distributed EV fleet model as: SoEi(k+1) = SoEi(k) + SoEin,i(k)nin,i(k) − SoEout,i(k)nout,i(k) + η ch Pc,i(k)ΔT Emax,ind , (6) with related constraints: 0≤SoEi(k) ≤ ncb,i(k),ncb,i∈ {0,1},(7) 0≤Pc,i(k) ≤ ncs,i(k)Pcmax,ind,ncs,i∈ [0,1],(8) where Nv from Eq. (1) is now set to 1 and thus omitted in Eq. (6), SoEin,i and SoEout,i are SoE values of the i th EV when it connects to and disconnects from the grid, respectively, and nin,i and nout,i are binary variables taking the value of 1 if connection/disconnection of i th EV takes place within the k th step, and 0, otherwise. The state variable SoEi is defined similarly to the definition of SoEagg in Eq. (2): SoEi(k) = Ec,i(k)/Emax,ind, where Ec,i equals zero if i th EV is disconnected. The variable nc from Eqs. (3) and (4) is replaced by ncb,i in Eq. (7) and by ncs,i in Eq. (8), where ncb,i represents the binary variable taking the value of 1 if the i th EV is connected within k th step (partially or fully), and 0, otherwise, while ncs,i represents a share of EV connection time within the k th step (e.g., ncs,i=0.1 means that a related EV was connected 10 % of time step duration ΔT). 2.3. Distributed electric vehicle fleet model for simulation study The model state equation (6), originally proposed in [11], is modified to strictly satisfy the lower SoE constraint in Eq. (7), 0 ≤SoEi(k), within the EV fleet simulation model: SoEi(k+1) = {SoEint,i(k),for nout,i(k) = 0, 0,for nout,i(k) = 1,(9) where SoEi at k +1 step takes an intermediate value SoEint,i if the EV is not disconnected at the k th step, while it equals 0, otherwise. The intermediate SoE value incorporates the SoE contributions brought by EV connection to the grid (SoEin,i) and charging with the power Pc,i (cf. Eq. (6)): SoEint,i(k) = SoEi(k) + SoEin,i(k)nin,i(k) + η ch Pc,i(k)ΔT Emax,ind .(10) The SoE on departure, SoEout,i,is updated in the (k +1) th step to SoEint,i(k)only if a new driving mission starts at the k th step (nout,i(k) = 1): SoEout,i(k+1) = {SoEout,i(k),for nout,i(k) = 0, SoEint,i(k),for nout,i(k) = 1.(11) On the other hand, the SoE of an EV arriving from a driving mission and connecting to the grid in the k th step, SoEin,i, i.e. when nin,i(k) = 1 holds, is calculated as a function of the SoE at previous departure SoE (i.e., SoEout,i(k)) and a travelled distance di(k)of that driving mission: SoEin,i(k) = {0,for nin,i(k) = 0, fSoE(SoEout,i(k),di(k)),for nin,i(k) = 1.(12) The upper constraints on individual charging powers are set to: Pc,max,i(k) = min(ncs,i(k)Pcmax,ind,1−SoEi(k) − SoEin,i(k)nin,i(k) η chΔTEmax,ind ), (13) where the first term within the operator min(.) corresponds to the upper constraint of Eq. (8), while the second one is to satisfy the upper SoE limit from Eq. (7) (derived from Eq. (10) with SoEint,i limited to 1; recall that the lower SoE limit from (7) is ensured through the modified state equation (9) and (10)). For the purpose of post-analysis, the SoE and charging power values of individual EVs from the distributed model can be aggregated for each time step k as: SoEagg(k) = ∑Nv i=1SoEi(k)/Nv and Pc,agg(k) = ∑Nv i=1Pc,i(k). The charging power can be supplied from the grid (Pg) or from the local renewable energy sources (RES; Pres), with the priority of charging being given to RES while covering the eventual power deficit from the grid: Pg(k) = {Pc,agg(k) − Pres(k),for Pc,agg(k) − Pres(k)>0, 0,otherwise.(14) 3. Offline charging management optimization This section presents a formulation of offline optimization problem for the charging management framework outlined in Fig. 1. The formulation differs for the cases of excluded and included production from RES. The offline optimization is aimed at providing a globally optimal solution to serve as a benchmark for verification of online charging strategies. 3.1. Problem formulation The main aim of EV fleet charging optimization is to minimize the total cost of energy drawn from the grid: Cbatt =∑ Nt−1 k=0 Cel(k)Pg(k)ΔT 1000 ,(15) where Cel(k)is the electricity unit price time profile (given in EUR/ kWh), and the term Pg(k)ΔT/1000 denotes the grid-supplied charging energy increment in the k th step (expressed in kWh). The problem is subject to SoE dynamics (1) and charging power inequality constraints (3)-(5) in the case of aggregate model; and (5)-(8) in the case of distributed model. Also, it is required that the final SoE values, i.e., SoEagg(Nt)in the case of aggregate model and SoEi(Nt),i∈ {1,2,⋯,Nv}, in the case of distributed model, are equal to a pre-determined target value SoEfinal, which is set to be equal to the initial SoE value: SoEfinal = SoEinit, to satisfy the charge sustaining condition. 3.2. No renewable energy sources included – Linear programming formulation In the case of no electricity production from RES, the above optimization problem is linear both in the cost function and constraints and can be represented in the general linear programming (LP) form [30]: min xcTx,(16) s.t.Ax ≤b,x∈Rnx,(17) B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 4 where cT is a cost vector and x is an optimization vector to be determined to minimize the total cost cTx, while satisfying linear inequality constraints contained within the matrix expression Ax ≤b. Considering the cost function (15), the vectors cT and x are defined as: cT= [Cel(0)Cel(1)⋯Cel(Nt−1)]ΔT/1000 (18) x=[Pg(0)Pg(1)⋯Pg(Nt−1)]T.(19) In the case of distributed model, the control vector to be optimized for each EV within the fleet, is: ui=[Pc,i(0)Pc,i(1)⋯Pc,i(Nt−1)]T,i∈ {1,2,⋯,Nv},(20) while in the case of aggregate model the control vector is: u=[Pc,agg(0)Pc,agg(1)⋯Pc,agg(Nt−1)]T,i∈ {1,2,⋯,Nv}.(21) Note that Pg(k) = Pc,agg(k) = ∑Nv i=1Pc,i(k)in the case of Pres(k) = 0,∀k (see Eq. (14)). 3.3. Renewable energy sources included – Mixed integer linear programming formulation In the case of having RES production in the optimal control problem formulation, the discontinuity in the form of conditional expression from Eq. (14) prevents the optimization to be solved with standard LP solvers. To incorporate such discontinuity, a MILP formulation is employed [31], which augments the LP formulation (16) and (17) with additional integer optimization variables, z∈Znz, used to decompose Eq. (14) into a set of logical equations: z(k) = 1→Pg(k) = − Pres(k) + ∑Nv i=1ui(k),(22) z(k) = 0→Pg(k) = 0,(23) where z= [z(0)z(1)⋯z(Nt−1)]T, and z(k)in k th step implies the cases for the grid power Pg(k): z(k) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1,for −Pres(k) + ∑Nv i=1ui(k) ≥ 0, 0,for −Pres(k) + ∑Nv i=1ui(k)<0. (24) Transformation of the logical relations (22)-(24) into equivalent set of inequality constraints appropriate for related MILP solvers is given in Appendix A. The above LP and MILP formulations are implemented within Matlab environment by using YALMIP optimization toolbox [32], and linprog and intlinprog solvers are used to solve the respective optimization problems. 3.4. Dynamic programming Solving the optimization problem using DP algorithm is beneficial for the general nonlinear/nonconvex EV fleet models and constraints [33], because it guarantees globally optimal solution for such a general problem [5]. However, it is limited only to the aggregate EV fleet model, having a single control (Pc,agg) and state variable (SoEagg), as its computational complexity prohibitively increases in the case of distributed model with N v state and N v control variables. The SoE constraints are accounted for within the DP formulation via soft constraints L(k)added to the cost function (15) as: J=∑ Nt−1 k=0 Cel(k)Pg(k)ΔT 1000 +L(k) ⏟⏞⏞⏟ F(k) ,(25) L(k) = Kg,1(SoEagg(k+1) − 1)H(SoEagg(k+1) − 1)+Kg,2(−SoEagg(k +1))H(−SoEagg(k+1))+Kg,3(SoEagg(k +1) − nc(k+1) Nv)H(SoEagg(k +1) − nc(k+1) Nv)+Kg,4H(SoEfinal −SoEagg(k+1))H(k−Nt+1), (26) where the function H(.) represents the Heaviside function defined as: H(z) = 0 for z<0 and H(z) = 1 for z≥0. Relative importance of the individual terms/constraints are given via related weighting factors K g,i , i =1,…,4, which are all set to high values to enforce constraint satisfaction if possible. The aggregate charging power constraints given by Eqs. (4) and (5) are implemented as hard constraints within the DP algorithm. 4. Online charging management This section describes an MPC framework for the real-time EV fleet charging, relying on presented LP/MILP/DP charging optimization algorithms. 4.1. Model predictive control The online EV fleet charging management optimization is based here on the receding horizon MPC framework (denoted as MPC-REC). The control variable optimization problem is solved within the MPC as in the case of offline optimization (Section 3), with the main difference that it is now run online over the receding horizon of length N p =96 which corresponds to one day period (24 h) for ΔT =0.25 h. The MPC cost function to be minimized is (cf. Eq. (15)): J=∑ Np−1 j=0 Cel(j|k)Pg(j|k)ΔT 1000 ,(27) where k denotes the current simulation time step, j is the time step on the prediction horizon (relative to the current step k;j=0,1,⋯,Np−1). Two variants of online EV fleet charging are considered: (i) singlelevel charging, where individual EV charging powers, Pc,i(j|k),∀i, are optimized directly over the prediction horizon, and (ii) hierarchical charging, where the aggregate charging power, Pc,agg(j|k), is optimized and then distributed over individual EVs (see Fig. 1). As in the case of offline optimization, the MPC optimization problem in the former approach incorporates EV fleet dynamics model (6) for prediction of SoE evolution over the prediction horizon, SoE and charging power constraints (5), (7), (8), while in the latter approach SoE dynamics (1) and constraints (3), (4), and (5) are incorporated. Since the final step is generally outside of the receding horizon, the requirement on the final SoE is omitted here (e.g., by setting related weighting factor to zero within the DP optimization, Kg,4=0, see Eq. (26)). The optimization provides a sequence of optimal charging power values Pc,i(j|k),∀i, (or Pc,agg(j|k), j=0,1,⋯,Np−1, in the case of aggregate model) and only the first element Pc,i(0|k),∀i, (Pc,agg(0|k)) is applied in the current, k th sampling step, while the remaining ones are discarded. Another MPC approach considered performs the optimization on a shrinking horizon (denoted as MPC-DIM), which gradually diminishes as time progresses towards the end time of a day (set to be in the early morning period when the transport system is typically at rest preparing for upcoming day activities). Thus, the time-varying length of MPC-DIM prediction horizon Np,dim(k)is set as: Np,dim(k) = Np−k+⌊k Np⌋Np,k= {0,⋯,Nt−1},(28) B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 5 where Np is the fixed horizon length (equal to Np =96), and ⌊x⌋ is a mathematical operator providing a nearest lower integer of a real number x. Note that a full prediction horizon of length Np rebuilds when a new day starts. MPC-DIM relies on the same optimization problem as MPC-REC, but with the final SoE condition included, i.e., Kg,4=108, as its final step is now contained within the prediction horizon. The MPC-DIM approach is deemed as a reasonable alternative option since the fleet driving schedules are planned offline one day ahead. Apart from that, the MPC-DIM is characterized by an improved computational efficiency since its prediction horizon length is shorter in average when compared to MPC-REC, and thus related optimization executes faster. 4.2. Preparation of model predictive control input distributions The following input time profiles of individual EVs denoted by the subscript i=1,2,⋯,Nv should be predicted over the prediction horizon j=0,1,⋯,Np−1: nin,i(j|k), nout,i(j|k), SoEin,i(j|k), and SoEout,i(j|k), to serve for the single-level MPC optimization, and as a basis for calculating the following input time profiles needed for the hierarchical aggregate-level MPC optimization: nin(j|k), nout(j|k), nc(j|k), SoEin,avg(j|k), and SoEout,avg(j|k). While the arrival and departing times of each EV, nin,i(j|k) and nout,i(j|k)may be predicted from the planned driving schedules, the SoE of the arriving EVs, SoEin,i(j|k), should be predicted by using a transport energy demand model (below denoted by fSoE(.)). To maximize the vehicle range and also to simplify the energy demand model, it may be assumed that EV batteries are always fully charged when disconnecting from the grid and departing, i.e., SoEout,i(j|k) = 1 when nout,i(j|k) = 1 [5,11]. However, it may happen that an EV is parked and connected to the grid for a relatively very short amount of time between two driving missions and cannot be fully charged under present charging power limit of Pcmax,ind. To satisfy the departure schedule, it disconnects from the charger before the battery is full, and eventually rely on fast charging on road (at depot or e-hub) if the energy charged is not high enough to cover the trip energy demand. Thus, SoEout,i profiles should be carefully prepared to have the maximal possible values of 1, if possible, while not violating the individual charging power limit. For this purpose, the distributed model (9)-(13) is evaluated over the prediction horizon (j=0,1,⋯,Np−1; in the recursive sense) for the scheduled profiles nin,i(j|k)and nout,i(j|k), known initial conditions: ncs,i(0|k), SoEi(0|k), SoEin,i(0|k), nin,i(0|k), and SoEout,i(0|k), and the consistently applied maximum charging power Pc,i(j|k) = Pcmax,ind,∀j. The SoE is saturated to 1 if being reached prior to vehicle departure. The obtained SoE values at departures are used as inputs for the transport demand model fSoE( • ) in Eq. (12) to predict the SoE at the next arrivals (i.e., return to depot) and connections to the grid (when nin,i(j|k) = 1). Additional time profiles needed for MPC optimization, are related to the electricity price Cel(j|k)and the RES power production Pres(j|k), which should be also predicted, e.g., by using machine learning techniques based on historical data and meteorological forecasts. 4.3. Distribution of aggregate charging power to individual vehicles In the case of hierarchical charging approach, the aggregate charging power Pc,agg(k), obtained by MPC in a k th time step, should be distributed to connected individual EVs. For this purpose, a simple rule-based algorithm is established which prioritizes to charge EVs with lower SoE and sooner departure time. The related procedure given below is iterative since saturation of individual charging power due to the upper limits (13) may inhibit one-shot aggregate power distribution. The procedure starts by calculating the lower and upper individual charging power limits, Pc,min,i(k)and Pc,max,i(k), where Pc,max,i(k)is given by Eq. (13), while Pc,min,i(k)is determined according to the requirement that each EV is targeted to have the maximum possible SoE (equal to 1) each time when disconnecting from the grid (leading to the maximum EV range). More specifically, Pc,min,i(k)is derived from Eq. (10) under the assumption that i th EV will be charged with the maximum power Pcmax,ind from the following (k +1) th time step until the end of connection time tc,i. Eq. (29) is solved for Pc,i(k)to get the minimum charging power Pc,min,i(k)in the k th step under which the i th EV battery can still be fully charged until departure: Pc,min0,i(k) = 1 ncs,i(k)ΔT(Emax,ind η ch (1−SoEi(k) − SoEin,i(k)nin,i(k))−(tc,i(k) −ncs,i(k)ΔT)Pcmax,ind ). (30) The upper charging power constraint (13) is set to have priority over the lower constraint (30), i.e., the maximum charging power constraint cannot be violated, while the SoE at departure can be lower than 1 if the battery cannot be fully charged due to short connection/parking time. To this end, the lower limit Pc,min,i(k)of each EV is saturated to Pc,max,i(k) as: Pc,min,i(k) = {Pc,min0,i(k),for Pc,min0,i(k) ≤ Pc,max,i(k), Pc,max,i(k),for Pc,min0,i(k)>Pc,max,i(k).(31) The individual charging power values are then initialized to their lower limit values: Pc,i(k) = {Pc,min,i(k),for Pc,min,i(k)>0, 0,otherwise.(32) They are rescaled by the factor Pc,agg,max/∑Nv i=1Pc,i(k)if their sum exceeds the allowed aggregate charging power Pc,agg,max given by Eq. (5) (i.e., if ∑Nv i=1Pc,i(k)>Pc,agg,max). The remained aggregate charging power is then calculated as: Pc,agg,r(k) = Pc,agg(k) − ∑ Nv i=1 Pc,i(k),(33) which is distributed over individual EVs according to shares pi(k), set to be proportional to the deviation of corresponding SoE from 1 (i.e., from being fully charged), and inversely proportional to the remaining connection time tc,i(k): pi(k) = 1−SoEi(k) − SoEin,i(k)nin,i(k) − η ch Pc,i(k)ΔT Emax,ind tc,i(k).(34) These shares are calculated only for EVs connected to the grid, ncb,i(k) = 1, for which the currently designated charging power values Pc,i(k)are 1=SoEi(k) + SoEin,i(k)nin,i(k) + η ch ncs,i(k)Pc,i(k)ΔT+(tc,i(k) − ncs,i(k)ΔT)Pcmax,ind Emax,ind .(29) B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 6 lower than the related maximum values Pc,max,i(k),Pc,i(k)<Pc,max,i(k)(i. e., those that can still accommodate additional charging power). For other EVs, it is set to zero, pi(k) = 0. Then, the calculated shares are normalized: pi(k) = ⎧ ⎪ ⎨ ⎪ ⎩ pi(k) ∑Nv i=1pi(k),if ∑Nv i=1pi(k)>0, 0,otherwise, ∀i,(35) and as such they are used for distributing the remained aggregate charging power Pc,agg,r(k): Pʹ c,i(k) = Pc,i(k) + pi(k)Pc,agg,r(k),∀i.(36) Pc,i(k) = min(Pʹ c,i(k),Pc,max,i(k)),∀i.(37) The distribution procedure represented by Eqs. (33)-(37) is iteratively repeated until the remained aggregate power Pc,agg,r(k)given by Eq. (33), which is yet to be distributed, is brought to zero, or all shares pi(k) become zero (∑Nv i=1pi(k) = 0 in Eq. (35), i.e., there are no EVs available for charging). The presented distribution algorithm can be applied both in an offline and online manner. In the offline case, the whole aggregate power sequence is obtained offline (e.g., by the DP optimization and the aggregate battery model) and then it is distributed over individual EVs step-by-step by using the distribution algorithm (no feedback present). In the online case, the distribution algorithm is performed after getting the optimal charging power Pc,agg(0|k)by executing the MPC algorithm in the actual, k th sampling step, and using it to determine the individual charging power values in the same sampling step (feedback is present through the MPC path). 4.4. Baseline (dumb) charging strategy A so-called dumb charging strategy is introduced to serve as a baseline for verification of the developed MPC charging strategies. Its idea is to charge connected EVs as soon as possible, without accounting for electricity price or production from RES. Thus, individual charging power values in each time step k are set to their maximum values Pc,max,i(k)given by Eq. (13) if not violating the upper limit on the aggregate charging power given by Eq. (5) (see the first condition below; note that charging power of non-connected EVs is zero); otherwise, they are set to values obtained by scaling down Pc,max,i(k)in the way that satisfies the aggregate power limit (second condition below): Pc,i(k) = ⎧ ⎪ ⎨ ⎪ ⎩ Pc,max,i(k),for ∑Nv i=1Pc,max,i(k) ≤ Pc,agg,max Pc,max,i(k)Pc,agg,max ∑Nv i=1Pc,max,i(k),otherwise. ,∀i. (38) 5. Results This section describes a case study used for demonstrating proposed charging management methods and provides related comparative results and insights. The implementation aspects of certain MPC approaches and optimization algorithms are analysed through respective average online execution times. 5.1. Case study description and parametrization of electric vehicle fleet models EV fleet models described in Section 2 are parameterized by using the data recorded for a delivery vehicle fleet of a local retail company [34]. The data were recorded for ten mid-size Diesel engine-propelled delivery trucks by using GPS/GPRS equipment over a three-month period (from September to November). The vehicles mission was to deliver cargo from a distribution centre (a depot; DC) to different sales centres. These trucks were virtually converted to extended range electric vehicles (EREV) with similar power and torque characteristics as in the real trucks [2,35]. EREVs (denoted as EVs hereafter for the sake of brevity) were used instead of pure battery electric vehicles (BEV) to overcome limited range of BEVs and, thus, to be able to cover all recorded driving missions (both shortand long-distance ones). It was assumed that (i) they are controlled in CD/CS manner (charge depleting/charge sustaining), meaning the vehicle operates in fullyelectric mode until its state of charge (SoC) drops to 30 %, when it switches to hybrid driving mode that keeps the SoC at 30 % [35], and (ii) their charging could take place only at the DC during their parking periods between two driving missions. The recorded GPS positions were used to reconstruct the time periods of the vehicles being located within the DC and, thus, available for charging. From this data, the following time profiles of EV fleet models from Section 2 are derived: nin, nout, nc, ncb,i, ncs,i, nin,i, nout,i, where i=1,2,⋯,Nv. In the presented case study, the first full-week data were considered, while the remaining data of the three-month period were employed in sensitivity analysis. The week was set to start at 5 a.m., when all the vehicles were parked within the DC. The fleet was found to perform a repetitive activity over workdays (from Monday to Friday) with the onroad peak occurring around 10 a.m., while a reduced activity was observed over weekend days [34]. To obtain the transport demand model, the selected EREV model has been simulated over a range of initial battery SoC values (further denoted as SoE, to be aligned with the presented EV fleet models) and synthetic driving cycles with different travelled distances d [36]. The obtained maps of SoE-at-destination (i.e., at return to DC; SoEin) and related fuel consumption Vf, shown in Fig. 2, are used for EV fleet model simulations (cf. Eq. (12)). Note that SoEin is around the CS value of 0.3 (30 %) for large travelled distances or low initial SoE, where it should be noted that the SoE is kept around this value at the cost of fuel consumption (see Fig. 2b). This transport model is used in combination with Eqs. (9)-(13) to derive SoE time profiles of individual EVs, SoEin,i(k)and SoEout,i(k), where SoEout,i(k)is brought to 1 when possible, i.e., unless limited by the individual charging power limit (see Subsection 4.2 for more details). These SoE profiles are prepared offline and, thus, not influenced by charging optimization. Other time profiles used relate to the two-tariff electricity price model (currently applied in Croatia, Fig. 4c) and electricity production from solar panels hypothetically installed on the DC roofs (Fig. 5c, [35]). 5.2. Results for case of no renewable energy sources considered Firstly, the aggregate battery model (1)-(5) is used as an EV fleet simulation model for conducting the offline DP optimization in the case of no electric power production from RES. The obtained aggregate charging power is then distributed to individual EVs, represented by the distributed vehicle fleet model (9)-(13), by using the heuristic distribution algorithm (Section 4). Fig. 3 shows the aggregate SoE and charging power time profiles prior to and after performing the power distribution (see Eqs. (29)-(37)). Evidently, the power distribution does not perturb the aggregate charging power profile significantly (cf. [11]). Certain discrepancies can be explained by inaccuracies of the aggregate battery model, which cannot fully capture the distributed model dynamics. The low discrepancy level is confirmed through relatively high correlation indices of the two power profiles and the two SoE profiles, which are determined by Matlab function corrcoef(.) and equal 0.78 and 0.91, respectively (the ideal/maximum value is 1). The discrepancy influences the charging cost increase by around 10 %. Different online charging management methods are then tested for the case of distributed EV fleet model (9)-(13) and the same case of no electric power production from RES. The results obtained offline B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 7 through single/distributed-level LP optimization (Section 3), denoted as LP-OFF-S (see the legend in Table 1 for description of abbreviations), represent the globally optimal benchmark. The MPC-obtained (postdistribution) aggregated SoE and charging power profiles, given in Fig. 4, closely align with those of the LP-OFF-S benchmark. MPC-DPDIM-H relying on diminishing prediction horizon provides somewhat different profiles, which is due to the additional constraint on the final SoE to be equal to 0.95 (at the end of each day). The SoE profile of the baseline (dumb) charging strategy differs significantly from other strategies’ profiles, as it often tends to its upper limit due to applying maximum power charging whenever possible (see Section 4). Its unawareness of electricity price is reflected in relatively high charging power in the periods of high electricity cost, as opposed to other approaches which push the charging power towards low-cost hours. The related numerical results are given in Table 1. To account for certain differences in the total electrical energy consumption for different charging strategies, specific charging costs (EUR/kWh) are contained, as well. The specific charging cost of the offline hierarchical DP charging (DP-OFF-H) is 3.9 % higher than that of single-level charging benchmark (LP-OFF-S), while the same charging cost excess metrics equal 2.5 % and 0.0 % for MPC (online) strategies MPC-DP-RECH and MPC-DP-DIM-H, respectively. The performance gain of MPC vs. offline DP approach can be attributed to a feedback effect incorporated by taking exact current SoE values within MPC in each sampling time step thus mitigating the aggregate model deficiencies, which is not present when performing distribution of DP-OFF aggregate charging power profile. The hierarchical MPC strategy relying on LP, MP-LP-RECH, provides similar performance as MPC-DP-REC-H. On the other hand, a single-level MPC relying on LP, MPC-LP-REC-S, delivers better performance, which is comparable to that of the benchmark LP-OFF-S. Its even slightly better performance in terms of specific cost (−0.3 %) can be explained by slight differences in the final SoE, the total fuel and energy consumptions. Finally, all methods are significantly better than the DUMB charging strategy, which ends up in around 10 % higher charging cost when compared to the benchmark LP-OFF-S. 5.3. Results for case of renewable energy sources considered When including the power production from RES, the offline DP optimization, DP-OFF, tends to shift charging closer to solar noon, when the RES production is around its peak (see Fig. 5). It is interesting to note that the (aggregate) battery is not fully charged at 5 a.m. unlike the case Fig. 2. Map-based transport demand model providing SoE-at-destination (SoE in ), i.e. when arriving to DC (a), and related fuel consumption (V f ) for driving missions of varying length when arriving to DC (b), and related fuel consumption (V f ) for driving missions of varying length d. Fig. 3. Comparative plots of aggregate SoE and charging power profiles obtained directly by DP-OFF (blue) and after applying distribution algorithm and aggregation (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 8 Fig. 4. Aggregated SoE (a), and charging power profiles (b) obtained by different charging approaches applied to distributed fleet model for (c) two-tariff electricity price model. Fig. 5. Comparative DP-OFF optimization results for cases of RES production and no RES production. B. ˇ Skugor et al. Energy Conversion and Management 342 (2025) 120119 9