PHEMCEP.cpp

Go to the documentation of this file.

/****************************************************************************/
// Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
// Copyright (C) 2013-2025 German Aerospace Center (DLR) and others.
// This program and the accompanying materials are made available under the
// terms of the Eclipse Public License 2.0 which is available at
// https://www.eclipse.org/legal/epl-2.0/
// This Source Code may also be made available under the following Secondary
// Licenses when the conditions for such availability set forth in the Eclipse
// Public License 2.0 are satisfied: GNU General Public License, version 2
// or later which is available at
// https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
/****************************************************************************/
// Helper class for PHEM Light, holds a specific CEP for a PHEM emission class
/****************************************************************************/
#include <config.h>
 
#include <cmath>
#include <string>
#include <utils/common/StdDefs.h>
#include <utils/common/StringBijection.h>
#include <utils/common/UtilExceptions.h>
#include "PHEMCEP.h"
 
 
// ===========================================================================
// method definitions
// ===========================================================================
PHEMCEP::PHEMCEP(bool heavyVehicle, SUMOEmissionClass emissionClass, const std::string& emissionClassIdentifier,
                 double vehicleMass, double vehicleLoading, double vehicleMassRot,
                 double crossArea, double cdValue,
                 double f0, double f1, double f2, double f3, double f4,
                 double ratedPower, double pNormV0, double pNormP0, double pNormV1, double pNormP1,
                 double axleRatio, double engineIdlingSpeed, double engineRatedSpeed, double effectiveWheelDiameter,
                 double idlingFC,
                 const std::string& vehicleFuelType,
                 const std::vector< std::vector<double> >& matrixFC,
                 const std::vector<std::string>& headerLinePollutants,
                 const std::vector< std::vector<double> >& matrixPollutants,
                 const std::vector< std::vector<double> >& matrixSpeedRotational,
                 const std::vector< std::vector<double> >& normedDragTable,
                 const std::vector<double>& idlingValuesPollutants) {
    _emissionClass = emissionClass;
    _resistanceF0 = f0;
    _resistanceF1 = f1;
    _resistanceF2 = f2;
    _resistanceF3 = f3;
    _resistanceF4 = f4;
    _cdValue = cdValue;
    _crossSectionalArea = crossArea;
    _massVehicle = vehicleMass;
    _vehicleLoading = vehicleLoading;
    _massRot = vehicleMassRot;
    _ratedPower = ratedPower;
    _vehicleFuelType = vehicleFuelType;
 
    _pNormV0 = pNormV0 / 3.6;
    _pNormP0 = pNormP0;
    _pNormV1 = pNormV1 / 3.6;
    _pNormP1 = pNormP1;
 
    _axleRatio = axleRatio;
    _engineIdlingSpeed = engineIdlingSpeed;
    _engineRatedSpeed = engineRatedSpeed;
    _effictiveWheelDiameter = effectiveWheelDiameter;
 
    _heavyVehicle = heavyVehicle;
    _idlingFC = idlingFC;
 
    std::vector<std::string> pollutantIdentifier;
    std::vector< std::vector<double> > pollutantMeasures;
    std::vector<std::vector<double> > normalizedPollutantMeasures;
 
    // init pollutant identifiers
    for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
        pollutantIdentifier.push_back(headerLinePollutants[i]);
    } // end for
 
    // get size of powerPatterns
    _sizeOfPatternFC = (int)matrixFC.size();
    _sizeOfPatternPollutants = (int)matrixPollutants.size();
 
    // initialize measures
    for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
        pollutantMeasures.push_back(std::vector<double>());
        normalizedPollutantMeasures.push_back(std::vector<double>());
    } // end for
 
    // looping through matrix and assigning values for speed rotational table
    _speedCurveRotational.clear();
    _speedPatternRotational.clear();
    _gearTransmissionCurve.clear();
    for (int i = 0; i < (int)matrixSpeedRotational.size(); i++) {
        if (matrixSpeedRotational[i].size() != 3) {
            throw InvalidArgument("Error loading vehicle file for: " + emissionClassIdentifier);
        }
 
        _speedPatternRotational.push_back(matrixSpeedRotational[i][0] / 3.6);
        _speedCurveRotational.push_back(matrixSpeedRotational[i][1]);
        _gearTransmissionCurve.push_back(matrixSpeedRotational[i][2]);
    } // end for
 
    // looping through matrix and assigning values for drag table
    _nNormTable.clear();
    _dragNormTable.clear();
    for (int i = 0; i < (int) normedDragTable.size(); i++) {
        if (normedDragTable[i].size() != 2) {
            return;
        }
 
        _nNormTable.push_back(normedDragTable[i][0]);
        _dragNormTable.push_back(normedDragTable[i][1]);
 
    } // end for
 
    // looping through matrix and assigning values for Fuel consumption
    _cepCurveFC.clear();
    _powerPatternFC.clear();
    _normalizedPowerPatternFC.clear();
    _normedCepCurveFC.clear();
    for (int i = 0; i < (int)matrixFC.size(); i++) {
        if (matrixFC[i].size() != 2) {
            throw InvalidArgument("Error loading vehicle file for: " + emissionClassIdentifier);
        }
 
        _powerPatternFC.push_back(matrixFC[i][0] * _ratedPower);
        _normalizedPowerPatternFC.push_back(matrixFC[i][0]);
        _cepCurveFC.push_back(matrixFC[i][1] * _ratedPower);
        _normedCepCurveFC.push_back(matrixFC[i][1]);
 
    } // end for
 
    _powerPatternPollutants.clear();
    double pollutantMultiplyer = 1;
 
    _drivingPower = _normalizingPower = CalcPower(NORMALIZING_SPEED, NORMALIZING_ACCELARATION, 0, vehicleLoading);
 
    // looping through matrix and assigning values for pollutants
 
    if (heavyVehicle) {
        _normalizingPower = _ratedPower;
        pollutantMultiplyer = _ratedPower;
        _normalizingType = RatedPower;
    } else {
        _normalizingPower = _drivingPower;
        _normalizingType = DrivingPower;
    } // end if
 
    const int headerCount = (int)headerLinePollutants.size();
    for (int i = 0; i < (int)matrixPollutants.size(); i++) {
        for (int j = 0; j < (int)matrixPollutants[i].size(); j++) {
            if ((int)matrixPollutants[i].size() != headerCount + 1) {
                return;
            }
 
            if (j == 0) {
                _normailzedPowerPatternPollutants.push_back(matrixPollutants[i][j]);
                _powerPatternPollutants.push_back(matrixPollutants[i][j] * _normalizingPower);
            } else {
                pollutantMeasures[j - 1].push_back(matrixPollutants[i][j] * pollutantMultiplyer);
                normalizedPollutantMeasures[j - 1].push_back(matrixPollutants[i][j]);
            } // end if
        } // end for
    } // end for
 
    for (int i = 0; i < (int) headerLinePollutants.size(); i++) {
        _cepCurvePollutants.insert(pollutantIdentifier[i], pollutantMeasures[i]);
        _normalizedCepCurvePollutants.insert(pollutantIdentifier[i], normalizedPollutantMeasures[i]);
        _idlingValuesPollutants.insert(pollutantIdentifier[i], idlingValuesPollutants[i] * pollutantMultiplyer);
    } // end for
 
    _idlingFC = idlingFC * _ratedPower;
 
} // end of Cep
 
 
PHEMCEP::~PHEMCEP() {
    // free power pattern
    _powerPatternFC.clear();
    _powerPatternPollutants.clear();
    _cepCurveFC.clear();
    _speedCurveRotational.clear();
    _speedPatternRotational.clear();
} // end of ~Cep
 
 
double
PHEMCEP::GetEmission(const std::string& pollutant, double power, double speed, bool normalized) const {
    std::vector<double> emissionCurve;
    std::vector<double> powerPattern;
 
    if (!normalized && fabs(speed) <= ZERO_SPEED_ACCURACY) {
        if (pollutant == "FC") {
            return _idlingFC;
        } else {
            return _idlingValuesPollutants.get(pollutant);
        }
    } // end if
 
    if (pollutant == "FC") {
        if (normalized) {
            emissionCurve = _normedCepCurveFC;
            powerPattern = _normalizedPowerPatternFC;
        } else {
            emissionCurve = _cepCurveFC;
            powerPattern = _powerPatternFC;
        }
    } else {
        if (!_cepCurvePollutants.hasString(pollutant)) {
            throw InvalidArgument("Emission pollutant " + pollutant + " not found!");
        }
 
        if (normalized) {
            emissionCurve = _normalizedCepCurvePollutants.get(pollutant);
            powerPattern = _normailzedPowerPatternPollutants;
        } else {
            emissionCurve = _cepCurvePollutants.get(pollutant);
            powerPattern = _powerPatternPollutants;
        }
 
    } // end if
 
 
 
    if (emissionCurve.size() == 0) {
        throw InvalidArgument("Empty emission curve for " + pollutant + " found!");
    }
 
    if (emissionCurve.size() == 1) {
        return emissionCurve[0];
    }
 
    // in case that the demanded power is smaller than the first entry (smallest) in the power pattern the first two entries are extrapolated
    if (power <= powerPattern.front()) {
        double calcEmission =  PHEMCEP::Interpolate(power, powerPattern[0], powerPattern[1], emissionCurve[0], emissionCurve[1]);
 
        if (calcEmission < 0) {
            return 0;
        } else {
            return calcEmission;
        }
 
    } // end if
 
    // if power bigger than all entries in power pattern the last two values are linearly extrapolated
    if (power >= powerPattern.back()) {
        return PHEMCEP::Interpolate(power, powerPattern[powerPattern.size() - 2], powerPattern.back(), emissionCurve[emissionCurve.size() - 2], emissionCurve.back());
    } // end if
 
    // bisection search to find correct position in power pattern
    int upperIndex;
    int lowerIndex;
 
    PHEMCEP::FindLowerUpperInPattern(lowerIndex, upperIndex, powerPattern, power);
 
    return PHEMCEP::Interpolate(power, powerPattern[lowerIndex], powerPattern[upperIndex], emissionCurve[lowerIndex], emissionCurve[upperIndex]);
 
} // end of GetEmission
 
 
double
PHEMCEP::Interpolate(double px, double p1, double p2, double e1, double e2) const {
    if (p2 == p1) {
        return e1;
    }
    return e1 + (px - p1) / (p2 - p1) * (e2 - e1);
} // end of Interpolate
 
 
double PHEMCEP::GetDecelCoast(double speed, double acc, double gradient, double /* vehicleLoading */) const {
    if (speed < SPEED_DCEL_MIN) {
        return speed / SPEED_DCEL_MIN * GetDecelCoast(SPEED_DCEL_MIN, acc, gradient, _vehicleLoading); // !!!vehicleLoading
    } // end if
 
    double rotCoeff = GetRotationalCoeffecient(speed);
 
    int upperIndex;
    int lowerIndex;
 
    double iGear = GetGearCoeffecient(speed);
 
    double iTot = iGear * _axleRatio;
 
    double n = (30 * speed * iTot) / ((_effictiveWheelDiameter / 2) * M_PI2);
    double nNorm = (n - _engineIdlingSpeed) / (_engineRatedSpeed - _engineIdlingSpeed);
 
    FindLowerUpperInPattern(lowerIndex, upperIndex, _nNormTable, nNorm);
 
    double fMot = 0;
 
    if (speed >= 10e-2) {
        fMot = (-GetDragCoeffecient(nNorm) * _ratedPower * 1000 / speed) / 0.9;
    } // end if
 
    double fRoll = (_resistanceF0
                    + _resistanceF1 * speed
                    + pow(_resistanceF2 * speed, 2)
                    + pow(_resistanceF3 * speed, 3)
                    + pow(_resistanceF4 * speed, 4)) * (_massVehicle + _vehicleLoading) * GRAVITY_CONST; // !!!vehicleLoading
 
    double fAir = _cdValue * _crossSectionalArea * 1.2 * 0.5 * pow(speed, 2);
 
    double fGrad = (_massVehicle + _vehicleLoading) * GRAVITY_CONST * gradient / 100; // !!!vehicleLoading
 
    return -(fMot + fRoll + fAir + fGrad) / ((_massVehicle + _vehicleLoading) * rotCoeff); // !!!vehicleLoading
} // end of GetDecelCoast
 
 
double
PHEMCEP::GetRotationalCoeffecient(double speed) const {
    int upperIndex;
    int lowerIndex;
 
    PHEMCEP::FindLowerUpperInPattern(lowerIndex, upperIndex, _speedPatternRotational, speed);
 
    return PHEMCEP::Interpolate(speed,
                                _speedPatternRotational[lowerIndex],
                                _speedPatternRotational[upperIndex],
                                _speedCurveRotational[lowerIndex],
                                _speedCurveRotational[upperIndex]);
} // end of GetRotationalCoeffecient
 
double PHEMCEP::GetGearCoeffecient(double speed) const {
    int upperIndex;
    int lowerIndex;
 
    FindLowerUpperInPattern(lowerIndex, upperIndex, _gearTransmissionCurve, speed);
 
    return Interpolate(speed,
                       _speedPatternRotational[lowerIndex],
                       _speedPatternRotational[upperIndex],
                       _gearTransmissionCurve[lowerIndex],
                       _gearTransmissionCurve[upperIndex]);
} // end of GetGearCoefficient
 
double PHEMCEP::GetDragCoeffecient(double nNorm)  const {
    int upperIndex;
    int lowerIndex;
 
    FindLowerUpperInPattern(lowerIndex, upperIndex, _dragNormTable, nNorm);
 
    return Interpolate(nNorm,
                       _nNormTable[lowerIndex],
                       _nNormTable[upperIndex],
                       _dragNormTable[lowerIndex],
                       _dragNormTable[upperIndex]);
} // end of GetGearCoefficient
 
void PHEMCEP::FindLowerUpperInPattern(int& lowerIndex, int& upperIndex, const std::vector<double>& pattern, double value) const {
    if (value <= pattern.front()) {
        lowerIndex = 0;
        upperIndex = 0;
        return;
 
    } // end if
 
    if (value >= pattern.back()) {
        lowerIndex = (int)pattern.size() - 1;
        upperIndex = (int)pattern.size() - 1;
        return;
    } // end if
 
    // bisection search to find correct position in power pattern
    int middleIndex = ((int)pattern.size() - 1) / 2;
    upperIndex = (int)pattern.size() - 1;
    lowerIndex = 0;
 
    while (upperIndex - lowerIndex > 1) {
        if (pattern[middleIndex] == value) {
            lowerIndex = middleIndex;
            upperIndex = middleIndex;
            return;
        } else if (pattern[middleIndex] < value) {
            lowerIndex = middleIndex;
            middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
        } else {
            upperIndex = middleIndex;
            middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
        } // end if
    } // end while
 
    if (pattern[lowerIndex] <= value && value < pattern[upperIndex]) {
        return;
    } else {
        throw ProcessError("Error during calculation of position in pattern!");
    }
} // end of FindLowerUpperInPattern
 
 
double
PHEMCEP::CalcPower(double v, double a, double slope, double /* vehicleLoading */) const {
    const double rotFactor = GetRotationalCoeffecient(v);
    double power = (_massVehicle + _vehicleLoading) * GRAVITY_CONST * (_resistanceF0 + _resistanceF1 * v + _resistanceF4 * pow(v, 4)) * v;
    power += (_crossSectionalArea * _cdValue * AIR_DENSITY_CONST / 2) * pow(v, 3);
    power += (_massVehicle * rotFactor + _massRot + _vehicleLoading) * a * v;
    power += (_massVehicle + _vehicleLoading) * slope * 0.01 * v;
    return power / 950.;
}
 
 
double
PHEMCEP::GetMaxAccel(double v, double a, double gradient, double /* vehicleLoading */) const {
    UNUSED_PARAMETER(a);
    double rotFactor = GetRotationalCoeffecient(v);
    const double pMaxForAcc = GetPMaxNorm(v) * _ratedPower - PHEMCEP::CalcPower(v, 0, gradient, _vehicleLoading); // !!!vehicleLoading
    return (pMaxForAcc * 1000) / ((_massVehicle * rotFactor + _massRot + _vehicleLoading) * v); // !!!vehicleLoading
}
 
 
 
double PHEMCEP::GetPMaxNorm(double speed) const {
    // Linear function between v0 and v1, constant elsewhere
    if (speed <= _pNormV0) {
        return _pNormP0;
    } else if (speed >= _pNormV1) {
        return _pNormP1;
    } else {
        return PHEMCEP::Interpolate(speed, _pNormV0, _pNormV1, _pNormP0, _pNormP1);
    }
} // end of GetPMaxNorm
 
 
/****************************************************************************/