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