**Energy Storage Systems for Electric Vehicles**

Editor

**Erik Schaltz**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Erik Schaltz Aalborg University Denmark

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: https://www.mdpi.com/journal/energies/special issues/ storage EV).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-962-1 (Hbk) ISBN 978-3-03936-963-8 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

### **Contents**




### **About the Editor**

**Erik Schaltz** (Associate Professor) obtained his MSc and PhD degrees in electrical engineering from the Department of Energy Technology, Aalborg University, Aalborg, Denmark, in 2005 and 2010, respectively. From 2009 to 2012, he was an assistant professor, and since 2012, he has been an associate professor. Both positions are at the Department of Energy Technology, Aalborg University. In this Department, he is the program leader of the research program 'E-mobility and Industrial Drives', and the vice program leader of 'Battery Storage Systems'. He has been the main supervisor in four completed PhD projects, and has participated in more than 15 national and international research projects. He has been a guest and associate editor in several journals related to batteries and e-mobility, and he has authored or co-authored more than 100 publications. His research interests include usage of power electronics, electric machines, fuel cells, batteries, ultracapacitors, etc., in electric and hybrid electric vehicles. In addition, he is also focused on battery state-estimation, management (electric and thermal), and the modelling (electric, thermal and lifetime) of battery cells and packs.

### **Preface to "Energy Storage Systems for Electric Vehicles"**

This book is based on the papers published in the Special Issue of *Energies* on 'Energy Storage Systems for Electric Vehicles'. Energy storage systems for electric vehicles is a multidisciplinary field, and this book therefore covers a wide range of topics, e.g. electric engineering, mechanical engineering, control engineering, environmental engineering, material science, etc. Energy storage systems is a key technology for future transportation systems and it is my hope that this book will give a good overview of the relevant research disciplines applied for energy storage systems for electric vehicles.

I would like to thank all the authors for their contributions, and the *Energies* Editorial Office for the professional management of the Special Issue. I am very grateful for the support of the authors and the Editorial Office, as it would not have been possible to create this book without them.

> **Erik Schaltz** *Editor*

### *Article* **An Energy Consumption Model for Designing an AGV Energy Storage System with a PEMFC Stack**

**Roman Niestrój 1, Tomasz Rogala <sup>2</sup> and Wojciech Skarka 2,\***


Received: 2 May 2020; Accepted: 1 July 2020; Published: 3 July 2020

**Abstract:** This article presents a methodology for building an AGV (automated guided vehicle) power supply system simulation model with a polymer electrolyte membrane fuel cell stack (PEMFC). The model focuses on selecting the correct parameters for the hybrid energy buffering system to ensure proper operating parameters of the vehicle, i.e., minimizing vehicle downtime. The AGV uses 2 × 1.18 kW electric motors and is a development version of a battery-powered vehicle in which the battery has been replaced with a hybrid power system using a 300 W PEMFC. The research and development of the new power system were initiated by the AGV manufacturer. The model-based design (MBD) methodology is used in the design and construction of a complete simulation model for the system, which consists of the fuel cell system, energy processing, a storage system, and an energy demand models. The energy demand model has been developed based on measurements from the existing AGV, and the remaining parts of the model are based on simulation models tuned to the characteristics obtained for the individual subsystems or from commonly available data. A parametric model is created with the possibility for development and determination by simulation of either the final system or from the parameters of the individual models' elements (components of the designed system). The presented methodology can be used to develop alternative versions of the system, in particular the selection of the correct size of supercapacitors and batteries which depend on the energy demand profile and the development of the DC/DC converter and controllers. Additionally, the varying topology of the whole system was also analyzed. Minimization of downtime has been presented as one of many possible uses of the presented model.

**Keywords:** fuel cell; automated guided vehicle; hybrid energy storage system; model-based design; waveforms modeling; autoregressive models of nonstationary signals

#### **1. Introduction**

The use of electric drives in various types of vehicles is becoming increasingly popular. The growing use of such drives is due to the many advantages of electric motors compared to internal combustion engines. This is particularly evident in closed areas in internal transport where automated guided vehicles (AGVs) are heavily utilized. High torque, quiet operation, and zero-emissions are just some of the advantages over other primitive solutions. However, these vehicles have operational problems such as insufficient work duration and limited ranges. This is caused by the relatively low energy density of the energy sources used in these vehicles. The development dynamics of the basic energy sources used in AGVs, such as lithium-ion batteries, does not indicate that this problem will be solved in the short-term (within the next decade). For this reason, designers are searching for other energy sources that provide significantly higher levels of energy density while having the same advantages as modern batteries. One of the proposed solutions is to use hydrogen fuel cell stacks to power

AGVs. However, the power supply system itself, based on a fuel cell (FC), is much more complex than that of battery power. Usually, the fuel cell is supplemented with a hybrid energy storage system constituting an energy buffer that eliminates the imperfections associated with using hydrogen fuel cell stacks. This is due to the operational characteristics of the fuel cell that must produce electricity after commissioning; thus, the efficiency of electricity generation varies significantly depending on the load on the fuel cell. It is particularly unfavorable to operate the fuel cell at a very low/high load or with high dynamics of load change, which significantly reduces the efficiency of this device. Typically, an energy buffer comprises a battery and a set of supercapacitors with appropriately selected parameters. Control of the operation of the hydrogen fuel cell, integration of the appropriate battery size, buffering problems, multidirectional energy conversion, adaptation of electrical energy to various parameters, and hydrogen storage and supply all have specific characteristics and require appropriate adjustment of the power supply system to meet these energy demand characteristics. This means designing an optimal power supply system using a hydrogen fuel cell is a complex task.

As part of the work, the design of a hydrogen fuel cell-based power supply system for an existing AGV (Formica-1, AIUT Ltd., Gliwice, Poland) powered with a lithium-ion battery was undertaken, with an effort to preserve the vehicle's operational characteristics, minimize any structural changes, and significantly increase the vehicle's operating time.

Due to the limited number of commercially available FC's capacities, the fuel cell selection is usually based on average demand power. The authors note that the main problem in designing the entire power system based on FC is the correct selection of the energy buffer. Therefore, particular attention is paid to the selection of a hybrid energy storage system because the correct choice of this system allows one to adjust the characteristics of the entire system to one's needs, whilst minimizing the fuel cell's power.

The justification for using an energy buffer with an FC is to temper large fluctuations in power demand and to accumulate energy from regenerative braking. The energy buffer, in this case, corrects FC deficiencies as the FC is not able to rapidly increase its power output, has a limited peak power, and is not able to absorb braking energy. The nature of FC's work dedicates them rather to independent work in stationary applications. For traction applications, an energy buffer is needed that is tailored exactly to the nature of the energy demand.

To solve the selection problem for hybrid energy storage system in an AGV powered by a polymer electrolyte membrane fuel cell stack (PEMFC) outlined in this section, we urge the reader to refer to the literature review regarding the model-based design for methodologies used, FC modeling, for discussion on the components utilized, and for an overview of FC-based power systems (Section 2). Section 3 describes the assumptions of the general methodology for designing the entire system, and the assumptions, modeling methods, and model bases for the individual subsystems of the entire AGV, in particular, the energy demand at various operational states, the hybrid vehicle power systems model such as the FC, DC/DC converter model, as well as supercapacitors and batteries. For the system's application (Section 4), details of AGV development research involving the change of the power supply system from a battery system to a system based on FC is described. Section 4 highlights the identification of the vehicle's energy demand at various operating states, the model for this demand, the use of a power system model with the energy demand model for optimizating the newly developed power system based on a previously selected FC, and the selection of the structure and parameters of the buffering system energy. Optimizations were carried out through simulation experiments using developed models. The last section provides a detailed discussion of the results from earlier studies.

#### **2. State-of-the-Art**

Model-based design (MBD) methodology is often used to design complex mechatronic systems [1,2]. The methodology consists of building computer simulation models of the designed system and simulating its operation. The use of such a methodology is particularly beneficial in the design of systems requiring the cooperation of specialists from various fields and systems. In our case, the system crosses fields that include mechanical, electrical, and chemical sciences, which address difficult to describe phenomena and require personalized and very specific system control. A typical example of such a system is a drive and power supply system for vehicles, such as a power supply system consisting of a hydrogen fuel cell, battery, supercapacitors, and the respective control systems and energy flow processes. The purpose of using MBD methodology is to initially plan how to design a system, its operation, and control its parameters, all whilst meeting the criteria and fulfilling its desired functions. It is also possible to determine appropriate or optimal technical parameters of the individual subsystems, such as the technical parameters and the control method. Usually, the complete model includes not only the designed system, the power supply system, and the drive system, but also the entire facility on which the designed system is built, as well as external conditions affecting the operation of the whole. For vehicles, this is usually a power system model, the propulsion system model, and the entire vehicle model, and often includes the route model and the conditions they encounter. Depending on the situation, the complexity of the model should be adjusted to obtain satisfactory results [2,3]. If we have a prototype or a copy of the system available, we can determine an appropriate model using experimental tests, but if the system is in the concept phase we must build a model, e.g., a model based on the sum of the general theoretical phenomena occurring in the system. Likewise, the model for the power supply or propulsion system itself is much more complex, whilst the model of the routes and the whole vehicle is simplified or considers the relevant data to enable simulation. For hydrogen-powered electric vehicles, the most important and substantive input model is the hydrogen fuel cell itself, which forms the entire power supply system as the whole vehicle is driven and powered by such a system. Usually, choosing the correct system parameters makes the most sense when the vehicle is traveling along a fixed route or a finite and known set of routes where the load and driving conditions are set or predictable. With this knowledge, one can accurately determine the features of the power system. This is the case with certain types of vehicles such as AGVs or racing vehicles. However, if the load conditions and route conditions are unknown or difficult to predict, the task is much more difficult, and the results will not be as definitive as expected.

This section describes the problems associated with modeling vehicle system components and is was divided into two parts: The first concerns the hydrogen cells themselves and the second deals with the remaining elements of the energy conversion system. In these subsections, the authors focus primarily on the energy buffer, but FCs are also analyzed because it is the operational features of the FCs that have a significant impact on the selection of the energy buffer. Another element that affects the form and characteristics of the energy buffer is the fluctuating nature of energy demand and the need to recover energy under vehicle braking. Correct and detailed modeling of these sections of the whole system (and not only the energy buffering system) is, therefore, a condition for completing simulation experiments from which the energy buffering system will be selected.

#### *2.1. Modeling of the Fuel Cell Stack*

A hydrogen fuel cell is an electrochemical device that converts chemical energy via an electrolytic reaction directly into electricity. For modeling purposes, it is not necessary to know all physicochemical conditions related to energy production in the FC.

There are several classes of simulation models in the literature, which can be divided into three sub-groups; electrical, chemical, and experimental. Electric fuel cell models are used to compute power systems. This model treats the fuel cell as an element of an electrical circuit and does not include phenomena underlying electrical production. Phenomena such as particle diffusion, mass transport, and thermodynamic transformations are addressed in a chemical model.

The commonly used generic simulation model using a MATLAB/Simulink system includes two types of models: Simplified and detailed. Such models include the calculation of the irreversibles that affect the voltage drop of a cell during operation relative to the theoretical voltage resulting from a chemical reaction, which in turn results in changes in the energy characteristics. This is influenced by the following types of irreversibles: Activation losses, fuel crossover and internal currents, ohmic losses, mass transport, and concentration losses. The origins and descriptions of these irreversibles, as well as the modeling method, have been previously discussed [4,5]. When creating a simplified model, two points from the activation region and two points from the ohmic region are utilized from the polarization curve. However, for the construction of a detailed model, further data are required, such as the number of cells, nominal Low Heating Value (LHV) stack efficiency, nominal operating temperature, nominal air flow rate, absolute supply pressure, and the nominal composition of fuel and air; these are typically provided in the fuel cell manufacturer's documentation. When it comes to modeling fuel cell dynamics, current step and interrupt tests must be completed for a given cell. The necessary parameters to construct this part of the model are then determined from these tests or can be obtained directly from the manufacturer's data, because they depend on the fuel cell itself. If such tests cannot be performed, the data can be assumed from a recommended range [4]. Occasionally, the fuel cell manufacturer does not provide basic technical data for the FC, and in this case more tests on the system are required to determine a full range of parameters. Other parameters obtained in tests depend on the whole system in which the FC works and its load, and they are specific for a particular configuration of the system.

Other fuel cell equivalent circuit models for passive mode testing and dynamic mode design have been compared in [6]. This comparison includes the following dynamic models: Larmie [7], Dicks–Larminie [5], Yu-Yuvarajan [8], Choi [9], and shows that complex models are not always effective for practical applications. These four dynamic models are used to simulate the power-generating cell, whilst the passive equivalent circuit model represents the fuel cell which is not producing electric power. These models represent the response to an external electric stimulus to determine the condition of the fuel cell. Additionally, in [6], Page [10], and Garnier [11], passive models are compared. The work [6] does not present any relationship between passive mode test responses and dynamic mode performance.

Not all of the fuel cell's irreversibles are relevant under normal operating conditions. While commissioning and rated conditions are the most common conditions, overloading is not a common condition. Some systems do not function under FC operation with such overloading conditions at all. Therefore, irreversibles that affect work under such unusual conditions are not considered or modeled at all. However, sometimes this is needed, and irreversible mass transport and concentration losses must be modeled. A model for mass transport losses in the form of a theoretical model is presented in [12] and in the form of an empirical model in [4]. This model was developed to simulate transport phenomena in a proton exchange membrane fuel cell (PEMFC).

The hydrogen fuel cell is complex and expensive, and in systems with high dynamics of power demand where it is required to supplement such a cell with additional elements such as startup batteries, buffers, inverters/converters, then the whole system needs to be modified to handle a specific load. Testing such a system can be completed using a computer simulation model presented later in the article, but it is also possible to create a physical simulator which is a cheap alternative to testing. Such a solution built based on a programmable DC power supply, control interface, and software written in LabVIEW has been proposed for testing the entire system and acts as a guide in the development of power conditioning equipment [13] with the ability to work in steady-state and transient modes.

#### *2.2. Modeling of the System Using a Fuel Cell Stack*

To generate energy in FCs, it is necessary to use a hydrogen tank together with a hydrogen pressure reduction system and a control valve mediated by a controller to regulate the amount of hydrogen supplied on an ongoing basis. Oxygen is usually supplied from the air through a fan system to the fuel cell. It is also possible to supply oxygen from a high-pressure tank similar to the whole hydrogen supply system.

For large FCs (larger than 10 kW), the installation of the FC itself becomes very complicated and maintaining balanced operational parameters becomes a problem. These issues are the subject of separate research, and balance of the plant (BOP) [14] and incorrect configuration and selection of inappropriate operating parameters of individual elements of the FC system can lead to insufficient cell performance and rapid degradation of the cell. A simple solution to the complexity of the installation on the FC preparation side can be made by using an FC configuration based on a dead-end anode (DEA) structure. This type of installation, unlike the flow-through anode (FTA) configuration, significantly simplifies the need to prepare hydrogen and guarantees the appropriate humidity of the cell, ensuring close to 100% hydrogen use by controlling the (normally closed) purge valve [15,16]. This configuration is popular for low power FCs, but also developed for higher power applications. DEA installations operating differently to FTA are not managed by a regulated control valve and have to be purged periodically by the purge valve [17].

Since the fuel cell itself is an energy generator operating under specific parameters, usually this produced energy must be adapted to the purposes of the energy demand characteristics. If the power take-off is not variable, this system may be simpler, but with high variability of energy demand, it is necessary to consider the electric converter/inverter and energy buffer supplied via batteries or supercapacitors. Supervising the work of these devices can take place at various levels, most often at the basic level through ongoing control of the parameters of individual devices, and frequently at the strategic level by adapting the operating parameters not only to the ongoing demand but also to future demand.

Modeling the power supply load is a separate problem. A power supply load model takes the form of a specific load profile based on the behavior of the powered system and optimized with measurements taken during the experiment or by considering physical phenomena, e.g., a model outlining the dynamics of a moving object. The choice to develop this model depends on the design phase. If one has a prototype or physical copy of the system required to be powered, one can choose the first solution, but if one only has the concept or accurate documentation, the second solution is needed. Interesting solutions can be found in various works for modeling system fragments or the entire system oriented at determining specific parameters. An example model for a complete power supply system for hybrid vehicles is described in [18]. The modeled system consists of a hydrogen power supply, DC/DC converter, battery, inverter, electric motor, and the vehicle body. A complete model of the system was developed based on the experimental data. The model was then used to develop a power management control algorithm for fuel cell hybrid vehicles using a stochastic programing technique. This approach requires a complete system which can be subjected to a series of experiments. Another approach is to use model-based design (MBD), where a model is created at the design and concept stage, and numerical simulation experiments outline various potential solutions and determine the impact of various parameters on the system's performance (sensitivity analysis), or to formally optimize the system or its components [1,2].

Improved modeling of the Proton Exchange Membrane (PEM) fuel cell power stack for electric vehicles in which a separate oxygen tank supply system was used to improve performance is presented in [19]. Simulation calculations were oriented towards finding an optimal control strategy for the pressure that facilitates the output power according to the power demand of the load.

In addition to the holistic approach to modeling the hydrogen fuel cell system, researchers are interested in individual elements of the system. Furthermore, the hydrogen cell itself, with a series of tanks, controllers, control valves, and fans supplying air, may include power electronics which process and adapt energy to meet demand from energy buffering units, including various types of batteries and supercapacitors. Regulators and controllers are indispensable to these units and operate at various levels, and often operate with a complex strategy for a given application.

Selecting the power electronics for the FC's energy conversion system is quite a difficult task. The situation is additionally complicated by the fact that the energy-receiving system requires the conversion of energy to different voltages, types of currents, and their power simultaneously. We chose to only focus on work completed on general modeling of energy flow and power losses, not energy-electronic phenomena or their modeling. Therefore, only models for average value converters were assessed, and those cooperating with basic energy buffers and thus implementing alternative

strategies of constant current and constant output voltage were analyzed. There are no universal solutions to control of energy flow because the characteristics of the energy demand received from FC are application dependent. The selection of elements of the entire FC system is of interest to many scientists. An essential element in the system is the boost converter. A simple model of the cell as an electrical circuit has been previously described [20]. The model, taking into account a portion of the irreversibles appropriate to the nature of the application, is used to select the suitable type of DC/DC boost converter and to select the parameters of the energy storage constituting the energy buffer which compensates for rapid changes in energy demand. Various connection options (behind or before the converter) of the supercapacitors are also discussed.

For energy buffers, there has been significant progress in the development of the latest types of batteries. The multitude of solutions is not conducive when making optimal decisions, especially at the development stage of the system. Therefore, simple battery models using the most popular battery types are used. The basic battery models are lithium-ion, lead-acid, nickel-cadmium, and nickel–metal hydride [21], and their various parameters are also defined, including charge and discharge, temperature effects, and ageing. This enables the modeling of various connections cells in series and/or in parallel [22–25].

A "Theoretical Modeling Methods for Thermal Management of Batteries" review has been previously completed [26]. In addition to typical models, various new approaches are presented, e.g., in [27,28].

In [27] a novel lumped electrothermal circuit of a single battery cell was presented, including the extraction procedure of the parameters of the single-cell from experimental evidence and a simulation environment, given in SystemC-WMS for the simulation of a battery pack.

In [28] a new open-circuit voltage (OCV) model is proposed. The new model can simulate the OCV curves of a lithium iron magnesium phosphate (LiFeMgPO4) battery at different temperatures. It also considers both charging and discharging. The most remarkable feature from the different models, in addition to the proposed OCV model, is their integration into a single hybrid electrical model. A lumped thermal model is implemented to simulate the temperature development in the battery cell. The synthesized electro-thermal battery cell model is extended to model a battery pack of an actual electric vehicle.

Typically, the problem of choosing a buffer system includes what type of energy buffers will be selected, the size of the buffer, and the features of individual parts (batteries, supercapacitors). Buffer hybridization is a common solution which involves a combination of a supercapacitor with a battery and is outlined in [29]. Various configurations and sets of supercapacitors and batteries together with DC/DC converters are discussed in several papers [30–32]. The correct selection of the buffer parameters and the topology of this system allows one to overcome most of the FC's weaknesses. Selecting the optimal parameters and topology for these subsystems in the FC is important, as the FC is strongly dependent on the energy demand characteristics in the system [30,32].

Modeling supercapacitors (SC) requires consideration of the electrical, self-discharge, and thermal behavior. A comprehensive review of the modeling techniques is described [33–35]

The equivalent mathematical model derived from the electrical model, which was used to simulate the voltage response of the supercapacitor, is presented in [33].

The review presented in [33] discusses SC modeling, state estimation, and their industrial applications, intending to summarize recent research progress and stimulate innovative thoughts for SC control/management. For the SC modeling, state-of-the-art models for electrical, self-discharge, and thermal behavior are systematically reviewed, where the electrochemical, equivalent circuit, intelligent, and fractional-order models describing the electrical behavior simulation are highlighted. For SC state estimation, methods for state-of-charge (SOC) estimation and state-of-health (SOH) monitoring are covered, together with an underlying analysis of the ageing mechanism and its influencing factors.

The models which are described in the literature have various advantages and disadvantages, ranging from the ease of use down to the complexity of characterization and parameter identification. Work presented in [35] presents a comprehensive review and compares these models, specifically focusing on the models that predict the electrical characteristics of double-layer capacitors (DLC), showing the strengths and weaknesses of different available models and their various areas for improvement.

Experience in implementing the various applications of the hydrogen fuel cell system is very helpful when designing a complete system. One can find many interesting descriptions of applications with different degrees of maturity and covering both stationary and mobile applications in ground, water, and aerial vehicles. Research has described the various aspects of the whole system and its hybridization [36–39], current energy management and energy management strategy [40–44], energy control and processing [45,46], optimization of power systems based on fuel cells for matching operational parameters [47,48], power transmission in hydrogen cell-powered propulsion systems [49], and general aspects of the development of hydrogen cell-based systems [50,51].

#### **3. Model of Energy Transfer in the System**

A general methodology for building an energy transfer model enabling simulation experiments when designing a hybrid power supply (HPS) system based on a hydrogen cell stack for an AGV is shown in Figure 1.

**Figure 1.** A general methodology for building a model enabling simulation experiments for designing a hybrid energy storage system with a hydrogen cell stack for an automated guided vehicle (AGV).

Conducting simulation experiments requires the definition of the HPS system in the AGV. Since these vehicles are designed for close repetitive transport operations over long periods and have known operating conditions, i.e., speed and load, one can adapt the HPS system to individual needs, such as the demand for instantaneous power during a specific operating condition. At this stage, the criterion for assessing the designed HPS system should also be determined. For the next step, it is necessary to measure and identify the instantaneous power demand by the AGV at different operating conditions. These measurements should include the power demand for expected operating conditions over the planned route. From this, work can be completed on the data preprocessing, modeling, and validation of the models representing the power demand. These models are identified based on data from instantaneous power measurements at various operating conditions. Based on a set of such models, it is possible to simulate the power demand for the new AGV route and other operating conditions. A detailed discussion on this subject is presented in Section 3.1.

Simultaneously, by the defining power demand models, it is possible to create component models of the HPS system. It should be noted that these models of HPS can be identified based on additional measurements or characteristics provided by the manufacturers. More information about creating and identifying component models of an HPS including a hydrogen cell, models of storage components, and other auxiliary components, are described in Section 3.2. To validate the hybrid power supply system model, it is possible to provide the load in the form of played-back real values of instantaneous power demand and creating comparisons, e.g., concerning the current power supply system installed in AGV.

The proposed methodology described herein and in further sections of this manuscript allowed us to design a customized HPS system for operating conditions over a preplanned route. This can be achieved by conducting simulation experiments to find the optimal solution or a set of possible solutions which satisfy defined criteria. The optimal criteria can refer to finding the optimal battery or supercapacitor capacity for the HPS or other objectives. More information on this subject is discussed in Sections 4.2 and 4.3.

#### *3.1. A Generic Model for Instantaneous Power Demand*

This section describes a generic procedure for building a model to compute instantaneous power demand. This generic model is used to estimate the instantaneous power demand under the AGV's different operating conditions during working duty cycles. The model results are used as a load for the hybrid hydrogen power supply system model discussed in Section 3.2. The use of both models makes it possible to perform different simulation experiments, which allows one to examine different configurations of a power supply system with varying parameters. The generic model for instantaneous power demand is the first part of this model.

The presented methodology for building an instantaneous power demand model, ultimately to develop a new hybrid vehicle power supply system, depends on the available data sources. Two possible options defining the data source availability can be distinguished:


In the further sections of this script, only variant II is considered. This approach requires one to define the following operating conditions:


For the aforementioned values, under a combination of operating conditions, a bank of autoregressive models has been applied. These models are representations of signals which, for selected operating conditions, represent instantaneous power demand for the selected type of vehicle. The main task of the models, in detail, is to:


#### 3.1.1. Models for Stationary Conditions

The autoregressive model of the signal [53–55] of instantaneous power demand is given by the formula:

$$M(k):\ y(k) = E\{y\} + \frac{\epsilon(k) \times Var\{y\}}{A\{q\_p^{-1}\}} \tag{1}$$

where *y* is the instantaneous power demand, is the noise which follows a Gaussian distribution, *A q*−<sup>1</sup> *n* is a polynomial of the n order represented by *A q*−<sup>1</sup> *n* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>a</sup>*<sup>1</sup> <sup>×</sup> *<sup>q</sup>*−<sup>1</sup> <sup>+</sup> *<sup>a</sup>*<sup>2</sup> <sup>×</sup> *<sup>q</sup>*−<sup>2</sup> <sup>+</sup> *an* <sup>×</sup> *<sup>q</sup>*−*n*, and *E* - *y* , *Var*- *y* are the expected value and variance of the instantaneous power demand. The expected value and variance can be additionally represented by other linear or quadratic functions of *f*(*V*, *L*). To account for dynamic changes in the instantaneous power demand, the model can be represented in the frequency domain [55] using the following formula:

$$M(j\omega): P\_y(e^{j\omega}) = \frac{\epsilon(k) \times Var\{Y\}}{\left|A\{e^{j\omega}\}\right|^2} = E\{Y\}^2 + \frac{\epsilon(k) \times Var\{Y\}}{|1 + \sum\_{k=1}^p a(k) \times e^{-j\omega k}|^2} \tag{2}$$

where *Py ej*<sup>ω</sup> represents the power spectral density of the modeled signal. The above model can be applied under stationary operating conditions.

#### 3.1.2. Models for Nonstationary Conditions

Similarly to stationary conditions, a signal model can be built for the nonstationary conditions [55,56]. This applies to parameters such as acceleration, braking, and emergency braking, etc. The model to apply for this case has the following formula:

$$M(k):\ y(k) = E\{y\_d\} + \frac{\epsilon(k) \times Var\{y\_d\} \times En\{y\_d\}}{A\{q\_p^{-1}\}} + Tr\{y\} \tag{3}$$

where *Tr*- *y* is the linear model of the acceleration, deceleration ramp, etc. This part of the model can be determined using a least-squares criterion. *En*- *yd* is the envelope model established for the signal after removing the trend from the nonstationary signal *yd*. The model of the envelope can be evaluated for the following signal:

$$En\{y\} = \frac{|z[k]| + |z[k-1]| + \dots + |z[n-N]|}{N-1} \tag{4}$$

where *z*[*k*] is a module of an analytical signal obtained using a Hilbert transform [57] and *<sup>N</sup>* is a length of the moving average filter. The envelope signal can be represented by a regressive model given by Equation (1). If the envelope is flat and monotonical then a linear model can be used.

After removing the trend and by eliminating the second-order nonstationarity resulting from the variable variance, the frequency assessment of the model presented in the Equation (2) can then also be used.

#### 3.1.3. Model Validation

The validation of the model describing the route section under selected operating conditions can be calculated by using the following measures:

• Using a training data set to develop the model and validation data *yval*(*k*) , the following measures of model compliance can be determined:

$$Err\_{Energy} = \frac{\left| \int\_0^t M(k)dt - \int\_0^t y\_{val}(k)dt \right|}{\int\_0^t y\_{val}(k)dt} 100\% \tag{5}$$

where *<sup>t</sup>* <sup>0</sup> *<sup>M</sup>*(*k*)*dt* is the energy computed for the signal model, and *<sup>t</sup>* <sup>0</sup> *yval*(*k*)*dt* is the energy of the validation signal.

• The second measures (as a functional feature) of model compliance are executed with the use of relative error of power in the frequency domain:

$$Err\_{FreqStart} = \left| P\_y \left( e^{j\omega} \right) - P\_{\mathcal{Y}\_{val}} \left( e^{j\omega} \right) \right| \tag{6}$$

where *Py ej*<sup>ω</sup> , *Pyval ej*<sup>ω</sup> are the power spectral densities of the model obtained as an output of the model and the power spectral density of validation data, respectively. The measure determined here is a functional assessment in the frequency domain and it determines the difference in signal power for the frequency components. The selection of the model order is determined, based on the similarity of the power spectral density characteristics, to reflect the dynamics of the signal changes by the signal model.

#### 3.1.4. Combining Models

After validating the individual models representing the signal from instantaneous power demand, a selected scenario can be built which represents the AGV route. Usually, this route is planned and the AGV moves along the route under established operating condition parameters such as speed, load, etc.

Before creating a power demand model for a selected vehicle scenario, it was necessary to divide the scenario into appropriate route sections for which appropriate models would be assigned to generate the instantaneous power demand signals.

An important element when building full waveforms for the entire scenario was the points where the signals of the partial models would be combined. To combine waveforms of the individual models, it is possible to use the following window (Equation (7)), which is a modified version of the window previously shown in [58], the length of which can correspond to the length of the modeled waveforms:

$$w(y) = \begin{cases} \begin{pmatrix} -\frac{1}{2} \times d\_l + \frac{1}{2} \end{pmatrix} \times \begin{cases} 1 + \cos\left(\frac{2\Pi}{r} \times \left[\mathbf{x} - \frac{r}{2}\right]\right) \end{pmatrix} + d\_l, & 0 \le y < \frac{r}{2} \\\ 1, & \frac{r}{2} \le y < 1 - \frac{r}{2} \\\ \begin{pmatrix} -\frac{1}{2} \times d\_l + \frac{1}{2} \end{pmatrix} \times \begin{pmatrix} 1 + \cos\left(\frac{2\Pi}{r} \times \left[\mathbf{x} - 1 + \frac{r}{2}\right]\right) \end{pmatrix} + d\_l, & 1 - \frac{r}{2} \le y < 1 \end{cases} \tag{7}$$

#### *3.2. Hybrid Power Supply System Model for the AGV*

The model of the hybrid power supply system for the AGV was developed in the MATLAB/Simulink environment partly using the Simscape Electrical library components. This model is a numerical tool supporting the selection of elements for the hybrid power supply system. The block diagram of a hybrid power supply system is shown in Figure 2. The numerical model was built based on this block diagram (Supplementary Materials). This model could be used to optimize the parameters of the power supply system after a specific operation scenario for the AGV is chosen (length and diversity of the route, load, driving dynamics) and after assuming the optimization criteria (for example, minimizing the capacity of the main energy store).

**Figure 2.** A block diagram describing the hybrid power supply system for the AGV.

The electrical energy source in the hybrid power supply system was the fuel cell stack fueled by hydrogen. It was assumed that hydrogen was stored in a metal hydrides tank equipped with a pressure regulator [59]. The flow of hydrogen through the fuel cell stack was regulated by a proportional control valve. The control signal for this valve was generated using a hydrogen flow regulator. This regulator was a component of the fuel cell stack controller. The controller additionally protected the stack against operation from moving outside the safe operating range of the electrical and thermal parameters. In addition, the fuel cell stack controller contained the SCU (short circuit unit), which periodically short-circuited the stack and improved its performance [60]. Due to the operation of the SCU, it was necessary to install an auxiliary supercapacitor in the system, which maintained the supply voltage for the duration of the stack short-circuit, and additionally provided an energy buffer for rapid changes in the load current of the stack when the stack was not able to impulsively provide adequate power due to limitations imposed by its own dynamics and the hydrogen fueling system dynamics.

Electrical energy from the fuel cell stack was supplied to the main AGV power busbars through a DC/DC Constant Current - Constant Voltage (CCCV) converter working at a Constant Current (CC) or Constant Voltage (CV) output, where the output current setpoint for CC mode could be invariable or could be set by the stack load power regulator, which was part of the converter control system. The method for determining the output current setpoint depended on the configuration of the hybrid power supply system used and the method of its optimization.

A lithium-ion battery or supercapacitor could act as the main electrical energy storage. The reasons for using electrical energy storage (as energy buffer) together with the fuel cell stack, in traction applications, were the large fluctuations in the power demand and the need to accumulate energy from regenerative braking. The energy storage, in this case, complemented the deficiencies of the fuel cell stack, meaning the stack was not able to increase the output power impulsively, had limited peak power, and was not able to absorb braking energy. The nature of the fuel cell stack was rather dedicated to independent work in stationary applications. In traction applications, an additional energy storage device was necessary [30–32].

There was a management system between the main power busbars and the main energy storage, the primary role of which was to protect the energy storage against operation outside the safe range of electrical and thermal parameters. The management system also allowed for pre-charging of the main energy storage with energy from the fuel cell stack after starting the hybrid power supply system, which was needed when the supercapacitor acted as the energy storage. It was assumed that the energy storage could also be charged from an external energy source, depending on the adopted configuration of the hybrid power supply system and the scenario of the AGV operation.

While developing the numerical model for the hybrid power supply system, assumptions were considered from the practical conditions or were the result of previous preliminary analyses. The initial selection of the fuel cell stack was guided by the average power demand of the AGV and from economic criteria. The cheapest fuel cell stack was selected that would meet the AGV requirements according to preliminary estimates. It was assumed that a horizon fuel cell stack, type H-300, with 300 W power, a rated voltage 36 V, and rated current of 8.3 A would be used [61]. This stack consisted of 60 PEM fuel cells connected in series, low-temperature operation, powered by hydrogen from the pressure tank and oxygen obtained from atmospheric air. The nominal efficiency of the H-300 stack was 40%. This was a low power fuel cell stack that had a very simple "balance of plant" structure. The stack was equipped with three fans that provided cooling to the stack with a suitable amount of the air. The fuel cell stack was equipped with a factory controller that regulated the rotation speed of the fans by supplying them using the Pulse Width Modulation (PWM) method, and which controlled the hydrogen two-state valves. This fuel cell stack with factory controller functioned as a dead-end anode stack [15,16] without external humidification and hydrogen recirculation. It was assumed that the functionality of this controller could be extended to meet the needs of the power supply system under development by controlling the proportional hydrogen valve for flow-through anode operation [17], which was included in the numerical model.

The presented numerical model was primarily used to determine the flow of electrical energy in a hybrid power supply system, so several simplifications were assumed when developing this model. It was assumed that the fuel cell operated at a constant temperature and the airflow from which oxygen was extracted was always sufficient, regardless of the power load of the cell. The assumption regarding airflow was also fulfilled for the modeled fuel cell in the absence of external restrictions, which has been previously determined [62], where it was stated that even with the smallest used fan efficiency the cell worked with an air excess coefficient of ~20. Both thermal phenomena occurring in the hydrogen tank and the hydrogen release dynamics from the metal hydride storage were not taken into consideration. It was assumed that the hydrogen in the fueling system always had sufficient pressure to achieve the required hydrogen flow. Additionally, thermal phenomena in other elements of the power supply system were deemed to be negligible, assuming that they worked in optimal and constant thermal conditions. The phenomena related to the pulse operation of power electronic devices in the DC/DC converter were also not taken into account together with any ageing of the lithium-ion battery.

It was assumed that an external energy source was required to start the hybrid power supply system, ensuring the power needed to start the fuel cell stack and the stack controller, especially when the main energy storage was discharged. A low-capacity start-up battery could be used as an auxiliary energy source, which, if necessary, could be charged from an external source and, after starting the

power supply system, could be recharged from the main power busbars. The energy needed to start the power supply system was small; however, the starter battery model was omitted for simplicity.

Modeling of the AGV drive system had been simplified to just model the instantaneous power demand, while the demand for the power of the components of the drive system (inverters, motors) during vehicle movement and related to the operation of the vehicle's control, safety, and signaling systems also had to be taken into account. The instantaneous power demand model for a selected AGV operation scenario was created by submitting multiple data samples obtained during measurements made by a real AGV with different load states and with different operating states, both during steady driving and in dynamic states (acceleration, braking). The data samples were recorded for an AGV powered by a standard (factory) lithium-ion battery that was charged from an external source at the end of the operation. Then the data were subjected to filtering and processing as described in Section 3.1. It was assumed that the instantaneous power demand for a vehicle powered by a standard battery and in a vehicle powered by a hybrid power supply system with a fuel cell stack under the same operating conditions and load conditions was the same. In connection with the adopted method of modeling the AGV drive system, the phenomena associated with switching power electronic devices in the inverters of the vehicle's drive nodes were excluded from the research.

Optimization of the structure and parameters of the hybrid power supply system could be carried out considering various criteria by setting selected parameters for the numerical model and analyzing the obtained waveforms, both utilizing experiments performed by trial and error and by automatic optimization algorithms. Usually, the parameters of the fuel cell stack were assumed at the beginning of the optimization process because the choice of the stack was not very flexible and the rated powers of the available stacks were highly graduated. The choice of energy storage was more flexible, so the parameters of this storage device could be optimized. During the simulation, the ongoing analysis of the selected waveforms of electrical quantities were carried out in terms of exceeding the defined criteria (critical values). This analysis is conducted regardless of the applied optimization method in the numerical model. If such an exceedance occurred during the simulation, then the simulation would be stopped and the model would return an error code that determined which criterion had been violated. A total of fifteen different criteria were defined in the numerical model for the various components of the hybrid power supply system. These criteria are:


These criteria resulted from the catalogue of real element parameters of the hybrid power supply system and the conditions imposed by the elements of the AGV drive system (e.g., for inverters: The minimum and maximum supply voltage). Not all the criteria needed to be active at the same time. The selection of active criteria depended on which power supply parameters were unknown in the design aid process and which were imposed as project assumptions. For example, if the required minimum DC/DC converter power rating was unknown, then the criteria related to the power overload of the converter was turned off. If a specific DC/DC converter type needed to be used in the design, then in this situation the parameters of this converter should have been treated as project assumptions and the appropriate criteria values in the model were to be set, following the datasheet of the converter. In addition, the model for the hydrogen fueling system analyzed the hydrogen consumption during

the simulation and returned the appropriate error code if the hydrogen tank was emptied. In this situation, the simulation was also stopped.

The numerical model of the hybrid power supply system defined the allowable voltage range and allowable state-of-charge (SOC) range of the main energy storage. Exceeding the voltage or state of charge for energy storage was not treated as a critical error and did not stop the simulation. However, it affected the way the energy storage worked, which was signaled in the model by the appropriate status signals. If the minimum voltage or the minimum state of charge was exceeded during discharge, the energy storage could only be charged. If with such limited use of energy storage, there was an increased power demand from the AGV model, the voltage of the main power busbars would fall below the criterion value. Similarly, if the maximum voltage or maximum charge was exceeded during charging, the energy storage could only be discharged. If under this condition, the AGV model attempted to achieve a return of braking energy to the energy storage, then the voltage of the main power busbars would rise above the criterion value. Exceeding the criterion values of the main power busbars voltage was treated as a critical error and stopped the simulation by returning an appropriate error code. In this situation, the error code had to be analyzed together with the main energy storage status to detect the reason for stopping the simulation.

The simulation model developed in the MATLAB/Simulink environment was built according to the block diagram shown in Figure 2. In addition to the blocks outlined in Figure 2, it also contained elements that allowed one to record the simulation results in the MATLAB for automatic optimization, and it also contained elements that allowed an ongoing view of waveforms, important parameters, error and status signals for the trial and error experiments.

To model the fuel cell stack, a block from the Simscape Electrical library was used, which is described in detail in [4]; the addition of concentration or mass transport losses in accordance with the method presented in [5] was applied. The losses of concentration or mass transport Δ*V*trans are described by the equation:

$$
\Delta V\_{\text{trans}} = m \times \exp(n \times I\_{\text{FC}}) \tag{8}
$$

where the coefficients *m* and *n* are selected experimentally and *I*FC is the stack load current. To tune the model for the fuel cell stack's activation area and load losses (ohmic losses), the results of measurements completed on the real H-300 stack and the genetic algorithm were used. During measurements this stack operated as a dead-end anode with the factory controller. In addition, the concentration losses model was experimentally tuned to obtain the appropriate stack voltage drop when overloaded. The thresholds for stack voltage and current were taken into account, and when they reached the stack were disconnected from the load by the factory stack controller.

The power of the fuel cell stack's own needs ("balance of plant") was modeled as being linearly dependent on the stack load power. The H-300 stack balance of plant was very simple (containing only fans, a controller, and hydrogen valves). However, it would be possible to model the balance of plant for a more sophisticated system, if the power demand characteristics of the components were available.

The fuel cell stack controller model included a hydrogen flow regulator that generated the *FFR*(ref) control signal for the hydrogen proportional control valve, which determined the flow through the anode of the stack. The principle of proportional control for this regulator was derived from the equations of the fuel cell stack model used in MATLAB presented in publications [4] and [5]. This regulator calculated the hydrogen flow needed to meet the hydrogen needs of the fuel cell stack at a given load current and a given hydrogen utilization. With a set number of cells in the stack, stack temperature, pressure and purity of hydrogen, the control principle is described as follows:

$$FFR\_{\text{(ref)}} = C\_{FFR} \times \frac{I\_{\text{FC(avg)}}}{I\_{\text{HF2\%}(\text{ref})}} \tag{9}$$

where the value of the coefficient *CFFR* can be determined using the relationships given in [4] or [5]. *U*H2%(ref) is the percentage setpoint of hydrogen utilization and the input quantity is the average

current *I*FC(avg) of the fuel cell stack. This regulator ensured the hydrogen flow when the load current in the stack increased dynamically, which in turn ensured the rapid opening of the hydrogen control valve and prevented a voltage drop in the stack. Due to the strong averaging of the stack load current at the regulator input and the dynamics of the control valve (which was modeled by first-order inertia), the setpoint of hydrogen utilization by the stack should have been slightly less than the nominal hydrogen utilization to ensure proper fueling of hydrogen in fast transient states. The nominal hydrogen utilization could be calculated using the stack's rated parameters and relationships, as given in [4]. For an H-300 stack, it was 83%. When starting the hybrid power supply system and its associated transient states, the flow regulator ensured a sufficiently high initial hydrogen flow. The stack controller model contained a stack power demand model (power of its own needs), implemented as an array of values with interpolation that models "balance of plant". This power demand was included in the load model of the main power busbars.

The characteristics of an H-300 stack for nominal hydrogen utilization, obtained by the numerical model and tuned based on the results of the measurements are presented in Figure 3.

**Figure 3.** The characteristics of the H-300 fuel cell stack obtained from a numerical model at a temperature of 40 degrees Celsius, absolute hydrogen pressure of 1.5 bar, nominal hydrogen utilization of 83%. The stack's rated parameters are 36 V, 8.3 A, 300 W, with nominal efficiency of 40%. (**a**) The current-voltage characteristic, (**b**) the fuel flow rate vs. stack load power, (**c**) the stack load (gross) power and available (net) power vs. the stack load current, and (**d**) stack efficiency and system efficiency (stack efficiency taking into account "balance of plant").

The "auxiliary supercapacitor" block in Figure 2 also contains the controller that charges the auxiliary supercapacitor in a precise manner during the power supply system start-up to the required minimum voltage and then connects it to the output busbars of the fuel cell stack.

The DC/DC converter model was an average value model that considered the efficiency characteristics implemented as an array of values with interpolation and a no-load current. Additionally, it included the characteristics of the output power limitation as a function of the converter supply voltage. The output power limitation could be used interchangeably with the power threshold detection (and error code) depending on the purpose of the simulation test. The setpoint of the output current in CC mode could be constant or it could come from the regulation of the fuel cell stack load power. It was a Proportional Integral (PI) type, anti-windup regulator.

The model of the main energy storage management system, depending on the state of charge and voltage of the energy storage, allowed for its normal operation (such as charging and discharging) or to operate with restrictions (only discharging or only charging). This allowed a pre-charge of the energy storage after starting the power supply system if this function was needed.

The main energy storage model contained models of supercapacitor or lithium-ion battery, alternatively selected.

The model of the main power busbar loading system is included in the "AGV" block shown in Figure 2, which loads the power supply system with the power required by the AGV. Additionally, the power for the fuel cell's own need is represented by the "auxiliary DC/DC converter" block in the same diagram. The power required by an AGV is shown in the value tables, containing samples of the power demand while driving and samples of the vehicle's own needs.

The hybrid power supply system model included control signals that enforced the appropriate order of switching on its elements during start-up, thus mapping the operation of the real system.

#### **4. Optimization Process Use Case**

#### *4.1. Automated Guided Vehicle (AGV)*

An automated guided vehicle is designed for the transport of goods, materials, and semi-finished products as part of internal transport carried out in closed production or warehouse halls. The vehicle is designed to travel at ground level and can transport goods directly by itself by placing a loaded pallet on the upper loading surface of the vehicle or by pulling an attached transport trolley. The vehicle moves independently throughout the hall, performing tasks independently without human assistance in accordance with its pre-planned action and along a planned route. Usually, the vehicle travels along fixed routes according to a fixed schedule adapted in conjunction with the production cycle. The reproducible nature of the travel route and loads is important for matching the planned hydrogen fuel cell stack-based power supply system to the application. The vehicle monitors the surroundings via a sensor system to avoid collisions with them. The vehicle is powered by a lithium-ion battery placed in an easily accessible and replaceable cassette and the drive consists of two electric motors. A low-power AGV (Formica-1, AIUT Ltd., Gliwice, Poland) was used in this research.

#### *4.2. Instantaneous Power Demand Model—Route Scenario*

#### 4.2.1. Identification Experiment

Identification of the instantaneous power demand model whose output is the input of the hybrid power supply system model requires proper planning of the identification experiment. The first step of these activities was to develop a common test plan for different operating conditions that take into account various stationary and nonstationary operations carried out on the real AGV.

The experiment was completed for different operating conditions at different route sections. The experiments are listed in Table 1. Due to the autonomous operation of the AGV control and the stochastic nature of the interaction between the vehicle surface and the AGV, the selected experiments were repeated several times and the average results obtained in this way were used for testing the signal models.


**Table 1.** List of conducted experiments using AGV.

During the conducted experiments, the following values were recorded: The voltage and the current returned by the batteries, the current values recorded on the main drive, and the current value on the stabilizing converter. Additionally, measurements of the resistance of the drive that was not directly measured were made. Due to these measurements, it was possible to record the instantaneous power demand. A schematic of the measuring system is shown in Figure 4. The data were recorded using an oscilloscope and with a sampling frequency of 100 or 50 kHz, depending on the duration of the selected route section.

**Figure 4.** Block diagram of the low-power AGV drive system including the oscilloscope probes used to measure the power demands.

Restrictions on the safety and control of AGV are specified in the standard [52], including various responsibilities imposed on manufacturers and users. Due to the above reasons, the AGV was equipped with a logger system to record or monitor selected parameters during operating conditions around the route.

Selected logger data was used to observe the operating conditions. The data gathered concerned the rotational speeds of integrated Tekno TO-62 drives (left and right drive nodes according to Figure 4) equipped with an induction motor (nominal power 1.18 kW), the mechanical transmission with gear ratio 8.12 with a maximum continuous wheel torque of 25 Nm, and the power with a nominal voltage of 33 V. This element was also equipped with a 48 VDC nominal brake and a 5000 pulses speed encoder. The data were recorded using the AGV's inbuilt logger with a sampling frequency of ~2.5 Hz and were not synchronized with the instantaneous power demand signals recorded with the use of an

oscilloscope. The recorded speed data were used to identify the operating conditions associated with the route section covered and its identification. This is a necessary part of the proposed approach, in particular, which is forced by conducting measurements in situ conditions where synchronization of measurements with the logger data (operating condition) was not possible. Figure 5 shows selected waveforms, speed signals from the logger, and the auxiliary computed signals.

**Figure 5.** Example waveforms achieved from the logger and the additionally computed auxiliary signals.

To synchronize the measurements and thereby identify individual sections of the route for which signal models can be developed, the data were preprocessed by determining the auxiliary signals which were used to enhance the recognition of different operating conditions (some examples are shown in Figure 5), using resampling methods and identifying common starting points for both sources of data. For the obtained segments of the labeled measurement data related to route sections and selected operating conditions, models for stationary and nonstationary conditions were identified defining the banks of models.

Figure 6 shows the selected labeled measurement data based on previously determined data labels from the logger and auxiliary data. Based on the data labels, it was possible to segment the data and create a bank of signal models representing the instantaneous power demand for selected operating conditions.

**Figure 6.** An example of labeling of the measurement results (in this case, the measurement of the current on the left inverter) based on operating conditions (in this case, based on the speed of straight route both backwards and forwards).

#### 4.2.2. Instantaneous Power Demand Model for a Selected Scenario

The route scenario presented in Figure 7, developed for the AGV, consists of a section of the slalom route with the load (marked in red) and the rest of route unloaded.

**Figure 7.** An example scenario of the AGV route for which the instantaneous power demand model is being built.

For the presented scenario, a signal of instantaneous power demand for a section of the route without load, shown in Figure 7, was modeled with the use of a set of models. For this section of the route the following models were prepared:


The needed lengths (number of samples) of the individual waveforms computed by the models was determined based on information about the time necessary to achieve the required speed (in the case of braking and accelerating). The model output did not compute any velocity, as this value could be read from the inverse of the average power demand versus the average velocity which had been identified based on the collected data sets presented in Table 1. This was determined using the linear approximation *Pinst*−*ave* <sup>=</sup> *<sup>C</sup>*<sup>1</sup> <sup>×</sup> *vave* <sup>+</sup> *<sup>C</sup>*2, where *<sup>C</sup>*<sup>1</sup> is 351.4 Ws <sup>m</sup> and *C*<sup>2</sup> is 279.3 W (valid for the average velocities *vave* between 0.3 and 0.8 m/s). The required number of samples for a constant speed period could be determined from the required length of the route and sampling frequency.

The waveforms were generated for the considered route section shown in Figure 7 by using previously listed models. The errors of the individual models are presented in Table 2. An example assessment of the selected model (with constant speed) for the instantaneous power signal distribution in the frequency domain is presented in Figure 8. The calculated errors were obtained from the test measurement data. Next, the individual waveforms generated with signal models were combined using the window indicated in Equation (7). An example of the joined data from two models is shown in Figure 9.

**Table 2.** List of models and their relative errors for the considered scenario.

**Figure 8.** A periodogram estimate of Power Spectral Density (PSD) of instantaneous power for a model with twentieth order and the PSD measurement data for a straight profile at a constant speed.

**Figure 9.** The combination of the two waveforms generated from the signal models for nonstationary and stationary signals (**a**) without using a prepared window; (**b**) with the use of window; and (**c**) with a used window to combine the two signals from the models.

The relative energy error for the modeled route was mainly (excluding the influence of windowing) a weighted average of errors for the individual models used to determine the instantaneous power demand, and the weights of this average resulted from the fraction used of the individual signals in the whole combined waveform.

The example presented in this section shows the possibility of modeling the instantaneous power demand using a well-known class of autoregressive models of signals. The proposed approach requires a simulation experiment by recording the instantaneous power demand for various operating conditions. The advantage of the presented method is the lack of interference from the AGV software, including its control system where this information is often unavailable due to company intellectual property issues, and there is no need to create a dynamic vehicle model.

#### *4.3. An Example of Using the Model to Optimize the Hybrid Power Supply System*

To demonstrate the practical use of the numerical model for an AGV hybrid power supply system, a short model route was designed, as outlined in Figure 7. The AGV moved with a load of 1.2 tons along a model route (marked in red) and then moved along a route without a load (marked in blue). The loading and unloading points and control points are marked in green, where the vehicle stopped for a maximum of a few seconds. When driving without a load, the vehicle accelerated and braked more rapidly than when driving with a load.

The demand for power (*P*AGV) for the AGV during the model route was determined by the measurement results obtained for a real AGV, using the processing methods described in Section 3. An example of the waveform of the power demand while driving is shown in Figure 10. The results of measurements for the power demand when the vehicle was stopped were used to model the AGV stoppage.

**Figure 10.** The power demand for an AGV driving on a model route; brown circle—loading point, blue circle—unloading point, green circles—control points (the points marked with circles coincide with points in Figure 7).

It was assumed that the model cycle for AGV operation included: Waiting time for the first drive after starting the power supply system of 30 s, five drives along the model route, a standstill after each drive, and waiting time for switching off after the driving cycles of 10 s. An example of the AGV's power demand waveform during the operation cycle is shown in Figure 11. The standstill time after driving was one of the parameters that changed during the optimization process and, in this case, was 255 s. The total electric energy consumption during the entire operation cycle was 131.2 Wh, with an average power demand of 236 W.

**Figure 11.** The AGV's demand for power during the model operation cycle.

The optimization aimed to minimize the energy storage capacity and the duration of stops between drives, assuming that all the energy needed to power the AGV came from the fuel cell stack, i.e., the state of charge of the energy storage should have been the same after an entire operation cycle as at its beginning. In the model hybrid power supply system, none of the fifteen electrical parameter criteria could be violated. It was essential that, during the simulated vehicle operation between the main energy storage voltage and the threshold (criterion) values of the supply voltage of the load system, a safety margin of ~3 V was maintained. These threshold values were 30 and 60 V, respectively. Additional parameters that were tuned in the optimization process and which had an essential impact on the results obtained were the allowable range of the energy storage voltage (in particular the energy pre-charge storage voltage), the output voltage of the DC/DC converter in CV mode, and the output current of the DC/DC converter in CC mode. An important result obtained from the model was the hydrogen consumption for the assumed operation cycle of the AGV, which allowed one to choose the required capacity of the hydrogen tank.

The preliminary simulation tests were carried out for the assumed operation cycle, assuming that the main energy storage was a LiFePO4 battery with a capacity of 10 Ah and a rated voltage of 48 V, achieved from a real test bench. This is a low-cost battery that can provide a large enough impulse discharge current, with an expected value of ~3 C without degradation. This battery is built of sixteen 10 Ah prismatic cells connected in series. The energy storage model was tuned using optimization methods and the datasheet from the battery cells. It was assumed that the battery was pre-charged before starting the AGV power supply system (the initial state of charge was 50%). The selected results of preliminary simulation tests are shown in Figure 12. The results for the turned off SCU are presented so that transients that are associated with the operation of the SCU do not impair their readability.

**Figure 12.** The results of preliminary simulation tests obtained assuming that the main energy storage was a LiFePO4 battery with a capacity of 10 Ah; (**a**) the total demand for power; (**b**) the load power of the fuel cell stack; (**c**) the hydrogen flow rate; (**d**) the output current of the DC/DC converter; (**e**) the voltage of the battery and (**f**) the state of charge of the battery.

Figure 12a shows the total demand for power *P*req, including *P*AGV power for the AGV and *P*Aux power for the hybrid power supply system's own needs. Figure 12b shows the load power of the fuel cell stack *P*FC. It can be seen that the stack was utilized optimally and correctly (without overloading) throughout the entire operating cycle of the AGV stack and was loaded with power close to the rated power.

Figure 12c shows the hydrogen flow rate *FFR* obtained at the control valve output. Integrating this waveform, after converting it to the standard liters, it can be determined that ~141 L of hydrogen were consumed during the entire AGV cycle of operation, which corresponds to 423 Wh hydrogen energy, assuming that the change in enthalpy of formation was equal to the lower heat value [5]. When considering the energy consumption of the AGV from Figure 12a (~144 Wh), the efficiency of the hybrid power supply system was 34%. Figure 12d shows the current *I*DC/DC output waveform of the DC/DC converter. It can be seen that this converter works permanently in CC mode and the setpoint of the output current is constant and equal to 5 A. Figure 12e shows the *U*Batt voltage waveform of the 10 Ah battery. This voltage was practically constant, which resulted from the relatively rigid discharge characteristics and the slight changes in the state of charge *SOC*Batt% for this battery, shown in Figure 12f. A practically constant voltage of the battery at a constant output current of the DC/DC converter caused the fuel cell stack to be loaded with constant power throughout the entire operation cycle. The standstill time needed to restore the battery state of charge to the condition before driving was 257 s.

The battery used in the preliminary simulation tests had too large a capacity for the energy demand for the selected scenario of AGV operation, which was uneconomical. In subsequent tests, the battery capacity was reduced to a value of 0.2 Ah; this still ensured the correct operation of the AGV. The capacity of the battery was chosen so that its SOC varied from 20% to 80% during the operation of the AGV. It was assumed that the battery was pre-charged to an SOC of 20%, with an SOC of at least 70% required to start the vehicle. Therefore, when the power supply system was turned on, the battery was pre-charged by the fuel cell stack. The results did not change significantly. A slightly poorer utilization of the fuel cell stack was obtained while the AGV was driving. This was due to greater voltage drops in the smaller capacity battery, which, with a constant output current of the DC/DC converter, resulted in a decrease in the stack load power. The consumption of hydrogen increased to ~148 standard liters due to the initial charging of the battery, which absorbed 7.3 standard liters of hydrogen and lasted ~90 s. After omitting the supercapacitor pre-charge energy, the system efficiency was similar to previously, at ~34%. The used battery had a capacity of only 0.2 Ah, which was practically impossible due to too low current values for batteries with such a small capacity. Simulation tests using a selected scenario for vehicle operation and a 0.2 Ah battery were not of practical importance but were used to present the issue and how to use the model as a design aid. Therefore, during these tests, no criterion values for battery current were determined. A similar battery operation regime with real, higher capacity could be obtained for the real scenario AGV operation with higher energy demand. The use of battery capacity in such a work regime seems optimal, but it should also be noted that a battery working continuously in such a regime can quickly degrade.

Due to the possibility of quick battery degradation, further simulation tests were completed using a supercapacitor as the main energy storage. The most important advantage of supercapacitors, outlined in [30,32], is their use as an energy buffer for a fuel cell stack in traction applications due to their higher power density compared to batteries. Supercapacitors have higher efficiency and a higher number of charge and discharge cycles without degradation compared to the battery. By optimization, the criterion for the minimum capacity of the supercapacitor was chosen. The capacity of the selected SC maintained the voltage of the main power busbars over their required operating range whilst preserving a safety margin from the criterion values. The results of these simulation tests are given in Figure 13. It was assumed that the supercapacitor was not pre-charged, so its pre-charging was implemented by the fuel cell stack on start-up of the power supply system.

**Figure 13.** Selected results from the simulation tests obtained using a supercapacitor with a capacity of 35 F as the main energy storage; (**a**) the total demand for power and (**b**) the load power of the fuel cell stack.

The main difference between the operation of the power supply system with a LiFePO4 battery and a supercapacitor is due to the different discharge characteristics of these energy storage devices. The voltage *U*SC of the supercapacitor changes significantly more when discharged than the voltage *U*Batt of the LiFePO4 batteries (at least in the SOC range between 20% and 80%). There is a constant output current of the DC/DC converter in CC mode, resulting in a much worse utilization of the fuel cell stack due to significantly lower stack load power at a low supercapacitor voltage (Figure 13b). The presented results were obtained using a real supercapacitor model which consisted of 44 component supercapacitors with a capacity of 385 F each in a 2P22S connection system (two connected in parallel, 22 in series), which gave a resultant capacity of 35 F and a rated voltage of 61.6 V.

Using a supercapacitor, the hydrogen consumption was now 251 standard liters, of which 19.4 standard liters were required for pre-charging of the supercapacitor. The pre-charge time was approximately 340 s. The energy consumption of the AGV was 243.7 Wh (Figure 13a) and the hydrogen energy used for the AGV operation was ~695 Wh. The obtained power supply system efficiency was similar to that previously shown in the battery case, but the vehicle's standstill time after driving required to charge the supercapacitor to its pre-driving condition was 670 s. Such a long standstill time was due to the low power utilization of the fuel cell stack at low supercapacitor voltage when the stack provided slightly more power than that of the AGV's own needs.

Further simulation tests were completed with the fuel cell load power regulator turned on, which affected the setpoint of the DC/DC converter output current under the CC mode. As a result of the re-optimization, the capacity of the supercapacitor in the main energy storage was reduced to 24.5 F, obtained by connecting the 2P22S component supercapacitors with a capacity of 270 F each. The results of the simulation tests are shown in Figure 14. It can be seen that the utilization of the fuel cell stack was again optimal and correct (Figure 14b). At the same time, a much shorter standstill time (240 s) after driving is needed to charge the supercapacitor to its pre-driving condition. The hydrogen consumption was 150 standard liters, of which ~14.5 standard liters are required for the initial charge of the supercapacitor, which lasts ~200 s. The energy consumption of the AGV is 139.6 Wh (Figure 14a) and hydrogen energy used for AGV operation was ~406.5 Wh. Again, a system efficiency of ~34% was obtained, but the required standstill was much shorter than for a 35 F supercapacitor. Improvement in the use of the fuel cell stack compared to the previous simulation was obtained as the stack load power regulator increased the setpoint of the DC/DC converter output current in the CC mode at low supercapacitor voltage, just enough not to overload the stack. The fast transients shown in Figure 14b–d result from when the supercapacitor was charged to a certain maximum voltage, and the DC/DC converter then went into CV mode. In this situation, the load power of the fuel cell stack dropped sharply, and after the converter returned to CC mode, it increased again rapidly.

**Figure 14.** The results of simulation tests obtained assuming that the main energy storage is a supercapacitor with a capacity of 24.5 F and the fuel cell stack load power regulator is turned on; (**a**) the total demand for power; (**b**) the load power of the fuel cell stack; (**c**) the hydrogen flow rate; (**d**) the output current of the DC/DC converter; (**e**) the voltage of the supercapacitor and (**f**) the state of charge of the supercapacitor.

It can be seen that the presented simulation studies achieved the optimization goals and aid the design of the hybrid power supply system. The given results were valid for the assumed scenario of the AGV operation. In subsequent simulation tests, how the SCU operation affects the obtained optimization results was checked.

The SCU short-circuited the fuel cell stack for 100 ms every 10 s. However, transients lasted longer than 100 ms and were associated, among other things, with the need to recharge the auxiliary supercapacitor which was partially discharged when supporting the DC/DC converter supply voltage during stack short-circuit. The obtained results were accurate and similar to those presented in Figure 14, the main difference was that the standstill time after driving was extended to 250 s to charge the supercapacitor back to its condition before driving. Hydrogen consumption increased slightly to 155 standard liters and the power supply system efficiency decreased by 0.7%.

#### **5. Discussion**

The energy transfer numeric model for hybrid power supply system was the main objective of the article. The model consisted of the elements of supply system such as: Power converter, energy buffer, FC, and control devices. For the load of the supply system, the instantaneous power demand model was used. This model represents a generic instantaneous power demand for a given route for stationary and nonstationary conditions for an AGV under experimentally determined selected operating conditions. The main reason for developing the generic instantaneous power demand model is its simplicity and the fact that it does not require any additional information about the subsystems of AGVs and any supplementary information from the AGV manufacturers. It should be emphasized that the method proposed here could be improved by using further models that consider nonstationary conditions, i.e., that include nonstationarity of frequency components.

Section 4.3 describes the optimization process of the hybrid power supply system parameters for a model route and an AGV operation scenario. A similar process of optimization and design-aid can be completed using the real vehicle operation cycle by obtaining the route parameters and driving scenario from the vehicle manufacturer or vehicle user. In a situation where the AGV cannot make stops after driving, which allow for recharging of the main energy storage as described in Section 4.3, it is possible to minimize the capacity of the energy storage using a numerical model with the appropriate utilization of a fuel cell stack, under the assumption that the energy storage will be recharged using an external source after the entire operation cycle of the vehicle. In this situation, a compromise can be made between the capacity of the main energy storage (lithium-ion battery) and the hydrogen consumption, to consequently determine the capacity of the hydrogen tank. It is possible to use a stack load power regulator to intentionally reduce stack utilization and not exceed the assumed hydrogen consumption.

Further development of the presented numerical model, in the part related to the production of electrical energy, may include issues such as the dynamics of hydrogen release from the metal hydride storage under various operating conditions, and the dynamics of the cell response to a change in the hydrogen flow at the control valve output by considering the dynamics of the hydrogen distribution inside the cell. In the section of the model in which the power demand for the AGV is modeled, a dynamic model of the vehicle can be used which requires the drive torque for the specific route conditions, vehicle load, and the traction parameters (speed and acceleration). The hybrid power supply system power demand can be calculated based on the required drive torque in the dynamic models for the vehicle's drive nodes. Further development of the numerical model requires conducting additional tests on the fuel cell stack together with the hydrogen tank and examination of the real AGV to collect additional data and verify the extended numerical model.

The overall assessment of the proposed solution was carried out quantitatively for selected model elements (models of instantaneous power demand and the tuned hydrogen cell model). For other elements of the model, the assessment is qualitative as it is dependent on the specific instance of the AGV.

#### **6. Conclusions**

The research aim was to develop a model of a hybrid power supply system with a fuel cell stack for designing an energy storage system. The power supply system model is a numerical tool supporting the design and optimization of the power supply system following the MBD methodology. The article presents an example of the process of energy storage optimization for the AGV hybrid power supply system, which implements an example cycle of operation. Data for the AGV power demand model were obtained from measurements carried out on a real factory battery-powered AGV. These data were processed, and allowed the extract models for standard route fragments that can be interpolated under various load conditions. Based on these fragmentary models, it is possible to develop a power demand model for any route and optimize the hybrid power supply system for this route.

The conclusions from the generic instantaneous power demand model are:


Conclusions related to the model of the hybrid power supply system:


**Supplementary Materials:** The following are available online at http://www.mdpi.com/1996-1073/13/13/3435/s1: Preview of simulation model, graphical abstract, figures.

**Author Contributions:** Conceptualization, W.S.; data curation, T.R.; formal analysis, R.N., T.R., and W.S.; funding acquisition, W.S.; investigation, R.N., T.R., and W.S.; methodology, R.N., T.R., and W.S.; software, R.N. and T.R.; supervision, W.S.; validation, R.N. and T.R.; writing—original draft, R.N., T.R., and W.S.; writing—review and editing, R.N., T.R., and W.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** Publication supported by the Rector's professor's grant. Silesian University of Technology grant number 10/060/RGP18/0102.

**Acknowledgments:** The authors would like to express their gratitude to AIUT Ltd., Gliwice, Poland for providing the research objects and materials for the construction of this research and to Barbara Skarka for support in developing the English version of the article.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Battery-Aware Electric Truck Delivery Route Exploration †**

#### **Donkyu Baek 1, Yukai Chen 2,\*, Naehyuck Chang 3, Enrico Macii <sup>4</sup> and Massimo Poncino <sup>2</sup>**


Received: 13 March 2020; Accepted: 16 April 2020; Published: 24 April 2020

**Abstract:** The energy-optimal routing of Electric Vehicles (EVs) in the context of parcel delivery is more complicated than for conventional Internal Combustion Engine (ICE) vehicles, in which the total travel distance is the most critical metric. The total energy consumption of EV delivery strongly depends on the order of delivery because of transported parcel weight changing over time, which directly affects the battery efficiency. Therefore, it is not suitable to find an optimal routing solution with traditional routing algorithms such as the Traveling Salesman Problem (TSP), which use a static quantity (e.g., distance) as a metric. In this paper, we explore appropriate metrics considering the varying transported parcel total weight and achieve a solution for the least-energy delivery problem using EVs. We implement an electric truck simulator based on EV powertrain model and nonlinear battery model. We evaluate different metrics to assess their quality on small size instances for which the optimal solution can be computed exhaustively. A greedy algorithm using the empirically best metric (namely, distance × residual weight) provides significant reductions (up to 33%) with respect to a common-sense heaviest first package delivery route determined using a metric suggested by the battery properties. This algorithm also outperforms the state-of-the-art TSP heuristic algorithms, which consumes up to 12.46% more energy and 8.6 times more runtime. We also estimate how the proposed algorithms work well on real roads interconnecting cities located at different altitudes as a case study.

**Keywords:** Electric Truck Simulator; Electric Vehicle (EV); Vehicle Routing Problem (VRP); Traveling Salesman Problem (TSP); least-energy routing algorithm; energy efficiency; EV batteries; metric evaluation

#### **1. Introduction**

Electric Vehicles (EVs) are currently a tiny fraction of the car market, which is dominated by Internal Combustion Engine Vehicles (ICEVs); however, the growth of the EV market over last ten years is remarkable, and they are expected to progressively replace ICEVs in the next 20 years.

EVs have high energy efficiency and are sustainable transportation, since their electric motor has high dynamic performance and they are more environmentally-friendly than ICEVs on the market today; even when evaluating the emissions generated during electricity production for the charge of EVs, their overall Greenhouse Gas (GHG) emission is up to 58% that lower than the emissions of average mid-size passenger ICEVs [1]. Moreover, the impact on climate change by the production of electricity and operation of EVs is up to 30% less compared with ICEVs when considering the average generation of electricity in Europe. In recent years, the landscape of EVs is widening and expands to domains such as electric racing cars, electric buses, and electric trucks.

In particular, Tesla announced electric trucks will replace existing ICE trucks in the future [2]. The electric truck can accelerate more quickly than conventional diesel trucks because of the characteristic of the electric motor: high torque at low Rotations Per Minute (RPM) with high efficient. Furthermore, 98% of the kinetic energy can be recovered with electric energy during regenerative braking, which makes the electric truck more energy efficient. According to the announcement by Tesla, electric truck owners can save more than \$200,000 over a million miles compared to fuel cost of a conventional truck.

The optimal energy-efficient delivery route can further improve the energy efficient of electric trucks. The typical delivery scenario is defined as an electric truck loads all packages for customers at a depot, visits each customer to deliver their package, and then returns to the depot without payload. For a conventional ICEV, the "cost" of a path is strongly driven by the distance (even if weight also matters) and the problem nicely fits into the well-known Traveling Salesman Problem (TSP) using distance as a metric. However, when considering EVs, the solution is not as straightforward; although distance obviously matters, the total energy consumption also strongly depends on the order of delivery as the efficiency of the EV is affected by the total (vehicle + payload) weight. Since the impact of the package weight on the battery SOC of electric trucks is much higher than that of the ICE trucks due to the characteristics of battery, the energy-optimal delivery method for electric trucks should be newly considered. As a matter of fact, one key characteristic of a battery is that it is progressively less efficient in delivering its energy as its State Of Charge (SOC) decreases [3,4]. A fully charged lithium-ion battery in EVs has better performance to deliver a high power demand than when it is partially discharged [5]. As the power consumption of the electrical motor strongly depends on the total weight, apparently, if we deliver the heaviest package first, the overall vehicle weight is reduced the most after unloading this package and following such order would be optimal [5]. On the other hand, it is obvious that the delivery distance should be considered; if we deliver the heaviest package first and this package has a very long distance from the depot, the battery will be discharged by carrying the heaviest weight for a long time, which might lead to a non-optimal delivery route if following the rule of deliver heaviest package first.

One first difference with respect to a plain ICEVs delivery is therefore in terms of metrics: for EVs, some combination of weight and distance should be considered. However, the most significant difference (and complication) lies in the fact that the calculation of the optimal energy path cannot be done incrementally, as the energy cost of a path is "dynamic", i.e., it depends on the previous choices as a consequence of the dependence on the residual weight, which means the weight and distance are the variables during the delivery process.

In this paper, we propose an overall electric trucks delivery simulation framework implemented by SystemC and SystemC-AMS for the least-energy electric truck delivery routing problem. We first implement an electric truck simulator with a powertrain model and a non-linear battery model of the Tesla Semi [6] that adopts the methodology introduced in [3] to trace the SOC evolution during the package delivery. From the simulation results, we show that a conventional metric for TSP, total delivery distance, as with any other "static" metric, does not minimize total energy consumption for an EV delivery. Since only an exhaustive exploration of all path guarantees to find the optimal path, we evaluate different static metrics (functions of weight and distance) on small graph instances to assess their quality; then, using the best metric derived in this calibration phase, we show how a greedy algorithm using that metric can provide significant reductions (up to 33%) with respect to the common-sense heaviest first package delivery.

#### *Energies* **2020**, *13*, 2096

The work is an extension of our previous conference paper [7], the main contribution of this work is to devise a heuristic algorithm to determine the energy-optimal routing of an electric truck. This main results encompasses a number of elements that we summarize as follows:


The paper is organized as follows. Section 2 introduces the motivation for the least-energy electric truck delivery routing problem. A comparison between the traditional shortest distance route and heaviest first route presented in this section to illustrate the limitation of distance-based routing method for the electric truck delivery routing problem. Section 3 describes the EV powertrain model and battery model used in our simulation framework and the typical vehicle routing problem. Section 4 shows our analysis of the routing problem; we discuss new metric candidates considering both of delivery distance and package weight. Then, new heuristic methods are proposed including greedy algorithm and TSP methods. Section 5 firstly shows how to implement the powertrain model of an electric truck and the related battery model in our simulation; then, the model validation follows. After that, we compare the energy consumption for the electric truck delivery problem by the metrics and approximate algorithms. To evaluate our proposal in the real delivery environment, we performed a case study for the delivery in a province located in Italy. Finally, Section 6 gives the conclusion of our work.

#### **2. Motivation**

In order to clearly illustrate the motivation of this work, we built an example to indicate how it is not possible to derive an energy-optimal delivery policy simply using a single "static" metric. Figure 1a shows a simple three-destinational delivery task from a depot (D) with a rough mapping on the plane, and the distance matrix between any pair of destinations. We assumed symmetry between node pairs to guarantee the generality.

**Figure 1.** Example to illustrate the motivation of energy-optimal delivery route exploration cannot only depend on single static metric: a delivery task (**a**) and delivery routes for two different delivery weights (**b**,**c**).

To evaluate the energy cost of a delivery path, we use a time diagram that plots the evolution of the total transported weight over delivery distance. We use weight as a proxy of electric truck power consumption, since it is proportional to the weight of the vehicle plus the total payload. Notice that it is clearly a simplification and does not take into account non-linearities of a battery, but even such approximation can help reveal the point we are making.

Distance is used as a proxy of time; we assume a constant speed for the deliveries. Again, this is an approximation of the real setting, where speed can be extremely varying due to the driving habits and the road traffic. When using a real battery model in the loop, however, the real speed profile of the vehicle can be accounted for. Therefore, evolution of weight over distance is a proxy of power over time, and the area of one such curve is then an estimate of the energy consumed for that delivery route.

Figure 1b shows two such delivery routes for a case in which the weights are *W*<sup>1</sup> = 10, *W*<sup>2</sup> = 20, *W*<sup>3</sup> = 30 and vehicle weight is *Wv* = 40. The dotted red curve represents a route for which packages are delivered in heaviest-first order (*D* → 3 → 2 → 1 → *D*), whereas the solid blue line denotes a route with the minimum total distance (*D* → 1 → 3 → 2 → *D*). In this specific case, the "min. total distance" policy works best (smaller area under the curve) and, by exhaustive exploration of the 3!=6 combinations of deliveries, it can be shown to yield the best value of the metric. Notice that the (*D* → 3 → 2 → 1 → *D*) route is not just slightly sub-optimal, and it ranks fourth out of the six routes, *D* → 3 → 1 → 2 → *D* being the worst one.

Figure 1c shows two other waveforms for the same delivery task in which the weights change to *W*<sup>1</sup> = 30, *W*<sup>2</sup> = 20, *W*<sup>3</sup> = 10, that is, in which the heaviest package corresponds to the closest destination (node 1). In this case, the red dotted profile of the "heaviest-first" yields the best value of the metric, while the "min. total distance" yields a slightly worse value. Notice that the blue solid line corresponds to the same order of delivery (yet with a different cost) as in Figure 1b as the distance has not changed in the two examples.

Notice also that (due the symmetric distances) there are two paths (i.e., *D* → 1 → 3 → 2 → *D* and *D* → 2 → 3 → 1 → *D*) with the same distance but with different "energy" cost, the cost of the path *D* → 2 → 3 → 1 → *D* being larger than the other one shown in Figure 1c. Therefore, an algorithm that picks edge simply based on distance could even get farther from the optimal solution.

In this straightforward example, although there are several approximations during the power consumption estimation, it shows the main two critical points raised by our work. Firstly, no simple single static metric can solve ideally the problem of searching the energy-optimal delivery route. Secondly, due to the state-dependent characteristic of the cost function, only the brutal exhaustive exploration can find the optimal solution, but, since this is only feasible for very small instances, we need to find a provably good static metric that can be used in a heuristic algorithm. According to the above finding from the motivating example, this static metric should combine both *weight and distance* as twinned factors affecting the energy consumption of the EV.

#### **3. Background and Related Work**

#### *3.1. Powertrain Model in EV*

As is common, the vehicle powertrain model is from the vehicle dynamics. There are four forces acting on a vehicle driving on a road with *θ* road slope, as shown in Figure 2: rolling resistance *FR*, gradient resistance *FG*, inertia resistance *FI*, and aerodynamic resistance *FA*.

The powertrain model by vehicle dynamics (*Pdyna*) is described as

$$P\_{dyna} = T\omega = F\frac{d\mathbf{s}}{dt} = (F\_R + F\_G + F\_I + F\_A)\upsilon$$

$$F\_R \propto \mathbb{C}\_{rr} \mathcal{W}, F\_G \propto \mathbb{W}sin\theta, F\_I \propto ma, \text{and } F\_A \propto \frac{1}{2}\rho \mathbb{C}\_d A \upsilon^2 \tag{1}$$

$$P\_{dyna} \approx (a + \beta \sin \theta + \gamma a + \delta \upsilon^2) m\upsilon$$

where *Crr* is the rolling coefficient, *W* is vehicle weight, *θ* is road slope, *m* is vehicle mass, *v* is vehicle speed, *a* is vehicle acceleration, *Cd* is drag coefficient, and *A* is the vehicle facial area [8]. The coefficients

*α*, *β*, *γ*, and *δ* of the dynamic power *Pdyna* are coefficients of rolling resistance, gradient, inertia, and aerodynamic, respectively.

**Figure 2.** Forces acting on a truck.

In addition to the forces, there are several losses on a rotating motor: a copper loss is proportional to the square of motor torque and iron and friction loss are related to motor RPM. In addition, there is a constant loss while the vehicle is operating. The powertrain model of EV (*PEV*) includes the motor losses in addition to *Pdyna*:

$$P\_{EV} = P\_{dyna} + \mathcal{C}\_0 + \mathcal{C}\_1 \upsilon + \mathcal{C}\_2 \upsilon^2 + \mathcal{C}\_3 T^2. \tag{2}$$

where *C*0, *C*1, *C*2, and *C*<sup>3</sup> are the coefficients for constant loss, iron and friction loss, copper loss, and drivetrain loss, respectively.

EVs and hybrid vehicles mostly use regenerative braking during deceleration. The regenerative braking converts kinetic energy to electric energy from generation process. The harvested energy is related to the electromagnetic flux inside of the motor, which is proportional to the motor RPM. Therefore, the regenerative braking model can be simplified as

$$P\_{\text{regen}} = \epsilon \, T\upsilon + \zeta. \tag{3}$$

where  and *ζ* are regenerative braking coefficients to model the regenerative power as a function of the current velocity.

#### *3.2. Battery Model in EV*

The EV is normally powered by a battery pack that includes a large number of Lithium battery cells connected in parallel and series. The battery pack of Tesla Model 3 with long range version comprises 4416 2170-size lithium-ion cells of 4800 mAh nominal capacity with 46p96s arrangement, and 800 km range. Tesla Semi truck has a 750 kWh battery pack and weighs about 5.1 tons. To build the battery pack model directly is a non-trivial task, therefore, the battery pack can be built by composing the individual battery cell model, and the model must be able to accurately account for the varied load current and SOC variations of the usable battery capacity to capture the non-ideal discharge characteristics of battery. We select a circuit-equivalent battery model that can model the effects of load current magnitude and dynamics on real-time battery usable capacity [9], and then use this model to compose the whole battery pack model.

Figure 3 depicts the circuit-equivalent model of one battery cell adopted in this work. The left-hand part for modeling the battery lifetime (usable capacity) and the the right-hand part represents the transient battery voltage. Notice that the left-hand part also account for current

magnitude and load frequency dependency on actual battery capacity. The left-hand part comprises a capacitor C representing the nominal storage capacity of the battery in *Ah* and a current generator representing the battery current *Ibatt* requested by the load. As the available capacity of the battery is affected by the load current variations distribution, there are two voltage generators on the left part: one represents the dependency of the battery capacity on the current values, while the other generator models the dependency on load current frequency. Both decrease the voltage at node SOC, representing the SOC, for larger current magnitudes and frequencies. The right-part has a variable voltage generator affected by the SOC of battery, the internal resistance is also influenced by the current SOC of battery, and two pairs of RC express the instantaneous battery voltage.

The methodology to extract the relationship between SOC and internal resistance, capacitance, and open circuit voltage is presented in [10], and the implementation of adding two different dependencies of load current on the left-hand side of model is introduced in [9].

**Figure 3.** Adopted circuit-equivalent model for one single battery cell.

Based on this single cell model, we derive the battery pack model by ideally scaling all electrical parameters according to the series/parallel configuration in the pack, which leads to faster simulation runs and a higher flexibility in the modeling of large battery pack, so that not all cells of a large pack have to be modeled and simulated individually. Notice that the implementation of battery pack model is somehow ideal (e.g., cell mismatches are not considered), while this is still more accurate than a linear battery model that neglects state-dependent battery characteristic.

Given this model, we can track the energy consumed by the EV by applying to the model the drawn power (as current and voltage waveforms) corresponding to the electrical motor consumption on a given leg of the route. In the most general case, there is a non-ideal power conversion step between the electrical motor and the battery. In this case, it suffices to scale the motor current and voltage according to the converter efficiency *η* < 1, which can be any complex function of the motor parameters, i.e., *Pbatt* = *Pmotor* · *η*. We assume the converter efficiency is a fixed valued in this work as in [3].

#### *3.3. Work Flow of Proposed Methodology*

Figure 4 shows the conceptual flow of our methodology for the estimation of the operation range of EVs. Three descriptive datasets are required as inputs.

Vehicle data includes (1) motor information: motor efficiency by motor torque and RPM, operation range of motor torque and RPM, (2) vehicle information: weight, drive train and body shape, and (3) other electrical systems. Facial area of the vehicle affects aerodynamic resistance. Route information include road distance between cities, road slope, and traffic on each road. Also payload by a delivery task is given. Route information and vehicle data is used in the vehicle model, which generates instantaneous power demand (V(t) and I(t)) during delivery.

We implement a battery model from a given battery specification: nominal capacity, voltage-to-SOC curve, impedance, and structure of the battery pack. Then we combine the vehicle model and battery model together to conduct simulation, the power demand of given delivery tasks

over time derived from vehicle model is pass into battery model, finally the battery model computes the residual SOC of battery pack.

**Figure 4.** Overall concept of the proposed methodology.

#### *3.4. Vehicle Routing Problem*

The vehicle routing problem is formulated as a graph *G*(*V*, *E*, *C*) where *V* = {*v*0, ··· , *vN*} is the set of vertices including N destinations and a depot, *E* = {*eij*|*i*,*j*∈*V*} is the set of edges between two vertices *vi* and *vj*, and *C* = {*cij*|*i*,*j*∈*V*} is the cost related to each edge *eij*. Vertex *v*<sup>0</sup> is the depot, while the remaining vertices in *V* represent customers that need to be served. The TSP consists in finding a route based at the depot, such that each of the vertices is visited exactly once while minimizing the overall routing cost.

The formulation of the vehicle routing problem is generalized as the TSP. Because TSP is known as NP-hard, several approximation algorithms are proposed during last several decades. Christofides designed an approximation algorithm for TSP using the Minimum Spanning Tree (MST) algorithm, which obtains approximated results less than 1.5 times the optimal solution [11]. From the general TSP, there are several variants of TSP to consider various constraints and delivery requirements. There is a variant of TSP considering a set of potential customers living near secured customers [12]. The salesman finds the shortest path to cover all potential customers within a certain distance from the path. A fleet of delivery vehicles characterized by different capacities and costs is an important variant of TSP [13,14]. There is a set of customers and a set of different types of vehicles. Each vehicle has different capacities in terms of the number of customers and operation cost; the goal is to find a set of routes for each vehicle minimizing total delivery cost. Some studies consider the number of customers that each vehicle should be responsible for, but they do not consider the vehicle weight changing with each delivery.

Recently, the vehicle routing problem with pick-up and delivery considers the situation in which packages have to be picked-up from one of customers and delivered to another location [15,16]. During the pick-up and delivery process, visiting each pickup and delivery places occurs exactly once and total package weight during the delivery should not exceed the capacity. This problem considers the weight of each package; however, it does not consider the energy consumption that changes after unloading each package.

There is a paper minimizing energy in the vehicle routing problem [17]. This paper solves the vehicle routing problem using integer linear programming to minimize the product of distance and weight of each arc. However, there is no result validation using energy simulation. Therefore, there is no analysis of energy consumption in the points of view of package weight distribution and routing methods in the real road condition.

#### **4. Analysis of Routing Algorithms**

The examples in Figure 1 suggest that a reasonable metric to track the energy spent on a delivery path could be to use a quantity that is correlated to the product of *distance* × *weight*. More precisely, let *i* and *j* be two vertices of the delivery graph, *dij* the distance between them, and *Wi* and *Wj* the weights to be delivered, respectively, at *i* and *j*. This metric should be proportional to *dij* × (*Wcurrent* − *Wi*), assuming the edge is traveled in the direction *i* → *j*. *Wcurrent* is the current weight of the vehicle when reaching node *i* (Figure 5a); this metric, which we call *D*x*Wr* hereafter (i.e., distance times residual weight), uses the information contained in the solid oval circle.

**Figure 5.** Features of the Target Metric: Generic instance between two nodes *i* and *j* (**a**); and a greedy decision based on this metric (**b**).

By accruing this quantity over the delivery sequence, we would therefore be able to measure the area below the curves in the (weight, distance) space of Figure 1 as a proxy of the total energy spent.

Therefore, as the delivery problem is an instance of TSP, intuitively one could be tempted to run some TSP heuristic algorithm using the above cost in place of distance as in traditional TSP instances. Although approximate, (besides also neglecting battery non-idealities), this strategy will leverage well-consolidated heuristics for the solving the TSP and could be relatively efficient. Unfortunately, there is a subtlety in this argument. TSP formulations do assume the use of a "static" metric, i.e., *whose value does not depend on the currently built solution*, such as distance. As a matter of fact, state-of-the-art TSP heuristics rely on the calculation of the MST algorithm as a pre-processing. This is because *the cost of a MST is the simplest lower bound for the TSP*: it can be shown that the removal of one edge from any Hamiltonian cycle (i.e., a solution of the TSP) yields a spanning tree [11].

Algorithms to compute the MST systematically grow the tree by *greedily picking edges* in increasing order of the cost function, which clearly implies the need of a "static" metric, otherwise it would not be possible to guarantee global optimality by using local optimal choices (e.g., shortest edge) at each step. It is therefore immediately clear that the above cost function *dij* × (*Wcurrent* − *Wi*) cannot be used in a conventional TSP algorithm and in particular based on MST, for two reasons. First, MST runs on an undirected graph and there is no intuition about in what direction the edge is traversed. Secondly, and most importantly, MST algorithms take a "local" greedy decision independent of the specific previous decisions; this is at the basis of greedy algorithms, in which the global optimum is a sequence of locally optimal choices. Therefore, even replacing the traditional distance-based metric in the MST with *D*x*Wr*, it would simply not work, as shown in Figure 5b) (notice that for simplicity we assume the graph is directed to emphasize the direction of the delivery path).

Consider the decision to be taken at node *i*, at which we arrive with a given value of *Wcurrent*. Without loss of generality, let us assume that there are only two possible choices at *i*, i.e., nodes *j* and *k*. A MST algorithm based on *D*x*Wr*, since *Wcurrent* represents a "state" information of the current path at *i* and as such is a fixed value, will then greedily pick the edge with smallest value min(*dij*, *dik*) to decide about whether to grow the MST along *j* or *k*. In other words, as *Wcurrent* does not change, min((*Wcurrent* − *Wi*) ∗ *dij*,(*Wcurrent* − *Wi*) ∗ *dik*) and min(*dij*, *dik*) result in the same choice.

Therefore, as already shown by the example of Section 2, the energy-optimal (or more precisely, *D*x*Wr*-optimal) solution can only be obtained by enumeration of all the possible paths. A greedy, local metric used for a MST-based TSP algorithm would then result in strongly sub-optimal solutions even if we should adapt an MST algorithm to use *D*x*Wr*.

#### *Algorithms Used in Our Analysis*

Based on the previous discussion, an exhaustive exploration of all paths to collect the one with the smallest cumulative *D*x*Wr* is the only way to achieve the optimal solution. This algorithm has obviously factorial complexity (*O*(*n*!) for *n* vertices), which is clearly applicable only to small instances and can be used for evaluating the quality of different approximate algorithms.

Concerning TSP-based algorithms, the computational complexity of the TSP heuristic with the best approximation, i.e., Christofides' algorithm, is *O*(*n*3), assuming that the graph is fully connected [11]. Although polynomial, cubic complexity can already be significant as the number of instances grow in the order to a few tens. Therefore, given the approximations of a TSP-based solution (the intrinsic approximation of the algorithm plus that of the metric described above), it makes sense to devise a simpler and *greedy* algorithm that *builds up the cycle as a path, one edge at a time, starting from the depot*, possibly using different greedy metrics. The greedy heuristic would clearly be linear in the number of nodes. Should some of the greedy heuristics be roughly as approximate as the TSP heuristic, it would at least guarantee that it can handle larger problem instances.

This choice would allow one to use the *D*x*Wr* metric that more closely tracks the energy value; by forming a path, in fact, we can calculate the equivalent of *Wcurrent* for the path being built. Moreover, since we start from the depot node, edges have an implicit direction and *D*x*Wr* can be calculated correctly. The approximation lies obviously in the fact that the greedy solution is not optimal, and, unlike the TSP heuristic, the approximation cannot be bounded. Christofides' algorithm, for example, can be shown to yield a solution that is no more than 3/2 of the optimal cost.

In our analysis, we therefore compare three classes of algorithms to solve the optimal routing problem:


They are listed in Table 1 where algorithms belonging to the three above categories are separated by a double line in the table. Each algorithm is labeled with an abbreviation for ease of reference in the simulation results.


**Table 1.** Exhaustive exploration of paths: list of algorithms.

The *TSPDxW* heuristic of the last line requires a further explanation. Given that, as shown above, using *D*x*Wr* in a TSP based algorithm would be immaterial as it will coincide with a distance-based metric, we devised an alternative metric that mimics *D*x*Wr* but that is suitable for a TSP-based heuristic. This metric, which we call *D*x*W*, assumes: (i) a 50% chance of traveling the edge in each direction; and (ii) approximates *Wcurrent* with the total weight (vehicle + payload). This results in a quantity that does not depend on the currently built solution, and uses also the information about the destination node (as shown in the dashed oval of Figure 5a), i.e.,

$$d\_{i\bar{j}} \times (\mathcal{W}\_{total} - (\mathcal{W}\_{\bar{i}} + \mathcal{W}\_{\bar{j}})/2).$$

Notice that a plain "*Shortest first*" and "*Smallest DxWr*" would be the same algorithm. Therefore, to resolve this problem, for "*Smallest DxWr first*" algorithm, we used *dij* ∗ (*Wcurrent* − *Wj*), which considers both: (1) shortest distance first with *dij*; and (2) heaviest package first by abstracting weight of the destinations *Wj* from *Wcurrent*.

#### **5. Simulation Results**

#### *5.1. Simulation Setup*

#### 5.1.1. Powertrain Model

We implemented a powertrain model of a Tesla Semi truck from the vehicle specification based on the presentation by Elon Musk; this is currently the only source of information for the specs as Tesla is preparing to release the Semi in 2020 or later [6,18]. The powertrain consists of four Model 3 electric motors; each motor is a three-phase AC permanent magnet electric motor with maximum power of 192 kW from 4700 to 9000 RPM, and maximum torque is 410 Nm below 4500 RPM [19,20]. We estimated curb weight of Semi as the sum of typical weight of class 8 truck and battery pack weight [21].

We first implemented a vehicle model in ADVISOR (ADvanced VehIcle SimulatOR) [22] by using the above vehicle specification. Then, we extracted the coefficients of EV powertrain model with a number of ADVISOR simulations, as described in [23]. Table 2 summarizes the model coefficients of Tesla Semi. We compared the results computed between the derived power model and ADVISOR to validate the model we used. Figure 6 shows the difference between the estimation of power consumption by the ADVISOR vehicle simulator and the powertrain models we derived; the normalized root-mean-square error is 4.93%.

**Table 2.** Powertrain model coefficients for Tesla Semi truck.


**Figure 6.** Powertrain model validation results compared with ADVISOR.

#### 5.1.2. Battery Pack Model

There is no exact specification of battery pack in Tesla Semi truck until now; therefore, we assumed that each electric motor is connected to one battery pack of Model 3 in our experiments. Each battery pack is composed of four modules that are connected in series; each module consists of Panasonic NCR18650B 3400 mAh Lithium battery cells arranged in a 46p24s configuration [24].

Table 3 summarizes the physical electrical parameters of each cell, each module, and the whole battery pack.



We built our battery single cell model based on the measurement data by adopting the method described in [9]. We assumed such 7104 battery cells in the pack to be ideally balanced in the following experiments, and then built battery pack model, as indicated in Section 3.2. Concerning the regenerative braking phase, we assumed that regenerative charging efficiency is 20% in our simulation, i.e., 20% of the kinetic energy is converted to electric energy and transferred into the battery pack.

#### 5.1.3. Simulation Framework

We adopted the simulation framework proposed in [3], which targets the modeling and simulation of energy flow in the EV. The simulation framework is built by SystemC and SystemC-AMS, which are the extension of C/C++ with libraries to describe HW constructs and analog/mixed-signal subsystems. SystemC-AMS provides different abstraction levels to cover a wide variety of domains, three different Model of Computations (MoCs) supported by SystemC-AMS that allow the simulation framework to integrate circuit equivalent battery model and empirical powertrain model simultaneously. Another main advantage of the SystemC-AMS simulation framework is the fast simulation speed while keeping the same accuracy with regard to state-of-the-art tools such as Matlab/Simulink, with speedups up to two orders of magnitude and a high level of accuracy. Such quick estimation of energy consumption of EV and the battery lifetime give the opportunity to conduct exhaustive exploration in a short time for different delivery routes in our following experiments.

#### *5.2. Simulation Results*

#### 5.2.1. Comparison Against Exact Results

In this section, we compare the energy consumption of various delivery strategies based on different policies for a set of small-sized (4–7) instances for which an exhaustive exploration of all the possible delivery paths is feasible.

For each number of destinations, we randomly generated 50 instances with different distributions of locations of the depot and of the destinations by uniformly distributing them in a 30 km × 30 km area. We selected the area for the delivery so that all the delivery sequences can be completed without exhausting the battery energy before returning to a depot. For each of the 50 instances, package weights for each destination were chosen as uniformly distributed from 0.1 to 3 ton.

For each problem instance (destination and weight distribution), we calculated by exhaustive exploration the route yielding the smallest value of energy for a number of metrics. Energy was calculated by conducting simulation described in Section 5.1. In this test example, a constant speed of 76 km/h was assumed, which was the average truck speed in metropolitan area interstates in the US in 2015.

Figure 7a shows the energy consumption of the optimal route, averaged over the 50 instances, for problems with 4–7 destinations and for the set of algorithms described in Table 1.

The leftmost blue bars represent the optimal routes yielding the minimum energy consumption among all routes, obtained by computing the actual energy consumption per each segment using the battery model. This is the reference value against which the other results were compared. All the other bars refer to solutions (i.e., routes) returned by the above algorithms and evaluated using the battery model. The objective of the simulation was to check how the greedy algorithms (TSP-based or not) differ from the optimal solution and how the error increases with increased problem sizes. Bars in the plot are in the same order as in Table 1. For ease of reading, bars referring to path-based algorithms (first part of the table) are shown as solid bars, whereas those referring to greedy or TSP algorithms are shown as patterned bars.

Concerning path-based algorithms, we can notice how the *MinDxWr* metric (third bar from left) tracks very well the true energy value, much better than distance alone (second bar from left). Concerning approximation algorithms, as a first general comment, we can see that all algorithms overestimate the actual energy consumption. Then, we immediately observe that weight alone (Heaviest-first) tracks quite poorly the actual consumed energy, somehow contradicting the intuition suggested by the battery property; the distance from the reference is already > 20% even for the four-destination instance. Although the actual error may differ depending on the weight distribution (as shown in Figure 1), the results are the average 50 different runs, thus we can safely assume this is not a good metric. Notice also that the tracking error increases with larger instances.

Another observation is that a traditional TSP with distance metric (second bar from the right) performs reasonably only for the smallest instance; that average error increases quickly and is already around 18% for seven destinations. Therefore, we can also rule out this algorithm from the list.

The remaining ones (*Smallest DxWr first, Shortest first*, and *TSPDxW*) have errors below 10%, with the greedy algorithms being below 5% and scaling better with problem size than *TSPDxW*. Especially, *Smallest DxWr first* algorithm saves energy consumption from 6.58% to 12.46% compared with traditional *TSPD* algorithm. The gap between the best one and *TSPD* increases by the number of destinations.

Figure 7b shows the worst case error among the 50 instances for the same set of algorithms. Results are consistent with average error, with the maximum error being significantly larger than the average one. The greedy algorithms have show again the best results, in terms of both error and scalability. The *Smallest DxWr* is the only algorithm with worst-case error around 20% (as opposed to about 30–35% of the others) for the seven-node case.

(**a**) Average energy consumption comparison.

(**b**) Worst energy consumption comparison.

#### 5.2.2. Comparison by Package Weight Distribution

In this section, we discuss how the error increases with the increased package weight. Figure 8 shows the energy consumption of different metrics by package weight distribution. Each column means different weight distribution of delivery package: 0.1–0.5, 0.5–1.0, 1.0–1.5, 1.5–2.0, and 2.5–3.0 ton, respectively. In this comparison, one depot and seven destinations were uniformly distributed in a 30 km × 30 km area, and we extracted average energy consumption of 50 instances by different metrics. In order not to provide unnecessary information, we do not consider *MinDxWr* and *Heaviest first* metrics here because these metrics show either sufficiently accurate or irrelevant results from the previous section.

**Figure 8.** Energy consumption by package weight distribution.

When the weight of the package is less than 0.5 ton, the sum of all package weight can be ignored in most instances compared with total weight of the electric truck. Thus, all results by different metrics show less than 10% error. As the range of the package weight increases, however, the effect of the package weight on the overall energy consumption increases. Therefore, metrics considering distance only show worse results, namely the exhaustive exploration with distance (*minD*) and the Traditional heuristic TSP (*TSPD*). On the other hand, heuristic TSP with *D*x*W* (*TSPDxW*) and two other greedy heuristics show less than 10% error for all weight ranges.

#### 5.2.3. Application to Larger-Scale Instances

We generated a number of instances with 10, 20, 30, 50, and 100 destinations; for each problem size, we generated 20 random instances and collected the average value of energy and execution time. In all cases, weights were scaled so that the delivery task could be completed.

Figure 9 compares the absolute energy values for the three competitive algorithms resulting from the previous section: one TSP with the proposed metric (*TSPDxW*) and the the two greedy heuristics (*Shortest first* and *Smallest DxWr first*). From the results in Figure 7a we know that all approximations are overestimations, thus we can assume that lower values of energy imply higher accuracy. In Figure 9, the *Smallest DxWr first* shows the best results: its energy consumption is 10% smaller than the TSP heuristic and 4% smaller than the shortest first one, for the larger 100-node instance.

**Figure 9.** Greedy and approximate TSP algorithms on large problem instances.

Figure 10 shows the slowdown of the TSP heuristic with respect to the *Smallest DxWr first* algorithm. The TSP execution time is obviously independent of the metric used (*D* vs. *D*x*W*). The TSP heuristic is significantly slower than the greedy method; the slowdown increases for increasing problem sizes, reaching 8.6 times for the 100 destination case.

**Figure 10.** Slowdown of TSP heuristics vs. greedy algorithm.

#### *5.3. Case Study: Routing Problem in Real Roads*

In the experiments presented in Section 5.2, the distribution of the locations were synthetically generated on a plane. In this section, we show the routing algorithms in a real case, consisting of a set of locations taken from a map and for which actual distances, road slope, and road traffic between destinations are taken into account. We generated 50 instances, in which package weights for each destination were chosen as uniformly distributed from 0.1 to 3 ton.

#### 5.3.1. Extraction of Road and Elevation Information

We used Piedmont region in Italy as the delivery destinations. There are 10 destinations in Figure 11 including a depot. The destinations of the delivery are limited to towns or cities in the province.

**Figure 11.** Destinations for deliveries on real roads.

To travel between destinations, we cannot drive on a straight path, but we must use given roads. There are several route options connecting destinations to each other. Among them, the most recommended one by Google Maps is picked. We extracted the distance of each route and related average driving time.

Table 4 shows the altitude information of each city. We extracted average road slope from road distance and altitude difference between two cities. The average altitude difference among 10 cities is 84 m.

The practical distance of the route, road slope, and driving time are different by direction. Therefore, we implemented several matrices containing road distance, road slope, and driving time. The average distance of the roads is 44 km, and the average time for driving the roads is 51 min. The driving time was used to calculate average electric truck velocity. Therefore, the energy consumption for the delivery was obtained from the practical road distance, velocity, and road slope.


#### 5.3.2. Simulation Results

Figure 12a shows the average energy consumption of the optimal route in the case study described in Figure 11 on the 50 instances with routing algorithms described in Table 1. When the number of destinations was four, we randomly picked one depot and four destinations among the 10 cities in Figure 11.

(**b**) Worst energy consumption comparison on real roads.

**Figure 12.** Energy consumption for different metrics (exhaustive exploration) on real roads.

Similar to Figure 7a, the leftmost blue bars (MinE) represent energy consumption by the energy-optimal routes among all routes. This is the reference value to compare with the other results. Bars in the plot are in the same order as in Table 1.

As confirmed by the results in Figure 7a, the path-based algorithms using the metric *MinDxWr* (third bar from the left) track very well the true energy value, much better than distance alone (second bar from left), in all simulation results with different number of destinations. Concerning approximation algorithms, weight alone (Heaviest-first) tracks quite poorly the actual consumed energy as in Section 5.2. The average error is up to 47.9%.

On the other hand, shortest first and *Smallest DxWr first* metrics track the actual consumed energy very well. The average error is less than 10% for all numbers of destinations. Especially, *Smallest DxWr first* shows the best results for all different number of destinations. A traditional TSP with a distance metric (second bar from the right) tracks well only for the smallest instance as in Figure 7a; average error increases up to 15.0%. TSP with a metric *D*x*W* shows better results than distance metric.

Figure 12b shows the worst case error among the 50 instances for the same set of algorithms in Figure 12a. The results are consistent with average error, but significantly larger. The best algorithms in Figure 12a (*Smallest DxWr first*, *Shortest first*, and *TSPDxW*) again show the best results, in terms of both error and scalability. *TSPD* and *TSPDxW* methods are not better in some cases (e.g., six destinations) because triangular inequality is not applicable, and the distance between cities is asymmetric.

#### **6. Conclusions**

The total energy consumption for parcel delivery with an electric truck strongly depends on the order of delivery because battery efficiency is affected by how the transported weight changes over time. However, it is impossible to consider the transported weight changes with the traditional routing algorithms, which use "static" metrics such as distance. In this paper, we demonstrate that the functions of weight and distance as metrics provide significant energy reductions with respect to the traditional routing algorithms. A traditional TSP with distance metric performs reasonably only for the smallest instance; average error increases quickly. On the other hand, a greedy algorithm

minimizing driving distance by residual weight (*D*x*Wr*) shows better results with almost 10× faster runtime. We also performed a comparison in the real delivery case, where curved and sloped roads connect cities. In this case study, the greedy algorithm also shows better results than traditional TSP. In addition, the package weight affects the result of routing algorithm. As we increase the package weight distribution, the gain by the proposed method increases.

**Author Contributions:** Conceptualization, N.C. and M.P.; Writing—original draft, D.B., Y.C. and M.P.; and Writing—review and editing, N.C. and E.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. NRF- 2018R1A2B3007894).

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

#### **References**


c 2020 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Torque Coordination Control of an Electro-Hydraulic Composite Brake System During Mode Switching Based on Braking Intention**

#### **Yang Yang 1,2,\*, Yundong He 2, Zhong Yang 3, Chunyun Fu 1,2 and Zhipeng Cong <sup>2</sup>**


Received: 3 March 2020; Accepted: 17 April 2020; Published: 19 April 2020

**Abstract:** The electro-hydraulic composite braking system of a pure electric vehicle can select different braking modes according to braking conditions. However, the differences in dynamic response characteristics between the motor braking system (MBS) and hydraulic braking system (HBS) cause total braking torque to fluctuate significantly during mode switching, resulting in jerking of the vehicle and affecting ride comfort. In this paper, torque coordination control during mode switching is studied for a four-wheel-drive pure electric vehicle with a dual motor. After the dynamic analysis of braking, a braking force distribution control strategy is developed based on the I-curve, and the boundary conditions of mode switching are determined. A novel combined pressure control algorithm, which contains a PID (proportional-integral-derivative) and fuzzy controller, is used to control the brake pressure of each wheel cylinder, to realize precise control of the hydraulic brake torque. Then, a novel torque coordination control strategy is proposed based on brake pedal stroke and its change rate, to modify the target hydraulic braking torque and reflect the driver's braking intention. Meanwhile, motor braking torque is used to compensate for the insufficient braking torque caused by HBS, so as to realize a smooth transition between the braking modes. Simulation results show that the proposed coordination control strategy can effectively reduce torque fluctuation and vehicle jerk during mode switching.

**Keywords:** electric vehicles; electro-hydraulic braking; braking intention; mode switching; torque coordinated control

#### **1. Introduction**

The electro-hydraulic composite braking system of an electric vehicle (EV) consists of the motor braking system (MBS) and the hydraulic braking system (HBS), which realize the pure electric, pure hydraulic, and hybrid braking modes. The composite braking system converts the kinetic energy of the vehicle into electric energy and ensures braking stability and braking efficiency during braking [1–4]. The braking modes of the electro-hydraulic composite system switch between each other as braking conditions vary. However, the MBS and HBS dynamic response characteristics are not consistent, which leads to total braking torque fluctuations during mode switching, thus affecting braking safety and ride comfort. Therefore, it is of great significance to study the braking torque coordination control during braking mode switching.

Current research on the electro-hydraulic composite braking system mainly focuses on the distribution of braking forces and recovery of braking energy. For example, for the problem of braking force distribution in different modes, Sun et al. [5] established the optimal distribution coefficient response surface by optimizing the distribution coefficient of hydraulic braking torque and regenerative braking torque offline, which improved braking stability and energy recovery efficiency during braking. Shi et al. [6] designed a regenerative braking system that can achieve braking energy recovery during emergency braking. Considering the tire, hydraulic, and motor losses, Sun et al. [7] proposed an on-line control strategy for electro-hydraulic composite braking force, which improved the regenerative braking power. For the problem of braking torque coordination control during mode switching, Okano et al. [8] adopted a filtering algorithm to assign the MBS response to high frequency braking torque and the HBS response to low frequency braking torque, making full use of the dynamic characteristics of both systems. He et al. [9] designed a combined controller for torque disturbance in mode switching, incorporating linear quadratic optimal and sliding mode controllers—the former controller is used for anti-interference, and the latter is used to compensate the performance index offset of the nonlinear part—and achieved a good match of the target speed during mode switching. Considering the influence of the half axle elasticity and backlash nonlinearity of the transmission system on the control performance and dynamic characteristics of the MBS, Lv et al. [10,11] proposed an active control algorithm based on a hierarchical structure to realize the clearance compensation for the transmission system, and Zhang et al. [12] proposed a method of backlash sliding mode compensation and an elastic double closed-loop PID compensation for the control of a permanent magnet synchronous motor; the approaches of both research groups effectively compensated for the influence of the transmission system on the control performance of a permanent magnet synchronous motor. According to whether the HBS provides braking torque, Yang et al. [13] proposed to reduce the torque fluctuation by controlling the change rate of the clutch engagement torque and motor braking torque, and by modifying the target braking torque at different stages during mode switching. Yu et al. [14,15] proposed a double closed-loop feedback control and motor braking torque modifying method, based on the differences in the characteristics of the MBS and the HBS, using the motor braking torque to modify the hydraulic braking torque, thus reducing vehicle jerk during mode switching. Although the aforementioned research has improved braking torque coordination control during mode switching, leading to better control of the electro-hydraulic composite braking system and reduced vehicle jerking, they did not reflect the driver's braking intention during mode switching; that is, the motor and hydraulic braking torque cannot be adjusted reasonably according to whether the driver pays more attention to brake safety or ride comfort.

In this paper, the problem of the total braking torque fluctuation and the jerk of the complete vehicle is addressed for the electro-hydraulic braking system of a four-wheel-drive pure electric vehicle with a dual motor. Firstly, the dynamic characteristics of the vehicle during braking are analyzed. Braking force distribution control strategies that take into account braking safety and regenerative braking energy recovery are established based on the vehicle state parameter constraints. Then, a combined control method, which contains a PID controller and a fuzzy controller, is used to control the brake pressure in each wheel cylinder. Finally, the fuzzy control rules based on the brake pedal stroke and its change rate are designed to modify the target hydraulic braking torque and to reflect the driver's braking intention; that is, the target hydraulic braking torque is modified according to whether the driver pays more attention to brake safety or ride comfort. At the same time, the rapid response of the motor braking torque is used to compensate for the insufficient braking torque caused by the slow response of the HBS, so as to realize the smooth transition of the braking mode, which enhances braking safety and ride comfort during mode switching.

#### **2. Vehicle Dynamics Model and Braking Force Distribution Control Strategies**

#### *2.1. The Electric Vehicle Structure*

The pure electric vehicle analyzed in this paper is a front–rear centralized dual-motor driving system with two three-phase permanent magnet synchronous motors (PMSM), as shown in Figure 1. The driving forces are transmitted to the wheels from the motors through final drives I and II. The vehicle-state variables, such as vehicle speed, braking pedal displacement, and battery state of charge (SOC), are collected by an electronic control system and transmitted to the vehicle control unit (VCU) via a controller area network (CAN) bus, and the braking torques required of the MBS and the HBS are determined by the vehicle control unit. The motor braking torque and hydraulic braking torque are controlled by the motor control unit (MCU) and hydraulic control unit (HCU), respectively. The basic parameters of the vehicle are shown in Table 1.

**Figure 1.** Structure of a four-wheel-drive pure electric vehicle.

**Table 1.** Parameters of the vehicle.


#### *2.2. Dynamics Analysis of Braking*

When a vehicle is braking, and the air resistance moment, rolling resistance moment, and moment of inertial generated by the rotating mass are ignored, the normal acting force of the ground on the front wheel [7] is

$$F\_{Z1} = \frac{G\left(L\_r + zH\_\mathcal{S}\right)}{L} \tag{1}$$

The normal force of the ground acting on the rear wheel is

$$F\_{Z2} = \frac{G\left(L\_f - zH\_\mathcal{S}\right)}{L} \tag{2}$$

where *FZ*<sup>1</sup> and *FZ*<sup>2</sup> are the normal acting force of the ground on the front and rear wheels; *L* is wheelbase; *Hg* denote the height of the center of mass of vehicle; *Lf* and *Lr* are the distance from the center of mass to front and rear axles; *G* is the weight of vehicle; and *z* represents the braking strength, and it is the ratio of the vehicle deceleration to the gravitational acceleration.

The dynamic equation of wheel rotation during braking is

$$J\_t \dot{\phi}\_{\mathcal{W}} = F\_{Xb} r - (T\_{\mathcal{W}} + T\_{\mathcal{h}}) \tag{3}$$

where *Jt* represent the moment of inertia of transmission equivalent to the wheel; *FXb* denote brake force of ground; *r* is the tire radius; *Tm* is the braking torque acting on the wheel by PMSM I or PMSM II; *Th* is hydraulic braking torque; and ω*<sup>w</sup>* is angular velocity of the wheel.

Dynamics equation of the vehicle transmission system during braking is

$$k\_1 T\_{m1} i\_{O1} + k\_2 T\_{m2} i\_{O2} + k\_3 T\_h = \frac{G}{g} \frac{dv}{dt} r \tag{4}$$

where *v* denote vehicle speed; *g* is acceleration of gravity; according to the working state of PMSM I, PMSM II, and HBS, *k*1, *k*2, and *k*<sup>3</sup> is 0 or 1; and *iO*<sup>1</sup> and *iO*<sup>2</sup> are the gear ratio of the final drive I.

#### *2.3. The Braking Force Distribution Control Strategies*

The distribution of braking force is determined by the braking state of the vehicle and must meet the requirements of the brake regulations. Due to deceleration during braking of the vehicle, the stable braking strength provided by the motors varies between 0.086 and 0.2. According to the theory of braking stability, the braking force distribution curve under the I-Curve can ensure the stability of the vehicle; that is, to maintain the ability of straight-line driving of the vehicle during braking [16–18]. Therefore, the maximum braking strength provided by the motors in the pure electric braking mode is determined to be 0.17. The dynamic distribution control strategies based on the I-Curve are shown in Figure 2.

**Figure 2.** Braking force distribution curve of the front and rear axles.

The braking strength at points *A*, *B*, *C*, *D*, and *E* in Figure 2 are *z*(*A*), *z*(*B*), *z*(*C*), *z*(*D*), and *z*(*E*). I-Curve is the ideal braking force distribution over the front and rear axles. *Fb f* and *Fbr* denote the braking force of the front and rear axles. *Fm f* \_max and *Fmr*\_max denote the maximum braking forces provided by the PMSM I and the PMSM II. The fixed braking force distribution coefficients of the front and rear axles are β<sup>1</sup> and β2. μ is the road adhesion coefficient.

The minimum speed of the vehicle at which the motors maintain a stable braking torque is *v*min; the maximum speed of the vehicle at which the motors can perform regenerative braking is *v*max; and the maximum state of charge in which the battery can be charged is *SOCh*. During braking, the speed of the vehicle and the *SOC* of the battery should obey the restriction that *v*min ≤ *v* ≤ *v*max and *SOC* ≤ *SOCh*, respectively. The motor and hydraulic braking force distribution based on the braking strength required by driver is as follows:

When *SOC* < *SOCh* and *v*min ≤ *v* ≤ *v*max:

(1) 0 < *z* ≤ *z*(*A*), in order to guarantee the braking stability of the vehicle during braking, the PMSM I is given priority to provide braking force, the braking force is provided by PMSM I individually in this condition.

(2) *z*(*A*) < *z* ≤ *z*(*B*), the braking force is provided by PMSM I and PMSM II simultaneously, and PMSM I provides the maximum braking force, the remaining force is provided by PMSM II.

(3) *z*(*B*) < *z* ≤ *z*(*C*), the HBS starts to provide braking force in this condition, and the braking force of the front and rear axles are distributed according to the fixed braking force distribution coefficient β1. PMSM I maintains the maximum braking force and the insufficient braking force of the front axle is provided by the HBS. The braking force of PMSM II continues to increase with braking strength, until the maximum braking force of PMSM II is reached.

(4) *z*(*C*) < *z* ≤ *z*(*E*), the required braking force is provided by the motor braking system and hydraulic braking system simultaneously in this condition. Both PMSM I and PMSM II are working at the maximum braking force that can be supplied, and the insufficient braking force of the front and rear axles are provided by the hydraulic braking system.

When *SOC* ≥ *SOCh* or *v* < *v*min or *v* > *v*max, or *z* ≥ *z*(*E*), the braking force is provided only by the HBS, and the braking force at the front and rear axles are distributed according to the fixed braking force distribution coefficient β<sup>1</sup> and β2.

#### **3. Modeling and Characteristics Analysis of Braking Systems**

In order to obtain the dynamic response characteristics of the MBS and the HBS, a simulation model was established in MATLAB/Simulink (2016b, MathWorks, MA, USA) based on the mathematical models of the motor and the hydraulic brake system. In order to accurately control the hydraulic braking torque, a combined control method, which contains a PID controller and a fuzzy controller, was designed to control the brake pressure in each wheel cylinder.

#### *3.1. The Modeling of PMSM*

The main parameters of the PMSM used in this paper are shown in Table 2. The PMSM is a strong complex-coupling, high-order, and multivariable nonlinear system [19,20]. In order to realize the vector control of the motor, the mathematical model of the PMSM in the two-phase rotating reference frame (*d*-*q* axis) was used to establish the simulation model. The three-phase PMSM in the *d*-*q* axis can be described as [21]

$$
\begin{pmatrix} u\_d \\ \
u\_q \end{pmatrix} = \begin{pmatrix} R\_\delta + L\_d \frac{d}{dt} & -\omega \nu\_m L\_q \\ \omega \nu\_m L\_d & R\_\delta + L\_q \frac{d}{dt} \end{pmatrix} \begin{pmatrix} i\_d \\ i\_q \end{pmatrix} + \begin{pmatrix} 0 \\ \omega \nu\_m \psi\_f \end{pmatrix} \tag{5}
$$

The electromagnetic torque equation of the PMSM in *d*-*q* axis is

$$T\_m = 1.5 p\_n \left[ \psi\_f i\_q + \left( L\_d - L\_q \right) i\_d i\_q \right] \tag{6}$$

where *ud* and *uq* are the armature voltage components in the *d*-*q* axis, respectively; *id* and *iq* denote the armature current in the *d*-*q* axis, respectively; *Ld* and *Lq* represent the equivalent armature inductance in the *d*-*q* axis respectively; ψ*<sup>f</sup>* is the rotor flux corresponding to the permanent magnet; *Rs* represent the stator resistance; ω*<sup>m</sup>* denote the rotational angular velocity of *d*-*q* axis; and *pn* represents the pole pairs of motor.

**Table 2.** Key parameters of the motor.


#### *3.2. The Modeling of Hydraulic Components*

The hydraulic braking system mainly includes system components, such as the brake master cylinder, the wheel cylinder, and the brake pedal simulator, as well as the control components, such as the high-speed on–off valve [22,23]. The structure of the hydraulic braking system used in this article is shown in Figure 3. The pressure of the wheel cylinder is adjusted by the high-speed on–off valve to meet the requirement of the hydraulic braking torque, while the hydraulic braking torque is directly determined by the pressure in the wheel cylinder; i.e., these components reflect the dynamic characteristics of the HBS. Therefore, the mathematical models of the high-speed on–off valves and the brake wheel cylinders are mainly described.

**Figure 3.** The hydraulic braking system (HBS) structure diagram.

3.2.1. The Modeling of the High-Speed On–Off Valve

The pressure of a brake wheel cylinder is controlled by a pair of high-speed on–off valves: the inlet valve that is normally opened and the outlet valve that is normally closed. The structure of the outlet valve and the force analysis of the valve core are shown in the Figure 4. The key parameters of the high-speed on–off valves selected in this paper are summarized in Table 3.

**Figure 4.** The structure diagram of the outlet valve and the force analysis of the core.

According to the control signals of the vehicle controller, the pressure in the wheel cylinder is adjusted by the hydraulic controller through the combined control of the inlet and outlet valves. Through the force analysis of the valve core in Figure 4, the kinetic equation of the valve core can be obtained [24]:

$$\begin{cases} \frac{d\mathbf{v}}{dt} = \frac{1}{m\_f} \left[ F\_m(\mathbf{x}, \mathbf{i}) - K(\mathbf{x} + \mathbf{G}\_0) - F\_P(\mathbf{x}) - b\mathbf{v}\_f - F\_f - F\_j \right] \\\\ \mathbf{v}\_f = \frac{d\mathbf{x}}{dt} \end{cases} \tag{7}$$

where *Fm* is the electro-magnetic force; *Ff* denotes the frictional force, it takes an estimated value of 0.01 N in this study; *Fj* represents the jet force; *FP* is the flow force of the core assembly; *mf* is the mass of the core assembly; *K* is the stiffness of the return spring; *b* and *x* are the velocity viscosity coefficient and the displacement of the core, respectively; and *G*<sup>0</sup> represents the return spring's pre-compression force.

**Table 3.** Main parameters of the high-speed switch valve.


#### 3.2.2. The Modeling of the Wheel Cylinder

During hydraulic braking, the dynamic characteristics of the wheel cylinder piston can be expressed as a spring–mass–damper system, and its dynamic equation [25] is

$$p\_{\overline{w}}A\_p - F\_{k0} = m\_p \ddot{\mathbf{x}}\_p + C\_p \dot{\mathbf{x}}\_p + k\_p \mathbf{x}\_p \tag{8}$$

where *Fk*<sup>0</sup> represents pre-tightening force; *pw* denotes the wheel cylinder pressure; *Ap* is the effective action area of the piston; *mp* denotes the mass of piston; *xp* represents the displacement of piston; *Cp* is the damping of the brake; and *kp* is the equivalent stiffness.

The relationship between the pressure in the wheel cylinder and the hydraulic braking torque in front and rear axles can be expressed as [13]

$$\begin{cases} \quad T\_{hf} = 2p\_w \frac{\pi D\_f^2}{4} R\_f K\_f \\ \quad T\_{hr} = 2p\_w \frac{\pi D\_r^2}{4} R\_r K\_r \end{cases} \tag{9}$$

where *Th f* and *Thr* represent the hydraulic braking torque of the front and rear axles, respectively; *Kf* and *Kr* are the brake factors of the front and rear axles; *Df* and *Dr* denote the diameter of the front and rear wheel cylinders; and *Rf* and *Rr* represent the effective radius of the front and rear wheel brake discs.

The key parameters of the front and rear brakes is shown in Table 4.

**Table 4.** Key parameters of front and rear wheel brakes.


#### 3.2.3. Design of a Combined Controller for Hydraulic Braking Torque

The HBS should have good control performance to achieve a fast and accurate response to the hydraulic braking force. A cooperative control strategy is used to improve the response speed and control the accuracy of the hydraulic braking torque. When the pressure error Δ*p*(*t*) between the target pressure and the tracking pressure is larger than the value of the threshold, the PID controller controls the wheel cylinder pressure to achieve a rapid adjustment. Conversely, when the pressure error Δ*p*(*t*) is smaller than the value of the threshold, the fuzzy controller is used to stabilize the wheel cylinder pressure near the target pressure and reduce the hydraulic pressure fluctuation [26,27]. The control schematic is shown in Figure 5. Wheel cylinder control is determined by the logical threshold controller according to the pressure error Δ*p*(*t*) , and then either the PID or fuzzy controller determines the duty ratio of the high-speed on–off valve. The inlet and outlet valves of the wheel cylinder are directly controlled by PWM (pulse width modulation). Thus, the tracking pressure of the wheel cylinder follows the target pressure. A threshold value of 0.5 bar was set for Δ*p*(*t*) as simulation results indicated that both the response speed and the control accuracy of the pressure in the wheel cylinder were optimal at this level.

**Figure 5.** The pressure control algorithm of the wheel cylinder.

When the pressure error Δ*p*(*t*) is larger than the value of the threshold, the PID controller is used to adjust the pressure in the wheel cylinder, so that the tracking pressure in the wheel cylinder rapidly responds to the target pressure. The mathematical model of the PID controller [24] is

$$D(t) = k\_p \Delta p(t) + k\_i \int\_{t\_0}^{t\_1} \Delta p(t)dt + k\_d \frac{d\Delta p(t)}{dt} \tag{10}$$

where *D*(*t*)is the duty ratio of the PWM signal; *kp*, *ki*, and *kd* are the proportional, integral, and differential coefficients of the PID controller, respectively; *t*<sup>0</sup> and *t*<sup>1</sup> are the time when the HBS starts and ends to work.

When the pressure error Δ*p*(*t*) is smaller than the value of the threshold, the fuzzy controller is used to adjust the pressure of the wheel cylinder. The opening degree of the inlet and outlet valves are determined by the target pressure *p*(*t*) and the pressure error Δ*p*(*t*). When the target pressure *p*(*t*) is small (S') and the pressure error Δ*p*(*t*) is negative (N), the tracking pressure needs to be reduced, so that the duty ratio of the inlet valve *Din* takes a smaller value (S') and the duty ratio of the outlet valve *Dout* takes a larger value (VL'). When the target pressure *p*(*t*) is large (VL') and the pressure error Δ*p*(*t*) is positive (P), the tracking pressure needs to be increased, so that the duty ratio of the inlet valve *Din* takes a larger value (VL') and the duty ratio of the outlet valve *Dout* takes a smaller value (S'). The fuzzy controller rules for the degree of opening of the high-speed on–off valve are shown in Tables 5 and 6.

 In Table 5, N, Z, and P denote less than zero, equal to zero, and greater than zero, respectively; whereas S', MS', M', ML', L', and VL' represent small, small medium, medium, medium large, large, and very large, respectively.

The pressure variation of a wheel cylinder controlled by the combined controller under the input of a sinusoidal signal is shown in Figure 6. The tracking pressure in the wheel cylinder is precisely controlled, and the pressure error between the target pressure and tracking pressure is very small.


**Table 5.** The fuzzy control rules of the inlet valve.


**Table 6.** The fuzzy control rules of the outlet valve.

**Figure 6.** The response of the wheel cylinder with a sinusoidal input.

#### *3.3. Dynamic Characteristics Analysis of the MBS and the HBS*

Using the mathematical models described above, a dynamic model of the electro-hydraulic composite braking system was built in MATLAB/Simulink, and a separate HBS physical model was built in the Simulink sub-module, Simscape. The dynamic responses of the MBS and HBS are shown in Figure 7. Under the input of the same demand braking torque, the response time of the MBS is *tm*, the response time of the HBS is *th*, and the response time difference is Δ*tmh*. Compared with the HBS, the dynamic response of the MBS is fast, and the braking torque rise time is short, but there is a certain amount of overshoot. There are two main reasons for the differences in dynamic characteristics: In the initial stage of the hydraulic braking torque response, the HBS needs high-pressure brake fluid to fill the circuit and liquid chamber, and during the rising period of the hydraulic braking torque, there are viscosity resistance, flow force, and orifice compensation of the hydraulic braking system.

Because of the differences in dynamic response characteristics between the MBS and the HBS, the total braking torque fluctuates significantly during mode switching, which cannot meet the braking torque required by the driver, and may also lead to an increase in the jerk of the complete vehicle and the false trigger of the ABS braking, as shown in the simulation part of this study. So it is necessary to

coordinate the motor braking torque and hydraulic braking torque in the process of mode switching, to ensure the stability of the braking torque during braking.

**Figure 7.** Dynamic response characteristics of the HBS and the motor braking system (MBS).

#### **4. The Torque Coordinated Control Strategy of Mode Switching**

#### *4.1. The Condition Analysis of Mode Switching*

By controlling the working state of the motor and switching of the HBS coupling valve, the electro-hydraulic composite braking system of the electric vehicle can realize various braking modes. According to the braking force distribution control strategies of this paper, the working state of each component of the braking system of each braking mode are shown in Table 7. Based on the response differences between the MBS and HBS, and setting aside the discontinuous change due to brake pedal action, the processes responsible for the torque fluctuation during mode switching mainly exist in the working conditions that the braking torque step changed. Therefore, the coordinated control strategy is mainly applied to the braking conditions in which the braking torque of the MBS or HBS step changes.


**Table 7.** The working state of each component in the different braking modes.

In Table 7, "•" represents MBS working or HBS working, and "-" represents MBS not in operation or HBS not in operation.

#### *4.2. The Design of the Coordination Controller*

The dynamic coordination control strategy of brake mode switching developed in this paper is shown in the Figure 8. Firstly, the target hydraulic braking torque, *Th*\_*req*, and the target motor braking torque, *Tm*\_*req*, are preliminarily distributed based on the vehicle state parameters by the braking force distribution controller. Secondly, the target hydraulic braking torque is modified through fuzzy control rules based on the pedal opening and its change rate, to reflect the driver's braking intention. Then, the target motor braking torque is corrected by the actual hydraulic braking torque (*Th*) output by the HBS, and the rapid response of the MBS is used to compensate for the hydraulic braking torque, to achieve a smooth transition of the braking mode. Finally, the HBS and the MBS are controlled by the hydraulic and motor braking controller, respectively, to respond to the modified target braking torque *T h*\_*req* and *T <sup>m</sup>*\_*req*, and then output the actual hydraulic braking torque, *Th*, and actual motor braking torque, *Tm*, to decelerate the vehicle.

**Figure 8.** Dynamic coordination control structure.

#### 4.2.1. The Modification of Target Hydraulic Braking Torque

The target change rates of target hydraulic braking torque during braking are

$$k\_{\rm l\\_req} = \frac{d}{dt} T\_{\rm h\\_req} \tag{11}$$

where *kh*\_*req* represents the target change rates of hydraulic braking torque.

During the mode switching, the upper limit of the change rate of the hydraulic braking torque, *kh*, is determined by the fuzzy controller, which is designed based on the brake pedal stroke, α, and its change rate, . α. Then the upper limit of the change rate, *kh*, is compared with the target change rate *kh*\_*req*, to determine the modified target change rate.

$$k'\_h = \min(k\_{h\\_rcq\\_},\_{h})\tag{12}$$

The increment of the modified target hydraulic braking torque, Δ*Th*, can be obtained by the integration of the modified target change rate. Therefore, the modified target braking torque of the HBS is

$$T'\_{h\\_req} = T'\_h + \int \min(k\_{h\\_req}, k\_h) dt\tag{13}$$

where *T <sup>h</sup>*\_*req* represents the target braking torque of the HBS modified by the coordination controller; and *T <sup>h</sup>* is the initial braking torque at the moment when the mode is switched.

During mode switching, the driver's braking intention is reflected by the brake pedal stroke and the brake pedal stroke change rate. Then the fuzzy controller outputs the upper limit of the change rate of the target hydraulic braking torque, so as to realize the modification of the target hydraulic braking torque. The fuzzy subsets of brake pedal stroke, brake pedal stroke change rate, and the upper limit of the change rate of the hydraulic braking torque are {VS, S, MS, M, ML, L, VL}; therefore, the input and output of the fuzzy controller can be described as

$$\begin{aligned} \{\alpha\} &= \{\text{VS, S, MS, M, ML, VL, VL} \\ \{\dot{\alpha}\} &= \{\text{VS, S, MS, M, ML, VL, VL} \\ \{k\_{\text{lt}}\} &= \{\text{VS, S, MS, M, ML, VL, VL} \end{aligned} \tag{14}$$

where VS, S, MS, M, ML, L, and VL represent very small, small medium, medium, medium large, large, and very large, respectively.

The fuzzy control rules are shown in Table 8, and the membership function of input and output variables of fuzzy controller are shown in Figures 9–11. The fuzzy control rules are formulated based on the following experiences:

Criterion 1: If <sup>α</sup> and . α are S, then *kh* is MS. In this case, the brake pedal stroke and its change rate are small; it can be considered that the driver pays more attention to the ride comfort during mode switching, and the upper limit of the change rate of the hydraulic braking torque takes a medium-small value.

Criterion 2: If <sup>α</sup> is S and . α is L, then *kh* is ML. In this case, the braking pedal opening is small and its change rate is large, which indicates that the driver pays more attention to the braking safety during mode switching. Therefore, the upper limit of the change rate of the hydraulic braking torque takes a medium-large value.

Criterion 3: If <sup>α</sup> is VL and . α is S, then *kh* is L. In this case, the braking pedal opening is very large and its change rate is small, indicating that the driver pays attention to both braking safety and ride comfort during mode switching, so the upper limit of the change rate of the hydraulic braking torque takes a large value.

**Table 8.** Fuzzy control rules of target hydraulic braking torque change rate.


**Figure 9.** Pedal opening degree membership function.

**Figure 10.** Pedal opening change rate membership function.

**Figure 11.** The upper limit of the change rate of the hydraulic braking torque membership function.

According to the modified target hydraulic braking torque *T <sup>h</sup>*\_*req*, the hydraulic controller determines the duty ratio of high-speed on–off valve, and directly controls the inlet valve and outlet valve of wheel cylinder through PWM modulation, so as to make the tracking pressure of wheel cylinders follow the modified target pressure changes.

#### 4.2.2. The Modification of Target Motor Braking Torque

According to the analysis in Section 3.3 of this paper, the response time of the MBS is shorter than that of the HBS, thus the rapid response of MBS can be used to compensate for the insufficient braking torque caused by the slow response of HBS, so as to solve the fluctuation of total braking torque and the jerk of the complete vehicle during brake mode switching. Therefore, during mode switching, the MBS need to provide the target motor brake torque *Tm*\_*req*, which is determined by the braking force distribution controller and additionally provide the difference between the target hydraulic braking torque *Th*\_*req* and the current actual hydraulic braking torque *Th*; that is,

$$T\_{m\\_req}' = T\_{m\\_req} + T\_{h\\_req} - T\_h \tag{15}$$

The sum of target hydraulic braking torque *Th*\_*req* and target motor braking torque *Tm*\_*req* is the total braking torque *Tb*\_*req* required by the driver, so Equation (15) can be rewritten as

$$T\_{m\\_req}^{\prime} = T\_{b\\_req} - T\_h \tag{16}$$

The motor control parameters *id* and *iq* are output by the motor controller according to the modified target motor braking torque *T <sup>m</sup>*\_*req*, so that the MBS outputs the actual motor braking torque *Tm* to act on the vehicle.

According to the braking force distribution control strategy of this paper, the HBS starts to provide braking torque when the braking strength required by driver is greater than *z*(*B*). If the braking torque required by driver changes, in order to maintain the coordinated compensation ability of the MBS to the hydraulic braking torque during mode switching, the target braking torque *Tm*\_*req* allocated by the braking force distribution controller to the MBS should be smaller than the maximum braking torque *Tm*\_*max* that can be provided by the MBS. When the braking torque required by the driver continues to increase, and the required braking strength satisfy *z*(*B*) < *z* ≤ *z*(*C*), the target braking torque *Tm*\_*req* determined by braking force distribution controller should increase to *Tm*\_*max* gradually. If the braking torque required by driver remains unchanged, the maximum braking torque that can be provided by the MBS remains at *Tm*\_*max* at the braking torque distribution stage.

#### **5. Simulation and Analysis**

The forward simulation model of a pure electric vehicle was established in MATLAB/Simulink in this study, as shown in Figure 12. This simulation model includes the PMSM I and PMSM II models, the HBS model, the battery model, the vehicle longitudinal resistance model, and the controller model. Because the switching between different braking modes is similar, a specific mode switching can be selected for verification and analysis. The same conclusion can be obtained by the simulation of other switching processes with the proposed torque coordinated control method.

**Figure 12.** The simulation model of the electric vehicle.

The conditions of constant and variable braking strength are simulated to verify the effectiveness of the dynamic coordination control strategy. For the switching from a constant braking strength, the mode switching between the pure electric and pure hydraulic braking modes were selected for the simulation test, and the motor braking switching to a hybrid braking mode was selected for simulation verification in the mode switching of variable braking strength. In the switching from a constant braking strength, the motor and hydraulic braking torque with coordination will respond to the modified target braking torque output by the coordination controller, while the motor and hydraulic braking torque without coordination will respond to the target braking torque that is assigned by the braking force distributor. Under the condition of a variable braking strength, the braking torque and the jerk of the complete vehicle, focusing on safety and ride comfort, were compared, to verify whether the coordinated control strategy reflects the driver's braking intention.

#### *5.1. Simulation and Verification of the Constant Braking Strength of Mode Switching*

#### 5.1.1. Switch from Pure Hydraulic to Pure Electric Braking Mode

The speed variation condition shown in Figure 13a is designed to verify the effectiveness of the coordinated control algorithm during pure hydraulic switching to the pure electric braking mode. In this braking condition, the initial speed of the vehicle is 110 km/h, the initial *SOC* of the battery is 0.6, and the road adhesion coefficient is 0.8. The braking strength required by the driver is increased from 0 to 0.15 within 0 to 0.5 s, and then remains constant.

At the start of braking, the speed of the vehicle is too high for the motors to perform regenerative braking; hence, braking torque is provided by the HBS. After 1.5 s, when the speed of the vehicle has decreased sufficiently, the braking mode switches to pure electric, and braking torque is provided only by the MBS. Because of the compensation of the motor braking torque to the hydraulic braking torque, the *SOC* with coordination is slightly higher than the *SOC* without coordination at the start of braking; moreover, the slow increase of the motor braking torque with coordination during mode switching resulted in a lower *SOC* with coordination than that without coordination, as shown in Figure 13a. As shown in Figure 13b, during the increase of braking strength, the response speed of the HBS is slow due to the orifice compensation and viscosity resistance within the HBS, hence the total braking torque without coordination cannot quickly respond to the target total braking torque required by the driver. The coordinated total braking torque follows the target total braking torque well because the MBS can compensate for the insufficient braking torque caused by the slow response of the HBS.

During mode switching, Figure 13c demonstrates that the motor braking torque without coordination increased rapidly, while the slower responding HBS was still providing a high braking torque. According to Figure 13d, the rapid increase in total braking torque without coordination resulted in a 31.29 m/s<sup>3</sup> jerk of the complete vehicle. In contrast, the motor braking torque with coordination increases as the hydraulic braking torque decreases, and the total braking torque is maintained at a constant level as far as possible. The maximum jerk of the vehicle was 5.91 m/s3; thus, the ride comfort of the vehicle is improved during mode switching.

**Figure 13.** Simulation results of the pure hydraulic switch to a pure motor braking mode under a constant braking strength: (**a**) The variation in vehicle speed and *SOC*, (**b**) the variation in total braking torque, (**c**) the variation in motor and hydraulic braking torque, and (**d**) the variation in vehicle jerk.

#### 5.1.2. Switch from Pure Electric to Pure Hydraulic Braking Mode

The speed variation condition shown in Figure 14a is designed to verify the effectiveness of the coordinated control algorithm during pure electric switching to the pure hydraulic braking mode. In this braking condition, the initial speed of the vehicle is 25 km/h, the initial *SOC* of the battery is 0.6, and the road adhesion coefficient is 0.8. The braking strength required by driver increased from 0 to 0.1 within 0 to 0.5 s, and then remains unchanged.

With the deceleration of the vehicle, the vehicle speed was reduced to 20 km/h at 1.67s, which is too slow for the MBS to maintain a stable regenerative braking torque; hence, the braking mode is switched from pure electric to pure hydraulic. Due to the coordination of motor braking torque to hydraulic braking torque during mode switching, the *SOC* with coordination is higher than that

without coordination, as shown in Figure 14a. The simulation results of Figure 14b,c show that the mode switching without coordination cannot provide the braking torque required by the driver because of the rapid withdrawal of the motor braking torque, giving a wrong braking felling to the driver, resulting in a mis-operation by the driver. With coordination, the motor braking torque decreases as the hydraulic braking torque increases, creating a smooth transition, so that the total braking torque changes gradually and the speed of vehicle is steadily reduced. According to Figure 14d, the maximum jerk of the complete vehicle during mode switching with and without coordination was -3.14 m/s3 and -21.18 m/s3, respectively. Thus, the coordination control strategy improves the vehicle's safety and ride comfort during braking mode switching.

**Figure 14.** Simulation results of the pure hydraulic switch to a pure motor braking mode under a constant braking strength: (**a**) The variation in vehicle speed and *SOC*, (**b**) the variation in total braking torque, (**c**) the variation in motor and hydraulic braking torque, and (**d**) the variation in vehicle jerk.

#### *5.2. Simulation and Verification of Variable Braking Strength of Mode Switching*

Under the condition of a variable braking strength, the initial speed of vehicle is 60 km/h, the initial SOC of battery is 0.6, and the road adhesion coefficient is 0.8. The braking strength required by the driver is increased from 0 to 0.1 within 0 to 0.5 s, and then remains unchanged. At 1.5 s, the driver depressed the brake pedal with two different brake pedal stroke change rates, and the braking strength gradually increases to 0.3. The variation in the braking strength is shown in Figure 15a.

The driver depressed the brake pedal to the same opening with two different brake pedal stroke change rates; therefore, both the switching processes are switched from pure electric to the hybrid braking mode. Because the change rate in brake pedal stroke is different, one switching process focuses on brake safety while the other focuses on ride comfort. The simulation results in Figure 15 are all obtained by the coordinated control strategy proposed in this paper. It can be seen from Figure 15a that since the variation in each braking torque is almost the same for both the safety-focused and comfort-focused intentions, the *SOC* changes of the two braking intention are almost the same, but the SOC with a safety-focused intention is slightly higher than that of the comfort-focused intention.

As shown in Figure 15b, c, since the motor and hydraulic braking torque is modified by the torque coordination controller according to the different braking intentions of the driver, the coordinated total braking torque, motor, and hydraulic braking torque that focus on ride comfort are changed more gently than that focusing on braking safety. As shown in Figure 15d, the jerk of the vehicle during mode switching that focuses on ride comfort is <sup>−</sup>4.84 m/s3, while the mode switching that focuses on braking safety is <sup>−</sup>9.97m/s3. The simulation results show that the coordinated control strategy of this paper can modify the motor and hydraulic braking torque reasonably according to the driver's braking intention, so as to take into account the braking safety and ride comfort during braking mode switching.

**Figure 15.** Simulation results of the pure motor switch to a hybrid braking mode under variable braking strengths: (**a**) The variation in the braking strength and *SOC*, (**b**) the variation in total braking torque, (**c**) the variation in motor and hydraulic braking torque, and (**d**) the variation in vehicle jerk.

#### **6. Conclusions**

The configuration of a four-wheel-drive pure electric vehicle with a dual motor is considered in this paper, and a braking torque dynamic coordinated control strategy, based on the braking intention of the driver, is proposed. The dynamic coordinated control strategy effectively reduces the torque fluctuation and the jerk of the complete vehicle during mode switching.

A controller combining a PID controller and a fuzzy controller is used to adjust the pressure in the wheel cylinder of the HBS to achieve precise control of the hydraulic braking torque. The brake pedal stroke and its change rate are used to reflect the braking intention of the driver, and to modify the target hydraulic braking torque based on a fuzzy control algorithm according to whether the driver pays more attention to brake safety or ride comfort. At the same time, the rapid response of the MBS is used to compensate for the insufficient braking torque caused by the slow response of the HBS, so as to ensure the ride comfort and stability of the braking during mode switching.

Based on the mathematical model of the electro-hydraulic composite braking system, a simulation verification platform is constructed. The typical mode switching of the constant and variable braking strength conditions are simulated to verify the effectiveness of the dynamic coordinated control strategy. The simulation results show that the torque coordination control strategy described in this paper can not only modify the motor and hydraulic braking torque according to the driver's braking intention, but also significantly reduce the braking torque fluctuation and the jerk of the complete vehicle, thereby improving ride comfort and safety. The influence of the proposed control strategy on energy recovery mainly depends on the type of mode switching condition.

Intelligent transportation system (ITS) and vehicle-to-vehicle communication (V2V) are future development trends; therefore, in future research, the authors will combine ITS to make active predictions of braking mode. In addition, energy consumption will be considered by reducing the frequency of mode switching.

**Author Contributions:** Conceptualization, Y.Y. and Y.H.; methodology, Y.Y. and Y.H.; software, Y.H. and Z.C.; validation, Y.H., Z.C. and C.F.; formal analysis, Y.Y.; investigation, Y.H.; resources, Y.Y. and Y.H.; data curation, Z.Y.; writing—original draft preparation, Y.H.; writing—review and editing, Y.Y., Y.H. and C.F.; visualization, Y.H.; supervision, Y.Y., C.F. and Z.Y.; project administration, Y.Y., C.F. and Z.Y.; funding acquisition, Y.Y. and C.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (grant No. 2018YFB0106100), and the National Natural Science Foundation of China (grant No.51575063).

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
