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