Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2001-2025 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 : /****************************************************************************/
14 : /// @file HelpersEnergy.cpp
15 : /// @author Daniel Krajzewicz
16 : /// @author Jakob Erdmann
17 : /// @author Michael Behrisch
18 : /// @author Mirko Barthauer
19 : /// @date Mon, 10.05.2004
20 : ///
21 : // Helper methods for HBEFA-based emission computation
22 : /****************************************************************************/
23 : #include <config.h>
24 :
25 : #include <complex> // due to sqrt of complex
26 : #include <utils/common/MsgHandler.h> // due to WRITE_WARNING
27 : #include <utils/common/SUMOTime.h>
28 : #include <utils/common/ToString.h>
29 : #include <utils/common/PolySolver.h>
30 : #include "HelpersEnergy.h"
31 :
32 :
33 : // ===========================================================================
34 : // method definitions
35 : // ===========================================================================
36 46284 : HelpersEnergy::HelpersEnergy():
37 46284 : PollutantsInterface::Helper("Energy", ENERGY_BASE, ENERGY_BASE)
38 46284 : { }
39 :
40 :
41 : double
42 2444882 : HelpersEnergy::compute(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
43 2444882 : if (e != PollutantsInterface::ELEC) {
44 : return 0.;
45 : }
46 698342 : if (param == nullptr) {
47 0 : param = EnergyParams::getDefault();
48 : }
49 698342 : if (param->isOff()) {
50 : return 0.;
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 621147 : const double lastV = v - ACCEL2SPEED(a);
56 621147 : const double mass = param->getTotalMass(myDefaultMass, 0.);
57 :
58 : // calculate power needed for potential energy difference
59 621147 : double power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
60 :
61 : // power needed for kinetic energy difference of vehicle
62 621147 : power += 0.5 * mass * (v * v - lastV * lastV) / TS;
63 :
64 : // add power needed for rotational energy diff of internal rotating elements
65 621147 : power += 0.5 * param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, myDefaultRotatingMass) * (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 621147 : power += 0.5 * 1.2041 * param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, myDefaultFrontSurfaceArea) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, myDefaultAirDragCoefficient) * 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 621147 : power += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, myDefaultRollDragCoefficient) * GRAVITY * mass * v;
82 :
83 : // Energy loss through friction by radial force [Ws]
84 : // If angle of vehicle was changed
85 621147 : const double angleDiff = param->getAngleDiff();
86 621147 : if (angleDiff != 0.) {
87 : // Compute new radio
88 6207 : 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 6207 : if (radius < 0.0001) {
92 : radius = 0.0001;
93 6044 : } 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 6207 : power += param->getDoubleOptional(SUMO_ATTR_RADIALDRAGCOEFFICIENT, myDefaultRadialDragCoefficient) * mass * v * v / radius / TS;
100 : }
101 :
102 : // E_Bat = E_kin_pot + EnergyLoss;
103 621147 : if (power > 0) {
104 : // Assumption: Efficiency of myPropulsionEfficiency when accelerating
105 377871 : power /= param->getDoubleOptional(SUMO_ATTR_PROPULSIONEFFICIENCY, myDefaultPropulsionEfficiency);
106 : } else {
107 : // Assumption: Efficiency of myRecuperationEfficiency when recuperating
108 243276 : power *= param->getDoubleOptional(SUMO_ATTR_RECUPERATIONEFFICIENCY, myDefaultRecuperationEfficiency);
109 243276 : if (a != 0) {
110 : // Fiori, Chiara & Ahn, Kyoungho & Rakha, Hesham. (2016).
111 : // Power-based electric vehicle energy consumption model: Model
112 : // development and validation. Applied Energy. 168. 257-268.
113 : // 10.1016/j.apenergy.2016.01.097.
114 : //
115 : // Insaf Sagaama, Amine Kchiche, Wassim Trojet, Farouk Kamoun
116 : // Improving The Accuracy of The Energy Consumption Model for
117 : // Electric Vehicle in SUMO Considering The Ambient Temperature
118 : // Effects
119 20383 : power *= (1 / exp(param->getDoubleOptional(SUMO_ATTR_RECUPERATIONEFFICIENCY_BY_DECELERATION, myDefaultRecuperationEfficiencyByDeceleration) / fabs(a)));
120 : }
121 : }
122 :
123 : // EnergyLoss,constantConsumers
124 : // Energy loss through constant loads (e.g. A/C) [W]
125 621147 : power += param->getDoubleOptional(SUMO_ATTR_CONSTANTPOWERINTAKE, myDefaultConstantPowerIntake);
126 :
127 : // convert from [W], to [Wh/s] (3600s / 1h):
128 621147 : return power / 3600.;
129 : }
130 :
131 :
132 : double
133 70 : HelpersEnergy::acceleration(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const EnergyParams* param) const {
134 70 : if (e != PollutantsInterface::ELEC) {
135 : return 0.;
136 : }
137 :
138 70 : if (param == nullptr) {
139 0 : 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 70 : const double mass = param->getTotalMass(myDefaultMass, 0.);
161 70 : const double rotMass = param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, myDefaultRotatingMass);
162 : double const1, const2, const3;
163 : double Prest;
164 : int numX;
165 : double x1, x2, x3;
166 :
167 : // Power delivered by the drive
168 70 : if (P > 0) {
169 : // Assumption: Efficiency of myPropulsionEfficiency when accelerating
170 70 : Prest = P * 3600 * param->getDoubleOptional(SUMO_ATTR_PROPULSIONEFFICIENCY, myDefaultPropulsionEfficiency);
171 : } else {
172 : // Assumption: Efficiency of myRecuperationEfficiency when recuperating
173 0 : Prest = P * 3600 / param->getDoubleOptional(SUMO_ATTR_RECUPERATIONEFFICIENCY, myDefaultRecuperationEfficiency);
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 70 : Prest -= mass * GRAVITY * sin(DEG2RAD(slope)) * (v);
180 70 : 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->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * v;
184 : // i.e. in terms of 'lastV' and 'a': power += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * (lastV + TS * a);
185 70 : Prest -= param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, myDefaultRollDragCoefficient) * GRAVITY * mass * v;
186 70 : const1 += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, myDefaultRollDragCoefficient) * 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->getDoubleOptional(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 70 : const1 += 0.5 * mass * (2 * v);
196 70 : 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 70 : const1 += 0.5 * rotMass * (2 * v);
202 70 : 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->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA) * param->getDoubleOptional(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 70 : const double drag = 0.5 * 1.2041 * param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, myDefaultFrontSurfaceArea) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, myDefaultAirDragCoefficient);
208 70 : Prest -= drag * (v * v * v);
209 70 : const1 += drag * (3 * v * v * TS);
210 70 : const2 += drag * (3 * v * TS * TS);
211 70 : const3 = drag * (TS * TS * TS);
212 :
213 : // Prest = const1*a + const2*a^2 + const3*a^3
214 : // Solve cubic equation in `a` for real roots, and return the number of roots in `numX`
215 : // and the roots in `x1`, `x2`, and `x3`.
216 70 : std::tie(numX, x1, x2, x3) = PolySolver::cubicSolve(const3, const2, const1, -Prest);
217 :
218 :
219 70 : switch (numX) {
220 : case 1:
221 : return x1;
222 0 : case 2:
223 : return MAX2(x1, x2);
224 70 : case 3:
225 : return MAX3(x1, x2, x3);
226 0 : default:
227 0 : WRITE_ERROR(TL("An acceleration given by the power was not found."));
228 0 : return 0;
229 : }
230 :
231 : }
232 :
233 :
234 : /****************************************************************************/
|
