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 55625 : HelpersPHEMlight5::HelpersPHEMlight5() :
39 : HelpersPHEMlight("PHEMlight5", PHEMLIGHT5_BASE, -1),
40 111250 : myIndex(PHEMLIGHT5_BASE) {
41 55625 : }
42 :
43 :
44 55625 : HelpersPHEMlight5::~HelpersPHEMlight5() {
45 55653 : for (const auto& cep : myCEPs) {
46 28 : delete cep.second;
47 : }
48 111250 : }
49 :
50 :
51 : SUMOEmissionClass
52 34 : HelpersPHEMlight5::getClassByName(const std::string& eClass, const SUMOVehicleClass vc) {
53 34 : if (eClass == "unknown" && !myEmissionClassStrings.hasString("unknown")) {
54 0 : myEmissionClassStrings.addAlias("unknown", getClassByName("PC_EU4_G", vc));
55 : }
56 34 : if (eClass == "default" && !myEmissionClassStrings.hasString("default")) {
57 0 : myEmissionClassStrings.addAlias("default", getClassByName("PC_EU4_G", vc));
58 : }
59 : if (myEmissionClassStrings.hasString(eClass)) {
60 6 : return myEmissionClassStrings.get(eClass);
61 : }
62 28 : if (eClass.size() < 6) {
63 0 : throw InvalidArgument("Unknown emission class '" + eClass + "'.");
64 : }
65 28 : const OptionsCont& oc = OptionsCont::getOptions();
66 28 : myVolumetricFuel = oc.getBool("emissions.volumetric-fuel");
67 : std::vector<std::string> phemPath;
68 28 : phemPath.push_back(oc.getString("phemlight-path") + "/");
69 28 : if (getenv("PHEMLIGHT_PATH") != nullptr) {
70 0 : phemPath.push_back(std::string(getenv("PHEMLIGHT_PATH")) + "/");
71 : }
72 28 : if (getenv("SUMO_HOME") != nullptr) {
73 56 : phemPath.push_back(std::string(getenv("SUMO_HOME")) + "/data/emissions/PHEMlight5/");
74 : }
75 84 : 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 28 : myHelper.setCommentPrefix("c");
95 28 : myHelper.setPHEMDataV("V5");
96 28 : myHelper.setclass(eClass);
97 28 : if (!myCEPHandler.GetCEP(phemPath, &myHelper, myCorrection)) {
98 0 : throw InvalidArgument("File for PHEMlight5 emission class " + eClass + " not found.\n" + myHelper.getErrMsg());
99 : }
100 28 : PHEMlightdllV5::CEP* const currCep = myCEPHandler.getCEPS().find(myHelper.getgClass())->second;
101 28 : int index = myIndex++;
102 28 : if (currCep->getHeavyVehicle()) {
103 0 : index |= PollutantsInterface::HEAVY_BIT;
104 : }
105 56 : myEmissionClassStrings.insert(eClass, index);
106 28 : myCEPs[index] = currCep;
107 28 : myEmissionClassStrings.addAlias(StringUtils::to_lower_case(eClass), index);
108 28 : return index;
109 28 : }
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 75824 : HelpersPHEMlight5::getEmission(PHEMlightdllV5::CEP* currCep, const std::string& e, const double p, const double v, const double drivingPower, const double ratedPower) const {
131 75824 : return currCep->GetEmission(e, p, v, &myHelper, drivingPower, ratedPower);
132 : }
133 :
134 :
135 : double
136 202272 : HelpersPHEMlight5::calcPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const {
137 : // copy of CEP::CalcPower
138 202272 : const double power = calcWheelPower(currCep, v, a, slope, param) / PHEMlightdllV5::Constants::_DRIVE_TRAIN_EFFICIENCY;
139 202272 : if (!(currCep->getCalcType() == "HEV" || currCep->getCalcType() == "BEV")) {
140 202272 : return power + param->getDoubleOptional(SUMO_ATTR_CONSTANTPOWERINTAKE, currCep->getAuxPower() * 1000.) / 1000.;
141 : }
142 : return power;
143 : }
144 :
145 :
146 : double
147 202272 : HelpersPHEMlight5::calcWheelPower(PHEMlightdllV5::CEP* currCep, const double v, const double a, const double slope, const EnergyParams* param) const {
148 : // copy of CEP::CalcWheelPower
149 202272 : const double rotFactor = currCep->GetRotationalCoeffecient(v);
150 202272 : const double mass = param->getDoubleOptional(SUMO_ATTR_MASS, currCep->getVehicleMass());
151 202272 : const double massRot = param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, currCep->getVehicleMassRot());
152 202272 : const double load = param->getDoubleOptional(SUMO_ATTR_LOADING, currCep->getVehicleLoading());
153 202272 : const double cw = param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, currCep->getCrossSectionalArea()) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, currCep->getCWValue());
154 202272 : const double rf0 = param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, currCep->getResistanceF0());
155 :
156 202272 : double power = (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * currCep->getResistance(v, rf0) * v;
157 202272 : power += (cw * PHEMlightdllV5::Constants::AIR_DENSITY_CONST / 2) * std::pow(v, 3);
158 202272 : power += (mass * rotFactor + massRot + load) * a * v;
159 202272 : power += (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * slope * 0.01 * v;
160 202272 : return power / 1000.;
161 : }
162 :
163 :
164 : double
165 75866 : HelpersPHEMlight5::getModifiedAccel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const {
166 75866 : PHEMlightdllV5::CEP* currCep = myCEPs.count(c) == 0 ? nullptr : myCEPs.find(c)->second;
167 75866 : if (currCep != nullptr) {
168 75866 : if (v == 0.) {
169 : return 0.;
170 : }
171 : // this is a copy of CEP::GetMaxAccel
172 60060 : const double rotFactor = currCep->GetRotationalCoeffecient(v);
173 60060 : const double mass = param->getDoubleOptional(SUMO_ATTR_MASS, currCep->getVehicleMass());
174 60060 : const double massRot = param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, currCep->getVehicleMassRot());
175 60060 : const double load = param->getDoubleOptional(SUMO_ATTR_LOADING, currCep->getVehicleLoading());
176 60060 : const double ratedPower = param->getDoubleOptional(SUMO_ATTR_MAXIMUMPOWER, currCep->getRatedPower());
177 60060 : const double pMaxForAcc = currCep->GetPMaxNorm(v) * ratedPower - calcPower(currCep, v, 0, slope, param);
178 60060 : 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 108376 : HelpersPHEMlight5::getCoastingDecel(const SUMOEmissionClass c, const double v, const double a, const double slope, const EnergyParams* param) const {
187 108376 : PHEMlightdllV5::CEP* const currCep = myCEPs.find(c)->second;
188 : // this is a copy of CEP::GetDecelCoast
189 108376 : if (v < PHEMlightdllV5::Constants::SPEED_DCEL_MIN) {
190 23078 : return v / PHEMlightdllV5::Constants::SPEED_DCEL_MIN * getCoastingDecel(c, PHEMlightdllV5::Constants::SPEED_DCEL_MIN, a, slope, param);
191 : }
192 85298 : const double rotFactor = currCep->GetRotationalCoeffecient(v);
193 85298 : const double mass = param->getDoubleOptional(SUMO_ATTR_MASS, currCep->getVehicleMass());
194 85298 : const double load = param->getDoubleOptional(SUMO_ATTR_LOADING, currCep->getVehicleLoading());
195 85298 : const double cw = param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, currCep->getCrossSectionalArea()) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, currCep->getCWValue());
196 85298 : const double ratedPower = param->getDoubleOptional(SUMO_ATTR_MAXIMUMPOWER, currCep->getRatedPower());
197 85298 : const double wheelRadius = param->getDoubleOptional(SUMO_ATTR_WHEELRADIUS, currCep->getWheelRadius());
198 85298 : const double rf0 = param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, currCep->getResistanceF0());
199 :
200 85298 : const double fRoll = currCep->getResistance(v, rf0) * (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST;
201 85298 : const double fAir = cw * PHEMlightdllV5::Constants::AIR_DENSITY_CONST * 0.5 * std::pow(v, 2);
202 85298 : const double fGrad = (mass + load) * PHEMlightdllV5::Constants::GRAVITY_CONST * slope / 100;
203 :
204 85298 : return -(currCep->getFMot(v, ratedPower, wheelRadius) + fRoll + fAir + fGrad) / ((mass + load) * rotFactor);
205 : }
206 :
207 :
208 : double
209 75866 : HelpersPHEMlight5::compute(const SUMOEmissionClass c, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
210 75866 : if (param != nullptr && param->isEngineOff()) {
211 : return 0.;
212 : }
213 : const double corrSpeed = MAX2(0.0, v);
214 : assert(myCEPs.count(c) == 1);
215 75866 : PHEMlightdllV5::CEP* const currCep = myCEPs.find(c)->second;
216 75866 : const double corrAcc = getModifiedAccel(c, corrSpeed, a, slope, param);
217 75866 : const bool isBEV = currCep->getFuelType() == PHEMlightdllV5::Constants::strBEV;
218 75866 : const bool isHybrid = currCep->getCalcType() == PHEMlightdllV5::Constants::strHybrid;
219 75866 : const double power_raw = calcPower(currCep, corrSpeed, corrAcc, slope, param);
220 75866 : const double ratedPower = param->getDoubleOptional(SUMO_ATTR_MAXIMUMPOWER, currCep->getRatedPower());
221 75866 : const double power = isHybrid ? calcWheelPower(currCep, corrSpeed, corrAcc, slope, param) : currCep->CalcEngPower(power_raw, ratedPower);
222 :
223 75866 : if (!isBEV && corrAcc < getCoastingDecel(c, corrSpeed, corrAcc, slope, param) &&
224 10192 : 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 66346 : const double drivingPower = calcPower(currCep, PHEMlightdllV5::Constants::NORMALIZING_SPEED, PHEMlightdllV5::Constants::NORMALIZING_ACCELARATION, 0, param);
229 66346 : switch (e) {
230 : case PollutantsInterface::CO:
231 9478 : return getEmission(currCep, "CO", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
232 9478 : case PollutantsInterface::CO2:
233 18956 : return currCep->GetCO2Emission(getEmission(currCep, "FC", power, corrSpeed, drivingPower, ratedPower),
234 : getEmission(currCep, "CO", power, corrSpeed, drivingPower, ratedPower),
235 9478 : getEmission(currCep, "HC", power, corrSpeed, drivingPower, ratedPower), &myHelper) / SECONDS_PER_HOUR * 1000.;
236 : case PollutantsInterface::HC:
237 9478 : return getEmission(currCep, "HC", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
238 : case PollutantsInterface::NO_X:
239 9478 : return getEmission(currCep, "NOx", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
240 : case PollutantsInterface::PM_X:
241 9478 : return getEmission(currCep, "PM", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.;
242 9478 : case PollutantsInterface::FUEL: {
243 9478 : 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 9478 : 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 9478 : if (isBEV) {
250 : return 0.;
251 : }
252 9478 : return getEmission(currCep, "FC", power, corrSpeed, drivingPower, ratedPower) / SECONDS_PER_HOUR * 1000.; // still in mg even if myVolumetricFuel is set!
253 : }
254 9478 : case PollutantsInterface::ELEC:
255 9478 : 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 : /****************************************************************************/
|
