Eclipse SUMO - Simulation of Urban MObility
Loading...
Searching...
No Matches
HelpersEnergy.cpp
Go to the documentation of this file.
1/****************************************************************************/
2// Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3// Copyright (C) 2001-2024 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
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// ===========================================================================
40 PollutantsInterface::Helper("Energy", ENERGY_BASE, ENERGY_BASE)
41{ }
42
43
44double
45HelpersEnergy::compute(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
47 return 0.;
48 }
49 if (param == nullptr) {
51 }
52 //@ToDo: All formulas below work with the logic of the euler update (refs #860).
53 // Approximation order could be improved. Refs. #2592.
54
55 const double lastV = v - ACCEL2SPEED(a);
56 const double mass = param->getDouble(SUMO_ATTR_MASS) + param->getDoubleOptional(SUMO_ATTR_LOADING, 0.);
57
58 // calculate power needed for potential energy difference
59 double power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
60
61 // power needed for kinetic energy difference of vehicle
62 power += 0.5 * mass * (v * v - lastV * lastV) / TS;
63
64 // add power needed for rotational energy diff of internal rotating elements
65 power += 0.5 * param->getDouble(SUMO_ATTR_ROTATINGMASS) * (v * v - lastV * lastV) / TS;
66
67 // power needed for Energy loss through Air resistance [Ws]
68 // Calculate energy losses:
69 // EnergyLoss,Air = 1/2 * rho_air [kg/m^3] * myFrontSurfaceArea [m^2] * myAirDragCoefficient [-] * v_Veh^2 [m/s] * s [m]
70 // ... with rho_air [kg/m^3] = 1,2041 kg/m^3 (at T = 20C)
71 // ... with s [m] = v_Veh [m/s] * TS [s]
72 // divided by TS to get power instead of energy
73 power += 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * v * v * v;
74
75 // power needed for Energy loss through Roll resistance [Ws]
76 // ... (fabs(veh.getSpeed())>=0.01) = 0, if vehicle isn't moving
77 // EnergyLoss,Tire = c_R [-] * F_N [N] * s [m]
78 // ... with c_R = ~0.012 (car tire on asphalt)
79 // ... with F_N [N] = myMass [kg] * g [m/s^2]
80 // divided by TS to get power instead of energy
81 power += param->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * v;
82
83 // Energy loss through friction by radial force [Ws]
84 // If angle of vehicle was changed
85 const double angleDiff = param->getDouble(SUMO_ATTR_ANGLE);
86 if (angleDiff != 0.) {
87 // Compute new radio
88 double radius = SPEED2DIST(v) / fabs(angleDiff);
89
90 // Check if radius is in the interval [0.0001 - 10000] (To avoid overflow and division by zero)
91 if (radius < 0.0001) {
92 radius = 0.0001;
93 } else if (radius > 10000) {
94 radius = 10000;
95 }
96 // EnergyLoss,internalFrictionRadialForce = c [m] * F_rad [N];
97 // Energy loss through friction by radial force [Ws]
98 // divided by TS to get power instead of energy
99 power += param->getDouble(SUMO_ATTR_RADIALDRAGCOEFFICIENT) * mass * v * v / radius * v;
100 }
101
102 // EnergyLoss,constantConsumers
103 // Energy loss through constant loads (e.g. A/C) [W]
105
106 // E_Bat = E_kin_pot + EnergyLoss;
107 if (power > 0) {
108 // Assumption: Efficiency of myPropulsionEfficiency when accelerating
110 } else {
111 // Assumption: Efficiency of myRecuperationEfficiency when recuperating
113 if (a != 0) {
114 // Fiori, Chiara & Ahn, Kyoungho & Rakha, Hesham. (2016).
115 // Power-based electric vehicle energy consumption model: Model
116 // development and validation. Applied Energy. 168. 257-268.
117 // 10.1016/j.apenergy.2016.01.097.
118 //
119 // Insaf Sagaama, Amine Kchiche, Wassim Trojet, Farouk Kamoun
120 // Improving The Accuracy of The Energy Consumption Model for
121 // Electric Vehicle in SUMO Considering The Ambient Temperature
122 // Effects
123 power *= (1 / exp(param->getDouble(SUMO_ATTR_RECUPERATIONEFFICIENCY_BY_DECELERATION) / fabs(a)));
124 }
125 }
126
127 // convert from [W], to [Wh/s] (3600s / 1h):
128 return power / 3600.;
129}
130
131
132double
133HelpersEnergy::acceleration(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const EnergyParams* param) const {
134 if (e != PollutantsInterface::ELEC) {
135 return 0.;
136 }
137
138 if (param == nullptr) {
139 param = EnergyParams::getDefault();
140 }
141
142 // Inverse formula for the function compute()
143
144 // Computes the achievable acceleration using the given speed and drive power in [Wh/s]
145 // It does not consider friction losses by radial force and the acceleration-dependent
146 // recuperation efficiency (eta_recuperation is considered constant)
147 //
148 // The power `Prest` used for acceleration is given by a cubic polynomial in acceleration,
149 // i.e.the equation
150 //
151 // Prest = const1*a + const2*a^2 + const3*a^3
152 //
153 // and we need to find `a` given `Prest`.
154 //
155 // The solutions of `a(P)` is stored in variables `x1`, `x2`, and `x3` returned by
156 // the method `PolySolver::cubicSolve()`, see below.
157 //
158 // Power used for accelerating, `Prest`, is the total used power minus power wasted by running resistances.
159
160 const double mass = param->getDouble(SUMO_ATTR_MASS) + param->getDoubleOptional(SUMO_ATTR_LOADING, 0.);
161 const double rotMass = param->getDouble(SUMO_ATTR_ROTATINGMASS);
162 double const1, const2, const3;
163 double Prest;
164 int numX;
165 double x1, x2, x3;
166
167 // Power delivered by the drive
168 if (P > 0) {
169 // Assumption: Efficiency of myPropulsionEfficiency when accelerating
170 Prest = P * 3600 * param->getDouble(SUMO_ATTR_PROPULSIONEFFICIENCY);
171 } else {
172 // Assumption: Efficiency of myRecuperationEfficiency when recuperating
173 Prest = P * 3600 / param->getDouble(SUMO_ATTR_RECUPERATIONEFFICIENCY);
174 }
175
176 // calculate power drop due to a potential energy difference
177 // in inverse original formula: power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
178 // i.e. in terms of 'lastV' and 'a': power = mass * GRAVITY * sin(DEG2RAD(slope)) * (lastV + TS * a);
179 Prest -= mass * GRAVITY * sin(DEG2RAD(slope)) * (v);
180 const1 = mass * GRAVITY * sin(DEG2RAD(slope)) * (TS);
181
182 // update coefficients of a(P) equation considering power loss through Roll resistance
183 // in inverse original formula: power += param->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * v;
184 // i.e. in terms of 'lastV' and 'a': power += param->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * (lastV + TS * a);
185 Prest -= param->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * v;
186 const1 += param->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * (TS);
187
188 //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
189 //Power loss through constant loads (e.g. A/C) [W]
190 //Prest -= param->getDouble(SUMO_ATTR_CONSTANTPOWERINTAKE) / TS;
191
192 // update coefficients of a(P) equation considering kinetic energy difference of vehicle
193 // in inverse original formula: power += 0.5 * mass * (v * v - lastV * lastV) / TS;
194 // i.e. in terms of 'lastV' and 'a': power += 0.5 * mass * (2 * lastV * a + TS * a * a);
195 const1 += 0.5 * mass * (2 * v);
196 const2 = 0.5 * mass * (TS);
197
198 // update coefficients of a(P) equation considering rotational energy diff of internal rotating elements
199 // in inverse original formula: power += 0.5 * rotMass * (v * v - lastV * lastV) / TS;
200 // i.e. in terms of 'lastV' and 'a': power += 0.5 * rotMass * (2 * lastV * a + TS * a * a);
201 const1 += 0.5 * rotMass * (2 * v);
202 const2 += 0.5 * rotMass * (TS);
203
204 // update coefficients of a(P) equation considering energy loss through Air resistance
205 // in inverse original formula: power += 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * v * v * v;
206 // i.e. in terms of 'lastV' and 'a': power += 0.5 * rotMass * (lastV^3 + 3* lastV^2 * TS * a + 3* lastV * TS^2 *a^2 + TS^3 * a^3);
207 Prest -= 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (v * v * v);
208 const1 += 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (3 * v * v * TS);
209 const2 += 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (3 * v * TS * TS);
210 const3 = 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (TS * TS * TS);
211
212 // Prest = const1*a + const2*a^2 + const3*a^3
213 // Solve cubic equation in `a` for real roots, and return the number of roots in `numX`
214 // and the roots in `x1`, `x2`, and `x3`.
215 std::tie(numX, x1, x2, x3) = PolySolver::cubicSolve(const3, const2, const1, -Prest);
216
217
218 switch (numX) {
219 case 1:
220 return x1;
221 case 2:
222 return MAX2(x1, x2);
223 case 3:
224 return MAX3(x1, x2, x3);
225 default:
226 WRITE_ERROR(TL("An acceleration given by the power was not found."));
227 return 0;
228 }
229
230}
231
232
233/****************************************************************************/
#define DEG2RAD(x)
Definition GeomHelper.h:35
#define GRAVITY
Definition GeomHelper.h:37
#define WRITE_ERROR(msg)
Definition MsgHandler.h:304
#define TL(string)
Definition MsgHandler.h:315
#define SPEED2DIST(x)
Definition SUMOTime.h:45
#define ACCEL2SPEED(x)
Definition SUMOTime.h:51
#define TS
Definition SUMOTime.h:42
int SUMOEmissionClass
@ SUMO_ATTR_MASS
@ SUMO_ATTR_ROLLDRAGCOEFFICIENT
Roll Drag coefficient.
@ SUMO_ATTR_CONSTANTPOWERINTAKE
Constant Power Intake.
@ SUMO_ATTR_LOADING
additional mass loaded on the vehicle
@ 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_RADIALDRAGCOEFFICIENT
Radial drag coefficient.
@ SUMO_ATTR_ROTATINGMASS
Mass equivalent of rotating elements.
@ SUMO_ATTR_PROPULSIONEFFICIENCY
Propulsion efficiency.
@ SUMO_ATTR_FRONTSURFACEAREA
Front surface area.
T MAX2(T a, T b)
Definition StdDefs.h:82
T MAX3(T a, T b, T c)
Definition StdDefs.h:96
An upper class for objects with additional parameters.
double getDouble(SumoXMLAttr attr) const
double getDoubleOptional(SumoXMLAttr attr, const double def) const
static const EnergyParams * getDefault()
double compute(const SUMOEmissionClass c, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams *param) const
Computes the emitted pollutant amount using the given speed and acceleration.
double acceleration(const SUMOEmissionClass c, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const EnergyParams *param) const
Computes the achievable acceleration using the given speed and amount of consumed electric power.
HelpersEnergy()
Constructor (initializes myEmissionClassStrings)
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.