Line data Source code
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 : /****************************************************************************/
14 : /// @file HelpersEnergy.cpp
15 : /// @author Daniel Krajzewicz
16 : /// @author Jakob Erdmann
17 : /// @author Michael Behrisch
18 : /// @date Mon, 10.05.2004
19 : ///
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>
26 : #include <utils/common/PolySolver.h>
27 : #include "HelpersEnergy.h"
28 :
29 : //due to WRITE_WARNING
30 : #include <utils/common/MsgHandler.h>
31 :
32 : //due to sqrt of complex
33 : #include <complex>
34 :
35 :
36 : // ===========================================================================
37 : // method definitions
38 : // ===========================================================================
39 55625 : HelpersEnergy::HelpersEnergy():
40 55625 : PollutantsInterface::Helper("Energy", ENERGY_BASE, ENERGY_BASE)
41 55625 : { }
42 :
43 :
44 : double
45 2322679 : HelpersEnergy::compute(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
46 2322679 : if (e != PollutantsInterface::ELEC) {
47 : return 0.;
48 : }
49 648859 : if (param == nullptr) {
50 0 : param = EnergyParams::getDefault();
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 648859 : const double lastV = v - ACCEL2SPEED(a);
56 648859 : const double mass = param->getDouble(SUMO_ATTR_MASS) + param->getDoubleOptional(SUMO_ATTR_LOADING, 0.);
57 :
58 : // calculate power needed for potential energy difference
59 648859 : double power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
60 :
61 : // power needed for kinetic energy difference of vehicle
62 648859 : power += 0.5 * mass * (v * v - lastV * lastV) / TS;
63 :
64 : // add power needed for rotational energy diff of internal rotating elements
65 648859 : 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 648859 : 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 648859 : 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 648859 : const double angleDiff = param->getDouble(SUMO_ATTR_ANGLE);
86 648859 : if (angleDiff != 0.) {
87 : // Compute new radio
88 5067 : 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 5067 : if (radius < 0.0001) {
92 : radius = 0.0001;
93 5058 : } 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 5067 : 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]
104 648859 : power += param->getDouble(SUMO_ATTR_CONSTANTPOWERINTAKE);
105 :
106 : // E_Bat = E_kin_pot + EnergyLoss;
107 648859 : if (power > 0) {
108 : // Assumption: Efficiency of myPropulsionEfficiency when accelerating
109 594286 : power /= param->getDouble(SUMO_ATTR_PROPULSIONEFFICIENCY);
110 : } else {
111 : // Assumption: Efficiency of myRecuperationEfficiency when recuperating
112 54573 : power *= param->getDouble(SUMO_ATTR_RECUPERATIONEFFICIENCY);
113 54573 : 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 19249 : 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 648859 : return power / 3600.;
129 : }
130 :
131 :
132 : double
133 71 : HelpersEnergy::acceleration(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const EnergyParams* param) const {
134 71 : if (e != PollutantsInterface::ELEC) {
135 : return 0.;
136 : }
137 :
138 71 : 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 71 : const double mass = param->getDouble(SUMO_ATTR_MASS) + param->getDoubleOptional(SUMO_ATTR_LOADING, 0.);
161 71 : 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 71 : if (P > 0) {
169 : // Assumption: Efficiency of myPropulsionEfficiency when accelerating
170 71 : Prest = P * 3600 * param->getDouble(SUMO_ATTR_PROPULSIONEFFICIENCY);
171 : } else {
172 : // Assumption: Efficiency of myRecuperationEfficiency when recuperating
173 0 : 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 71 : Prest -= mass * GRAVITY * sin(DEG2RAD(slope)) * (v);
180 71 : 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 71 : Prest -= param->getDouble(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * v;
186 71 : 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 71 : const1 += 0.5 * mass * (2 * v);
196 71 : 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 71 : const1 += 0.5 * rotMass * (2 * v);
202 71 : 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 71 : Prest -= 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (v * v * v);
208 71 : const1 += 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (3 * v * v * TS);
209 71 : const2 += 0.5 * 1.2041 * param->getDouble(SUMO_ATTR_FRONTSURFACEAREA) * param->getDouble(SUMO_ATTR_AIRDRAGCOEFFICIENT) * (3 * v * TS * TS);
210 71 : 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 71 : std::tie(numX, x1, x2, x3) = PolySolver::cubicSolve(const3, const2, const1, -Prest);
216 :
217 :
218 71 : switch (numX) {
219 : case 1:
220 : return x1;
221 0 : case 2:
222 : return MAX2(x1, x2);
223 71 : case 3:
224 : return MAX3(x1, x2, x3);
225 0 : default:
226 0 : WRITE_ERROR(TL("An acceleration given by the power was not found."));
227 0 : return 0;
228 : }
229 :
230 : }
231 :
232 :
233 : /****************************************************************************/
|
