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       45820 : HelpersEnergy::HelpersEnergy():
      40       45820 :     PollutantsInterface::Helper("Energy", ENERGY_BASE, ENERGY_BASE)
      41       45820 : { }
      42             : 
      43             : 
      44             : double
      45     1345756 : HelpersEnergy::compute(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
      46     1345756 :     if (e != PollutantsInterface::ELEC) {
      47             :         return 0.;
      48             :     }
      49      363532 :     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      363532 :     const double lastV = v - ACCEL2SPEED(a);
      56      363532 :     const double mass = param->getDouble(SUMO_ATTR_MASS) + param->getDoubleOptional(SUMO_ATTR_LOADING, 0.);
      57             : 
      58             :     // calculate power needed for potential energy difference
      59      363532 :     double power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
      60             : 
      61             :     // power needed for kinetic energy difference of vehicle
      62      363532 :     power += 0.5 * mass * (v * v - lastV * lastV) / TS;
      63             : 
      64             :     // add power needed for rotational energy diff of internal rotating elements
      65      363532 :     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      363532 :     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      363532 :     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      363532 :     const double angleDiff = param->getDouble(SUMO_ATTR_ANGLE);
      86      363532 :     if (angleDiff != 0.) {
      87             :         // Compute new radio
      88        4188 :         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        4188 :         if (radius < 0.0001) {
      92             :             radius = 0.0001;
      93        4181 :         } 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        4188 :         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      363532 :     power += param->getDouble(SUMO_ATTR_CONSTANTPOWERINTAKE);
     105             : 
     106             :     // E_Bat = E_kin_pot + EnergyLoss;
     107      363532 :     if (power > 0) {
     108             :         // Assumption: Efficiency of myPropulsionEfficiency when accelerating
     109      351497 :         power /= param->getDouble(SUMO_ATTR_PROPULSIONEFFICIENCY);
     110             :     } else {
     111             :         // Assumption: Efficiency of myRecuperationEfficiency when recuperating
     112       12035 :         power *= param->getDouble(SUMO_ATTR_RECUPERATIONEFFICIENCY);
     113       12035 :         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       12035 :             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      363532 :     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             : /****************************************************************************/