Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2013-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 HelpersPHEMlight5.cpp
15 : /// @author Daniel Krajzewicz
16 : /// @author Michael Behrisch
17 : /// @author Nikolaus Furian
18 : /// @date Sat, 20.04.2013
19 : ///
20 : // Helper methods for PHEMlight-based emission computation
21 : /****************************************************************************/
22 : #include <config.h>
23 :
24 : #include <limits>
25 : #include <cmath>
26 : #include <foreign/PHEMlight/V5/cpp/Constants.h>
27 : #include <foreign/PHEMlight/V5/cpp/Correction.h>
28 : #include <utils/common/StringUtils.h>
29 : #include <utils/geom/GeomHelper.h>
30 : #include <utils/options/OptionsCont.h>
31 :
32 : #include "EnergyParams.h"
33 : #include "HelpersPHEMlight5.h"
34 :
35 :
36 : // ===========================================================================
37 : // method definitions
38 : // ===========================================================================
39 57416 : HelpersPHEMlight5::HelpersPHEMlight5() :
40 : HelpersPHEMlight("PHEMlight5", PHEMLIGHT5_BASE, -1),
41 114832 : myIndex(PHEMLIGHT5_BASE) {
42 57416 : }
43 :
44 :
45 57416 : HelpersPHEMlight5::~HelpersPHEMlight5() {
46 57450 : for (const auto& cep : myCEPs) {
47 34 : delete cep.second;
48 : }
49 114832 : }
50 :
51 :
52 : SUMOEmissionClass
53 44 : HelpersPHEMlight5::getClassByName(const std::string& eClass, const SUMOVehicleClass vc) {
54 44 : if (eClass == "unknown" && !myEmissionClassStrings.hasString("unknown")) {
55 0 : myEmissionClassStrings.addAlias("unknown", getClassByName("PC_EU4_G", vc));
56 : }
57 44 : if (eClass == "default" && !myEmissionClassStrings.hasString("default")) {
58 0 : myEmissionClassStrings.addAlias("default", getClassByName("PC_EU4_G", vc));
59 : }
60 : if (myEmissionClassStrings.hasString(eClass)) {
61 10 : return myEmissionClassStrings.get(eClass);
62 : }
63 34 : if (eClass.size() < 6) {
64 0 : throw InvalidArgument("Unknown emission class '" + eClass + "'.");
65 : }
66 34 : const OptionsCont& oc = OptionsCont::getOptions();
67 34 : myVolumetricFuel = oc.getBool("emissions.volumetric-fuel");
68 : std::vector<std::string> phemPath;
69 34 : phemPath.push_back(oc.getString("phemlight-path") + "/");
70 34 : if (getenv("PHEMLIGHT_PATH") != nullptr) {
71 0 : phemPath.push_back(std::string(getenv("PHEMLIGHT_PATH")) + "/");
72 : }
73 34 : if (getenv("SUMO_HOME") != nullptr) {
74 68 : phemPath.push_back(std::string(getenv("SUMO_HOME")) + "/data/emissions/PHEMlight5/");
75 : }
76 102 : if (myCorrection == nullptr && (!oc.isDefault("phemlight-year") || !oc.isDefault("phemlight-temperature"))) {
77 4 : myCorrection = new PHEMlightdllV5::Correction(phemPath);
78 8 : if (!oc.isDefault("phemlight-year")) {
79 4 : myCorrection->setYear(oc.getInt("phemlight-year"));
80 : std::string err;
81 2 : if (!myCorrection->ReadDet(err)) {
82 0 : throw InvalidArgument("Error reading PHEMlight5 deterioration data.\n" + err);
83 : }
84 2 : myCorrection->setUseDet(true);
85 : }
86 8 : if (!oc.isDefault("phemlight-temperature")) {
87 4 : myCorrection->setAmbTemp(oc.getFloat("phemlight-temperature"));
88 : std::string err;
89 2 : if (!myCorrection->ReadTNOx(err)) {
90 0 : throw InvalidArgument("Error reading PHEMlight5 deterioration data.\n" + err);
91 : }
92 2 : myCorrection->setUseTNOx(true);
93 : }
94 : }
95 34 : myHelper.setCommentPrefix("c");
96 34 : myHelper.setPHEMDataV("V5");
97 34 : myHelper.setclass(eClass);
98 34 : if (!myCEPHandler.GetCEP(phemPath, &myHelper, myCorrection)) {
99 0 : throw InvalidArgument("File for PHEMlight5 emission class " + eClass + " not found.\n" + myHelper.getErrMsg());
100 : }
101 34 : PHEMlightdllV5::CEP* const currCep = myCEPHandler.getCEPS().find(myHelper.getgClass())->second;
102 34 : int index = myIndex++;
103 34 : if (currCep->getHeavyVehicle()) {
104 0 : index |= PollutantsInterface::HEAVY_BIT;
105 : }
106 68 : myEmissionClassStrings.insert(eClass, index);
107 34 : myCEPs[index] = currCep;
108 34 : myEmissionClassStrings.addAlias(StringUtils::to_lower_case(eClass), index);
109 34 : return index;
110 34 : }
111 :
112 :
113 : std::string
114 0 : HelpersPHEMlight5::getFuel(const SUMOEmissionClass c) const {
115 0 : const std::string name = myEmissionClassStrings.getString(c);
116 0 : std::string fuel = "Gasoline";
117 0 : if (name.find("_D_") != std::string::npos) {
118 : fuel = "Diesel";
119 : }
120 0 : if (name.find("_BEV_") != std::string::npos) {
121 : fuel = "Electricity";
122 : }
123 0 : if (name.find("_HEV") != std::string::npos) {
124 0 : fuel = "Hybrid" + fuel;
125 : }
126 0 : return fuel;
127 : }
128 :
129 :
130 : double
131 23 : HelpersPHEMlight5::getWeight(const SUMOEmissionClass c) const {
132 23 : return myCEPs.find(c)->second->getVehicleMass();
133 : }
134 :
135 :
136 : double
137 76292 : HelpersPHEMlight5::getEmission(PHEMlightdllV5::CEP* currCep, const std::string& e, const double p, const double v, const double drivingPower, const double ratedPower) const {
138 76292 : return currCep->GetEmission(e, p, v, &myHelper, drivingPower, ratedPower);
139 : }
140 :
141 :
142 : double
143 203583 : HelpersPHEMlight5::calcPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const {
144 : // copy of CEP::CalcPower
145 203583 : const double power = calcWheelPower(currCep, v, a, slope, param) / PHEMlightdllV5::Constants::_DRIVE_TRAIN_EFFICIENCY;
146 203583 : if (!(currCep->getCalcType() == "HEV" || currCep->getCalcType() == "BEV")) {
147 203583 : return power + param->getDoubleOptional(SUMO_ATTR_CONSTANTPOWERINTAKE, currCep->getAuxPower() * 1000.) / 1000.;
148 : }
149 : return power;
150 : }
151 :
152 :
153 : double
154 203583 : HelpersPHEMlight5::calcWheelPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const {
155 : // copy of CEP::CalcWheelPower
156 203583 : const double rotFactor = currCep->GetRotationalCoeffecient(v);
157 203583 : const double mass = param->getDoubleOptional(SUMO_ATTR_MASS, currCep->getVehicleMass());
158 203583 : const double massRot = param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, currCep->getVehicleMassRot());
159 203583 : const double load = param->getDoubleOptional(SUMO_ATTR_LOADING, currCep->getVehicleLoading()) + param->getTransportableMass();
160 203583 : const double cw = param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, currCep->getCrossSectionalArea()) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, currCep->getCWValue());
161 203583 : const double rf0 = param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, currCep->getResistanceF0());
162 :
163 203583 : double power = (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * currCep->getResistance(v, rf0) * v;
164 203583 : power += (cw * PHEMlightdllV5::Constants::AIR_DENSITY_CONST / 2) * std::pow(v, 3);
165 203583 : power += (mass * rotFactor + massRot + load) * a * v;
166 203583 : power += (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * slope * 0.01 * v;
167 203583 : return power / 1000.;
168 : }
169 :
170 :
171 : double
172 76338 : HelpersPHEMlight5::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const {
173 76338 : PHEMlightdllV5::CEP* currCep = myCEPs.count(c) == 0 ? nullptr : myCEPs.find(c)->second;
174 76338 : if (currCep != nullptr) {
175 76338 : if (v == 0.) {
176 : return 0.;
177 : }
178 : // this is a copy of CEP::GetMaxAccel
179 60455 : const double rotFactor = currCep->GetRotationalCoeffecient(v);
180 60455 : const double mass = param->getDoubleOptional(SUMO_ATTR_MASS, currCep->getVehicleMass());
181 60455 : const double massRot = param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, currCep->getVehicleMassRot());
182 60455 : const double load = param->getDoubleOptional(SUMO_ATTR_LOADING, currCep->getVehicleLoading()) + param->getTransportableMass();
183 60455 : const double ratedPower = param->getDoubleOptional(SUMO_ATTR_MAXIMUMPOWER, 1000. * currCep->getRatedPower()) / 1000.;
184 60455 : const double pMaxForAcc = currCep->GetPMaxNorm(v) * ratedPower - calcPower(currCep, v, 0, slope, param);
185 60455 : const double maxAcc = (pMaxForAcc * 1000) / ((mass * rotFactor + massRot + load) * v);
186 : return MIN2(a, maxAcc);
187 : }
188 : return a;
189 : }
190 :
191 :
192 : double
193 109012 : HelpersPHEMlight5::getCoastingDecel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const {
194 109012 : PHEMlightdllV5::CEP* const currCep = myCEPs.find(c)->second;
195 : // this is a copy of CEP::GetDecelCoast with an adapted slope force calculation
196 109012 : if (v < PHEMlightdllV5::Constants::SPEED_DCEL_MIN) {
197 23214 : return v / PHEMlightdllV5::Constants::SPEED_DCEL_MIN * getCoastingDecel(c, PHEMlightdllV5::Constants::SPEED_DCEL_MIN, a, slope, param);
198 : }
199 85798 : const double rotFactor = currCep->GetRotationalCoeffecient(v);
200 85798 : const double mass = param->getDoubleOptional(SUMO_ATTR_MASS, currCep->getVehicleMass());
201 85798 : const double load = param->getDoubleOptional(SUMO_ATTR_LOADING, currCep->getVehicleLoading()) + param->getTransportableMass();
202 85798 : const double cw = param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, currCep->getCrossSectionalArea()) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, currCep->getCWValue());
203 85798 : const double ratedPower = param->getDoubleOptional(SUMO_ATTR_MAXIMUMPOWER, 1000. * currCep->getRatedPower()) / 1000.;
204 85798 : const double wheelRadius = param->getDoubleOptional(SUMO_ATTR_WHEELRADIUS, currCep->getWheelRadius());
205 85798 : const double rf0 = param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, currCep->getResistanceF0());
206 :
207 85798 : const double fRoll = currCep->getResistance(v, rf0) * (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST;
208 85798 : const double fAir = cw * PHEMlightdllV5::Constants::AIR_DENSITY_CONST * 0.5 * std::pow(v, 2);
209 85798 : const double fGrad = (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * sin(DEG2RAD(slope));
210 :
211 85798 : return -(currCep->getFMot(v, ratedPower, wheelRadius) + fRoll + fAir + fGrad) / ((mass + load) * rotFactor);
212 : }
213 :
214 :
215 : double
216 76338 : HelpersPHEMlight5::compute(const SUMOEmissionClass c, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
217 76338 : if (param != nullptr && param->isEngineOff()) {
218 : return 0.;
219 : }
220 : const double corrSpeed = MAX2(0.0, v);
221 : assert(myCEPs.count(c) == 1);
222 76338 : PHEMlightdllV5::CEP* const currCep = myCEPs.find(c)->second;
223 76338 : const double corrAcc = getModifiedAccel(c, corrSpeed, a, slope, param);
224 76338 : const bool isBEV = currCep->getFuelType() == PHEMlightdllV5::Constants::strBEV;
225 76338 : const bool isHybrid = currCep->getCalcType() == PHEMlightdllV5::Constants::strHybrid;
226 76338 : const double power_raw = calcPower(currCep, corrSpeed, corrAcc, slope, param);
227 76338 : const double ratedPower = param->getDoubleOptional(SUMO_ATTR_MAXIMUMPOWER, 1000. * currCep->getRatedPower()) / 1000.;
228 76338 : const double power = isHybrid ? calcWheelPower(currCep, corrSpeed, corrAcc, slope, param) : currCep->CalcEngPower(power_raw, ratedPower);
229 :
230 76338 : if (!isBEV && corrAcc < getCoastingDecel(c, corrSpeed, corrAcc, slope, param) &&
231 10227 : corrSpeed > PHEMlightdllV5::Constants::ZERO_SPEED_ACCURACY) {
232 : return 0.;
233 : }
234 : // TODO: this is probably only needed for non-heavy vehicles, so if execution speed becomes an issue this could be optimized out
235 66790 : const double drivingPower = calcPower(currCep, PHEMlightdllV5::Constants::NORMALIZING_SPEED, PHEMlightdllV5::Constants::NORMALIZING_ACCELARATION, 0, param);
236 66790 : switch (e) {
237 : case PollutantsInterface::CO:
238 9502 : return getEmission(currCep, "CO", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
239 9502 : case PollutantsInterface::CO2:
240 19004 : return currCep->GetCO2Emission(getEmission(currCep, "FC", power, corrSpeed, drivingPower, ratedPower),
241 : getEmission(currCep, "CO", power, corrSpeed, drivingPower, ratedPower),
242 9502 : getEmission(currCep, "HC", power, corrSpeed, drivingPower, ratedPower), &myHelper) / SECONDS_PER_HOUR * 1000.;
243 : case PollutantsInterface::HC:
244 9502 : return getEmission(currCep, "HC", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
245 : case PollutantsInterface::NO_X:
246 9502 : return getEmission(currCep, "NOx", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
247 : case PollutantsInterface::PM_X:
248 9502 : return getEmission(currCep, "PM", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
249 9778 : case PollutantsInterface::FUEL: {
250 9778 : if (myVolumetricFuel && currCep->getFuelType() == PHEMlightdllV5::Constants::strDiesel) { // divide by average diesel density of 836 g/l
251 0 : return getEmission(currCep, "FC", power, corrSpeed, drivingPower, ratedPower) / 836. / SECONDS_PER_HOUR * 1000.;
252 : }
253 9778 : if (myVolumetricFuel && currCep->getFuelType() == PHEMlightdllV5::Constants::strGasoline) { // divide by average gasoline density of 742 g/l
254 0 : return getEmission(currCep, "FC", power, corrSpeed, drivingPower, ratedPower) / 742. / SECONDS_PER_HOUR * 1000.;
255 : }
256 9778 : if (isBEV) {
257 : return 0.;
258 : }
259 9778 : return getEmission(currCep, "FC", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.; // still in mg even if myVolumetricFuel is set!
260 : }
261 9502 : case PollutantsInterface::ELEC:
262 9502 : if (isBEV) {
263 0 : const double auxPower = param->getDoubleOptional(SUMO_ATTR_CONSTANTPOWERINTAKE, currCep->getAuxPower() * 1000.) / 1000.;
264 0 : return (getEmission(currCep, "FC_el", power, corrSpeed, drivingPower, ratedPower) + auxPower) / SECONDS_PER_HOUR * 1000.;
265 : }
266 : return 0;
267 : }
268 : // should never get here
269 : return 0.;
270 : }
271 :
272 :
273 : /****************************************************************************/
|
