Optimal Management of Thermal Comfort and Driving Range in Electric Vehicles

The HVAC system represents the main auxiliary load in battery-powered electric vehicles (BEVs) and requires efficient control approaches that balance energy saving and thermal comfort. On the one hand, passengers always demand more comfort, but on the other hand the HVAC system consumption strongly impacts the vehicle’s driving range, which constitutes a major concern in BEVs. In this paper, a thermal comfort management approach that optimizes the thermal comfort while preserving the driving range during a trip is proposed. The electric vehicle is first modeled together with the HVAC and the passengers’ thermo-physiological behavior. Then, the thermal comfort management issue is formulated as an optimization problem solved by dynamic programing. Two representative test-cases of hot climates and traffic situations are simulated. In the first one, the energetic cost and ratio of improved comfort is quantified for different meteorological and traffic conditions. The second one highlights the traffic situation in which a trade-off between the driving speed and thermal comfort is important. A large number of weather and traffic situations are simulated and results show the efficiency of the proposed approach in minimizing energy consumption while maintaining a good comfort.


Introduction
Despite rapid evolution of battery performance and recharging infrastructure, the penetration of EVs in road transportation remains hindered by their limited driving range and by users' fear of running out of battery. A lot of research effort is focused the battery itself, but another area for improvement is to better anticipate the energy needs over the whole trip and manage them according to the energy available in the battery. The powertrain is the primary energy consumer, but in addtion the heating, ventilating and air conditioning system (HVAC) represents an important auxiliary load. The HVAC system is expected to provide good thermal comfort to the passengers, regardless of the surrounding context. In practice, its consumption depends on both the climatic conditions and the passengers' comfort requirements. It generally represents 20% of the total vehicle consumption, but it can reach up to 60% in urban areas and harsh conditions and can affect the vehicle driving range. Effective control approaches are therefore required in order to provide an acceptable balance between the passengers' comfort and the HVAC energy consumption, so as to preserve the vehicle autonomy management [21,22]. Some authors succeed at customizing the comfort for passengers, by integrating a learning module [23,24].
Although the above mentioned studies successfully reach a compromise between consumption and thermal comfort, they do not address the issue of the influence of the thermal comfort on the vehicle driving range. In EVs, this may be critical and the energy allocated to the HVAC should be constrained, if necessary in order for the driver to reach a charging station. This would require to work with a prediction horizon covering the whole trip. Another limitation of most published material is that they do multi-objective optimization through an aggregate cost functions, namely a weighted sum of a thermal comfort criterion and the HVAC consumption. Lastly, the thermal comfort criterion itself is an important issue. Very often, thermal comfort is simply modeled by a target temperature of 25 • C. However, this criterion alone can lead to very different sensations, depending on other parameters such as humidity, wind, or solar irradiance. Therefore, more sophisticated tools have been investigated in order to define and quantify thermal comfort, as explained in the next section.

Thermal Comfort
All the above mentioned studies model thermal comfort through a temperature set point. However, thermal comfort is a complex notion that cannot be reduced to a simple "comfortable temperature". Other physical quantities, such as humidity, are also important, but not only. Thermal comfort is defined by ASHRAE 55 standard [25] as "a subjective concept characterized by a sum of sensations, which produce a person's physical and mental wellbeing, condition for which a person would not prefer a different environment". The state of thermal comfort felt by a person is in close connection with their physical and mental conditions. Different thermal comfort indexes have been proposed over the years, the first two being the predicted mean vote (PMV) and the associated predicted percentage of dissatisfied (PPD) [26]. These two comfort indexes were developed for steady state conditions in buildings. In addition, they do not take into account the thermoregulatory response of the subject, which limit their ability to accurately assess the thermal comfort. Other indexes have been proposed to better assess thermal comfort, especially in vehicle cabin. The most important one is the equivalent temperature model, which is standardized in both ASHRAE [25] and ISO [27] standards. The equivalent temperature is a physical concept defined as "the temperature of an imaginary enclosure with a mean radiant temperature equal to the air temperature, still air, and in which a person would have the same heat exchange by convection and radiation as in the actual conditions" [28]. This definition holds if applied to the whole body or to each body part. In the first case, we will talk about the global body equivalent temperature, while in the second case we will talk about local equivalent temperature.
If one sticks to the above definitions, evaluating the equivalent temperature requires to computing two physical quantities: the heat exchange between the body and the surroundings on the one hand, and the skin temperature on the other hand. This can be done using a thermo-physiological model. Thermo-physiological models are mathematical descriptions of the complex physiological responses and heat transfer of the human body. These models are valuable tools for understanding the human thermoregulation reactions that maintain the core temperature around 37 • C, and for evaluating thermal comfort and sensations. The development of advanced thermoregulation models have started at the United States National Aeronautics and Space Administration (NASA) and the United States army for the sole purpose of assessing the effect of extreme environmental conditions on the human body [29]. In recent years, various thermo-physiological models have been developed to improve the prediction of thermal responses of buildings' occupants. The most referenced includes models developed by Gagge [30], Stolwijk [31], Tanabe [2], Fiala [32][33][34], the Berkeley Model [35,36] and ThermoSEM [37,38].
The joint use of a model of human thermo-physiology and a model of the thermal surrounding environment allows to calculate the equivalent temperature of all the body segments in given conditions. This equivalent temperature is then interpreted in terms of thermal comfort level, as expressed by a large number of subjects and represented by the mean thermal vote (MTV) scale [28]. Figure 1 Energies 2020, 13, 4471 5 of 31 shows the graphical representation of thermal comfort levels proposed in Nilsson's thesis [28]: for each body segment, three thermal comfort zones are defined as a function of the local equivalent temperature. The green zone represents the range of local equivalent temperatures that corresponds to the comfortable zone. The red and blue zones represent the "hot comfortable" and "cold comfortable" zones respectively, in which the subject feels either hot or cold but in a comfortable way. The bold dash line represents the equivalent temperature set points T eq,sp that corresponds to ideal comfort and that will be used latter in this work (Section 4.1). The comfortable temperature range depends on the subject's clothing. Figure 1a corresponds to an informal dress, whereas Figure 1b corresponds to a formal one, both in summer. Comfortable temperatures are higher in the first situation, because the subject is lightly dressed (short, tee-shirt). Other diagrams exists in the literature for different clothing, worn for instance in spring and winter.
Energies 2019, 12, x FOR PEER REVIEW 5 of 30 [28]: for each body segment, three thermal comfort zones are defined as a function of the local equivalent temperature. The green zone represents the range of local equivalent temperatures that corresponds to the comfortable zone. The red and blue zones represent the "hot comfortable" and "cold comfortable" zones respectively, in which the subject feels either hot or cold but in a comfortable way. The bold dash line represents the equivalent temperature set points , that corresponds to ideal comfort and that will be used latter in this work (Section 4.1). The comfortable temperature range depends on the subject's clothing. Figure 1a corresponds to an informal dress, whereas Figure 1b corresponds to a formal one, both in summer. Comfortable temperatures are higher in the first situation, because the subject is lightly dressed (short, tee-shirt). Other diagrams exists in the literature for different clothing, worn for instance in spring and winter.

System Description and Modeling
The studied system consists of four sub-systems: the HVAC system, the powertrain, the battery, which provides energy to the powertrain and to the HVAC system, and finally the human body through which thermal comfort is evaluated by the equivalent temperature index. The present work deals with thermal comfort in hot climates, and we will consider only the cooling function of the HVAC system. This section presents the four sub-systems and the associated models. At first reading, the reader does not need to go into the details of all the equations to understand the working principle of the system and the optimization problem stated in Section 4. The reader can come back to these equations later, in order to fully understand the results presented in Section 5. Figure 2 shows a schematic drawing of the HVAC system. It consists of three main components: one is the ventilation circuit (grey loop) that blows conditioned air in the cabin, and the two others are the cooling and heating systems (resp. blue and orange boxes) that exchange thermal power with the blown air in order to adjust its temperature and humidity.

HVAC System Description and Modeling
As depicted, the outside air flow is mixed with a certain amount of recirculated cabin air and blown by a fan into the ventilation system. The air is firstly cooled and dried up by yielding thermal energy to the cold loop of the cooling system through the evaporator, as detailed in Figure 3. At this stage, it is too cold to be directly blown into the cabin and it needs to be warmed up in the heating system. This can be done at a zero energetic cost by exchanging with the cooling system of the traction motor.

System Description and Modeling
The studied system consists of four sub-systems: the HVAC system, the powertrain, the battery, which provides energy to the powertrain and to the HVAC system, and finally the human body through which thermal comfort is evaluated by the equivalent temperature index. The present work deals with thermal comfort in hot climates, and we will consider only the cooling function of the HVAC system. This section presents the four sub-systems and the associated models. At first reading, the reader does not need to go into the details of all the equations to understand the working principle of the system and the optimization problem stated in Section 4. The reader can come back to these equations later, in order to fully understand the results presented in Section 5. Figure 2 shows a schematic drawing of the HVAC system. It consists of three main components: one is the ventilation circuit (grey loop) that blows conditioned air in the cabin, and the two others are the cooling and heating systems (resp. blue and orange boxes) that exchange thermal power with the blown air in order to adjust its temperature and humidity.    Using the notations defined in Figure 2, the ventilation sub-system can be modeled as follows. After mixing, the temperature and the specific humidity of the air (mass of water vapor present in a unit mass of moist air) are given by (1) and (2), where is the recycling ratio:

HVAC System Description and Modeling
= .
The air flows through the cooling and heating sub-system, after which the temperature and specific humidity are denoted by and . The HVAC system enables controling these quantities, as it will be described later on in this section.
The conditioned air is blown into the cabin, and the cabin temperature evolution is governed by (3), where is the cabin thermal capacitance (including walls and seats). is the thermal power exchanged through the windshield and walls by convection.
corresponds to the heat  As depicted, the outside air flow is mixed with a certain amount of recirculated cabin air and blown by a fan into the ventilation system. The air is firstly cooled and dried up by yielding thermal energy to the cold loop of the cooling system through the evaporator, as detailed in Figure 3. At this stage, it is too cold to be directly blown into the cabin and it needs to be warmed up in the heating system. This can be done at a zero energetic cost by exchanging with the cooling system of the traction motor.   Using the notations defined in Figure 2, the ventilation sub-system can be modeled as follows. After mixing, the temperature and the specific humidity of the air (mass of water vapor present in a unit mass of moist air) are given by (1) and (2), where is the recycling ratio: = .
The air flows through the cooling and heating sub-system, after which the temperature and specific humidity are denoted by and . The HVAC system enables controling these quantities, as it will be described later on in this section.     Using the notations defined in Figure 2, the ventilation sub-system can be modeled as follows. After mixing, the temperature and the specific humidity of the air (mass of water vapor present in a unit mass of moist air) are given by (1) and (2), where is the recycling ratio: = .
The air flows through the cooling and heating sub-system, after which the temperature and specific humidity are denoted by and . The HVAC system enables controling these quantities, as it will be described later on in this section.
The conditioned air is blown into the cabin, and the cabin temperature evolution is governed by  Using the notations defined in Figure 2, the ventilation sub-system can be modeled as follows. After mixing, the temperature T and the specific humidity φ of the air (mass of water vapor present in a unit mass of moist air) are given by (1) and (2), where β is the recycling ratio: Energies 2020, 13, 4471 7 of 31 The air flows through the cooling and heating sub-system, after which the temperature and specific humidity are denoted by T out and φ out . The HVAC system enables controling these quantities, as it will be described later on in this section.
The conditioned air is blown into the cabin, and the cabin temperature evolution is governed by (3), where C cab is the cabin thermal capacitance (including walls and seats).
. Q ext is the thermal power exchanged through the windshield and walls by convection. . Q f low corresponds to the heat exchange due to the air circulation.
. Q driver represents the heat flow rate produced by the driver. Its mathematical expression (22) will be presented in Section 4.1: The cabin specific humidity φ cab evolves according to the water mass balance Equation (4), where V cab , ρ air and . m air respectively denote the cabin volume, the air density and the air flow rate. The term . m sw represents the mass flow of water vapor produced by sudation. Its mathematical expression (23) will be presented in Section 4.1: The cabin wall temperature T w evolves according to Equation (5), where C wall is the cabin wall heat capacity and . Q wall is the heat exchange between the wall and the cabin air. The term . Q sun represents the thermal power received from solar radiation, and is expressed by the product of solar irradiance . q sun W/m 2 and total wall surface S w m 2 : .
Let us now focus on the cooling system, sketched in Figure 3. The cold loop is composed of four main entities: compressor (a), condenser (b), expansion valve (c) and evaporator (d). A refrigerant fluid circulates throughout these entities and undergoes a thermodynamic cycle during which it receives thermal energy from the ventilated air and yields it outside.
The thermodynamic cycle is depicted in Figure 4. The first stage is the compression, in which the compressor brings the refrigerant from a low pressure, gaseous state to a higher pressure and temperature state (transformation 1 → 2 ). Then, during the condensation stage, the fluid releases heat to the outside air, while turning into a high pressure saturated liquid state (transformation 2 → 5 ). The moto-ventilator group forces the circulation of outside air in order to help evacuating heat from the refrigerant. Next, the saturated liquid undergoes an isenthalpic expansion through the expansion valve: its pressure and temperature decrease and it turns into a cold two-phase fluid (transformation 5 → 6 ). Finally, during the evaporation phase, this two-phase fluid absorbs heat from the ventilation circuit air. It evaporates and exits the evaporator at low-pressure gas state (transformation 6 → 1 ). A crucial point to note is that the whole cold loop is controlled by the compressor rotational speed, denoted by N comp . The cold loop model consists in evaluating at each time step, the refrigerant enthalpy at each point of the thermodynamic cycle, and then the heat exchanges with the air and the power consumed by the compressor. Practically, the model computes the refrigerant evaporating and condensing temperature T 7 (evaporator) and T 4 (condenser) by solving the following heat balances Equations (7) and (8). In these equations, is the specific enthalpy (J/kg) at point i, calculated using the refrigerant temperature and pressure. The quantity . m re f is the mass flow rate of the refrigerant, controlled by the compressor rotational speed. The complete model of cold loop is detailed in Appendix A: .
When exiting the cooling system, the air is very dry and cold: only a few degrees Celsius. Hence, it needs to be warmed up before being blown into the cabin. The heating system consists of two heating resistors and an exchanger with the water cooling system of the electric machine. For the weather conditions considered here, the heating resistors are not needed and free heating is provided by the water cooling system of the electric machine. Again, an important point to note is that the heating is controlled by the ratio of air flow derived into the heat exchanger, denoted by α.
The next point deals with the energy consumption of the HVAC system. Three components need power feeding: the compressor (a), the fan (e), and the moto-ventilator group MVG (f). The compressor consumption P elec,comp is given by (10), where h 1 and h 2 denote the refrigerant enthalpy at the points 1 and 2 of the thermodynamic cycle represented in Figure 4, . m re f is the refrigerant mass flow rate (linked to the compressor rotational speed), η meca and η elec are the mechanical and electrical efficiencies of the compressor: The consumptions of the fan and the MVG depend on the external temperature according to tabulated data provided by Groupe PSA. The total electric consumption of the HVAC system P HVAC is the sum of the electric power consumed by the compressor, the fan and the MVG (11). P HVAC = P elec,comp + P f an + P MVG (11)

Powertrain Model
The second sub-system is the powertrain, modeled using a forward model in a very classical way. The electrical machine provides a torque T EM , which is transmitted to the wheels and converted there into a tractive force F traction . The resulting vehicle speed evolution is calculated according to (12), where m is the vehicle equivalent mass, accounting for all moving parts, and F road (t) denotes the sum of the external forces applied to the vehicle: aerodynamic drag, gravitational force, and rolling resistance [1]. F traction is negative during regenerative braking: The transmission chain between the electrical machine and the wheels is modelled by a fixed speed ratio, and a variable efficiency given by tabulated data as a function of torque and speed. The electrical machine and its control electronics are modelled by a measured losses map. Hence, the electrical machine power consumption P EM is given by Equation (13), where ω EM denotes the machine rotation speed, which is proportional to the vehicle speed: Energies 2020, 13, 4471 9 of 31 In the present study, we consider that the driver is modeled as a PI regulator [39] that controls the electrical machine torque in order to follow the speed profile of a given driving cycle. The actual speed profile v(t) and the resulting electrical consumption P EM (t) are calculated according to the equations afore mentioned.

Battery Model
The battery provides the energy to the electrical machine, the HVAC system and various auxiliaries, according to Equation (14): The battery is modeled by its open circuit voltage V bat and internal resistance R bat [40][41][42]. Formulas (15) and (16) give the current as a function of the power and the resulting battery state of charge SOC variation. Q 0 is the nominal battery capacity:

Thermo-Physiological Model
Human beings are homeotherms. This means that they can regulate their core (internal) body temperature with physiological and behavioral actions, within a certain range of ambient thermal conditions. The nature of the human body is to maintain a constant body core temperature around 36.5 • C at rest, and up to 41 • C during heavy exercise. Exposure to extreme environmental conditions can lead however to a poor regulation thus inducing hyperthermia and hypothermia.
The human body interaction with the environment and the regulation mechanism have been the focus of a number of research work. Many different thermo-physiological models have been developed. Thermo-physiological models simulate the human thermal responses in a complex way by describing the heat transfer phenomena inside the body, and also the heat exchange with the surrounding environment. Generally speaking, these types of models are able to take into account: (1) the body constitution (weight, height, fat percentage, etc.), thermoregulation and cardiovascular system; (2) environmental conditions (air temperature, air speed, relative humidity, and mean radiant temperature); (3) personal factors (activity level, clothing insulation, and water vapor permeability). One of the most important thermo-physiological models is the one developed by Tanabe in 2002 [2].
The Tanabe's model is a segmented physiological model. It represents an average man with a weight of 74.43 kg and a body surface area of 1.87 m 2 . As depicted in Figure 5, the body is divided into 16 segments (head, chest back, pelvis, left hand, right hand, left thigh, right thigh, left leg, right leg, left foot, and right foot), indexed by letter i in the model Equations (17)- (21). Each body segments is in turn divided into four layers (core, muscle, fat, and skin), indexed by the letter j. Each 2-uplet (i, j) constitutes a node of the model. The last node is the blood, which brings the model to a 65 node thermo-physiological model. A set of 65 differential equations is built by writing the heat balance at each node of the model. into 16 segments (head, chest back, pelvis, left hand, right hand, left thigh, right thigh, left leg, right leg, left foot, and right foot), indexed by letter in the model Equations (17)(18)(19)(20)(21). Each body segments is in turn divided into four layers (core, muscle, fat, and skin), indexed by the letter . Each 2-uplet ( , ) constitutes a node of the model. The last node is the blood, which brings the model to a 65 node thermo-physiological model. A set of 65 differential equations is built by writing the heat balance at each node of the model. The four layers of each of the sixteen body parts are modeled by the heat balance Equations (17)(18)(19)(20), where ( , ) and ( , ) denote the temperature and the heat capacity of the node ( , ), whereas and denote the temperature and heat capacity of the blood. Equation (21) represents the heat balance of the blood, which circulates through the whole body: The four layers of each of the sixteen body parts are modeled by the heat balance Equations (17)-(20), where T (i,j) and C (i,j) denote the temperature and the heat capacity of the node (i, j), whereas T b and C b denote the temperature and heat capacity of the blood. Equation (21) represents the heat balance of the blood, which circulates through the whole body: The term . Q (i,j) m models the metabolic heat flow produced at the node (i, j). It depends on the physical activity of the person. The conduction between adjacent layers is modeled by d , which corresponds to the heat flowing from the node (i, j) to the node (i, j + 1). The convection through the blood circulation is modeled by b , which represents the heat flowing from the node (i, j) to the blood. The heat exchange with the environment due to respiration is modeled by res , a term that exists only in the core layer of the chest segment (node (2, 1)). The term δ 2,i refer to the Kronecker symbol. Other heat transfers to the environment take place at the skin level (nodes (i, 4)), due to convection, radiation and perspiration evaporation, respectively modeled by the heat fluxes . If a node is in contact with the seat, the heat exchanges by convection, radiation, and evaporation are considered to be null. They should be replaced by the heat contact between the segment and the seat, but as our cabin model does not account for seats, the contact heat fluxes are not considered here.
Thermoregulation depends on perspiration. The thermo-physiological model is mainly used to compute the skin's temperature and heat exchanges by convection and radiation. Those quantities are needed to evaluate the thermal comfort of each body segment, and then the body's global comfort. The model also allows to quantify: (i) the total heat transfer between the driver's body and the cabin's air . Q driver , calculated by (22) and reported in (3), and (ii) the mass flow rate of water vapor due to sudation . m sw , given by (23), where λ H 2 O is the latent heat of water, and reported in (4): .
For more information about the complete model, we refer readers to the original paper of Tanabe [30].

Thermal Comfort Index
The next step is to model thermal comfort. In the present work, the local equivalent temperature is used to build a thermal comfort criterion, since it is more accurate and suitable to vehicular environment than other thermal comfort indexes such as a set point temperature. We recall that the equivalent temperature is defined as the temperature of an imaginary enclosure with still air, and a wall temperature equal to the air temperature, in which a person has the same heat exchange by convection and radiation as in the actual conditions. Figure 6 illustrates this concept: the left figure corresponds to the real environment and the right one to the imaginary environment. The posture, the activity level and the clothing are the same in both environments. The real environment is characterized by given T cab , T w and v air . In these conditions, the driver exchanges heat flow . Q real with the surroundings, due to the convection and radiation. The imaginary environment is defined by a uniform temperature T cab = T w = T eq and still air (v air = 0 m/s). In these conditions, the convection and radiation heat exchange between the driver and the surroundings is . Q eq . By definition of the equivalent temperature, . Q eq is equal to . Q real .
Energies 2019, 12, x FOR PEER REVIEW 11 of 30 characterized by given , and . In these conditions, the driver exchanges heat flow with the surroundings, due to the convection and radiation. The imaginary environment is defined by a uniform temperature ( = = ) and still air ( = 0 m/s). In these conditions, the convection and radiation heat exchange between the driver and the surroundings is . By definition of the equivalent temperature, is equal to . The concept of equivalent temperature can be applied at a local level: the driver's body is divided into 16 segments (i). For each of them, a local equivalent temperature ( ) is mathematically defined as the solution of the Equation (24): where ( ) is the sum of the heat exchanged by convection ( , ) and radiation ( , ) , both calculated by the thermo-physiological model. The heat exchange ( ) is thus, given by Equation (25): The concept of equivalent temperature can be applied at a local level: the driver's body is divided into 16 segments (i). For each of them, a local equivalent temperature T (i) eq is mathematically defined as the solution of the Equation (24): real is thus, given by Equation (25): The term . Q (i) eq represents the equivalent convection and radiation heat exchange for an air temperature and wall temperature equal to T (i) eq , and a null air speed.
. Q (i) eq is given by Equation (26): Q eq , the local equivalent temperature T (i) eq is obtained by solving Equation (26) for each body segment. These equivalent temperatures are used to build a global thermal comfort criterion in Section 4.1.

Optimal Energy Management
The objective of the HVAC management system is to control the cabin temperature, the relative humidity and the wall temperature, so that the driver feel comfortable at the lowest energy cost. On the other hand, the driver first priority is to reach his destination, and the battery state of charge may be too low in order to insure the ideal comfort during the whole trip. Hence, it may be necessary to limit the power provided to the air conditioning, in order to extend the vehicle driving range up to its final destination. This energy management problem can be described as an optimal control problem and we propose to use dynamic programing to solve it. In the rest of this section, we first introduce the cost function and then formulate the optimal energy management problem.

Cost Function
The global comfort criterion is built by comparing each local equivalent temperature to the local ideal comfort temperature T eq,sp,i . As a first approximation, we propose to use a global quadratic discomfort criterion L, defined for the whole body by Equation (27): where T eq,sp,i is the set point equivalent temperature for the body segment i, whereas T eq,max,i and T eq,min,i are the upper and lower limit values of the comfort zone. We use equivalent temperature data extracted from Nilsson thesis [28]. In addition to the discomfort criterion L, a mean global equivalent temperature for the whole body is defined by:

Problem Formulation
The role of the proposed energy management is to determine the best HVAC system control in order to reach the lowest thermal discomfort for a given trip. The trip is modeled by a driving cycle, which is known in advance. The powertrain model is used to calculate the instantaneous power required for the electric machine P EM (t), and then the total energy required for the traction over the considered trip. The difference between the energy embedded in the battery at the beginning of the trip and the energy needed for the traction gives the energy available for air conditioning, denoted by E HVAC max .
The control model of the whole system is built by assembling the models of the four sub-systems previously described in Section 3. Two state variables are defined: x = [T cab , φ cab , T wall ] , which corresponds to the thermodynamic quantities needed to estimate the passenger's thermal comfort, and E HVAC the energy consumed by the HVAC system. The four control variables are gathered in a vector, denoted by u = N comp , . m air , α, β The outside weather and the vehicle speed constitute the Using these notations, the overall system model can be formalized by Equations (29) and (30), where f is the thermodynamic evolution equation of the system: . .
Different operational constraints, such as the battery power limitation, complete the model but they are not all detailed here.
The optimization problem consists in determining the command u that minimizes the global discomfort (31) for the considered trip: The finite energy constraint is an inequality end constraint (32):

Problem Solving by Dynamic Programing
The optimization problem is discretized in time and energy and reformulated as a multi-stage decision-making process, which can be effectively solved by DP algorithms [43,44]. The system state trajectory is discretized in the 2d-space t k , E HVAC,i k=0,N;i=0,M , where t k and E HVAC, i respectively denote the discrete values of time t k = t 0 + k.∆t, and energy E HVAC,i = M.∆E HVAC . In the rest of the paper, the subscript k refers to the value of the corresponding quantity at time t k . Using this notation, the command to look for is u = {u k } k=0,N−1 and Equations (29)-(31) become (33)- (35): Furthermore, since we work in the discretized space t k , E HVAC,i , u k must be such that E HVAC,k+1 is equal to one of the discrete values E HVAC,i . DP is usually implemented starting from a given final state, with a backward exploration phase, followed by a forward reconstruction phase. In our case, the final state of the system is not known: we work with an inequality constraint on the final battery state of energy and no constraint at all on the final thermal parameters in the cabin. Hence, dynamic programming must be applied starting from the initial state, with a forward exploration phase followed by a backward reconstruction phase. This does not change the algorithm principle and allows to study the influence of the HVAC energy consumption on the global comfort criterion J, simply by starting the backward reconstruction phase from any final energy state.
At time t 0 , the energy and thermal state are E HVAC,0 and x 0 respectively. Let us assume that E HVAC,M corresponds to the highest energy level, E HVAC,0 = E HVAC,M . Dynamic programming is based on the cost-to-go 2d-table V k,i defined by Equation (36): V k,i represents the lowest cost (cumulated thermal discomfort) to go from the initial state The cost-to-go table and optimal command matrix u * k,i are recursively constructed thanks to (37) and (38), where u k,j,i represents the command to go from t k , E HVAC,j to t k+1 , E HVAC,i and ∆V k,i,j the corresponding cost. The determination of u k,j,i involves solving () to determine the compressor speed N comp . It involves also in the meantime, minimizing the thermal discomfort with respect to the other HVAC control variables . m air , α and β. An infinite cost is affected if the transition from t k , E HVAC,j to t k+1 , E HVAC,i is impossible: Of course, one needs to first calculate V 1,i before starting the recursive algorithm. In addition, several 2d-tables are used in order to keep track of the intermediate optimal commands and corresponding thermal states. These technical details are not described here.
After completing the forward exploration phase, the optimal command matrix u * k,i is used for the reconstruction phase and determining the optimal trajectories for any given consumed energy. The reconstruction phase is done backward, starting from a given energy at the final time t f . For a given t f , E HVAC,j , the optimal command matrix is used to determine optimal state at time t f −1 , and then at t f −2 , t f −3 and so on, until reaching the initial time t 0 . This process enables to reconstruct the optimal trajectories of T cab (t), RH cab (t) and T w (t) and other variables such as P bat (t), for any quantity of energy consumed by the HVAC over the whole trip.

Simulation Results
This section presents the application of the proposed approach and shows its effectiveness in managing thermal comfort on a given trip, for a large number of meteorological and traffic conditions. The first one highlights the tradeoffs that can be reached between the energy consumption of the HVAC system and the thermal comfort. Some preliminary results regarding the HVAC consumption for ideal comfort in different weather conditions will first be presented. In particular, the relationship between steady state cabin temperature and solar irradiance will be studied. Thereafter, we illustrate the possible tradeoff for a given energy consumption, in different weather and traffic situations. Finally, we analyze the state variables, the commands, and the thermoregulation mechanisms trajectories. The second test-case illustrates how the speed profile can be adjusted to reduce the total consumption and preserve the ideal thermal comfort. This is done by applying a k-homothety on the speed profile, so that the driver travels the same distance at a lower speed. Results show how the traction and HVAC system's energy consumption evolve as a function of k. The results gives also, in some situations, the possible energy savings.

Simulation Data and Technical Characteristics of the Vehicle
Each scenario corresponds to a trip completed in given weather and traffic conditions, with a certain type of clothing. As the HVAC system efficiency depends on the vehicle speed, four traffic conditions, representative of congested urban, fluid urban, road and highway environments are modeled, using the UL1, UF1, R2 and A2 INRETS driving cycles (Figure 7). The INRETS driving cycles are a set of ten driving cycles built with data logged around Lyon, France, in various driving environments. The main characteristics of these cycles are reported in Table 1. In order to build a trip with a duration around one hour, each type of cycle is repeated several times. When doing so, one obtains speed profiles with a duration of 1 h 06 min. certain type of clothing. As the HVAC system efficiency depends on the vehicle speed, four traffic conditions, representative of congested urban, fluid urban, road and highway environments are modeled, using the UL1, UF1, R2 and A2 INRETS driving cycles (Figure 7). The INRETS driving cycles are a set of ten driving cycles built with data logged around Lyon, France, in various driving environments. The main characteristics of these cycles are reported in Table 1. In order to build a trip with a duration around one hour, each type of cycle is repeated several times. When doing so, one obtains speed profiles with a duration of 1 h 06 min.  Weather scenarios with the same duration than the speed profiles have been defined. In these scenarios, the solar irradiance is kept constant, while the external temperature and humidity evolve, according to parametrized time profiles. The temperature profile ( ) is shown in Figure 8. The  Weather scenarios with the same duration than the speed profiles have been defined. In these scenarios, the solar irradiance is kept constant, while the external temperature and humidity evolve, according to parametrized time profiles. The temperature profile T ext (t) is shown in Figure 8. The temperature rises from T ext,0 at the beginning of the trip to T ext,0 + 2 • C at halfway, and then decreases back to T ext,0 . As for humidity, weather conditions are described in terms of relative humidity RH, whereas thermodynamic equation are written using the specific humidity φ. The relationship between these two quantities is given by (39), where P atm is the atmospheric pressure and P sat denotes the saturated vapor pressure of water at the considered temperature: Energies 2020, 13, 4471

of 31
The chosen scenario assumes a linear increase of the specific humidity over the whole trip, from , to , + 0.05 , where , is calculated according to the values of the initial relative humidity , and temperature , . Figure 8 shows an example of the time evolution of the temperature, relative humidity and solar power for the following scenario parameters: , = 35 °C, , = 55% and , = 1000 W/m . Two clothing levels are considered. The first one corresponds to a driver dressed in an informal way (t-shirt, short and shoes), whereas the second one corresponds to a formal dress code (shirt, underpants, trousers, socks and shoes). The comfort zones corresponding to each dress code are depicted in Figure 1.
The initial temperature and relative humidity in the cabin are assumed the same as outside (vehicle parked in the shade). For all the tested scenarios, we consider that driver is at rest, seated.
In summary, the benchmark built from the abovementioned scenarios corresponds to four given trips carried out in 36 weather conditions and two clothing levels. This combination brings the total number of simulated test cases at 4 × 36 × 2 = 288 The vehicle is assumed to be a small size car, with a nominal power around 60 kW and a 40 kWh Li-ion battery. A summary of simulation data is reported in Table 2. The results are presented in the following sections.  The chosen scenario assumes a linear increase of the specific humidity over the whole trip, from φ ext,0 to φ ext,0 + 0.05, where φ ext,0 is calculated according to the values of the initial relative humidity RH ext,0 and temperature T ext,0 . Figure 8 shows an example of the time evolution of the temperature, relative humidity and solar power for the following scenario parameters: T ext,0 = 35 • C, RH ext,0 = 55% and Two clothing levels are considered. The first one corresponds to a driver dressed in an informal way (t-shirt, short and shoes), whereas the second one corresponds to a formal dress code (shirt, underpants, trousers, socks and shoes). The comfort zones corresponding to each dress code are depicted in Figure 1.
The initial temperature and relative humidity in the cabin are assumed the same as outside (vehicle parked in the shade). For all the tested scenarios, we consider that driver is at rest, seated.
In summary, the benchmark built from the abovementioned scenarios corresponds to four given trips carried out in 36 weather conditions and two clothing levels. This combination brings the total number of simulated test cases at 4 × 36 × 2 = 288 The vehicle is assumed to be a small size car, with a nominal power around 60 kW and a 40 kWh Li-ion battery. A summary of simulation data is reported in Table 2. The results are presented in the following sections.

Ideal Comfort Results
This section reports some observations and analysis on three physical quantities: (i) the HVAC system consumption E * HVAC that provide ideal comfort, (ii) the ratio r * HVAC of E * HVAC with respect to the total energy (sum of E * HVAC and traction energy), (iii) the steady state cabin temperature T cab,steady . Figures 9 and 10 summarize the results of E * HVAC for a driver dressed in an informal and formal way, in a suburban traffic (i.e., INRETS R2 driving cycle). The graphs show the evolution of E * HVAC as a function of solar irradiance for different T ext,0 and RH ext,0 . We observe that E * HVAC is an increasing function of the solar irradiance, outside temperature and relative humidity except in some situations.
When the driver is dressed in an informal way, and that . q sun = 0 and T ext,0 ≥ 32 • C, the energy consumed during the trip does not depend on the outside temperature, nor on the relative humidity. The reason is that at night, the comfort temperature is slightly less than 28 • C, as it will be explained later in this section. This temperature can be maintained by operating the compressor at the minimal rotational regime, using the other control variables. Since the compressor is by far the main consumer in the HVAC system, E * HVAC is consequently almost the same for all scenarios where . q sun = 0 and T ext,0 ≥ 32 • C. The same explanation applies to Figure 10a. , . Figures 9 and 10 summarize the results of for a driver dressed in an informal and formal way, in a suburban traffic (i.e., INRETS R2 driving cycle). The graphs show the evolution of * as a function of solar irradiance for different , and , . We observe that * is an increasing function of the solar irradiance, outside temperature and relative humidity except in some situations. When the driver is dressed in an informal way, and that = 0 and , ≥ 32 °C, the energy consumed during the trip does not depend on the outside temperature, nor on the relative humidity. The reason is that at night, the comfort temperature is slightly less than 28 °C, as it will be explained later in this section. This temperature can be maintained by operating the compressor at the minimal rotational regime, using the other control variables. Since the compressor is by far the main consumer in the HVAC system, * is consequently almost the same for all scenarios where = 0 and , ≥ 32 °C. The same explanation applies to Figure 10a. In a typical hot summer day (i.e., , = 32 °C, , = 55%, = 1000 W/m ), with a driver dressed in an informal way, * reaches 1.3 kWh. This consumption increases if the driver is dressed in a formal way, and reaches 1.8 kWh for the aforementioned weather conditions.    . HVAC system consumption E * HVAC for a driver dressed in an informal way, in suburban traffic (INRETS R2), for four T ext (t) and three RH ext (t) characterized by: (a) RH ext,0 = 35%, (b) RH ext,0 = 55%, and (c) RH ext,0 = 75%. , formal way, in a suburban traffic (i.e., INRETS R2 driving cycle). The graphs show the evolution of * as a function of solar irradiance for different , and , . We observe that * is an increasing function of the solar irradiance, outside temperature and relative humidity except in some situations. When the driver is dressed in an informal way, and that = 0 and , ≥ 32 °C, the energy consumed during the trip does not depend on the outside temperature, nor on the relative humidity. The reason is that at night, the comfort temperature is slightly less than 28 °C, as it will be explained later in this section. This temperature can be maintained by operating the compressor at the minimal rotational regime, using the other control variables. Since the compressor is by far the main consumer in the HVAC system, * is consequently almost the same for all scenarios where = 0 and , ≥ 32 °C. The same explanation applies to Figure 10a. In a typical hot summer day (i.e., , = 32 °C, , = 55%, = 1000 W/m ), with a driver dressed in an informal way, * reaches 1.3 kWh. This consumption increases if the driver is dressed in a formal way, and reaches 1.8 kWh for the aforementioned weather conditions.  In a typical hot summer day (i.e., T ext,0 = 32 • C, RH ext,0 = 55%, . q sun = 1000 W/m 2 ), with a driver dressed in an informal way, E * HVAC reaches 1.3 kWh. This consumption increases if the driver is dressed in a formal way, and reaches 1.8 kWh for the aforementioned weather conditions. Figures 11 and 12 summarize the results of r * HVAC for a driver dressed in the two ways presented in the previous section, in a moderate humid weather (i.e., RH ext,0 = 55%). The graphs show the values of r * HVAC with respect to four driving cycles: INRETS A2, R2, UF1 and UL1, for different T ext,0 and . q sun,0 . The results illustrate the relative high value of r * HVAC in congested and fluid urban conditions. This is mainly due to the low traction energy at low speed. For instance, in a typical hot summer day, in congested urban condition, this ratio reaches 80%. Results also show that for the same driving cycle, this ratio increases with outside temperature and solar irradiance.
conditions. This is mainly due to the low traction energy at low speed. For instance, in a typical hot summer day, in congested urban condition, this ratio reaches 80%. Results also show that for the same driving cycle, this ratio increases with outside temperature and solar irradiance.  It is worth underlining the capacity of the equivalent temperature model to account for the cabin wall radiant temperature. For a given clothing and solar irradiance, the steady state wall and cabin temperature are computed. Then, they are averaged over the number of outside temperatures and relative humidity profiles, characterized by , and , . As an illustrative example, we report results for different solar irradiance levels ∈ 0 ; 500 ; 1000 and two dress codes, for the INRETS R2 driving cycle. Table 3 reports the corresponding mean steady state cabin temperature   conditions. This is mainly due to the low traction energy at low speed. For instance, in a typical hot summer day, in congested urban condition, this ratio reaches 80%. Results also show that for the same driving cycle, this ratio increases with outside temperature and solar irradiance.  It is worth underlining the capacity of the equivalent temperature model to account for the cabin wall radiant temperature. For a given clothing and solar irradiance, the steady state wall and cabin temperature are computed. Then, they are averaged over the number of outside temperatures and relative humidity profiles, characterized by , and , . As an illustrative example, we report results for different solar irradiance levels ∈ 0 ; 500 ; 1000 and two dress codes, for the INRETS R2 driving cycle. Table 3 reports the corresponding mean steady state cabin temperature   It is worth underlining the capacity of the equivalent temperature model to account for the cabin wall radiant temperature. For a given clothing and solar irradiance, the steady state wall and cabin temperature are computed. Then, they are averaged over the number of outside temperatures and relative humidity profiles, characterized by T ext,0 and RH ext,0 . As an illustrative example, we report results for different solar irradiance levels . q sun ∈ {0 ; 500 ; 1000} and two dress codes, for the INRETS R2 driving cycle. Table 3 reports the corresponding mean steady state cabin temperature T cab,steady , and the mean steady state wall temperature T wall,steady , for the same trip. In the absence of solar irradiance (at night), T wall,steady = 28.3 • C, the ideal thermal comfort is obtained for T cab,steady = 28.0 • C. In sunny conditions, all other parameters being equals, T wall,steady = 40.0 • C and ideal comfort requires a cooler ambient temperature T cab,steady = 24.1 • C and more energy. If the driver is dressed in a formal way, in sunny conditions, all other parameters being equal, T wall,steady = 38.0 • C and ideal comfort requires an even cooler ambient temperature T cab,steady = 22.1 • C, and hence consumes more energy.

Test Case 1: Thermal Discomfort vs. Energy Cost Trade-off
In this test-case, we investigate how the thermal comfort management algorithm adjusts thermal comfort according to the energy available for the HVAC system in different situations. The thermal comfort management algorithm computes the control actions that optimize the thermal comfort according to the energy available for the HVAC system. For the same available energy, these actions differ, depending on external conditions. In order to show that, we fix an available energy of 1.1 kWh, and analyze how thermal comfort is adjusted for different driving cycle and a given outside temperature. The same analysis will be done for different outside temperature profiles and a given driving cycle. All simulations have been performed for RH ext,0 = 55%, . q sun = 1000 W/m 2 and two clothing ensembles. Figure 13a shows the optimal global equivalent temperature T eq,rms as a function of the HVAC system consumption E HVAC for four INRETS driving cycle: INRETS A2, R2, UF1 and UL1, an outside temperature profile with T ext,0 = 32 • C, and a driver dressed with an informal clothing. Those curves correspond to Pareto frontiers of the E HVAC , T eq,rms bi-criteria optimization problem. For an available energy of 1.1 kWh, T eq,rms is deduced from the Pareto frontiers. The corresponding mean local equivalent temperatures T eq,i , over the last 15 min of the trip, are reported in Figure 13b in the comfort zones diagram. The same results are reported for a driver dressed with a formal clothing in Figure 13c,d. The same conclusions can be formulated from the results depicted in Figure 13c,d, where the driver is dressed in a formal way. However, we notice also that 1.1 kWh is not sufficient to achieve a good thermal comfort for a formal clothing, due to a higher thermal insulation. As consequence, the driver might feel discomfort in a INRETS A2 driving cycle, especially in lower body parts ( Figure  13d).
Let us now do the same analysis for a given driving cycle, and different outside temperature profiles. Figure 14a shows the Pareto frontiers for three outside temperature profiles characterized by , [°C] ∈ 32 ; 35 ; 38 , for an INRETS R2 driving cycle and a driver dressed with an informal clothing. For an available energy of 1.1 kWh, the optimal global equivalent temperature , is deduced from the Pareto frontiers. The corresponding mean local equivalent temperatures , ,, are reported in Figure 14b in the comfort zone diagram. The same results are reported in Figure 14c,d for a driver dressed with a formal clothing.
We notice that the outside temperature has a strong influence on thermal comfort, larger than the vehicle speed. For an available energy of 1.1 kWh, a good thermal comfort is achieved for an outside temperature of , = 29 °C. In this situation, the optimal global equivalent temperature is q sun = 1000 W/m 2 , and where the driver is dressed in: (a,b) informal way, (c,d) formal way.
As expected, the optimal global equivalent temperature is a decreasing function of the HVAC system consumption. For example, in suburban conditions, T eq,rms is improved by 2 • C at an additional energy cost of 430 Wh (1.1% of the battery capacity). We also notice that the thermal comfort consumes more energy in urban traffic (INRETS UL1 and UF1) than in highway (INRETS A2). This can be explained by two concomitant effects. The first one is due to a better efficiency of the HVAC at high vehicle speed v veh . In other words, for the same power consumed by the compressor, the heat exchange at the evaporator is more important. This is mainly due to a better heat evacuation at the condenser, since the air flow generated by the moto-ventilator group is an increasing function of vehicle speed. The second effect is due to the important convection heat exchange between the vehicle's walls and the exterior at high v veh . For the same cabin temperature, the wall temperature is lower, which induce less power consumption.
This vehicle speed effect, has consequences on thermal comfort management. For a given energy, the discomfort decreases with vehicle speed. This is shown in Figure 13a, where the difference between optimal global equivalent temperatures in a INRETS A2 and UL1 driving cycle, is about 0.75 • C, for the same consumed energy. Another representation of this effect is shown in Figure 13b. The mean local equivalent temperatures lies in hot comfortable zone for an INRETS UF1 and UL1 driving cycle. However, almost all T eq,i are close to the set point equivalent temperatures (comfortable zone), for an INRETS A2 for the same consumed energy of 1.1 kWh.
The same conclusions can be formulated from the results depicted in Figure 13c,d, where the driver is dressed in a formal way. However, we notice also that 1.1 kWh is not sufficient to achieve a good thermal comfort for a formal clothing, due to a higher thermal insulation. As consequence, the driver might feel discomfort in a INRETS A2 driving cycle, especially in lower body parts (Figure 13d).
Let us now do the same analysis for a given driving cycle, and different outside temperature profiles. Figure 14a shows the Pareto frontiers for three outside temperature profiles characterized by T ext,0 [ • C] ∈ {32 ; 35 ; 38}, for an INRETS R2 driving cycle and a driver dressed with an informal clothing. For an available energy of 1.1 kWh, the optimal global equivalent temperature T eq,rms is deduced from the Pareto frontiers. The corresponding mean local equivalent temperatures T eq,i , are reported in Figure 14b in the comfort zone diagram. The same results are reported in Figure 14c,d for a driver dressed with a formal clothing. In this harsh scenario, all , lie outside the comfort zone, which means that the driver experiences an uncomfortable hot thermal sensation in the last 15 min of the trip. The same conclusions can be formulated in the case of the driver dressed in a formal way (Figures 14c,d).
As stated before, the algorithm takes actions on the HVAC system in order to choose the optimal thermal comfort for a given available energy. Hence, it is worth analyzing these actions and their We notice that the outside temperature has a strong influence on thermal comfort, larger than the vehicle speed. For an available energy of 1.1 kWh, a good thermal comfort is achieved for an outside temperature of T ext,0 = 29 • C. In this situation, the optimal global equivalent temperature is about 28 • C and the corresponding local equivalent temperatures, are presented in Figure 14b, all lie in the comfortable green zone. However, in hotter weather, T eq,rms increases significantly and it reaches 35 • C for T ext,0 = 38 • C for the same available energy.
In this harsh scenario, all T eq,i lie outside the comfort zone, which means that the driver experiences an uncomfortable hot thermal sensation in the last 15 min of the trip. The same conclusions can be formulated in the case of the driver dressed in a formal way (Figure 14c,d).
As stated before, the algorithm takes actions on the HVAC system in order to choose the optimal thermal comfort for a given available energy. Hence, it is worth analyzing these actions and their consequences on the cabin physical quantities and the driver's thermo-physiological responses. Figure 15 shows the temporal evolution of cabin temperature and relative humidity that corresponds to colored squares plotted in Figure 14c. The yellow cabin temperature profile has the ideal shape since it leads to an ideal comfort. First, the cabin temperature drops at the beginning of the trip, when the hot cabin needs to be cooled down. This phase cannot be reduced below a certain time due to thermal inertia and the maximum power of the HVAC system. The temperature decreases afterwards until reaching 22 • C, which is the comfort temperature at . q sun = 1000 W/m 2 ( Table 3). For an outside temperature T ext,0 = 32 • C, the available energy is not sufficient to maintain ideal comfort over the whole trip. In this situation, after the cooling phase, the cabin temperature is maintained around 23.7 • C for 20 min, then it increases slowly until reaching 26 • C at the end of the trip. For even higher T ext,0 , the temperature at final time can reach 34 • C. For instance, if the available energy cannot ensure the ideal comfort, the compressor speed is set at the minimum regime before the end of the trip (red and purple curve). In this situation, the air flow and the heating ratio are slightly increased to adjust the thermal comfort.
The driver thermo-physiological response is represented by four thermoregulation mechanisms:  Figure 15 shows also the temporal evolution of two HVAC control variables: the compressor rotation speed (Figure 15c) and the air flow (Figure 15d). The compressor speed is set at its maximum value during the cooling phase, then it is adjusted according to the available energy, in order to provide the "reachable" thermal comfort.
For instance, if the available energy cannot ensure the ideal comfort, the compressor speed is set at the minimum regime before the end of the trip (red and purple curve). In this situation, the air flow and the heating ratio are slightly increased to adjust the thermal comfort.
The driver thermo-physiological response is represented by four thermoregulation mechanisms: sweating, shivering, vasodilatation and vasoconstriction. Figure 16 shows the temporal evolution of these mechanisms, which correspond to the colored squares in Figure 14c. Firstly, we observe in Figure 16b that shivering is null for all outside temperatures. This result is expected since simulations are performed exclusively in hot weather. In such conditions, the human body regulates its internal temperature through sweating and vasodilatation. The first mechanism allows to evacuate metabolism heat production, through evaporation. The second one, increases the diameter of veins to improve the heat exchange between the skin and the ambient environment. In a comfortable state, both mechanism should be disabled, and the body is in thermal equilibrium. Figure 16a,b illustrate these phenomena, especially when the outside temperature is T ext,0 = 29 • C and that all local equivalent temperatures lie in the comfort zone ( Figure 14d). This behavior of sweating and vasodilatation mechanism is also noticed at some extent for T ext,0 = 32 • C. Although the sweating mechanism is negligible, the vasodilation is still activated in order to regulate internal temperature. This could be explained by the fact that the local equivalent temperatures of lower body parts, lies in "hot but comfortable zone" (red line in Figure 14d). For higher outside temperature, the temporal evolution of both thermo-physiological mechanism is similar to the cabin temperature's one.

Test Case 2: Vehicle Speed vs. Thermal Comfort Trade-off
In the second test case, the driving cycle is modified, so that the driver travels the same distance at a lower speed. Less energy is required for the traction, but more is needed for the ideal thermal

Test Case 2: Vehicle Speed vs. Thermal Comfort Trade-off
In the second test case, the driving cycle is modified, so that the driver travels the same distance at a lower speed. Less energy is required for the traction, but more is needed for the ideal thermal comfort, since the trip lasts longer. The idea is to evaluate which effect is dominating, and if there is a possible compromise between the traction and the HVAC system consumption.
A scaling factor k < 1 is applied to the speed as follows: v new,cycle (t) = v cycle (t).k (40) By decreasing the values of k, one reduces the energy needed for traction, but one also increases the HVAC system energy needed for the ideal thermal comfort. The total energy may be an increasing or a decreasing function of k, or a function that have a minimum. Simulations have been conducted on a combination of six INRETS driving cycles: A1, R3, R2, R1, UF3 and UF1, and the following range of outside temperature, relative humidity and solar irradiance: , 50, 80} and . q sun W/m 2 ∈ {0, 500, 1000}. We assume that the driver is dressed in an informal way and for sake of simplicity, that the weather conditions are not time varying.
First, we will analyze the results for a specific weather condition, then an evaluation of potential energy gains will be assessed for all traffic and weather scenarios. Figure 17a shows the traction energy, the HVAC system energy for ideal comfort, and their sum, as a function of k for the following scenario: INRETS R2 driving cycle, T ext = 36 • C, RH ext = 50% and . q sun = 1000 W/m 2 . We notice that the total energy presents a minimum for = 0.75. The corresponding energy is 9 kWh, which represents a 0.45 kWh economy if the driver slows down by 25%. If the driver slows down only by 10% (i.e., = 0.9), the energy saving is about 0.35 kWh. Figure 17b shows the same type of results for the INRETS UF3 driving cycle. For this specific scenario, there is no substantial energy gain if the driver slows down, since the total energy is an increasing function of . Figure 18 shows the total energy as a function of factor, for the six driving cycles, and the aforementioned weather scenario. The obtained results indicate that for INRETS A2, R3 and R2, the dominant consumption is associated to the traction, and reducing the speed decreases the total consumption. On the contrary, for urban trips, the gain is not obvious, as a decrease of traction energy can be compensated by an increase of HVAC system consumption. We notice that the total energy presents a minimum for k = 0.75. The corresponding energy is 9 kWh, which represents a 0.45 kWh economy if the driver slows down by 25%. If the driver slows down only by 10% (i.e., k = 0.9), the energy saving is about 0.35 kWh. Figure 17b shows the same type of results for the INRETS UF3 driving cycle. For this specific scenario, there is no substantial energy gain if the driver slows down, since the total energy is an increasing function of k. Figure 18 shows the total energy as a function of k factor, for the six driving cycles, and the aforementioned weather scenario. The obtained results indicate that for INRETS A2, R3 and R2, the dominant consumption is associated to the traction, and reducing the speed decreases the total consumption. On the contrary, for urban trips, the gain is not obvious, as a decrease of traction energy can be compensated by an increase of HVAC system consumption. energy gain if the driver slows down, since the total energy is an increasing function of . Figure 18 shows the total energy as a function of factor, for the six driving cycles, and the aforementioned weather scenario. The obtained results indicate that for INRETS A2, R3 and R2, the dominant consumption is associated to the traction, and reducing the speed decreases the total consumption. On the contrary, for urban trips, the gain is not obvious, as a decrease of traction energy can be compensated by an increase of HVAC system consumption.    Figure 19 reports the energy gain for k = 0.9 as a function of

Conclusions and Perspectives
This paper introduced an optimal control approach, using dynamic programming principles, for thermal comfort and driving range management in BEVs. The passenger's thermal comfort was modeled through the equivalent temperature index, estimated thanks to a model of human thermophysiology. Simulations were conducted for different traffic and meteorological conditions. The main results are the following ones.
Using the equivalent temperature index and a detailed model of the HVAC enables to account

Conclusions and Perspectives
This paper introduced an optimal control approach, using dynamic programming principles, for thermal comfort and driving range management in BEVs. The passenger's thermal comfort was modeled through the equivalent temperature index, estimated thanks to a model of human thermo-physiology. Simulations were conducted for different traffic and meteorological conditions. The main results are the following ones.
Using the equivalent temperature index and a detailed model of the HVAC enables to account for the influence of parameters such as the wall temperature, or the vehicle speed, on the comfort temperature and humidity in the cabin. Adjusting the four HVAC control variables enables to reach a given thermal comfort in the cabin at the lowest energetic cost.
A large number of simulation was performed and the energetic cost of thermal comfort in various conditions was quantified. Its strong impact in urban traffic conditions was underlined.
The results showed the usefulness of using DP for thermal comfort management. Indeed, when there is enough energy for ideal comfort, the algorithm drives all local equivalent temperatures to the green zone of the comfort diagram. The thermoregulation mechanism are disabled, which reflects a "comfortable" state. In contrast, when the available energy is low, the DP algorithm decreases the thermal comfort to meet the energy constraint, by putting the driver in a "hot but comfortable" state.
Furthermore, energy savings can be realized in some traffic conditions by decreasing the vehicle speed. For instance, when the driver slows down by 10% in highway, at night conditions, energy gains can reach 2.2 kWh for a one-hour trip.
Future work will focus on developing a real-time control for thermal comfort management in electric vehicles, based on the results of the present work. On the other hand, more research is needed in order to apply the algorithm in cold climate scenarios.   Q f low , due to the air blown by the ventilation system, is quantified using the enthalpy difference between the ventilation outlet air h out and the cabin air h cab . Both enthalpies are functions of air temperature and relative humidity, and can be computed through tabulated data of air thermodynamic properties. The heat flow rate . Q f low is expressed by (A1): . Q w represents the heat flow rate exchanged with the cabin walls windows and windshield, with a constant heat transfer coefficient h c,w . We assume that the heat transfer is done only by convection. The term . Q w is expressed by (A2): . Q ext is the convective heat follow rate exchanged with the outside through vehicle body leakages. The heat transfer coefficient h c,ext is a logarithmic function of vehicle speed v veh . The heat flow rate . Q ext is expressed by (A3): The evaporator model is based on NTU method [42]. This method consist to model the heat flow between the air and the refrigerant as a function of the maximal heat flow between the two fluids, and an efficiency factor η. The maximal heat flow occurs at the inlet of the evaporator when temperature difference is maximal. For a dry air, the maximal heat flow . Q evap,d,max is expressed by (A4) For a humid air, one must take its relative humidity into account. In this situation, the maximal heat flow whereĥ air refers to enthalpy functions, which are extracted from the tabulated data of the air [45]. The actual heat flow is given by (A6) and (A7), where η evap,d and η evap,w are the efficiencies of the evaporator [46] when the air is respectively dry and humid: . Q evap,d = .
Q evap,d,max η evap,d Computing the actual heat flow . Q e at the evaporator allows next to determine the temperature and pressure of the refrigerant fluid, and more importantly, the thermodynamic properties of the outlet air. By combining (A8) and (9), one can resolve the heat balance Equation (7), and calculate the evaporating temperature of the refrigerant. Since it is in thermodynamic equilibrium, its pressure depends only on T 7 : P 7 =P re f (T 7 ) (A9) whereĥ air refers to pressure functions built from the tabulated thermodynamic properties of the refrigerant R1234yf [41,47]. For simplicity sakes, we assume that refrigerant pressure at the outlet of the evaporator is equal to its saturated pressure: The temperature T 1 of the refrigerant at the outlet of the evaporator is higher than its saturation temperature T 7 , since it continues to exchange the heat flow with the air. The temperature difference is called the superheat T sh : Appendix A.1.3. Air Thermal Properties at AC Outlet The actual heat flow . Q e exchanged at the evaporator serves also to evaluate the thermodynamic properties of the outlet air. Its enthalpy h AC is governed by the heat balance (A12) where h in is the enthalpy of the inlet air. This quantity depends on the inlet temperature T in and relative humidity RH in : The thermal properties of the air at the outlet of the air conditioning system, i.e., the outlet temperature T AC,out and the relative humidity RH AC,out depends on the presence of steam condensation. The quantities T AC,out and RH AC,out are functions of h AC,out and can be computed through the tabulated data of the air.

.4. Expansion Valve
At the inlet of the expansion valve, the refrigerant is in the liquid state, with a temperature T 5 . The relationship (A13) allows to determine the enthalpy of the refrigerant as a function of its temperature T 5 and the high pressure P 4 : h 5 =ĥ re f (T 5 , P 4 ) (A13) The model of the expansion valve neglects heat losses, i.e., constant enthalpy is assumed throughout the expansion valve: h 6 = h 5 (A14) Appendix A.1.5. Compressor The actual refrigerant mass flow rate is given by (A15): .
where V d refers to the refrigerant volume trapped in the cylinder per revolution (assumed to be constant) and ρ 1 ,its specific mass at the compressor inlet [48]. The volumetric efficiency η vol is a function of the ratio of the inlet and outlet pressure P 1 and P 2 , the rotational speed N comp , and the refrigerant trapped in the volume V d . The refrigerant specific enthalpy h 2 at the outlet of the compressor is given by (A16): where h 1 is the enthalpy of the refrigerant at the pressure P 1 and the temperature T 1 , given by (A17): h 1 =ĥ re f (P 1 , T 1 ) (A17) The difference h 2,is − h 1 is equivalent to the energy consumed if the thermodynamic transformation of the refrigerant was isentropic. The isentropic enthalpy h 2,is at outlet of the compressor is given by (A18): h 2,is =ĥ re f (P 2 , s 1 ) (A18) The quantity η is is the isentropic efficiency of the compressor, modeled by an affine function of . m re f , N comp and T 1 .
Appendix A.1.6. Condenser The condenser is modeled by using the same NTU method used to express the actual heat flow in the evaporator. The maximal heat flow in the condenser, given by (A19), is a function of temperature difference between the outside air and the refrigerant at the inlet of the heat exchanger: .
where . m air,MVG is the air mass flow rate propelled by the moto-ventilator group (f) (Figure 3). c p,air is the specific heat capacity of the air (J/ • K.kg).
The actual heat flow is deduced by multiplying . Q c,max by the condenser efficiency η cond : . Q cond = .
Q cond,max η cond (A20) Computing the actual heat flow at the condenser allows next to determine the temperature and pressure of the refrigerant fluid. By combining (A20) and (9) one can resolve the heat balance Equation (8), and then calculate the condensing temperature of the refrigerant T 4 . Since it is in thermodynamic equilibrium, its pressure depends only on T 4 : For simplicity, we assume that the refrigerant pressure at the outlet and inlet of the condenser is equal to its saturated pressure: P 2 = P 5 = P 4 (A22) The saturation of the refrigerant usually occurs before it exits the condenser. It continues to evacuate the heat to the outside air until reaching the point 5 in the refrigeration cycle. The temperature difference is called the subcooling, and allows to determine the refrigerant temperature at the exit of the condenser by (A23):

. Heating System Modeling
After being cooled by the AC system, a fraction of air α is heated by the electrical machine water cooling system. The heating system is also modeled by the − NTU method, where the heat flow rate . Q heat exchanged with the cooled air is expressed by the following Equation (A24): . Q heat = . m water c p,water (T water − T AC,out )η heat (A24) where . m water is the mass flow rate of the cooling water, c p,water is the specific heat capacity of the water, and T water is the water temperature.
The actual heat flow . Q heat allows to evaluate the thermodynamic properties of the heated air. Its temperature T heat at the outlet of the heating system, is governed by the heat balance (A25): Since the cooled air is mixed with the heated one, the temperature of the air at outlet of the HVAC system T out , is weighted by the ratio α. The heated air passes through the ventilation duct before it is blown into the cabin: Finally, the relative humidity RH out is computed using the tabulated data of the air.