Eclipse SUMO - Simulation of Urban MObility
HelpersMMPEVEM.cpp
Go to the documentation of this file.
1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2002-2022 German Aerospace Center (DLR) and others.
4 // This program and the accompanying materials are made available under the
5 // terms of the Eclipse Public License 2.0 which is available at
6 // https://www.eclipse.org/legal/epl-2.0/
7 // This Source Code may also be made available under the following Secondary
8 // Licenses when the conditions for such availability set forth in the Eclipse
9 // Public License 2.0 are satisfied: GNU General Public License, version 2
10 // or later which is available at
11 // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
12 // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
13 /****************************************************************************/
18 // The MMP's emission model for electric vehicles.
19 // If you use this model for academic research, you are highly encouraged to
20 // cite our paper "Accurate physics-based modeling of electric vehicle energy
21 // consumption in the SUMO traffic microsimulator"
22 // (DOI: 10.1109/ITSC48978.2021.9564463).
23 // Teaching and Research Area Mechatronics in Mobile Propulsion (MMP), RWTH Aachen
24 /****************************************************************************/
25 
26 
28 
29 #ifdef _MSC_VER
30 #define _USE_MATH_DEFINES
31 #endif
32 #include <cmath>
33 
35 #include <utils/common/SUMOTime.h>
37 #include <utils/geom/GeomHelper.h>
38 
39 
40 
41 
42 // Constants
43 const double EPS = 1e-6;
44 const double G = 9.81; // Gravitational acceleration [m/s^2]
45 const double RHO_AIR = 1.204; // Air density [kg/m^3]
46 
47 
48 
49 
85 bool calcPowerConsumption(double m, double r_wheel, double Theta, double c_rr,
86  double c_d, double A_front, double i_gear, double eta_gear, double M_max,
87  double P_max, double M_recup_max, double P_recup_max, double R_battery,
88  double U_battery_0, double P_const,
89  const CharacteristicMap& ref_powerLossMap, double dt, double v, double a,
90  double alpha, double& ref_powerConsumption) {
91  bool b_stateValid = true;
92 
93  // Mass factor [1]
94  const double e_i = 1.0 + Theta / (m * r_wheel * r_wheel);
95  // Average speed during the previous timestep
96  const double v_mean = v - 0.5 * a * dt;
97 
98  // Force required for the desired acceleration [N]
99  double F_a = m * a * e_i;
100  // Grade resistance [N]
101  double F_gr = m * G * std::sin(DEG2RAD(alpha));
102  // Rolling resistance [N]
103  double F_rr = m * G * std::cos(DEG2RAD(alpha)) * c_rr;
104  if (std::abs(v_mean) <= EPS) {
105  F_rr = 0;
106  }
107  // Drag [N]
108  double F_d = 0.5 * c_d * A_front * RHO_AIR * v_mean * v_mean;
109  // Tractive force [N]
110  const double F_tractive = F_a + F_gr + F_rr + F_d;
111 
112  // Speed of the motor [rpm]
113  const double n_motor = v_mean / (2 * M_PI * r_wheel) * 60 * i_gear;
114  // Angular velocity of the motor [1/s]
115  double omega_motor = 2 * M_PI * n_motor / 60;
116  if (omega_motor == 0) {
117  omega_motor = EPS;
118  }
119 
120  // Torque at the motor [Nm]. Please note that this model, like most real EVs,
121  // utilizes the EM to hold the vehicle on an incline rather than the brakes.
122  double M_motor = F_tractive * r_wheel / i_gear;
123  if (F_tractive < 0) {
124  M_motor *= eta_gear;
125  } else {
126  M_motor /= eta_gear;
127  }
128  // Engine power [W]
129  double P_motor = M_motor * omega_motor;
130 
131  // Cap torque or power if necessary
132  // Accelerating
133  if (M_motor >= 0) {
134  if (M_motor > M_max) {
135  M_motor = M_max;
136  P_motor = M_motor * omega_motor;
137  b_stateValid = false;
138  }
139  if (P_motor > P_max) {
140  P_motor = P_max;
141  M_motor = P_motor / omega_motor;
142  b_stateValid = false;
143  }
144  }
145  // Recuperating
146  else {
147  // Even when capping, the state is still valid here because the extra energy
148  // is assumed to go to the brakes
149  if (M_motor < -M_recup_max) {
150  M_motor = -M_recup_max;
151  P_motor = M_motor * omega_motor;
152  }
153  if (P_motor < -P_recup_max) {
154  P_motor = -P_recup_max;
155  M_motor = P_motor / omega_motor;
156  }
157  }
158 
159  // Power loss in the electric motor + inverter [W]
160  double P_loss_motor =
161  ref_powerLossMap.eval(std::vector<double> {n_motor, M_motor})[0];
162  if (std::isnan(P_loss_motor)) {
163  P_loss_motor = 0.0;
164  b_stateValid = false;
165  }
166 
167  // Power demand at the battery [W]
168  double P_battery = P_motor + P_loss_motor + P_const;
169  // Total power demand (including the power loss in the battery) [W]
170  double P_total = (U_battery_0 * U_battery_0) / (2.0 * R_battery)
171  - U_battery_0 * std::sqrt((U_battery_0 * U_battery_0
172  - 4.0 * R_battery * P_battery) / (4.0 * R_battery * R_battery));
173 
174  ref_powerConsumption = P_total;
175  return b_stateValid;
176 }
177 
178 
179 
180 
185  : PollutantsInterface::Helper("MMPEVEM", MMPEVEM_BASE, MMPEVEM_BASE + 1)
186 { }
187 
188 
189 
190 
191 
209  const PollutantsInterface::EmissionType e, const double v, const double a,
210  const double slope, const EnergyParams* ptr_energyParams) const {
211  if (e != PollutantsInterface::ELEC) {
212  return 0.0;
213  }
214 
215  // Extract all required parameters
216  // Vehicle mass
217  const double m = ptr_energyParams->getDouble(SUMO_ATTR_VEHICLEMASS);
218  // Wheel radius
219  const double r_wheel = ptr_energyParams->getDouble(SUMO_ATTR_WHEELRADIUS);
220  // Internal moment of inertia
221  const double Theta
222  = ptr_energyParams->getDouble(SUMO_ATTR_INTERNALMOMENTOFINERTIA);
223  // Rolling resistance coefficient
224  const double c_rr
225  = ptr_energyParams->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT);
226  // Air drag coefficient
227  const double c_d = ptr_energyParams->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT);
228  // Cross-sectional area of the front of the car
229  const double A_front
230  = ptr_energyParams->getDouble(SUMO_ATTR_FRONTSURFACEAREA);
231  // Gear ratio
232  const double i_gear = ptr_energyParams->getDouble(SUMO_ATTR_GEARRATIO);
233  // Gearbox efficiency
234  const double eta_gear = ptr_energyParams->getDouble(SUMO_ATTR_GEAREFFICIENCY);
235  // Maximum torque
236  const double M_max = ptr_energyParams->getDouble(SUMO_ATTR_MAXIMUMTORQUE);
237  // Maximum power
238  const double P_max = ptr_energyParams->getDouble(SUMO_ATTR_MAXIMUMPOWER);
239  // Maximum recuperation torque
240  const double M_recup_max
241  = ptr_energyParams->getDouble(SUMO_ATTR_MAXIMUMRECUPERATIONTORQUE);
242  // Maximum recuperation power
243  const double P_recup_max
244  = ptr_energyParams->getDouble(SUMO_ATTR_MAXIMUMRECUPERATIONPOWER);
245  // Internal battery resistance
246  const double R_battery
247  = ptr_energyParams->getDouble(SUMO_ATTR_INTERNALBATTERYRESISTANCE);
248  // Nominal battery voltage
249  const double U_battery_0
250  = ptr_energyParams->getDouble(SUMO_ATTR_NOMINALBATTERYVOLTAGE);
251  // Constant power consumption
252  const double P_const
253  = ptr_energyParams->getDouble(SUMO_ATTR_CONSTANTPOWERINTAKE);
254  // Power loss map
255  const CharacteristicMap& ref_powerLossMap
256  = ptr_energyParams->getCharacteristicMap(SUMO_ATTR_POWERLOSSMAP);
257 
258  double P = 0.0; // [W]
259  bool b_stateValid = calcPowerConsumption(m, r_wheel, Theta, c_rr, c_d,
260  A_front, i_gear, eta_gear, M_max, P_max, M_recup_max, P_recup_max,
261  R_battery, U_battery_0, P_const, ref_powerLossMap, TS, v, a, slope, P);
262  P /= 3600.0; // [Wh/s]
263 
264  if (b_stateValid) {
265  return P;
266  } else {
267  return std::nan("");
268  }
269 }
270 
#define DEG2RAD(x)
Definition: GeomHelper.h:35
bool calcPowerConsumption(double m, double r_wheel, double Theta, double c_rr, double c_d, double A_front, double i_gear, double eta_gear, double M_max, double P_max, double M_recup_max, double P_recup_max, double R_battery, double U_battery_0, double P_const, const CharacteristicMap &ref_powerLossMap, double dt, double v, double a, double alpha, double &ref_powerConsumption)
Compute the power consumption of an EV powertrain in a certain state.
const double RHO_AIR
const double G
const double EPS
#define TS
Definition: SUMOTime.h:40
int SUMOEmissionClass
@ SUMO_ATTR_GEAREFFICIENCY
Gear efficiency.
@ SUMO_ATTR_MAXIMUMPOWER
Maximum Power.
@ SUMO_ATTR_INTERNALBATTERYRESISTANCE
Internal battery resistance.
@ SUMO_ATTR_MAXIMUMTORQUE
Maximum torque.
@ SUMO_ATTR_ROLLDRAGCOEFFICIENT
Roll Drag coefficient.
@ SUMO_ATTR_CONSTANTPOWERINTAKE
Constant Power Intake.
@ SUMO_ATTR_WHEELRADIUS
@ SUMO_ATTR_AIRDRAGCOEFFICIENT
Air drag coefficient.
@ SUMO_ATTR_POWERLOSSMAP
A string encoding the power loss map.
@ SUMO_ATTR_MAXIMUMRECUPERATIONPOWER
Maximum recuperation power.
@ SUMO_ATTR_MAXIMUMRECUPERATIONTORQUE
Maximum recuperation torque.
@ SUMO_ATTR_VEHICLEMASS
Vehicle mass.
@ SUMO_ATTR_GEARRATIO
Gear ratio.
@ SUMO_ATTR_INTERNALMOMENTOFINERTIA
Internal moment of inertia.
@ SUMO_ATTR_NOMINALBATTERYVOLTAGE
Nominal battery voltage.
@ SUMO_ATTR_FRONTSURFACEAREA
Front surface area.
The purpose of this class is to store a characteristic map (German: Kennfeld) of arbitrary dimensions...
std::vector< double > eval(const std::vector< double > &ref_p, double eps=1e-6) const
Evaluate a point in the map using linear interpolation.
An upper class for objects with additional parameters.
Definition: EnergyParams.h:41
double getDouble(SumoXMLAttr attr) const
const CharacteristicMap & getCharacteristicMap(SumoXMLAttr attr) const
Return the CharacteristicMap that belongs to a given attribute.
HelpersMMPEVEM()
Constructor.
double compute(const SUMOEmissionClass, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams *ptr_energyParams) const
Compute the amount of emitted pollutants for an emission class in a given state.
Helper methods for PHEMlight-based emission computation.
EmissionType
Enumerating all emission types, including fuel.
#define M_PI
Definition: odrSpiral.cpp:40