Eclipse SUMO - Simulation of Urban MObility
HelpersEnergy.cpp
Go to the documentation of this file.
1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2001-2020 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 /****************************************************************************/
20 // Helper methods for HBEFA-based emission computation
21 /****************************************************************************/
22 #include <config.h>
23 
24 #include <utils/common/SUMOTime.h>
25 #include <utils/common/ToString.h>
27 #include "HelpersEnergy.h"
28 
29 //due to WRITE_WARNING
31 
32 //due to sqrt of complex
33 #include <complex>
34 
35 
36 // ===========================================================================
37 // method definitions
38 // ===========================================================================
39 HelpersEnergy::HelpersEnergy() : PollutantsInterface::Helper("Energy", ENERGY_BASE) {
40  // default values from
41  // Kurczveil, T., López, P.Á., & Schnieder, E. (2014). Implementation of an Energy Model and a Charging Infrastructure in SUMO.
53 }
54 
55 
56 double
57 HelpersEnergy::compute(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const std::map<int, double>* param) const {
58  if (e != PollutantsInterface::ELEC) {
59  return 0.;
60  }
61  if (param == nullptr) {
62  param = &myDefaultParameter;
63  }
64  //@ToDo: All formulas below work with the logic of the euler update (refs #860).
65  // Approximation order could be improved. Refs. #2592.
66 
67  const double lastV = v - ACCEL2SPEED(a);
68  const double mass = param->find(SUMO_ATTR_VEHICLEMASS)->second;
69 
70  // calculate power needed for potential energy difference
71  double power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
72 
73  // power needed for kinetic energy difference of vehicle
74  power += 0.5 * mass * (v * v - lastV * lastV) / TS;
75 
76  // add power needed for rotational energy diff of internal rotating elements
77  power += 0.5 * param->find(SUMO_ATTR_INTERNALMOMENTOFINERTIA)->second * (v * v - lastV * lastV) / TS;
78 
79  // power needed for Energy loss through Air resistance [Ws]
80  // Calculate energy losses:
81  // EnergyLoss,Air = 1/2 * rho_air [kg/m^3] * myFrontSurfaceArea [m^2] * myAirDragCoefficient [-] * v_Veh^2 [m/s] * s [m]
82  // ... with rho_air [kg/m^3] = 1,2041 kg/m^3 (at T = 20C)
83  // ... with s [m] = v_Veh [m/s] * TS [s]
84  // divided by TS to get power instead of energy
85  power += 0.5 * 1.2041 * param->find(SUMO_ATTR_FRONTSURFACEAREA)->second * param->find(SUMO_ATTR_AIRDRAGCOEFFICIENT)->second * v * v * v;
86 
87  // power needed for Energy loss through Roll resistance [Ws]
88  // ... (fabs(veh.getSpeed())>=0.01) = 0, if vehicle isn't moving
89  // EnergyLoss,Tire = c_R [-] * F_N [N] * s [m]
90  // ... with c_R = ~0.012 (car tire on asphalt)
91  // ... with F_N [N] = myMass [kg] * g [m/s^2]
92  // divided by TS to get power instead of energy
93  power += param->find(SUMO_ATTR_ROLLDRAGCOEFFICIENT)->second * GRAVITY * mass * v;
94 
95  // Energy loss through friction by radial force [Ws]
96  // If angle of vehicle was changed
97  const double angleDiff = param->find(SUMO_ATTR_ANGLE)->second;
98  if (angleDiff != 0.) {
99  // Compute new radio
100  double radius = SPEED2DIST(v) / fabs(angleDiff);
101 
102  // Check if radius is in the interval [0.0001 - 10000] (To avoid overflow and division by zero)
103  if (radius < 0.0001) {
104  radius = 0.0001;
105  } else if (radius > 10000) {
106  radius = 10000;
107  }
108  // EnergyLoss,internalFrictionRadialForce = c [m] * F_rad [N];
109  // Energy loss through friction by radial force [Ws]
110  // divided by TS to get power instead of energy
111  power += param->find(SUMO_ATTR_RADIALDRAGCOEFFICIENT)->second * mass * v * v / radius * v;
112  }
113 
114  // EnergyLoss,constantConsumers
115  // Energy loss through constant loads (e.g. A/C) [W]
116  power += param->find(SUMO_ATTR_CONSTANTPOWERINTAKE)->second;
117 
118  // E_Bat = E_kin_pot + EnergyLoss;
119  if (power > 0) {
120  // Assumption: Efficiency of myPropulsionEfficiency when accelerating
121  power /= param->find(SUMO_ATTR_PROPULSIONEFFICIENCY)->second;
122  } else {
123  // Assumption: Efficiency of myRecuperationEfficiency when recuperating
124  power *= param->find(SUMO_ATTR_RECUPERATIONEFFICIENCY)->second;
125  if (a != 0) {
126  // Fiori, Chiara & Ahn, Kyoungho & Rakha, Hesham. (2016).
127  // Power-based electric vehicle energy consumption model: Model
128  // development and validation. Applied Energy. 168. 257-268.
129  // 10.1016/j.apenergy.2016.01.097.
130  //
131  // Insaf Sagaama, Amine Kchiche, Wassim Trojet, Farouk Kamoun
132  // Improving The Accuracy of The Energy Consumption Model for
133  // Electric Vehicle in SUMO Considering The Ambient Temperature
134  // Effects
135  power *= (1 / exp(param->find(SUMO_ATTR_RECUPERATIONEFFICIENCY_BY_DECELERATION)->second / fabs(a)));
136  }
137  }
138 
139  // convert from [W], to [Wh/s] (3600s / 1h):
140  return power / 3600.;
141 }
142 
143 
144 double
145 HelpersEnergy::acceleration(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const std::map<int, double>* param) const {
146  if (e != PollutantsInterface::ELEC) {
147  return 0.;
148  }
149 
150  if (param == nullptr) {
151  param = &myDefaultParameter;
152  }
153 
154  // Inverse formula for the function compute()
155  // Computes the achievable acceleration using the given speed and amount of consumed electric power
156  // It does not consider loss through friction by radial force and the recuperation efficiency dependent on |a| (eta_reuperation is considered constant)
157  // Prest = const1*a + const2*a^2 + const3*a^3
158 
159  const double mass = param->find(SUMO_ATTR_VEHICLEMASS)->second;
160  double const1, const2, const3;
161  double Prest;
162  int numX;
163  double x1, x2, x3;
164 
165  if (P > 0) {
166  // Assumption: Efficiency of myPropulsionEfficiency when accelerating
167  Prest = P * 3600 * param->find(SUMO_ATTR_PROPULSIONEFFICIENCY)->second;
168  } else {
169  // Assumption: Efficiency of myRecuperationEfficiency when recuperating
170  Prest = P * 3600 / param->find(SUMO_ATTR_RECUPERATIONEFFICIENCY)->second;
171  }
172 
173  // calculate power drop due to a potential energy difference
174  Prest -= mass * GRAVITY * sin(DEG2RAD(slope)) * (v);
175  const1 = mass * GRAVITY * sin(DEG2RAD(slope)) * (TS);
176 
177  // Power loss through Roll resistance [W]
178  // ... (fabs(veh.getSpeed())>=0.01) = 0, if vehicle isn't moving
179  // EnergyLoss,Tire = c_R [-] * F_N [N] * s [m]
180  // ... with c_R = ~0.012 (car tire on asphalt)
181  // ... with F_N [N] = myMass [kg] * g [m/s^2]
182  Prest -= param->find(SUMO_ATTR_ROLLDRAGCOEFFICIENT)->second * GRAVITY * mass * v;
183  const1 += param->find(SUMO_ATTR_ROLLDRAGCOEFFICIENT)->second * GRAVITY * mass * (TS);
184 
185  //Constant loads are omitted. We assume P as the max limit for the main traction drive. Constant loads are often covered by an auxiliary drive
186  //Power loss through constant loads (e.g. A/C) [W]
187  //Prest -= param->find(SUMO_ATTR_CONSTANTPOWERINTAKE)->second / TS;
188 
189  // kinetic energy difference of vehicle
190  const1 += 0.5 * mass * (2 * v);
191  const2 = 0.5 * mass * (TS);
192 
193  // add rotational energy diff of internal rotating elements
194  const1 += param->find(SUMO_ATTR_INTERNALMOMENTOFINERTIA)->second * (2 * v);
195  const2 += param->find(SUMO_ATTR_INTERNALMOMENTOFINERTIA)->second * (TS);
196 
197  // Energy loss through Air resistance [Ws]
198  // Calculate energy losses:
199  // EnergyLoss,Air = 1/2 * rho_air [kg/m^3] * myFrontSurfaceArea [m^2] * myAirDragCoefficient [-] * v_Veh^2 [m/s] * s [m]
200  // ... with rho_air [kg/m^3] = 1,2041 kg/m^3 (at T = 20C)
201  // ... with s [m] = v_Veh [m/s] * TS [s]
202  Prest -= 0.5 * 1.2041 * param->find(SUMO_ATTR_FRONTSURFACEAREA)->second * param->find(SUMO_ATTR_AIRDRAGCOEFFICIENT)->second * (v * v * v);
203  const1 += 0.5 * 1.2041 * param->find(SUMO_ATTR_FRONTSURFACEAREA)->second * param->find(SUMO_ATTR_AIRDRAGCOEFFICIENT)->second * (3 * v * v * TS);
204  const2 += 0.5 * 1.2041 * param->find(SUMO_ATTR_FRONTSURFACEAREA)->second * param->find(SUMO_ATTR_AIRDRAGCOEFFICIENT)->second * (3 * v * TS * TS);
205  const3 = 0.5 * 1.2041 * param->find(SUMO_ATTR_FRONTSURFACEAREA)->second * param->find(SUMO_ATTR_AIRDRAGCOEFFICIENT)->second * (TS * TS * TS);
206 
207  // Prest = const1*a + const2*a^2 + const3*a^3
208  // solve cubic equation in a
209 
210  std::tie(numX, x1, x2, x3) = PolySolver::cubicSolve(const3, const2, const1, -Prest);
211 
212 
213  switch (numX) {
214  case 1:
215  return x1;
216  case 2:
217  return MAX2(x1, x2);
218  case 3:
219  return MAX3(x1, x2, x3);
220  default:
221  WRITE_ERROR("An acceleration given by the power was not found.");
222  return 0;
223  }
224 
225 }
226 
227 
228 /****************************************************************************/
#define DEG2RAD(x)
Definition: GeomHelper.h:35
#define GRAVITY
Definition: GeomHelper.h:37
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:284
#define SPEED2DIST(x)
Definition: SUMOTime.h:43
#define ACCEL2SPEED(x)
Definition: SUMOTime.h:49
#define TS
Definition: SUMOTime.h:40
int SUMOEmissionClass
@ SUMO_ATTR_ROLLDRAGCOEFFICIENT
Roll Drag coefficient.
@ SUMO_ATTR_CONSTANTPOWERINTAKE
Constant Power Intake.
@ SUMO_ATTR_RECUPERATIONEFFICIENCY_BY_DECELERATION
Recuperation efficiency (by deceleration)
@ SUMO_ATTR_RECUPERATIONEFFICIENCY
Recuperation efficiency (constant)
@ SUMO_ATTR_AIRDRAGCOEFFICIENT
Air drag coefficient.
@ SUMO_ATTR_ANGLE
@ SUMO_ATTR_VEHICLEMASS
Vehicle mass.
@ SUMO_ATTR_RADIALDRAGCOEFFICIENT
Radial drag coefficient.
@ SUMO_ATTR_PROPULSIONEFFICIENCY
Propulsion efficiency.
@ SUMO_ATTR_INTERNALMOMENTOFINERTIA
Internal moment of inertia.
@ SUMO_ATTR_FRONTSURFACEAREA
Front surface area.
T MAX2(T a, T b)
Definition: StdDefs.h:79
T MAX3(T a, T b, T c)
Definition: StdDefs.h:93
double compute(const SUMOEmissionClass c, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const std::map< int, double > *param) const
Computes the emitted pollutant amount using the given speed and acceleration.
std::map< int, double > myDefaultParameter
The default parameter.
Definition: HelpersEnergy.h:81
HelpersEnergy()
Constructor (initializes myEmissionClassStrings)
double acceleration(const SUMOEmissionClass c, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const std::map< int, double > *param) const
Computes the achievable acceleration using the given speed and amount of consumed electric power.
Helper methods for PHEMlight-based emission computation.
EmissionType
Enumerating all emission types, including fuel.
static std::tuple< int, double, double, double > cubicSolve(double a, double b, double c, double d)
Solver of cubic equation ax^3 + bx^2 + cx + d = 0.
Definition: PolySolver.cpp:77