refactoring phemlight helper to allow param insertion #14285

eclipse-sumo · Jan 31, 2024 · 48b570d · 48b570d
1 parent c252b4a
commit 48b570d
Show file tree

Hide file tree

Showing 7 changed files with 93 additions and 19 deletions.
diff --git a/src/foreign/PHEMlight/V5/cpp/CEP.cpp b/src/foreign/PHEMlight/V5/cpp/CEP.cpp
@@ -198,6 +198,10 @@ namespace PHEMlightdllV5 {
         privateCalcType = value;
     }
 
+    bool CEP::isHBEV() const {
+        return getFuelType() == Constants::strBEV || getFuelType() == Constants::strHybrid;
+    }
+
     const double& CEP::getRatedPower() const {
         return privateRatedPower;
     }

diff --git a/src/foreign/PHEMlight/V5/cpp/CEP.h b/src/foreign/PHEMlight/V5/cpp/CEP.h
@@ -63,13 +63,32 @@ namespace PHEMlightdllV5 {
         void setFuelType(const std::string& value);
         const std::string&  getCalcType() const;
         void setCalcType(const std::string& value);
+        bool isHBEV() const;
 
     public:
         const double&  getRatedPower() const;
         void setRatedPower(const double&  value);
         double getAuxPower() const {
             return _auxPower * getRatedPower();
         }
+        double getVehicleMass() const {
+            return _massVehicle;
+        }
+        double getVehicleLoading() const {
+            return _vehicleLoading;
+        }
+        double getVehicleMassRot() const {
+            return _vehicleMassRot;
+        }
+        double getCrossSectionalArea() const {
+            return _crossSectionalArea;
+        }
+        double getCWValue() const {
+            return _cWValue;
+        }
+        double getResistance(const double speed) const {
+            return _resistanceF0 + _resistanceF1 * speed + _resistanceF4 * std::pow(speed, 4);
+        }
 
     protected:
         double _massVehicle;
@@ -142,7 +161,6 @@ namespace PHEMlightdllV5 {
     public:
         double GetMaxAccel(double speed, double gradient, bool HBEV);
 
-    private:
         double GetPMaxNorm(double speed);
 
         //--------------------------------------------------------------------------------------------------

diff --git a/src/tools/TrajectoriesHandler.cpp b/src/tools/TrajectoriesHandler.cpp
@@ -130,7 +130,7 @@ TrajectoriesHandler::computeEmissions(const std::string id, const SUMOEmissionCl
         }
     }
     if (myAccelZeroCorrection) {
-        a = PollutantsInterface::getModifiedAccel(c, v, a, s);
+        a = PollutantsInterface::getModifiedAccel(c, v, a, s, params);
     }
     if (a == INVALID_VALUE) {
         throw ProcessError(TL("Acceleration information is missing; try running with --compute-a."));

diff --git a/src/utils/emissions/HelpersPHEMlight5.cpp b/src/utils/emissions/HelpersPHEMlight5.cpp
@@ -116,11 +116,45 @@ HelpersPHEMlight5::getEmission(PHEMlightdllV5::CEP* currCep, const std::string&
 
 
 double
-HelpersPHEMlight5::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope) const {
+HelpersPHEMlight5::calcPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const {
+    // copy of CEP::CalcPower
+    const double power = calcWheelPower(currCep, v, a, slope, param) / PHEMlightdllV5::Constants::_DRIVE_TRAIN_EFFICIENCY;
+    if (!currCep->isHBEV()) {
+        return power + currCep->getAuxPower();
+    }
+    return power;
+}
+
+
+double
+HelpersPHEMlight5::calcWheelPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* /* param */) const {
+    // copy of CEP::CalcWheelPower
+    const double rotFactor = currCep->GetRotationalCoeffecient(v);
+    const double mass = currCep->getVehicleMass();
+    const double massRot = currCep->getVehicleMassRot();
+    const double load = currCep->getVehicleLoading();
+    const double cw = currCep->getCrossSectionalArea() * currCep->getCWValue();
+
+    double power = (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * currCep->getResistance(v) * v;
+    power += (cw * PHEMlightdllV5::Constants::AIR_DENSITY_CONST / 2) * std::pow(v, 3);
+    power += (mass * rotFactor + massRot + load) * a * v;
+    power += (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * slope * 0.01 * v;
+    return power / 1000.;
+}
+
+
+double
+HelpersPHEMlight5::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const {
     PHEMlightdllV5::CEP* currCep = myCEPs.count(c) == 0 ? nullptr : myCEPs.find(c)->second;
     if (currCep != nullptr) {
-        const bool isHBEV = currCep->getFuelType() == PHEMlightdllV5::Constants::strBEV || currCep->getFuelType() == PHEMlightdllV5::Constants::strHybrid;
-        return v == 0.0 ? 0.0 : MIN2(a, currCep->GetMaxAccel(v, slope, isHBEV));
+        if (v == 0.) {
+            return 0.;
+        }
+        // this is a copy of CEP::GetMaxAccel
+        const double rotFactor = currCep->GetRotationalCoeffecient(v);
+        const double pMaxForAcc = currCep->GetPMaxNorm(v) * currCep->getRatedPower() - calcPower(currCep, v, 0, slope, param);
+        const double maxAcc = (pMaxForAcc * 1000) / ((currCep->getVehicleMass() * rotFactor + currCep->getVehicleMassRot() + currCep->getVehicleLoading()) * v);
+        return MIN2(a, maxAcc);
     }
     return a;
 }
@@ -140,17 +174,16 @@ HelpersPHEMlight5::compute(const SUMOEmissionClass c, const PollutantsInterface:
     const double corrSpeed = MAX2(0.0, v);
     assert(myCEPs.count(c) == 1);
     PHEMlightdllV5::CEP* const currCep = myCEPs.find(c)->second;
-    const double corrAcc = getModifiedAccel(c, corrSpeed, a, slope);
+    const double corrAcc = getModifiedAccel(c, corrSpeed, a, slope, param);
     const bool isBEV = currCep->getFuelType() == PHEMlightdllV5::Constants::strBEV;
     const bool isHybrid = currCep->getFuelType() == PHEMlightdllV5::Constants::strHybrid;
-    const double power_raw = currCep->CalcPower(corrSpeed, corrAcc, slope, isBEV || isHybrid);
-    const double power = isHybrid ? currCep->CalcWheelPower(corrSpeed, corrAcc, slope) : currCep->CalcEngPower(power_raw);
+    const double power_raw = calcPower(currCep, corrSpeed, corrAcc, slope, param);
+    const double power = isHybrid ? calcWheelPower(currCep, corrSpeed, corrAcc, slope, param) : currCep->CalcEngPower(power_raw);
 
     if (!isBEV && corrAcc < currCep->GetDecelCoast(corrSpeed, corrAcc, slope) &&
             corrSpeed > PHEMlightdllV5::Constants::ZERO_SPEED_ACCURACY) {
         return 0.;
     }
-    const std::string& fuelType = currCep->getFuelType();
     switch (e) {
         case PollutantsInterface::CO:
             return getEmission(currCep, "CO", power, corrSpeed) / SECONDS_PER_HOUR * 1000.;
@@ -165,19 +198,19 @@ HelpersPHEMlight5::compute(const SUMOEmissionClass c, const PollutantsInterface:
         case PollutantsInterface::PM_X:
             return getEmission(currCep, "PM", power, corrSpeed) / SECONDS_PER_HOUR * 1000.;
         case PollutantsInterface::FUEL: {
-            if (myVolumetricFuel && fuelType == PHEMlightdllV5::Constants::strDiesel) { // divide by average diesel density of 836 g/l
+            if (myVolumetricFuel && currCep->getFuelType() == PHEMlightdllV5::Constants::strDiesel) { // divide by average diesel density of 836 g/l
                 return getEmission(currCep, "FC", power, corrSpeed) / 836. / SECONDS_PER_HOUR * 1000.;
             }
-            if (myVolumetricFuel && fuelType == PHEMlightdllV5::Constants::strGasoline) { // divide by average gasoline density of 742 g/l
+            if (myVolumetricFuel && currCep->getFuelType() == PHEMlightdllV5::Constants::strGasoline) { // divide by average gasoline density of 742 g/l
                 return getEmission(currCep, "FC", power, corrSpeed) / 742. / SECONDS_PER_HOUR * 1000.;
             }
-            if (fuelType == PHEMlightdllV5::Constants::strBEV) {
+            if (isBEV) {
                 return 0.;
             }
             return getEmission(currCep, "FC", power, corrSpeed) / SECONDS_PER_HOUR * 1000.; // still in mg even if myVolumetricFuel is set!
         }
         case PollutantsInterface::ELEC:
-            if (fuelType == PHEMlightdllV5::Constants::strBEV) {
+            if (isBEV) {
                 return (getEmission(currCep, "FC_el", power, corrSpeed) + currCep->getAuxPower()) / SECONDS_PER_HOUR * 1000.;
             }
             return 0;

diff --git a/src/utils/emissions/HelpersPHEMlight5.h b/src/utils/emissions/HelpersPHEMlight5.h
@@ -75,7 +75,7 @@ class HelpersPHEMlight5 : public HelpersPHEMlight {
      * @param[in] slope The road's slope at vehicle's position [deg]
      * @return the modified acceleration
      */
-    double getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope) const;
+    double getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const;
 
     /** @brief Returns the maximum deceleration value (as a negative number), which can still be considered as non-braking.
      * @param[in] c the emission class
@@ -97,6 +97,24 @@ class HelpersPHEMlight5 : public HelpersPHEMlight {
     */
     double getEmission(PHEMlightdllV5::CEP* currCep, const std::string& e, const double p, const double v) const;
 
+    /** @brief Returns the total power needed.
+     * @param[in] currCep the emission class
+     * @param[in] v the speed value
+     * @param[in] a the acceleration value
+     * @param[in] slope The road's slope at vehicle's position [deg]
+     * @return the total power needed
+     */
+    double calcPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const;
+
+    /** @brief Returns the power without auxiliaries.
+     * @param[in] currCep the emission class
+     * @param[in] v the speed value
+     * @param[in] a the acceleration value
+     * @param[in] slope The road's slope at vehicle's position [deg]
+     * @return the power without auxiliaries
+     */
+    double calcWheelPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const;
+
     /// @brief the index of the next class
     int myIndex;
     PHEMlightdllV5::CEPHandler myCEPHandler;

diff --git a/src/utils/emissions/PollutantsInterface.cpp b/src/utils/emissions/PollutantsInterface.cpp
@@ -183,10 +183,11 @@ PollutantsInterface::Helper::compute(const SUMOEmissionClass c, const EmissionTy
 
 
 double
-PollutantsInterface::Helper::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope) const {
+PollutantsInterface::Helper::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const {
     UNUSED_PARAMETER(c);
     UNUSED_PARAMETER(v);
     UNUSED_PARAMETER(slope);
+    UNUSED_PARAMETER(param);
     return a;
 }
 
@@ -372,8 +373,8 @@ PollutantsInterface::computeDefault(const SUMOEmissionClass c, const EmissionTyp
 
 
 double
-PollutantsInterface::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope) {
-    return myHelpers[c >> 16]->getModifiedAccel(c, v, a, slope);
+PollutantsInterface::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) {
+    return myHelpers[c >> 16]->getModifiedAccel(c, v, a, slope, param);
 }
 
 

diff --git a/src/utils/emissions/PollutantsInterface.h b/src/utils/emissions/PollutantsInterface.h
@@ -197,7 +197,7 @@ class PollutantsInterface {
          * @param[in] slope The road's slope at vehicle's position [deg]
          * @return the modified acceleration
          */
-        virtual double getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope) const;
+        virtual double getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const;
 
         /** @brief Returns the maximum deceleration value (as a negative number), which can still be considered as non-braking.
          * Default implementation returns always zero.
@@ -353,7 +353,7 @@ class PollutantsInterface {
      * @param[in] slope The road's slope at vehicle's position [deg]
      * @return the modified acceleration
      */
-    static double getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope);
+    static double getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param);
 
     /** @brief Returns the coasting deceleration value, useful for comparing with external PHEMlight references.
      * @param[in] c the emission class