*Article* **Comparative Assessment of Supervisory Control Algorithms for a Plug-In Hybrid Electric Vehicle**

**Nikolaos Aletras , Stylianos Doulgeris , Zissis Samaras and Leonidas Ntziachristos \***


#### **Highlights:**

#### **What are the main findings?**


**Abstract:** The study examines alternative on-board energy management system (EMS) supervisory control algorithms for plug-in hybrid electric vehicles. The optimum fuel consumption was sought between an equivalent consumption minimization strategy (ECMS) algorithm and a back-engineered commercial rule-based (RB) one, under different operating conditions. The RB algorithm was first validated with experimental data. A method to assess different algorithms under identical states of charge variations, vehicle distance travelled, and wheel power demand criteria is first demonstrated. Implementing this method to evaluate the two algorithms leads to fuel consumption corrections of up to 8%, compared to applying no correction. We argue that such a correction should always be used in relevant studies. Overall, results show that the ECMS algorithm leads to lower fuel consumption than the RB one in most driving conditions. The difference maximizes at low average speeds (<40 km/h), where the RB leads to more frequent low load engine operation. The two algorithms lead to fuel consumption differences of 3.4% over the WLTC, while the maximum difference of 24.2% was observed for a driving cycle with low average speed (18.4 km/h). Further to fuel consumption performance optimization, the ECMS algorithm also appears superior in terms of adaptability to different driving cycles.

**Keywords:** fuel consumption optimization; energy management system; hybrid vehicle control

### **1. Introduction**

Global warming due to increasing emissions of greenhouse gases (GHG) appears today as the main environmental pressure [1]. Transport is one of the key sources of manmade carbon dioxide (CO2) emissions [1–3]. This has led authorities around the world to set targets and take measures to reduce these emissions [3,4]. The European Union (EU) has set a target of reducing CO2 levels from new passenger cars by 37.5% by 2030, compared to 2021 [5]. Therefore, solutions such as electrified vehicles are promoted by the automotive industry to meet these targets [6].

Hybrid Electric Vehicles (HEVs) and Plug-in Hybrid Electric Vehicles (PHEVs) are currently the most widespread options for electrified vehicles in the market. HEVs and PHEVs have two independent energy sources, namely the battery and the fuel tank. However, only PHEVs can be charged directly by grid power and can cover substantial ranges (e.g., the latest models appear to have an electrical range of 100 km or 65 miles [7]) with electric

**Citation:** Aletras, N.; Doulgeris, S.; Samaras, Z.; Ntziachristos, L. Comparative Assessment of Supervisory Control Algorithms for a Plug-In Hybrid Electric Vehicle. *Energies* **2023**, *16*, 1497. https:// doi.org/10.3390/en16031497

Academic Editor: Wojciech Cieslik

Received: 29 December 2022 Revised: 16 January 2023 Accepted: 30 January 2023 Published: 2 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

power alone [6]. The share and mix between battery and engine power are constantly being decided during operation by an on-board energy management system (EMS) [1,6,8,9]. The EMS performance has a significant impact on fuel consumption (FC), which is directly linked to CO2 emissions. Therefore, EMS supervisory algorithms can be optimized to further decrease CO2 emissions from PHEVs.

Wu et al. [10] showed that there are a variety of principles for EMS algorithms. The main categories can be distinguished into rule-based (RB), optimization-based (OB), and learning-based (LB) ones. RB algorithms rely on a fixed set of rules, without a priori knowledge of driving conditions. OB algorithms are further split to offline or online ones. In online OB algorithms, such as the Equivalent Consumption Minimization Strategy (ECMS) [10–13], an instantaneous optimization is conducted based on current vehicle operation. In offline OB algorithms, a cost function is optimized for the complete driving cycle. Dynamic Programming (DP) [14,15] and Pontryagin's Minimum Principle (PMP) [10,16] are some of the common offline OB algorithms. LB algorithms are capable of instantaneously and in real time controlling and learning the optimal power split operations. Reinforcement Learning [17,18] and Artificial Neural Networks are some principles that are used in LB algorithms implementation [10].

There are only limited works in the literature on how different algorithm categories compare to each other under different operation conditions. Actually, Torreglosa et al. [19] mentioned that the optimization algorithms presented in the literature are seldom compared against commercial RB strategies. In their study, they compared different EMSs for HEVs with RB strategies using FASTSim, an open-source tool that includes validated HEV RB models. That analysis showed that optimum EMSs may provide fuel consumption benefits of 5% to 10%, compared to commercial RB EMSs. Wu et al. [20] proposed an optimization-based strategy that appeared to reduce the fuel consumption of a 2010 Toyota Prius hybrid by 3.5–6%, compared to an RB algorithm that was earlier published by Kim et al. [21]. Hwang et al. [22] applied particle swarm optimization to improve the fuel economy of a power split hybrid, and showed up to 9.4% improvement compared to an RB algorithm. This limited previous work showed that there are further margins to improve fuel consumption over commercially applied algorithms.

In assessing the performance of different algorithms, one needs to make sure that the exact operation profile is followed over computer simulations or real-world experimentation with the various EMS approaches. Fuel consumption differences of only a few percentage units, such as those expected when varying the EMS, can be observed only due to slight deviations of the original speed profile in consecutive simulations, e.g., due to underpowering accelerations. Moreover, it needs to be ensured that fuel consumption improvement is assessed under the same state of charge levels (SOC) to avoid part of the difference only being due to variation in the battery depletion levels, e.g., over consecutive simulations. Although such conditions may be self-evident, these are seldom if at all demonstrated in published EMS algorithm comparison studies.

The article focuses on the comparative assessment of commercial RB and ECMS based algorithms for a plug-in hybrid vehicle powertrain. The purpose is to examine if further FC reduction can be achieved by introducing an enhanced EMS over a commercial one. A method is first presented to compare fuel consumption over identical battery state of charge (SOC) levels, vehicle distance traveled and wheel power demand. The proposed method suggests novel corrections for the fuel energy consumption of the compared EMS algorithms. More specifically, the method introduces correction terms for the deviations in the final SOC values, propulsion energy and benefits from regenerative braking between the compared EMS algorithms. We propose that such a method needs to be used in all similar studies of fuel consumption comparison.

#### **2. Methods**

#### *2.1. Back-Engineered EMS Algorithm*

A vehicle simulator of a parallel P2 PHEV [23] has been built in the AVL Cruise simulation platform. Its overall performance has been validated with actual experimental data collected by tests on an actual vehicle in the chassis dynamometer of Aristotle University. The vehicle's technical specifications are listed in Table 1.

**Table 1.** Parallel P2 PHEV Technical Specifications.


Table 2 shows the tests conducted in the lab on the particular vehicle to understand the performance of its stock EMS algorithm. A Worldwide harmonized Light vehicles Test Cycle (WLTC) [24] and an ERMES cycle [25] have been used for the tests. The different cycles are distinguished into cold and hot start ones, depending on whether the start engine coolant temperature was lower than 35 ◦C or higher than 70 ◦C, respectively. A single case with intermediate start temperature is identified as warm start in Table 2.

**Table 2.** Driving cycles and specifications use for experimental validation of the back-engineered algorithm.


These experiments have been used to back-engineer the rules of the heuristic controller on-board the commercial vehicle. In this paper, the RB algorithm is a specialization of the general methodology described by Doulgeris et al. [26], while the specific controller algorithm is described in detail by Doulgeris et al. [27].

Figure 1 shows the flowchart of the RB algorithm. The engine switches on if the power demand, vehicle speed or acceleration, and SOC level are above specific thresholds. The power demand threshold for engine start depends on the SOC level. The engine always shuts off when the power demand becomes negative.

Figure 2 shows the engine power output decided by the algorithm curve, depending on the gear engaged (x-axis) and the vehicle speed (parameter). If the engine meets the criteria for switch-on according to Figure 1, then the engine power output is determined by Figure 2 depending on current vehicle speed and gear.

**Figure 1.** Stock rule-based algorithm extracted from experimental evidence.

**Figure 2.** Engine power output decided by the RB algorithm depending on gear and vehicle speed.

#### *2.2. Alternative Algorithm Description*

An Equivalent Consumption Minimization Strategy (ECMS) algorithm has been developed in the current study, as an alternative to the back-engineered RB one. Both of the algorithms—ECMS and the back-engineered one—are applied to the same vehicle simulator platform of parallel P2 PHEV that has been built in the AVL Cruise simulator platform. The ECMS algorithm aims at optimizing a predefined FC cost function for given operation conditions. The general cost function for fuel consumption optimization of a hybrid vehicle is given in Equation (1), where J is a performance index that needs to be minimized. The integral term represents the total fuel consumption over a complete driving profile, as it integrates the instantaneous FC ( . mfuel) from an initial (t0) to a terminal time stamp (tf). FC depends on the normalized engine load function u that ranges between 0 (engine shut off) to 1 (operation at full power—Pe,max), according to Equation (2). The

first term in Equation (1) is used as a soft constraint for the value of the SOC at the end of the cycle (SOCf). With the use of the ϕ function, SOC deviations from the final target value (SOCtarget) are being penalized. The choice of the SOCtarget value depends on the examined study. Usually, in PHEVs applications, the charge depletion is permitted because the battery can be charged from the electric power grid. As a result, for PHEVs applications, the SOCtarget can be lower than the initial value of SOC [8,9].

$$\mathbf{J} = \boldsymbol{\varphi} \left( \mathbf{SOC}\_{\mathbf{f}}, \mathbf{SOC}\_{\text{target}} \right) + \int\_{\mathbf{t}\_0}^{\mathbf{t}\_\mathbf{f}} \dot{\mathbf{m}}\_{\text{fuel}}(\mathbf{u}(\mathbf{t}), \mathbf{t}) \, \mathbf{d} \tag{1}$$

$$\mathfrak{u} = \frac{\mathcal{P}\_{\mathrm{e}}}{\mathcal{P}\_{\mathrm{e}, \mathrm{max}}} \tag{2}$$

The optimization of Equation (1) leads to the determination of the u for every second of the driving cycle, which in turn gives the instantaneous engine power output by means of Equation (2).

The ECMS optimization is subject to the conditions of the set of Equations (3)–(6). Equation (3) presents the power balance for the powertrain of a P2 parallel hybrid vehicle [28]. The sum of the demanded power at wheels (Preq,wheels) and the mechanical power losses (Pmech,losses ) must be equal to the mechanical power output from the main power units (electric motor power—Pem and engine power—Pe). If the vehicle velocity is known, then the Preq,wheels and Pmech,losses can be determined by a vehicle power-based model. We have set up a vehicle model in AVL Cruise for this purpose. So, with the use of Equation (3), the power output of the electric motor is determined.

$$\mathbf{P\_{req,wbeels}} + \mathbf{P\_{mech,losses}} = \mathbf{P\_e} + \mathbf{P\_{em}} \tag{3}$$

$$P\_{\mathbf{e}} \le P\_{\mathbf{e}, \max} \tag{4}$$

$$\mathcal{P}\_{\rm em,min} \le \mathcal{P}\_{\rm em} \le \mathcal{P}\_{\rm em,max} \tag{5}$$

$$\text{SOC}\_{\text{min}} \le \text{SOC} \le \text{SOC}\_{\text{max}} \tag{6}$$

Equations (4)–(6) correspond to the physical constraints of the powertrain components. Equation (4) suggests that the engine cannot overcome its full load curve, represented by the maximum engine power output (Pe,max). The electric motor is also limited by its full load curve, depending on whether it works as a motor (Pem > 0) or generator (Pem < 0) (Equation (5)), with corresponding limits given by Pem,max and Pem,min. Finally, the battery SOC cannot exceed a range of maximum (SOCmax) and minimum (SOCmin) levels for the purpose of maintaining battery life (Equation (6)).

Optimizing Equation (1) within the set of Conditions (3)–(6) is only possible when the operation mission is known a priori. In the real-world, a priori knowledge of the driving profile application is known only over in-lab tests and not for on-road driving. For on-road operation, the ECMS will have to be locally optimized according to the present driving conditions. Such local optimization is achieved by means of Equation (7), where the integral term of Equation (1) has been eliminated. In Equation (7), the engine fuel rate . mfuel and the battery power flow expressed in terms of an equivalent fuel rate . meq (Equation (8)—where QLHV is the fuel's lower heating value) result in an equivalent total fuel mass rate . mtot by means of the equivalence factor s (Equation (9)). The latter comprises the constant term s0 and a penalization term p(SOC) that depends on SOC. The s0 term is used as the main weighting factor of the . meq inside the cost function. The p penalizes deviations of the current SOC values from the target. The usage of the p term is similar to the one of the ϕ term in Equation (1). The difference is that the penalization is made for the instantaneous values of SOC instead of the SOC value at the end of the cycle, because the optimization is only carried out locally. The battery power (Pbatt) in Equation (8) can be positive for power outflux from the battery (Equation (10a)) and negative when the

EM acts as a generator that charges the battery (Equation (10b)), with ηbatt and ηem being the battery and electric motor efficiencies, respectively.

$$
\dot{\mathbf{m}}\_{\text{tot}} = \dot{\mathbf{m}}\_{\text{fuel}} + \mathbf{s} \times \dot{\mathbf{m}}\_{\text{eq}} \tag{7}
$$

$$
\dot{m}\_{\text{eq}} = \frac{\mathcal{P}\_{\text{batt}}}{\mathcal{Q}\_{\text{LHV}}} \tag{8}
$$

$$\mathbf{s} = \mathbf{s}\_0 + \mathbf{p}(\text{SOC}) \tag{9}$$

$$\mathbf{P\_{batt}} = \left(\frac{\mathbf{P\_{em}}}{\eta\_{\text{batt}} \times \eta\_{\text{em}}} \, | \mathbf{P\_{em}} \ge 0\right) \tag{10a}$$

$$\mathbf{P\_{batt}} = (\mathbf{P\_{em}} \times \boldsymbol{\eta\_{batt}} \times \boldsymbol{\eta\_{em}} \ | \mathbf{P\_{em}} < 0) \tag{10b}$$

The algorithm basically decides on the engine operation variable u(t) in Equation (2) that leads to the lowest total equivalent fuel mass ( . mtot). This procedure is repeated in every second of the complete mission profile. The algorithm takes into account two conditions regarding the potential battery charge or discharge. The first one is that a potential battery charge will lead to an SOC surplus, which can be utilized in the future. The second condition is that a present battery discharge generates a requirement for a future battery charge in order to retain the battery SOC within certain limits. An optimal solution can be guaranteed if the s term in Equation (7) is adapted appropriately (Equation (9)). In this way, although the s-by-s optimization cannot achieve as good a performance as the global optimal solution, it still produces a practical optimization solution that can be integrated in EMS without knowledge of the forthcoming driving profile.

Three alternative expressions for p have been examined in the current work (Equations (11a)–(11c)). In Type A expression, p is proportional to the difference of current SOC over a constant reference value SOCref (Equation (12a)). Therefore, this expression tries to keep SOC at a value close to the reference one over the complete driving profile—and it is tuned by a proportional term (kp). In Type B, the SOCref value varies with travelled distance D (Equation (12b) [29]). The algorithm also tries to keep the current SOC close to the SOCref value, as in Type A. More specifically, in Type B, the SOCref value starts with an initial value (SOCi) and then the SOCref decreases proportionally with the D until it reaches the target SOC value (SOCtarget). In Type B, the total driving distance (Dfinal) must be either known or estimated. Some research articles mention that this type of linear SOC trajectory with distance seems to be close to the global optimal solution [30,31]. In Type C expression, a specific SOC window is used for determining p [8] (Equation (11c)). More specifically, the *p* value dependents on SOC, a target value for the SOC (SOCtarget), selected maximum and minimum SOC values and a selected superscript for the penalization function (a). With this expression, SOC is retained above a certain level in order to ensure the battery physical constraints (SOC > SOCmin) proactively with the adaptation of the equivalence factor. Moreover, this expression constrains battery charging during charge sustain operation until a rational level (e.g., SOCmax = 18%).

Attention is required in selecting the parameters for each expression to achieve feasible solutions. For example, in our effort for parameters tuning, we spotted that some parameter combinations led to extremely low SOC levels or even that the vehicle could not follow the speed profile. So, after a trial-and-error basis in order to achieve feasible solutions, the setup of the algorithm parameters is presented in Table 3.

$$\mathbf{p(SOC)} = \mathbf{kp} \times (\mathbf{SOC} - \mathbf{SOC\_{ref}}) \tag{11a}$$

$$\text{kp(SOC)} = \text{kp} \times (\text{SOC} - \text{SOC}\_{\text{ref}}(\text{D})) \tag{11b}$$

$$\text{p(SOC)} = \text{kp} \times \left(1 - \left(\frac{(\text{SOC} - \text{SOC}\_{\text{target}})}{0.5(\text{SOC}\_{\text{max}} - \text{SOC}\_{\text{min}})}\right)^{\text{a}}\right) \tag{11c}$$

$$\text{SOC}\_{\text{ref}} = \text{SOC}\_{\text{target}} \tag{12a}$$

$$\text{SOC}\_{\text{ref}}(\text{D}) = \text{SOC}\_{\text{i}} - \left[ \left( \text{SOC}\_{\text{i}} - \text{SOC}\_{\text{target}} \right) \times \left( \frac{\text{D}}{\text{D}\_{\text{final}}} \right) \right] \tag{12b}$$

**Table 3.** Parameters selected for the three ECMS versions.


The ECMS algorithm flowchart is illustrated in Figure 3. ECMS requires power demand, vehicle velocity, gear number and current SOC as input data. The first step is to select the numerical values for the physical limitations according Equations (4)–(6). After that, the algorithm calculates the equivalent fuel rate . mtot for the different candidate operating points by adjusting the equivalent consumption of the electrical motor according to the SOC, as described in Equations (9) and (11a)–(12b). Finally, the algorithm selects the case with the minimum fuel mass, which is then translated into specific torque outputs of the ICE and the electric motor.

**Figure 3.** ECMS Procedure Overview.

#### *2.3. Corrections for the Assessment of Different Algorithms*

The assessment of different EMS algorithms needs to be carried out on a fair basis. Each EMS algorithm leads to different decisions for engine and motor engagement (Equation (2)) that may slightly affect the speed of the vehicle due to power availability and gear change interference. The different speed profiles will, in turn, result in a slightly different demanded power profile for each algorithm. These differences could have been avoided by backward modeling, because this approach guarantees that the vehicle exactly follows the target speed. However, in our approach, the forward modeling is chosen because forward simulators are based on physical causality. With these simulators, online control strategies can be developed [8]. Moreover, the SOC difference between trip start and end may differ between various algorithms. Nevertheless, when assessing the impacts of different algorithms on FC, one needs to make sure that the distance, demanded energy, and SOC differences are identical in the various simulations.

A method to adjust for such potential differences is, therefore, introduced. Assuming a total energy consumption over a theoretical accurate driving profile in each simulation, variations of this profile will lead to energy differences, not because of EMS performance but because of distance and SOC variations in each simulation. A corrected fuel-equivalent energy consumption (CE) can, therefore, be estimated from the simulated one (SE) by

correcting for deviations in the SOC (ΔE SOC), propulsion energy (ΔEPROP) and contribution of regenerative braking (ΔEREG) in Equation (13). ΔE SOC correction is needed to make sure that all simulations result in identical final SOC. The ΔEPROP term corrects for slight differences in the driving profile (speed, acceleration and distance) of the simulations of a given driving sequence. Finally, ΔEREG adjusts the energy consumption when the simulated regenerative braking energy benefits are different from the ones calculated in the theoretically accurate driving profile. In that way, the energy differences due to driving profile variations at the braking phases can be corrected. The corrected energy correction of Equation (13) should be implemented in all relevant works where different optimization algorithms are being compared.

$$\text{CE} = \text{SE} + \Delta \text{E}\_{\text{SOC}} + \Delta \text{E}\_{\text{PROP}} - \Delta \text{E}\_{\text{REG}} \tag{13}$$

Equation (14) describes ΔESOC as the fuel energy delivery that covers the difference between the simulated and reference depleted energies from the battery (Ebat). In the denominator of Equation (14), the average product of the individual components' efficiencies has been considered for the time moments that the battery is charged from the ICE, with η<sup>e</sup> being the ICE efficiency. The calculation for the reference value of the depleted battery energy is presented in Equation (15). The value is calculated as a difference from a final SOC level, which in our case has been selected to be 14%, with Cbat and Vbat being the battery capacity and average battery voltage, respectively.

$$\Delta\text{E\_{SOC}} = \frac{\text{E\_{bat} - \text{E\_{bat,SOCf}}}}{\left(\overline{\left(\eta\_{\text{batt}} \times \eta\_{\text{em}} \times \eta\_{\text{e}}\right)} \left| \mathcal{P}\_{\text{em}} < 0 \wedge \mathcal{P}\_{\text{e}} > 0\right)}\tag{14}$$

$$\rm{E\_{bat,SOC\_f}} = (\rm{SOC\_i} - \rm{SOC\_f}) \times \rm{C\_{bat}} \times \overline{V}\_{\rm{bat}} \tag{15}$$

The ΔEPROP—Equation (16)—is the fuel energy that should be supplied to equalize the simulated energy demand at gearbox (EGB) with the one calculated from the theoretical speed profile (EGB,th). The ICE efficiency should be the average one during positive power demand at gearbox inlet (PGB). EGB,th—Equation (17)—is the energy demand for vehicle motion for positive tractive force at the wheels (Fth), with vth and ηtr being the force, theoretical vehicle speed and average transmission efficiency from gearbox inlet to vehicle wheels, respectively. Ftr,th consists of a polynomial function of vehicle speed (which corresponds to road loads) and the term for vehicle acceleration (Equation (18)—F2, F1, F0 are coast down test coefficients and Mv is the vehicle mass).

$$
\Delta \text{EPRP} = \left(\frac{\text{E}\_{\text{GB,th}} - \text{E}\_{\text{GB}}}{\overline{\eta\_{\text{e}}}} \left| \mathbf{P}\_{\text{GB}} > 0 \right.\tag{16}
$$

$$\mathbf{E}\_{\rm CB,th} = \left(\frac{\sum \mathbf{F}\_{\rm tr,th} \times \mathbf{v}\_{\rm th}}{\overline{\eta}\_{\rm tr}} \mid \mathbf{F}\_{\rm tr,th} \ge \mathbf{0}\right) \tag{17}$$

$$\mathbf{F\_{tr,th}} = \mathbf{F\_2} \times \mathbf{v\_{th}^2} + \mathbf{F\_1} \times \mathbf{v\_{th}} + \mathbf{F\_0} + \mathbf{M\_V} \times \frac{\mathbf{dv\_{th}}}{\mathbf{dt}} \tag{18}$$

ΔEREG—Equation (19)—expresses the potential fuel energy that can be saved if the simulated speed profile was identical to the theoretical one (EBat, th) during decelerations. The simulated battery energy influx is calculated for the time instances that the power demand at gearbox inlet is negative. To convert the battery energy influx difference to fuel energy, the same average product of efficiencies as the one in Equation (14) has been used. The EB−Bat,th—Equation (20)—is the potential energy of braking that can be recuperated. In this calculation, the negative energy influx from the wheels (EB,th) is calculated based on Equation (21). In the current study, the share of the total braking energy that can be recuperated (bREG) is assumed to be 40% of the total braking energy, while the rest is consumed at the mechanical brakes. The losses from the gearbox inlet up to the battery

have been taken into account by the average product of the EM and battery efficiencies during the time events that PGB is negative.

$$
\Delta \mathcal{E}\_{\rm REC} = \left( \mathcal{E}\_{\rm B-But,th} - \mathcal{E}\_{\rm Bat} \left| \mathcal{P}\_{\rm GB} < 0 \right\rangle \times \left( \overline{\left( \eta\_{\rm bath} \times \eta\_{\rm em} \times \eta\_{\rm e} \right)} \left| \mathcal{P}\_{\rm em} < 0 \wedge P\_{\rm e} > 0 \right\rangle} \right. \tag{19}
$$

$$\mathbf{E\_{B-But,th}} = \left(\mathbf{E\_{B,th}} \times \mathbf{b\_{REG}} \times \overline{\eta}\_{\rm tr} \times \overline{\left(\eta\_{\rm batt} \times \eta\_{\rm em}\right)} \, \middle| \, \mathbf{P\_{GB}} < 0\right) \tag{20}$$

$$\mathbf{E}\_{\rm B,th} = \left(\sum (\mathbf{F}\_{\rm tr,th} \times \mathbf{v}\_{\rm th}) \, | \, \mathbf{F}\_{\rm tr,th} < 0\right) \tag{21}$$

#### *2.4. Drive Cycles Used in the Simulations*

The performance of the Rule Based (RB) and the different types of ECMS Based (EB) algorithms are compared in different driving cycles, according to Table 4. Firstly, they are compared in WLTC driving cycle. After that, various driving cycles have been used in order to examine different conditions and identify the best algorithm in each case. For WLTC, ERMES and 10-15 Mode [32] cycles, the individual segments have been examined as well. Table 4 displays a summary of the characteristics of the different driving cycles that have been used. All simulations have been conducted with the same initial SOC of 20%, which allows a comparison of the algorithms during the most challenging charge sustain mode.

**Cycle Repetitions Average Speed [km/h] RPA [m/s2] Total Trip Length [km]** 10 Mode 6 16.0 0.198 12.0 WLTC Low 6 18.4 0.217 18.4 10-15 Mode 5 22.9 0.172 20.8 JC08 [33] 2 26.9 0.184 20.6 ERMES Urban 6 32.1 0.188 29.7 WLTC Medium 6 40.0 0.209 28.2 WLTCx2 2 46.1 0.160 46.1 WLTC High 6 56.0 0.137 42.6 ERMES 2 66.0 0.106 48.0 ERMES Extra Urban 6 69.9 0.120 44.1 WLTC Extra High 6 90.7 0.131 49.0

**Table 4.** Simulated Driving Cycles and characteristics.

#### **3. Results**

#### *3.1. Validation of the Rule Based Algorithm*

Figures 4 and 5 present examples for model accuracy in SOC and fuel consumption, respectively, for four of the seven conducted tests (Table 2). The initial SOC was the actual one at the beginning of each test. The simulated SOC profile satisfactorily follows the measured one during the course of each cycle. This can be reflected by the high correlation coefficient values for SOC (RSOC), which are higher than 69%. There is only one exception, for WLTC CS HOT1, where our back-engineered RB results in higher battery depletion in the 1000–1200 s range compared to the measured one. Apart from this, the simulated SOC follows the measured trend with a rather constant offset.

ERMES Motor 6 96.0 0.086 75.0

**Figure 4.** SOC measurement vs. simulation—examples.

**Figure 5.** FC measurement vs. simulation—examples.

Table 5 compares the total simulated and measured FC and final SOC (SOCf) values for the seven tests conducted. The absolute error in the simulation of total FC is lower than 7% for five of the seven tests. In the WLTC CD and WLTC CDCS tests, the FC estimation error

is higher, but over very low FC values that overall lead to satisfactory deviations (WLTC CD: 0.5 g/km; WLTC CDCS: 4.2 g/km). The correlation coefficient between instantaneous measured and simulated fuel consumption (Rfuel) is always higher than 70% except in WLTC CD (57.7%). In WLTC CD, the initial SOC is rather high (71.4%) so the engine engagement is limited leading to few points where the instantaneous fuel consumption is non-zero, and this leads to a drop of the regression coefficient. With regard to battery SOC, the simulated final SOC values are quite close to the measured ones for WLTC. For the ERMES cases, the simulated SOC varies more than the measured one in WLTC. Additionally, for the ERMES CS the RSOC value is negative. This means that the simulated battery charging events do not follow the measured ones for that test. It is worth mentioning that the parameters tuning for engine start/shut-down events (Figure 1) were based on the WLTC tests experimental data. As the ERMES cycle has different characteristics, we expected that the model could not have the same accuracy as in the WLTC tests. For all the other test cases, the RSOC is higher than 69% which implies satisfactory model performance.



### *3.2. Comparative Assessment of RB and EB Algorithms in WLTC*

Table 6 summarizes the simulation results for the three ECMS types (EBA, EBB, EBC) and the RB algorithm over WLTC. The FC values show correspondence to the simulation output (SFC) and the corrected ones (CFC), according to Equations (13)–(21). The table values clearly show that if no correction was introduced then one would come up with totally wrong conclusions regarding the relative performance of the different algorithms. For example, EBC leads to the highest FC difference over the RB (−3.4%), which is actually lower than the magnitude of the correction (−5.6%). Had the correction not been applied, the EBC would have actually resulted in +4.9% higher FC than the RB, i.e., this would have led to an entirely opposite conclusion than what is actually reached using the corrected values. This can be explained because the ΔESOC correction turns out negative for the EB and positive for the RB cases. The negative value for ΔESOC means that the final SOC is higher compared to the reference one in EB, and vice versa for RB. Additionally, Table 6 shows that ΔEPRO has a higher impact on EB cases. This indicates that the simulated speed profile in the RB simulation better follows the theoretical one than in EB simulations.

The net energy flows between the different powertrain components during WLTC are presented in Figure 6. The shown energy flows stand for the positive propulsion instances. The electrical energy from the battery (2.89 MJ) and the energy demand at the gearbox inlet (13.2 MJ) have been adjusted to the exact same values in order to ensure that the comparison of the different algorithms is on a fair basis. Furthermore, the shown fuel energy consumption is the corrected one, which has been calculated from Equation (13).


**Table 6.** Consumption and mean efficiency values of the main powertrain components for WLTC.

**Figure 6.** Net energy flows (kJ) for positive propulsion instances in WLTC for ECMS-Based Type A (EBA) algorithm, ECMS-Based Type B (EBB) algorithm, ECMS-Based Type C (EBC) algorithm and Rule-Based (RB) algorithm. ICE: Internal Combustion Engine, EM: Electric Motor, B: Battery and GB: Gearbox, FT: Fuel tank. Negative numbers correspond to energy loss.

The electric motor may function either as a propulsion device or as a generator for battery charging, depending on power flux direction. Assuming positive power flux from the battery to the gearbox, the energy flow through the EM ranges from 2121 kJ with the EBB to 1622 kJ with the RB algorithm. The difference in net energy flows is explained by the decisions of the different algorithms.

The lowest net EM energy flow is balanced by the maximum total energy outflux from the FT (30.6 MJ) with RB compared to EB. The simulation results in Table 6 show that the average ICE efficiency is quite close for all cases. Additionally, the demanded mechanical output energy from the ICE is higher in the RB case (11.59 MJ) compared to the EB (11.1–11.2 MJ). In RB, the electric motor contributes more for propulsion compared to the EB cases and results in higher battery discharge. This requires higher power demand from the ICE during hybrid mode to recharge the battery. That higher power demand for battery charging could have been saved had the EMS controller decided to directly command ICE propulsion with the assistance of the electric motor.

The above analysis showed that the description of the energy flows is very useful in order to understand the FC differences between RB and ECMS algorithms. For analysis convenience, in this paper a quantifiable metric for energy flow comparisons has been used, namely the average propulsion efficiency (ηPROP). The ηPROP is defined as the ratio of the demanded energy at the gearbox (EGB) over the sum of the two net energy flows from the energy sources during the positive propulsion phase (Equation (22)). The provided fuel energy is the corrected one, but only during the propulsion phases. Therefore, the

correction terms of ΔE SOC and ΔEPROP have been added to the simulated fuel energy. Additionally, the reference battery depleted energy has been considered in the denominator. Table 6 shows that the RB algorithm results in the lowest propulsion efficiency and that leads to the highest fuel consumption.

$$\overline{\eta}\_{\text{PROP}} = \frac{\text{E}\_{\text{GB}}}{\text{E}\_{\text{bat,SOC}\_f=14\%} + \text{SE} + \Delta\text{E}\_{\text{SOC}} + \Delta\text{E}\_{\text{PROP}}} \tag{22}$$

As a benchmark, the results of Table 6 that have been obtained with localized optimization are compared to the global optimum for the known profile of WLTC. To do so, we have properly parameterized the hybrid electric vehicle model of Sundstrom et al. [34], which uses a DP a solution as an EMS algorithm. The parameterization has been achieved using exactly the same values for the individual components (ICE, EM, battery, Gearbox, axles, vehicle resistances and weight) with the ones used in the AVL Cruise simulated model. The DP model delivered 640.9 g as global optimum CFC, which is 7.3% lower than the one from the best-performing EBC (640.9 g vs. 691.4 g).

#### *3.3. Comparison of RB and EB Algorithms for Other Cycles*

Figure 7 illustrates the relative CFC difference between EB and RB algorithms for the different cycles. In most cases, EB achieves lower fuel consumption compared to RB. The highest FC reduction is observed for WLTC Low, which is 22.1–24.2% depending on the EB type, while the highest EB increase over RB is 2.7% over the ERMES Extra Urban. EB outperforms RB in all cycles of low speed. Evidently, the restrictions of RB on engine switch on criteria (Figure 1) and power output (Figure 2) have a cost on FC reached, while the adaptability of ECMS (Figure 3) allows for a much better result to be obtained.

**Figure 7.** Fuel consumption change in ECMS algorithm types vs. Rule Based algorithm case for different cycles. EBA: ECMS Based Type A, EBB: ECMS Based Type B, EBC: ECMS Based Type C.

Table 7 better explains how EBC and RB perform for a low speed (WLTC Low) and a high speed (ERMES Extra Urban) cycle. In WLTC Low, EBC results in a higher propulsion efficiency by eight percentage units compared to the RB, and this leads to 22% FC reduction (RB: 582 g and EBC: 453 g). The improvement in propulsion is mainly linked to the mean ICE efficiency in this case. As Figure 2 shows, RB selects to operate the engine at lower load for vehicle speeds 0–30 km/h, which results in low engine efficiency. In reality, this selection may have to do with the need to retain low noise, vibration and harshness (NVH) levels in the vehicle cabin during low-speed driving. On the other hand, ECMS algorithm decisions for the engine load are dependent on the optimization of a cost function—Equation (7). As a result of its optimization-based functionality, ECMS leads to higher ICE Mean efficiency by 3.7 percentage units compared to RB in WLTC Low.


**Table 7.** Consumption and mean efficiency values of the main powertrain components for WLTC Low and ERMES Extra Urban in the case of Rule Based (RB) and ECMS Based—C (EBC) algorithms.

The FC with ECMS appears higher than RB only in the ERMES cycle, and in particular in ERMES Extra Urban. In order to explain this, we can look the propulsion efficiency values at Table 7. The RB algorithm has a higher propulsion efficiency by 0.6 percentage units compared to EBC. The ERMES has a much higher power requirement and stronger accelerations than WLTC, which was used to tune the parameters of the ECMS algorithm (Equations (11a)–(12b) and Table 3). A better tuning for high power cycles could lead to lower FC compared to RB.

#### **4. Discussion and Conclusions**

This study makes an assessment and performance analysis for two types of energy management algorithms, including a back-engineered stock RB algorithm and an ECMS one.

A novel methodology to assess the different algorithms on a fair basis has been developed, introducing corrections for distance travelled, final SOC level and energy propulsion differences between alternative simulations. Our analysis showed that the impact of such corrections on fuel consumption can exceed 5%. Hou et al. [35] also considered SOC corrections when comparing different EMS algorithms, and found the impact of these corrections to be no more than 1% in fuel consumption. However, they did not correct for propulsion energy demand and regenerative braking. Our analysis shows that all three corrections can have a measurable result on fuel consumption values when different EMS algorithms are compared, and these should not be neglected in relevant studies.

The energy flow analysis also showed that the EMS affects both the efficiency of each powertrain component and the net energy flow between components. Therefore, it is the combination of these two magnitudes that affects overall total fuel consumption. As a result, the efficiencies of the individual components can be quite close for the two different EMS algorithms, but the FC can differ for the same energy demand.

Another important finding is that different driving conditions affect the magnitude of FC reduction that can be achieved with an alternative algorithm. Regarding the WLTC cycle, it has been found that the potential FC reduction with the use of ECMS algorithm is 3.4% over RB. For benchmarking purposes, a DP global optimization algorithm has been used for comparison. This requires a priori knowledge of the driving profile to find the optimum

solution. The DP led to a 7.3% lower FC compared to the one from ECMS. Sun et al. [13] demonstrated that, on average, the FC achieved with DP can be 6.7% lower compared to an ECMS algorithm. ECMS algorithms, therefore, seem to offer good adaptability to driving conditions and sufficiently good final fuel consumption, which is not very distant from the global optimum.

With the use of ECMS, we found that fuel consumption in driving cycles with low average speeds (lower than 20 km/h) can be improved by up to 24.2% compared to a stock RB algorithm. The magnitude of improvement achieved was actually similar to the ones reported by Geng et al. [36], who proposed a combination of DP and ECMS to reduce FC by up to 19.9% in NEDC compared to RB. In another study, Hao et al. [37] used an adaptive ECMS to decrease the FC by 8–15% compared to an RB algorithm on a mild parallel hybrid vehicle. In our case, this large difference was partly because the RB forced the engine to operate at low loads when switched on in low average speeds. Although this may have been mandated in order to retain low NVH in the cabin at low speeds or due to pollutants emission control restriction, one needs to admit that an improvement margin of more than 20% is a significant incentive to invest more in powertrain efficiency control, especially in urban conditions.

This paper examines the energy optimization of the hybrid powertrain for propulsion. However, an EMS algorithm may consider additional parameters, such as consumption of auxiliaries for cabin comfort and the thermal management of the emission control devices. These parameters are not addressed in the current work, as they add complexity and extra investigation is required to achieve balance between real-time implementation and optimality. However, they can be addressed in a future work by extending the cost function expression to cover these terms and by extending the list of conditions that have been considered to the implemented the algorithms. Such conditions can also take into account NVH requirements that can be added as cost penalizing terms in relevant optimization.

In general, ECMS algorithms are simple to enforce in an ECU as they do not require a full set of guidance for all foreseeable conditions and could always be used as back-up, fail-safe algorithms. For example, in driving situations that have not been thoroughly considered when setting up an RB algorithm and which cause higher FC than what would be technically possible, an ECMS could step in instead. Moreover, an ECMS could be identical for different vehicle model variants, without needing to set exact limits for each version of the vehicle when, e.g., engine or motor sizes change while powertrain architecture stays the same. An ECMS algorithm can, therefore, have a more universal usage due to its adaptability features.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/en16031497/s1, Figure S1: Examined Cycles.

**Author Contributions:** Conceptualization, N.A., S.D. and L.N.; Methodology, N.A., S.D. and L.N.; Software, N.A. and S.D.; Writing—Original draft preparation, N.A.; Visualization, N.A.; Writing— Review and Editing, L.N.; Supervision, L.N. and Z.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research work was supported by the Hellenic Foundation for Research and Innovation (HFRI) under the HFRI Ph.D. Fellowship grant (Fellowship Numbers: 607 and 6653).

**Data Availability Statement:** The data presented in this study are available in this article and the Supplementary Materials.

**Acknowledgments:** The authors would like to acknowledge support of this work by the personnel of the Laboratory of Applied Thermodynamics for performing the experimental campaign of the vehicle. The research work was supported by the Hellenic Foundation for Research and Innovation (HFRI) under the HFRI Ph.D. Fellowship grant (Fellowship Numbers: 607 and 6653).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**


#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Ana Pavli´cevi´c \* and Saša Mujovi´c**

Faculty of Electrical Engineering, University of Montenegro, Džordža Vašingtona bb, 81000 Podgorica, Montenegro

**\*** Correspondence: anaz@ucg.ac.me; Tel.: +382-67559535

**Abstract:** Climate change at the global level has accelerated the energy transition around the world. With the aim of reducing CO2 emissions, the paradigm of using electric vehicles (EVs) has been globally accepted. The impact of EVs and their integration into the energy system is vital for accepting the increasing number of EVs. Considering the way the modern energy system functions, the role of EVs in the system may vary. A methodology for analyzing the impact of reactive power from public electric vehicle charging stations (EVCSs) on two main indicators of the distribution system is proposed as follows: globally, referring to active power losses, and locally, referring to transformer aging. This paper indicates that there is an optimal value of reactive power coming from EV chargers at EVCSs by which active energy losses and transformer aging are reduced. The proposed methodology is based on relevant models for calculating power flows and transformer aging and appropriately takes into consideration the stochastic nature of EV charging demand.

**Keywords:** electric vehicles; electric vehicle charging stations; active power losses; power distribution transformers; thermal aging; reactive power

#### **1. Introduction**

In July 2009, the leaders of the European Union and the G8 announced the aim to reduce greenhouse gas emissions to at least 80% below 1990s' levels by 2050. The automotive sector has a very important role to play in achieving this ambitious goal [1].

One of the best ways to promote the use of electric vehicles (EVs) is to increase the number of available public electric vehicle charging stations (EVCSs). Connecting a large EV fleet in a small area faces a number of challenges such as: overheating distribution transformers, increasing power losses, voltage instability and harmonic distortion [2,3]. Therefore, there is a need for modeling and mitigating the impacts of EVCSs on the distribution grid. Connecting EVCSs to the distribution system has a significant impact on the exploitation of the existing energy infrastructure. This especially refers to increasing the total load, peak load and changing the load profile of the distribution transformer where EVCSs are connected. Taking into account a large number of distribution transformers as well as the fact that in most cases they do not operate in parallel in low voltage networks, the prevention of failures or possible losses of life (LOL) is of great importance. Since charging EVs at EVCSs increases load on distribution transformers, it could cause overloading and therefore overheating of distribution transformers. Apart from overheating, there are other sources of transformer aging such as: fault currents caused by short-circuits or inrush currents, overvoltage caused by lightning and switching impulses, contamination of the oil that could shorten their lifetime and cause premature failure [4,5]. Bearing in mind that the subject of this paper is the influence of EVCS on the power system and on transformer aging under normal operating conditions, an overview of the literature related to the aging of transformers under such conditions is given.

Researches have studied the impact of distributed generation (DG) and energy storage (ES) on transformer aging [6,7]. It has been concluded that with a high penetration of

**Citation:** Pavli´cevi´c, A.; Mujovi´c, S. Impact of Reactive Power from Public Electric Vehicle Stations on Transformer Aging and Active Energy Losses. *Energies* **2022**, *15*, 7085. https://doi.org/10.3390/en15197085

Academic Editor: Wojciech Cieslik

Received: 16 August 2022 Accepted: 21 September 2022 Published: 27 September 2022

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

EVs, and no support from photovoltaic (PV) and ES, the transformers will dramatically age. Furthermore, researchers have dealt with the influence of the charging method on transformers' aging [8–13]. The main idea is to find an optimal charging scheme to minimize the impact of EVs' charging demand on distribution transformers. The [8] proposed model considers time-of-use rates in order to minimize energy consumption costs and avoid transformer overloading and LOL based on load and meteorological data. The proposed smart charging scheme together with PV with ES can prevent transformer overloading and LOL. PV generation can reduce the energy purchased from the grid, while ES can assist during peak hours. Three smart ways of charging (central, decentralized and hierarchical) as well as dump charging have been analysed [9]. Based on the obtained results, it has been concluded that hierarchical charging strategies are the most desirable in terms of LOL, while the centralized charging strategy (valley filling) has shown the greatest LOL of transformers. In one paper [10], a simple method is used to coordinate the charging process of PEVs to avoid the overloading of transformers by using the fact that idle time represents a significant percentage of EV parking time. In the paper, [11] reducing the impact on transformer aging by implementing the capability of delivering a variable charging rate has been proposed. Two charging algorithms have been used as strategic ways for reducing the overload of distribution transformers. The results suggest that for both uncontrolled and controlled charging, the load on the transformer will be reduced. Since the current commercial, residential charging technology is categorized, infrastructure changes are necessary in order to support this idea. Another paper [12] proposes a temperature-based smart charging strategy that reduces transformer aging without substantially reducing the frequency by which EVs obtain a full charge. These reductions are substantially larger at hot climate locations compared to cool climate ones. The results indicate that simple smart charging schemes, such as delaying charging until after midnight can actually increase, rather than decrease, transformer aging. The authors of [13] analyze the impacts of a price incentive-based demand response on neighburhood distribution transformer aging. The results indicate that the integration of EVs in residential premises may indeed cause accelerated aging of distribution transformers, while the need to investigate the effectiveness of dynamic pricing mechanisms is evident. In [14], the number of EVs has been optimized by calculating additional steady-state hottest-spot temperature rises by ensuring that distribution reliability requirements are met.

In order to reduce transformer aging, the reactive power potential from on-board EV chargers (OBCs) at EVCSs is proposed in this paper. OBCs in modern EVs are bidirectional, which allows them to work in all four quadrants [15]. Using reactive power support, as a part of the vehicle-to-grid concept, is already familiar and has been investigated as, e.g., part of on-line control for improving the voltage stability of a microgrid with PV, wind, battery storage, residential and commercial building feeders [3], the load control of active and reactive power [16], the compensation of reactive load of wind turbines near the station [17] and the revenue potential of providing reactive power service [18]. In this paper, the respective value of reactive power is obtained by small oversizing charging converters, which does not require an upgrade of the charging infrastructure. In relation to [15], where the influence of the oversizing of converters on voltage deviations, peak load and grid losses in the network has been analyzed, this paper's analysis is extended to include transformer aging. Besides the impact on transformer aging, also the influence of reactive power potential of EVCSs is analyzed with the aim of reducing active energy losses without additional devices and facilities. Locations of EVCSs are chosen from the aspect of increasing the efficiency of the network, which refers to reducing total active energy losses.

The reduction of losses in the power system in the presence of EVCSs is achieved in different ways. A large number of papers can be found in the literature that treat the problem of reducing energy losses in the presence of electric vehicles, e.g., optimal operation of vehicles with distributed generators [19–23], optimal reconfiguration of the distribution network [24–28], optimal operation of the battery storage, [29,30], optimal charging strategies [31–34] and optimal location and size of reactive power compensators [35,36].

The contribution of this paper is related to the impact analysis of the influence of reactive power on active power losses and transformer aging. In other words, the paper proposes a methodology that allows the determination of the optimal injection of reactive power of the EVCSs with the aim of reducing transformer aging and reducing active power losses. The proposed method is based on the well-known model for calculating transformers' aging, which is described in detail in [37]. It is applicable to an arbitrary network with an arbitrary number of EVCSs. Considering the stochastic nature of EV charging demand, the presented paper clearly indicates the benefits and limitations of reactive power from OBCs at EVCSs to energy losses reduction and transformers aging.

Numerical results, obtained in MATLAB and on the IEEE 33 bus radial distribution network, indicate that it is possible to achieve the optimal value of injecting reactive power from EV vehicles in order to enable reductions of transformers' aging and a reduction in energy losses in relation to the case where there is no reactive power injection or when reactive power injection is the maximum possible. This analysis may serve focus on local distributions in order to prevent shortening transformers' working life. Furthermore, the presented results may be useful for the analysis of global distribution system parameters such as energy losses and node voltages.

The rest of the paper is organized as follows. Section 2 is dedicated to the system components modeling. In Section 3, a description of methodology is given, while in the final Section 4, the most important conclusions and directions for further research have been stated.

#### **2. Modeling System Components**

#### *2.1. Load Modeling*

There is certain number of load models in the literature that can be used for various analyses. In addition, in recent times, new nodal classifications have been made and their characteristics have been reviewed [38]. Considering that a steady state analysis is performed in this paper, basic information about node loads have been used in this paper. Namely, loads in the system are modeled as PQ nodes. On the basis of the type of loads connected to the nodes busses, three types have been adopted: residential, commercial and industrial busses, as seen in Table 1 [19]. There are 18 residential busses, 5 commercial buses and 9 industrial buses. Load changing over time has been taken into account by introducing a load scaling factor. Load scaling factor changes over time as a percentage of nominal load, which varies depending on a load type, Figure 1, [39].

**Figure 1.** Typical buses load curve, [39].

**Table 1.** Grouping of buses data, [19].


#### *2.2. Public EVCS Modeling*

According to the way EVs are connected to the electric grid, there are two types of charging strategy: conductive EV charging and inductive EV charging. In conductive charging, there is a physical connection between EV and the electrical grid. In inductive charging, also known as wireless power transfer [40], there is no galvanic connection between the vehicle and the electrical grid. Since conductive charging is the most common method of charging, this type of charging has been analyzed. There are two types of conductive charging in EVCSs: AC charging and DC charging with different power levels: Level 1, Level 2 and Level 3 charging. In this paper, the AC Level 2 type of charging station has been analyzed because of its simple infrastructure, low price and availability for the most commercial types of charging in public places.

There is a certain number of models for grid support by EVCSs [41]. Bearing in mind that vehicles are mostly stationary during the day, their role may be different. The literature recognizes that EVCSs could most effectively serve in the following ways: to ensure a fast frequency response [42,43] for shifting electricity [44,45], improve power quality [46] and as energy storage [44].

Significant work has been done towards minimizing grid losses due to PHEV charging through efficient charging algorithms with different objective functions [47–49]. All mentioned papers, in one way or another, with different stimulation techniques, have an influence on drivers' behavior. Nevertheless, this paper adopts an approach that does not impair the comfort of EV users, which means that two cases have been observed, which are as follows: active power is consumed during the charging of EVs; while in the second case, while charging with active power, vehicles inject reactive power.

#### 2.2.1. Operation Region of EV Chargers

In order to obtain reactive power and not to jeopardize EVs drivers' habits, the concept of small oversizing of 5.3% of on-board chargers ratings has been used, [15]. With this, a significant amount of reactive power which is 32.9% of active power is obtained. What is important is that this does not require an upgrade of the charging infrastructure, and charging active power is the same in both cases. The operating region of on-board charger charging is region IV, Figure 2. It can be seen that with small oversizing of an EV charger, in Figure 2 shown as a larger circle, with the same active charging power, a representative value of reactive power is obtained.

The charging operation in the second and third quadrants of the P-Q diagram was not considered, since the management of the active power of the chargers affects the convenience of the consumers themselves as well as battery life. Furthermore, the operation of chargers in the first quadrant, wherein EV can also consume reactive power, is not suitable for the analyzed network and loads.

During the charging of EVs, we have analyzed the influence of reactive power from equal to zero to the maximum available. The national network operator in our country does not yet support the injection of active power or the injection and absorption of reactive power by the EVCSs. The operation of chargers in all four or two quadrants is discussed in the literature at the level of conceptual solutions for which there are numerical or real experiments. The results obtained both in this work and in the literature represent guidelines for changing the network requirements of new electricity consumers/producers. This paper demonstrates the benefits of reactive power from EV charger compared to the current regulation, which stipulates that they only consume active power. Finally, it is

important to note that the network codes for the LV network, after looking at the possible positive impact of the operation of network inverters of photovoltanic power plants, have been changed in Germany. For example, in Germany, LV grid-connected PV installations rated above 3.68 kVA have to follow a specified PF droop curve as a function of their instantaneous power output [15].

**Figure 2.** Operation region of EVs chargers.

Furthermore, it is very important to regulate the interaction between EVs the and distribution system. The design of the battery charger will be crucial in this effort to effectively control the power flows. One of the important requirements of an EV charger is the amount of current distortion that it draws from the grid [50]. Parameters such as total harmonic distortion (THD) and total demand distortion (TDD) can be used to evaluate the harmonic content of the charger. Additionally, the chargers should meet the individual harmonic limits as well [50–52]. The tasks of OBC system controller is to follow the charging power and reactive power commands controlled by the grid operator. In [50], the proposed design of OBC has shown analytically and experimentally that chargers are able to symmetrically operate at all four quadrants of the power plane. Furthermore, grid current requirements such as that the ac-dc converter input current THD should be limited to 5% have been satisfied.

Figure 3 shows the proposed application of on-board EV chargers used in this paper. The chargers form EVCS injecting reactive power at the point of common coupling (PCC) and decreases transformer overloading. In the case of the AC charging station, electric vehicle supply equipment (EVSE) serves for monitoring, management and communication with a vehicle during charging, while energy conversion from AC to DC power is suitable for charging battery performed via OBC.

In this paper, it is assumed that the power factor management of EVCSs is performed based on the previous day's consumption forecast of the analyzed network and EVCS demand. After processing the necessary data, the distribution network control centre sends a signal for the value of the power factor of the station. The power factor is constant during the 24 h period.

**Figure 3.** Proposed reactive power support diagram using OBC from EVCS.

2.2.2. Mobility Stochastic Behavior Model

There is a large number of papers dealing with the optimal structure of EVCSs, which refers to number of charging spots, types of charging and technology of chargers themselves [53–55]. In this paper, the charging demand of EVCSs is obtained using the number of EVs that arrive at the EVCS during 24 h as an input. Furthermore, the random nature of EV arrival time, electrical efficiency, battery capacity and daily miles driven are respected in order to obtain EVCS demand. In this paper, the stochastic nature of EV load is modeled using probability density functions (PDF) for two parameters: plug-in time and miles traveled before last charge.

• Behavior of start charging time

A PDF of vehicles' arrival time is obtained from [56], where the results of a large number of measurements from public charging stations have been analyzed. Based on actual data, the multimodality of distribution has been described using the Beta Mixture Model (BMM), which has proved to be appropriate for analysing EV data, [56].

A corresponding density function from the BMM model is obtained by summing up the beta distributions of different parameters and weight factors. The beta distribution of probability density is defined by the following equation:

$$\text{Beta}(\boldsymbol{\pi} \mid a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \boldsymbol{\pi}^{a-1} \boldsymbol{\pi}^{b-1},\tag{1}$$

where Γ(.) is gamma function, *a* and *b* are parameters which define the shape of distribution and *x* is a parameter from the interval between [0, 1].

The BMM is represented by the following sum:

$$\text{BMM}(\mathbf{x}) = \sum\_{m=1}^{M} \left( w\_{\text{ln}} \text{Beta}(\mathbf{x} \mid a\_{m\prime} b\_{m}) \right) \tag{2}$$

wherein the sum of weights are equal to one:

$$\sum\_{m=1}^{M} w\_n = 1\tag{3}$$

**Figure 4.** PDF for plug in time for: (**a**) Working day; (**b**) Weekend.



### • Behavior of distance travelled and initial state of charge

In order to determine the energy required to charge vehicles, it is necessary to predict the level of the charge of a battery (SOC) of a single EV. An important factor for this is to determine daily miles traveled. The lognormal distribution of the distance travelled function has been shown to be suitable for describing the mileage distribution functions, and is represented by way of [6]:

$$\log(d,\mu,\sigma) = \frac{1}{d \cdot \sqrt{2\pi\sigma^2}} e^{-\frac{(\ln d - \mu)^2}{2\sigma^2}}, \; d > 0 \tag{4}$$

The distribution parameters depend on a driver's habits in the area being analyzed. Statistics on the conventional vehicles are often taken into account due to insufficient data of EVs. This paper uses statistical data that 50% of drivers drive less than 25 miles during the day, and 80% of drivers travel less than 40 miles [57]. The corresponding lognormal distribution of daily traveled miles and coresponding CDF are presented at Figures 5 and 6, with mean value μ = 3.37 and standard deviation σ = 0.5 [58].

**Figure 5.** PDF of daily distance traveled of EV.

**Figure 6.** CDF of daily distance traveled of an EV.

Now, the initial state of the charge can be determined as shown below [6]:

$$\text{SOC}\_{i}(\mathbf{x}) = \left(1 - \frac{E\_{\text{cons}} \cdot d}{\mathbb{C}\_{b}}\right) \cdot 100\tag{5}$$

where *Cb* battery capacity in kWh, *Econs* is energy consumtion in kWh per miles and *d* estimated daily distance traveled in miles.

#### 2.2.3. Charging Power Model

Considering the latest statistics of dominating shares of EV types, [59,60] all EVs charged in this study have 71% probability of being a battery electric vehicle (BEV) and 29% probability of being a plug-in hybrid electric vehicle (PHEVs).

Electric car economy changes all the time, depending on things such as road gradient and car speed. In this paper, electrical efficiency and battery capacity are taken as random numbers between values given in Table 3 depending on the type of EV.

**Table 3.** Parameters of BMM model.


In order to sustain a longevity, it is known that whole battery's energy capacity is usually not fully utilized. It has been found that between 78–95% of the full batteries capacity is used for different EV models [60]. The required value of *SOCreq* for the user is taken to be 95% [6,61,62]:

$$E\_{req} = \frac{\left(SOC\_{req} - SOC\_i\right) \cdot C\_b}{\eta \cdot 100} \tag{6}$$

where *η* is charging efficiency 0.95. The duration of EV charging depends on the required energy and charging power of EV and is equal to:

$$T\_{ch} = \frac{E\_{req}}{P\_{EV}} , P\_{EV} = \begin{cases} 3.4 \text{kW} / \text{PHEV} \\ 7.2 \text{kW} / \text{BVV} \end{cases} \tag{7}$$

#### *2.3. Power Flow Analysis*

Power flow or load flow calculations are very fundamental for all power system analyses. It is a very important tool for power systems planning, design, operations, maintenance, optimization and control. At the planing stage, load flow analysis is used to determine: the location of distributed generators, the location of capacitors, economic scheduling, power quality improvements, network reconfiguration, power systems optimization and other applications. Furthermore, in order to plan future growing load demands, power flow analysis is necessary. At the operation stage, it is run to explore system stability and to improve efficiency. Usage of power flow analysis enables that maintenance plans can proceed, without violating power system security. Load flow studies are performed for the determination of the steady state operating condition of a power system. Input parameters of power flow analysis of distribution network are network topology and network parameters, as well as load and generator models. Power flow calculations are performed by iterative methods. The most common power flow calculation methods are Gauss–Seidel Method, Newton–Raphson Method and Fast Decoupled Method. These methods are not appropriate for distribution systems due to their special characteristics such as low line *X*/*R* ratios, unbalanced load, radial or weakly meshed network structure and large number of nodes, etc. These features make distribution systems power flow computation different to analyze, as compared to the transmission systems. There is a number of methods proposed in the literature for the solution of power flow problem in radial distribution networks [63]. Back–forward sweeping (BFS) iterative method is suitable for calculations of radial distribution system load flow, which is analyzed in this paper [64]. The algorithm is very robust and numerically efficient for convergence. It is applicable for a wide variation of distribution networks. The algorithm begins with an initial solution for all node voltages and performs basic steps until a convergence criterion is satisfied. In the backward step, currents of each branch are calculated. Bus voltages are updated in a forward sweep starting from the first branch and moving towards end branches. After the convergence criterion has been satisfied, node voltages, currents in branches and active and reactive power losses are determined. This iterative method is implemented in Matlab for calculating power losses in radial IEEE 33 bus network. Different load types together with EVCS demands are implemented in this model. In this paper, power flow analysis has been used from the planning aspect, as a part of the optimization process for determinating optimal positions for EVCSs in an analyzed network. Additionally, it is used for calculating power losses for different EVCS demands level.

#### **3. Methodology**

#### *3.1. Description*

A flowchart of the proposed methodology is shown in Figure 7. The methodology includes both global and local solution approaches. The global one refers to the modeling and analysis of the power network and the determination of active power losses in the network, and the local one to the determination of the aging of the selected distribution transformer. The proposed algorithm consists of four steps. In the first step, network data and nodal loads are entered. Additionally, active and reactive nodal power are defined depending on a type of load (industrial, commercial or residential) according to Figure 1. Then, in this step, the number of stations are entered, as well as the data needed to determine the EVCS demand. The stochastic nature of EV charging is taken into account and explained in detail in Section 2. The second step of the proposed methodology includes the position of EVCSs in the power system. More precisely, in this step, a matrix that includes the possible positions of EVCSs is formed. In this step, the criterion that only one station can be placed in each part of the network has been met. The proposed methodology takes into account the fact that it is possible to divide the network into an arbitrary number of parts, and in each part of the network an arbitrary number of stations can be found. The third step represents the calculation of the electrical parameters of the network. In this step, the method for calculating power flows according to the characteristics of the network is implemented. In the fourth step, the locations of EVCSs are selected based on the criterion of minimum network losses, defined by the following equation:

$$J = \min P\_l = \sum\_{i=1}^{B} P\_{li} \tag{8}$$

where *i*, *B*, *Pl* and *Pli*, are the index, total number of branches, total active power and active power losses of branch *i*, respectively.

In this paper, optimization is achieved by searching the obtained results, so that the focus of paper is not the development and improvement of optimization methods. Nevertheless, it is important to note that in the last few years, different methods have been developed and applied to account for optimal power flows in the energy system or for different types of examples that comply with different limitations. Some of them are detailed and tested in papers [65–68]. In the fifth step, the final value of reactive power injection of the EVCSs is selected. Namely, in this step, by choosing the final value of the reactive power injected by EVCSs, the influence of the reactive power on the losses of active energy and on transformer aging through which the station is connected to the grid is affected. In this paper, the GRA method is used for determining the final value of reactive power.

After obtaining the best positions for EVCSs, off-grid managing of the EVCSs has been proposed. Based on the predicted daily diagram of EVCSs and other consumption, the station operator evaluates the optimal daily reactive power injection from EVs, i.e., power factor of on-board chargers of EVs in EVCSs. The proposed methodology assume that power factor of OBC is constant during the 24 h period. Detailed post-optimization analysis of the reactive power injection from EVs at EVCSs to energy losses and transformer aging has been performed.

It is important to emphasize that the performed post-optimization analysis represents an incremental contribution in relation to the paper [15], since it is extended to include transformer aging analyzis. Based on the performed analyses, the correlation of the impact of reactive power on energy losses and the impact of reactive power on transformer aging can be clearly seen.

The aim of this paper is to explore the possibility of the impact of reactive power from oversized OBCs on active power losses and on transformer aging. With regard to experimental confirmation of the obtained results, at this stage it was primarily necessary to computerize a model and verify the simultaneous influence of active energy losses and

transformer aging. On the other hand, the impact of chargers on losses in the distribution system in real-time circumstances is very difficult. Furthermore, one of the bigger obstacles so far is the number of EVs, because in whole of Montenegro there are only 126 EVs [69] and these are vehicles that have chargers that are not pre-dimensional. Relatively small oversizing of the OBCs is the basis of the assumption of this paper.

**Figure 7.** Flowchart of proposed methodology.

#### *3.2. Sample Generation with Monte Carlo*

The Monte Carlo (MC) statistical method was used to obtain the daily EVCS demand. The method requires that parameteres of the physical system, in this case EVCS demand, are described as PDFs.

When these functions are known, the MC simulation continues with a random selection of values from functions. For this purpose, the inverse transformation method was used, which states that any probability distribution can be obtained from a uniform probability distribution, if an inverse cumulative probability distribution can be determined [70]. A flowchart of the generation of a EVCS demand profile is given in Figure 8.

To terminate the MC process criteria number of MC iterations is used. The process was set to repeat until reaching 5000 iterations.

**Figure 8.** Flowchart for EVCS demand.

#### *3.3. Transformer Aging Model*

The insulation system of mineral oil transformers is composed of thermally upgraded oil-impregnated paper. Cellulose deterioration is influenced by hydrolysis, oxidation and pyrolysis which are the consequences of water, oxygen and heat, respectively. Taking into consideration that exposure to moisture and oxygen in transformers are generally reduced, the most significant determining factor to insulation deterioration is the heat. Transformer heating is caused primarily by energy losses. The majority of losses are located in magnetic core (no-load losses) and the windings (load losses). No-load are hysteresis and eddy current losses, while winding losses are primarily due to DC losses and stray load loss due to the eddy currents induced in other structural parts of the transformer. All these losses cause heating in the corresponding parts of the transformer. The most critical part for the transformer isolation system is placed where temperature has maximum value, which is known as hot spot temperature [4].

According to the loading guides, IEEE guide for loading mineral-oil-immersed transformers [37], the hot-spot temperature in a transformer winding consists of three components:

$$
\theta\_H = \theta\_A + \Delta\theta\_{TO} + \Delta\theta\_H \tag{9}
$$

where *θH*, is the winding hottest-spot temperature, ◦C, *θA*, is the average ambient temperature during the load cycle to be studied, ◦C, Δ*θTO* is the top-oil rise over ambient temperature, ◦C, and Δ*θ<sup>H</sup>* is the winding hottest-spot rise over top-oil temperature, ◦C.

As is proposed in the IEEE guide, *θ<sup>A</sup>* is constant and equals to 30 ◦C. The top-oil rise over ambient temperature is given in the following equation:

$$
\Delta\theta\_{\rm TO} = (\Delta\theta\_{\rm TO,\mu} - \Delta\theta\_{\rm TO,i}) \cdot \left(1 - e^{-t/\tau\_{\rm TO}}\right) + \Delta\theta\_{\rm TO,i} \tag{10}
$$

where *τTO* top-oil time is constant, Δ*θTO,i* and Δ*θTO,u* are the initial and ultimate top-oil rises over ambient temperature, respectively, which equals to:

$$
\Delta\theta\_{TO,i} = \Delta\theta\_{TO,R} \cdot \left[\frac{\left(K\_i^2 \cdot R + 1\right)}{R + 1}\right]^n \tag{11}
$$

$$
\Delta\theta\_{TO,\mu} = \Delta\theta\_{TO,R} \cdot \left[\frac{\left(K\_{\mu}^2 \cdot R + 1\right)}{R + 1}\right]^n \tag{12}
$$

In the previous equations, Δ*θTO,R* is the top-oil temperature rise over ambient temperature at rated load, *ki* and *ku* represent the ratio of initial and ultimate load to rated load, per unit, *n* is an empirically derived exponent which depends on a cooling mode and describes effects of change in oil resistance to change in load, while *R* is the ratio of a load loss at a rated load to a no-load loss. The top-oil time constant at a rated kVA is given in the following equation. The time constant equals to:

$$\tau\_{TO} = \tau\_{TO,R} \frac{\left(\frac{\Delta\theta\_{TO,\mu}}{\Delta\theta\_{TO,r}}\right) - \left(\frac{\Delta\theta\_{TO,i}}{\Delta\theta\_{TO,r}}\right)}{\left(\frac{\Delta\theta\_{TO,\mu}}{\Delta\theta\_{TO,r}}\right)^{\frac{1}{n}} - \left(\frac{\Delta\theta\_{TO,i}}{\Delta\theta\_{TO,r}}\right)^{\frac{1}{n}}} \tag{13}$$

Transient winding hottest-spot temperature rise over top-oil temperature is equal to:

$$
\Delta\theta\_H = (\Delta\theta\_{H,\mu} - \Delta\theta\_{H,i}) \cdot \left(1 - e^{-t/\tau\_w}\right) + \Delta\theta\_{H,i} \tag{14}
$$

In Equation (14), Δ*θH,i* and Δ*θH,u* represent the initial and ultimate winding hottestspot rises over top-oil temperature in ◦C, which equals to:

$$
\Delta\theta\_{H,i} = \Delta\theta\_{H,\mathbb{R}} \cdot \mathbb{K}\_i^{2m} \tag{15}
$$

$$
\Delta\theta\_{H,\mu} = \Delta\theta\_{H,R} \cdot K\_{\mu}^{2m} \tag{16}
$$

where Δ*θH,R* is the winding hottest-spot rise over top-oil temperature at rated load. Parameter *m* describes changes in resistance and oil viscosity with changes in load.

Calculating the temperature of the hottest spot allows the determination of important parameters that describe transformers' aging, which are the accelerated transformer aging factor *FAA* and LOL.

The aging acceleration factor for a given load and hottest spot temperature can be obtained as shown in Equation (17):

$$F\_{AA} = e^{\left[\frac{15000}{383} - \frac{15000}{\theta\_H + 273}\right]} \tag{17}$$

If *FAA* > 1, the transformer is experiencing accelerated aging. According to [37], normal aging occurs at the reference hottest-spot temperature of 110 ◦C, where *FAA* = 1. The equivalent aging factor *FEQA* is further obtained as indicated below:

$$F\_{EQA} = \frac{\sum\_{j=1}^{N} F\_{AA,j} \Delta t\_j}{\sum\_{n=j}^{N} \Delta t\_j} \tag{18}$$

where *j* is the time intervals index, *N* is the total number of time intervals, Δ*t* is the time interval, *FAA*,*<sup>j</sup>* is the aging acceleration factor for *j*-th time interval. The percentage of shortening the isolation life of a transformer is calculated as shown in the relation below:

$$\%LOL = \frac{F\_{EQA} \cdot t \cdot 100}{\text{Normal insulation life}} \tag{19}$$

The value of the normal duration of insulation, [37], is taken to be 18,000 h or 20.55 years. Therefore, when the temperature of the hottest point is 110 ◦C for 24 h, the percentage of the daily shortening of the lifespan is equal to:

$$\%LOL = \frac{24 \cdot 100}{180,000} = 0.01333\% \tag{20}$$

#### *3.4. Transformer Loading with EVCS*

Adding EVCS demand at node *n* to the distribution transformer affects its loading, and therefore can cause further transformer aging. The available reactive power of inverters of EVs can be utilized for reactive power compensation in order to improve power loss reduction and prevent transformer aging. The total active load of the transformer at node *k*, *PTLk* , is equal to:

$$P\_{TL\_k} = P\_{L\_k} + P\_{EVCS\_k} \tag{21}$$

where *PLk* and *PEVCSk* re the sum of load without EVs and active power demand of EVCSs at node *k*. In case, during the charging of EVs, reactive power is injected, the total reactive load of transformer *QTLn* is equal to:

$$Q\_{TL\_k} = Q\_{L\_k} - Q\_{EVCS\_k} \tag{22}$$

where *QLk* represents the sum of the load without EVs and *QEVCSi* is reactive power injected from EVCSs obtained from EV charging inverters at node *k*. Injecting reactive power to the transformer at node *k* has an influence on reducing the apparent power of the total load and therefore the loading and temperature of the transformer as well.

#### *3.5. Grey Relational Analysis (GRA)*

The selection of the final solution is done using Grey Relational Analysis (GRA). The GRA is a multiple-attribute decision-making method and it is widely applied in electric power systems [71–74]. GRA is performed in several steps. In the first step, *n* alternatives sequences with *m* criteria are formed. *Yi* represents *i*-th alternative sequence, and value *yij* represents the value of attributes *j* of alternative *i*:

$$Y\_{\mathbf{i}} = \left( y\_{\mathbf{i}1\prime} y\_{\mathbf{i}2\prime} \dots y\_{\mathbf{i}j\prime} \dots y\_{\mathbf{i}n} \right) \tag{23}$$

The second step of the GRA procedure is normalization, which converts the data to values between [0,1]. The smaller the normalization used in this case, the better, since smaller values are prefered in the problem desribed in this paper. The comparability sequence is as follows:

$$X\_{\mathbf{i}} = \left(\mathbf{x}\_{\mathbf{i}1\prime}\mathbf{x}\_{\mathbf{i}2\prime}\dots\mathbf{x}\_{\mathbf{i}\prime}\dots\mathbf{x}\_{\mathbf{i}\prime}\right) \tag{24}$$

is obtain using the equation below:

$$\text{max}\_{\text{ij}} = \frac{\max\left\{ y\_{\text{ij}} \mathbf{i} = 1, 2, \dots, m \right\} - y\_{\text{ij}}}{\max\left\{ y\_{\text{ij}} \mathbf{i} = 1, 2, \dots, m \right\} - \min\left\{ y\_{\text{ij}} \mathbf{i} = 1, 2, \dots, m \right\}} \tag{25}$$

where *i* = 1, 2, . . . , *m*, and *j* = 1, 2, . . . , *n*.

Reference sequence *X0* is defined as (*x01*, *x01*, ... , *x01*, ... *x01*) = (1, 1, ... , 1, ... 1). The aim is to find an alternative whose comparability sequence is closest to the reference sequence. For this purpose, the Grey Relational Coefficient (GRC) between *xij* and *x0j*, is calculated using the equation:

$$\gamma(\mathbf{x}\_{0j}, \mathbf{x}\_{ij}) = \frac{\Delta\_{\rm min} + \oint \cdot \Delta\_{\rm max}}{\Delta\_{ij} + \oint \cdot \Delta\_{\rm max}} \tag{26}$$

where *i* = 1, 2, . . . , *m*, and *j* = 1, 2, . . . , *n*.

The values of Δij, Δ*min*, Δ*max*, ξ are defined as:

$$\Delta\_{\text{ij}} = \left| \left. \mathbf{x}\_{0j} - \mathbf{x}\_{ij} \right| \right. \tag{27}$$

$$
\Delta\_{\text{min}} = \text{Min}\{\Delta\_{\text{ij}} \colon i = 1, 2, \dots, m; j = 1, 2, \dots n\}, \tag{28}
$$

$$
\Delta\_{\text{max}} = \text{Max}\{\Delta\_{\text{ij}}, i = 1, 2, \dots, m; j = 1, 2, \dots n\}, \tag{29}
$$

$$
\zeta \epsilon \left[ 0, 1 \right] \tag{30}
$$

For this purpose, 0.5 is usually used as the distinguishing coefficient value ξ.

For the problem described in this paper, the total number of alternatives which are compared equals the number of optimal solutions obtained from Pareto, while the elements of the sequence are the values of two criteria: transformer daily aging and active power losses. Now, in order to determine the best solution, it is necessary to determine the Grey Relation Grade (GRG). GRG is calculated by using the following equation:

$$\Gamma(X\_0, X\_i) = \sum\_{j=1}^n w\_j \cdot \gamma(\chi\_{oj}, \chi\_{ij}) \tag{31}$$

where *wj* represents the weight factor for the corresponding index. The alternative with the highest GRG would be the best choice.

#### **4. Analysis and Results**

For the purpose of the analysis, the IEEE medium voltage radial distribution network with 33 bus has been modeled, Figure 9. The base voltage of the network is 12.66 kV, and the base power is 100 MVA. Line parameters are taken from [75]. Daily active power demand, without public stations, is 62,556 kWh. Total daily active energy losses, obtained from the power flow analysis amount to 2924 kWh. The lowest value of voltage in the system occurs during the 18th hour at bus 18 with the value of 0.9016 p.u.

#### *4.1. EVCSs Location Impact of Active Energy Losses*

The proposed methodology, Figure 7, enables us to perform an analysis of the impact of EVs on losses for an arbitrary number of EVCSs. In this paper, it is adopted that it is necessary to build four EVCSs. Based on the calculation of power flow, the best and worst cases in terms of energy losses are determined by the "brute force" algorithm. The results of the power flow calculations and algorithm are shown in Table 4 below. OBCs of EVs are assumed to operate with a unit power factor.

**Table 4.** Best and worst cases for EVCSs locations.


Table 4 shows a comparison between the parameters with or without charging stations and a comparison between the best and worst cases for charging station positions when energy losses are concerned.

There is a clear difference in losses depending on the location of EVCSs. In the best case, the increase in daily losses is 3% and in the worst case, it is 53% compared to the case without charging stations. These results show a significant influence of location of charging stations on network losses.

#### *4.2. Impact of Reactive Power from EVCSs on Energy Losses*

From the previously obtained results shown in Table 4, it can be seen that the obtained locations of EVCSs are in a row, so it can be assumed that they will not satisfy the needs of EVs that are far from such station. For this purpose, the spatial division of the network into four parts is proposed, so that each separate part includes eight nodes, see Figure 9 [76].

The proposed methodology, Figure 7, also allows an arbitrary division of the network with an arbitrary number of nodes.

After the division of the network, it is necessary to determine one location in each zone of the network where a station for electric vehicles would be placed. In this part of the paper, the impact of the power factor of on-board chargers on active energy losses is analyzed.

Accordingly, two scenarios are considered: the first scenario is that EVs are charged with a unit power factor; the second scenario is that vehicles are charged with the same active power as in the first case, and the maximum reactive power capacity from on-board chargers is used. As is mentioned in Section 2, it is assumed that on-board charges are oversized by 5.3%. The results are summarized in Table 5 below.


**Table 5.** Locations for EVCSs and network parameters for the two scenarios.

As can be seen in Table 5, for analyses of the nework and proposed network division, reactive power injection from on-board chargers does not have influence on EVCSs location. On the other hand, reactive power injection has an influence on daily acitive energy losess. Namely, it is shown that with a relatively small oversizing of the inverter, it is possible to reduce the loss of active power by about 3% without any interruption of an EV driver's habits or charging behavior.

Bearing in mind that voltage level of 12.33 kV, at which the *R*/*X* ratio is not large, produces the expected results obtained by these voltages. More precisely, the lowest voltage in the network in case the stations provide reactive power, is slightly higher compared to the scenario where the stations' power factor is 1.

#### *4.3. Impact of Reactive Power from EVCS on Transformer Aging*

The effect of reactive power injection on transformer aging is shown in Tables 6 and 7 below. The percentages in the tables refers to the gradual increase in power of EVCSs in relation to higher EV penetration. The results were obtained based on the formulas from Section 3.2. Two standardized distribution transformers 11/0.433 kV of different rated power, T1 (200 kVA) and T2 (250 kVA) have been compared in order to determine the appropriate power transformer for the analysed bus 12. Thermal parameters of transformers are taken from [33].

**Table 6.** LOL for different levels of reactive power for T1 for different additional penetration level.



**Table 7.** LOL for different levels of reactive power for the T2 for different additional penetration levels.

There is a large number of papers the optimal transformer sizing [77] and especially transformer sizing with the presence of electric vehicle charging [78]. In this paper, we only take these to values for rated power for purpose of analysis based on the overall load which is sum of EVCS load and the other load for the analysed bus 12.

In Table 6, LOL percentage for a workday and weekend are represented for different values of reactive power injection during the EV charging for T1.

Since the daily energy consumption from EVCS during the weekend, and mainly concentrated in the central part of the day, is about 18% lower than during the weekday [56], the values for LOL are expected to be lower than during the weekday.

From Table 6, it can be concluded that different values of reactive power injection by the chargers also have different values of LOL. Namely, on the basis of Table 6, it is shown that, with an increased penetration of 4% per working day, LOL values change between 0.0129% and 0.0145 %. From Equation (20), it can be concluded that if the LOL value is greater than 0.01333%, additional aging of the transformer occurs. As we can see, with increased penetration of 4% per working day, with some injection level of reactive power, there is no additional aging. If reactive power injection is zero, than LOL = 0.0145% which means (0.0148 × 180,000/100 = 26.1 h) 26.1 h over the 24 h period. That means 2.1 h of additional ageing every day. It can also be noted that for values of reactive power *Q* = 2/3*Qmax*, LOL has a minimum value for each considered level of EV penetration. Within this case, with a higher degree level of reactive power, LOL is reduced compared to the case when there is no injection of reactive power from the vehicle charger. For this case, in Figure 10, the reactive load of the transformer per hour without reactive power from the vehicle (in blue color) is shown, and the reactive load of the transformer with the reactive power injection from EVs (in yellow) is shown as well. As it can be seen, there is a compensation of reactive power in the largest portion of the day. This has resulted in reduced apparent power for almost the entire duration of the day, Figure 10.

There is similar situation with the increase of 7% for the weekend, with adequate reactive power support there will be no aging during the weekend. If there is no reactive power support there will be 27.9 h daily aging over the 24 h period. That means 3.9 h of additional ageing every day.

As is presented in Table 6 for the example of the EV penetration increase of 7% at a EVCS, the daily aging value is 36.6 h (for cos*ϕ* = 1) to 30.24 h (for *Q* = 2/3*Qmax*) representing a decrease of 17.37%. During the working day with 7% additional EV penetration, daily transformer LOL reaches 0.0202%. With reactive power injection of 2/3*Qmax*, there is a decrease of daily aging from 36.36 h to 30.24 h, which is a 17% decrease. Even though it is significant reduction, hot spot temperature reaches 140 ◦C, which represents a critical value at which gas bubbles appear [35].

From Table 7, it can be seen that transformer T2 with an increased penetration of 40% increases aging on the working day with 25 h aging over the 24 h period with no reactive power support. What is interesting is that LOL is higher at maximum reactive power support compared to where cos*ϕ* = 1 and it is 26.64 h.

**Figure 10.** Reactive and apparent power loading of T1 with (*Stl*) and without reactive power (*QEVCS*) from EVCSs (*Sl*), with 4% increase in EVCS demand.

As is shown in Table 7, with an adequate value of reactive power, the value of LOL is less than 0.0133%. Furthermore, as in the previous case, there is an improvement in decreasing the aging, using reactive power from EVs. The change in daily aging of the transformer and daily energy losses in the function of reactive power injection are presented in Figures 11 and 12, for transformer T2 with EV adiditional penetration of 40%. From the obtained characteristics, it can be seen that for some power injection there is aditional aging (above 24 h). So, it can be concluded that there is an optimal value of reactive power support that can have a positive influence on both transformer aging and energy losses.

**Figure 11.** Change in active daily losses in function of reactive power from public EVCSs.

In order to better understand the mutual influence of reactive power on transformer aging and active energy losses, solutions that satisfy the Pareto principle of non-dominance are shown in Figure 13. The presented solutions have been obtained by changing the reactive power range from 0 to *Qmax* with a step of 0.05*Qmax*.

**Figure 12.** Change in daily transformer aging in function of reactive power factor of the EVCS.

**Figure 13.** Pareto front optimal solutions.

In order to choose the final solution, a post-optimized analysis based on GRA has been performed. In this paper, the total number of sequences compared equals to the number of optimal solutions obtained by analysed Pareto front solutions. The elements of the sequence are the values of two criteria functions. The results from the GRA analysis are shown in Table 8, the while the optimal solution is pointed out in Figure 13.



A set of Pareto front optimal solutions, with points, is shown in Figure 12. For the selected solutions, it is shown that for the analyzed range of reactive power injection, the reduction of active energy losses leads to a faster aging of the transformer.

The solution involves reactive power injection of 0.8*Qmax*, which means that the losses of active energy amount to 3533.4 kWh and that the daily aging of the transformer is 24.012 h. In the end, the adopted solution enables the reduction of active daily energy loss in the amount of 3.1%. Furthermore, with these values of reactive power injection, the reduction of transformer LOL amounts to 4% compared to the case when an EV operates with the unit power factor and 10.3% when an EV injects maximum available reactive power *Qmax*.

#### **5. Conclusions**

In this paper, the potential of reactive power from OBCs in EVCSs to reduce losses of active energy and the aging of transformers is analyzed in detail. For the analyzed network and the number of EVCSs, it is shown explicitly that the injection of reactive power contributes to a reduc-tion in the active energy losses, but contributes significantly more to the reduction of transformer aging. Namely, for the analyzed IEEE 33 bus distribution network and obtained EVCS demand, it is shown that if on-board EV chargers in EVCSs operate in the capacitive mode and with maximum possible reactive power, the losses of active energy are reduced by 3% compared to the scenario when stations operate with the unit power factor. On the other hand, the impact of reactive injection power on the aging of transformers is significantly more pronounced.

The main contribution of this paper is that it has been shown that a relatively small oversizing of the OBCs enables a significant reduction in transformer aging in addition to the reduction of active energy losses in the system. The proposed methodology is based on generally known models for calculating power flows and the widely used model for transformer aging. It is important to emphasize that the proposed methodology can be applied to any network and any number of EVCSs. Furthermore, the proposed methodology takes into account the stochastic nature of the EVs, using appropriate PDFs and MC simulations.

It is important to point out that the contribution of reactive power from EVs at EVCSs depends on the transformer rating and the EVCS demand. In particular, it was shown that if the EVCS demand is 4% higher than planned, it is possible to prevent additional aging of a 200 kVA transformer by injecting reactive power in the range from 0.154*Qmax* to *Qmax*. Furthermore, for a transformer of 250 kVA, for a 40% higher EVCS demand, additional aging of a transformer could be prevented by injecting reactive power in the range from 0.096*Qmax* to 0.799*Qmax*.

For different levels of the EVCS demand, there are different values of improvement in LOL. For example, for the increase in the EVCS demand for 7% 200 kVA transformers, there is an 18% improvement in LOL % in comparison to the case where there is no reactive power injecting from EVs. Moreover, with a 60% increase in the EVCS demand for 250 kVA transformers, for optimal reactive power injection, there is a 23.33% improvement of LOL compared to a case where there is a maximum reactive power injection and a 10.3% improvement when an EV operates with the unit power factor.

**Author Contributions:** All the authors have contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This Ministry of Science and Technological Development, Montenegro.

**Acknowledgments:** The authors are thankful to the anonymous reviewers for comments that helped us improve the paper content.

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

#### **References**

