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 PHEMCEP.cpp
15 : /// @author Nikolaus Furian
16 : /// @author Daniel Krajzewicz
17 : /// @author Jakob Erdmann
18 : /// @author Michael Behrisch
19 : /// @author Marek Heinrich
20 : /// @date Thu, 13.06.2013
21 : ///
22 : // Helper class for PHEM Light, holds a specific CEP for a PHEM emission class
23 : /****************************************************************************/
24 : #include <config.h>
25 :
26 : #include <cmath>
27 : #include <string>
28 : #include <utils/common/StdDefs.h>
29 : #include <utils/common/StringBijection.h>
30 : #include <utils/common/UtilExceptions.h>
31 : #include "PHEMCEP.h"
32 :
33 :
34 : // ===========================================================================
35 : // method definitions
36 : // ===========================================================================
37 0 : PHEMCEP::PHEMCEP(bool heavyVehicle, SUMOEmissionClass emissionClass, const std::string& emissionClassIdentifier,
38 : double vehicleMass, double vehicleLoading, double vehicleMassRot,
39 : double crossArea, double cdValue,
40 : double f0, double f1, double f2, double f3, double f4,
41 : double ratedPower, double pNormV0, double pNormP0, double pNormV1, double pNormP1,
42 : double axleRatio, double engineIdlingSpeed, double engineRatedSpeed, double effectiveWheelDiameter,
43 : double idlingFC,
44 : const std::string& vehicleFuelType,
45 : const std::vector< std::vector<double> >& matrixFC,
46 : const std::vector<std::string>& headerLinePollutants,
47 : const std::vector< std::vector<double> >& matrixPollutants,
48 : const std::vector< std::vector<double> >& matrixSpeedRotational,
49 : const std::vector< std::vector<double> >& normedDragTable,
50 0 : const std::vector<double>& idlingValuesPollutants) {
51 0 : _emissionClass = emissionClass;
52 0 : _resistanceF0 = f0;
53 0 : _resistanceF1 = f1;
54 0 : _resistanceF2 = f2;
55 0 : _resistanceF3 = f3;
56 0 : _resistanceF4 = f4;
57 0 : _cdValue = cdValue;
58 0 : _crossSectionalArea = crossArea;
59 0 : _massVehicle = vehicleMass;
60 0 : _vehicleLoading = vehicleLoading;
61 0 : _massRot = vehicleMassRot;
62 0 : _ratedPower = ratedPower;
63 : _vehicleFuelType = vehicleFuelType;
64 :
65 0 : _pNormV0 = pNormV0 / 3.6;
66 0 : _pNormP0 = pNormP0;
67 0 : _pNormV1 = pNormV1 / 3.6;
68 0 : _pNormP1 = pNormP1;
69 :
70 0 : _axleRatio = axleRatio;
71 0 : _engineIdlingSpeed = engineIdlingSpeed;
72 0 : _engineRatedSpeed = engineRatedSpeed;
73 0 : _effictiveWheelDiameter = effectiveWheelDiameter;
74 :
75 0 : _heavyVehicle = heavyVehicle;
76 0 : _idlingFC = idlingFC;
77 :
78 : std::vector<std::string> pollutantIdentifier;
79 : std::vector< std::vector<double> > pollutantMeasures;
80 : std::vector<std::vector<double> > normalizedPollutantMeasures;
81 :
82 : // init pollutant identifiers
83 0 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
84 0 : pollutantIdentifier.push_back(headerLinePollutants[i]);
85 : } // end for
86 :
87 : // get size of powerPatterns
88 0 : _sizeOfPatternFC = (int)matrixFC.size();
89 0 : _sizeOfPatternPollutants = (int)matrixPollutants.size();
90 :
91 : // initialize measures
92 0 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
93 0 : pollutantMeasures.push_back(std::vector<double>());
94 0 : normalizedPollutantMeasures.push_back(std::vector<double>());
95 : } // end for
96 :
97 : // looping through matrix and assigning values for speed rotational table
98 : _speedCurveRotational.clear();
99 : _speedPatternRotational.clear();
100 : _gearTransmissionCurve.clear();
101 0 : for (int i = 0; i < (int)matrixSpeedRotational.size(); i++) {
102 0 : if (matrixSpeedRotational[i].size() != 3) {
103 0 : throw InvalidArgument("Error loading vehicle file for: " + emissionClassIdentifier);
104 : }
105 :
106 0 : _speedPatternRotational.push_back(matrixSpeedRotational[i][0] / 3.6);
107 0 : _speedCurveRotational.push_back(matrixSpeedRotational[i][1]);
108 0 : _gearTransmissionCurve.push_back(matrixSpeedRotational[i][2]);
109 : } // end for
110 :
111 : // looping through matrix and assigning values for drag table
112 : _nNormTable.clear();
113 : _dragNormTable.clear();
114 0 : for (int i = 0; i < (int) normedDragTable.size(); i++) {
115 0 : if (normedDragTable[i].size() != 2) {
116 : return;
117 : }
118 :
119 0 : _nNormTable.push_back(normedDragTable[i][0]);
120 0 : _dragNormTable.push_back(normedDragTable[i][1]);
121 :
122 : } // end for
123 :
124 : // looping through matrix and assigning values for Fuel consumption
125 : _cepCurveFC.clear();
126 : _powerPatternFC.clear();
127 : _normalizedPowerPatternFC.clear();
128 : _normedCepCurveFC.clear();
129 0 : for (int i = 0; i < (int)matrixFC.size(); i++) {
130 0 : if (matrixFC[i].size() != 2) {
131 0 : throw InvalidArgument("Error loading vehicle file for: " + emissionClassIdentifier);
132 : }
133 :
134 0 : _powerPatternFC.push_back(matrixFC[i][0] * _ratedPower);
135 0 : _normalizedPowerPatternFC.push_back(matrixFC[i][0]);
136 0 : _cepCurveFC.push_back(matrixFC[i][1] * _ratedPower);
137 0 : _normedCepCurveFC.push_back(matrixFC[i][1]);
138 :
139 : } // end for
140 :
141 : _powerPatternPollutants.clear();
142 : double pollutantMultiplyer = 1;
143 :
144 0 : _drivingPower = _normalizingPower = CalcPower(NORMALIZING_SPEED, NORMALIZING_ACCELARATION, 0, vehicleLoading);
145 :
146 : // looping through matrix and assigning values for pollutants
147 :
148 0 : if (heavyVehicle) {
149 0 : _normalizingPower = _ratedPower;
150 : pollutantMultiplyer = _ratedPower;
151 0 : _normalizingType = RatedPower;
152 : } else {
153 : _normalizingPower = _drivingPower;
154 0 : _normalizingType = DrivingPower;
155 : } // end if
156 :
157 0 : const int headerCount = (int)headerLinePollutants.size();
158 0 : for (int i = 0; i < (int)matrixPollutants.size(); i++) {
159 0 : for (int j = 0; j < (int)matrixPollutants[i].size(); j++) {
160 0 : if ((int)matrixPollutants[i].size() != headerCount + 1) {
161 : return;
162 : }
163 :
164 0 : if (j == 0) {
165 0 : _normailzedPowerPatternPollutants.push_back(matrixPollutants[i][j]);
166 0 : _powerPatternPollutants.push_back(matrixPollutants[i][j] * _normalizingPower);
167 : } else {
168 0 : pollutantMeasures[j - 1].push_back(matrixPollutants[i][j] * pollutantMultiplyer);
169 0 : normalizedPollutantMeasures[j - 1].push_back(matrixPollutants[i][j]);
170 : } // end if
171 : } // end for
172 : } // end for
173 :
174 0 : for (int i = 0; i < (int) headerLinePollutants.size(); i++) {
175 0 : _cepCurvePollutants.insert(pollutantIdentifier[i], pollutantMeasures[i]);
176 0 : _normalizedCepCurvePollutants.insert(pollutantIdentifier[i], normalizedPollutantMeasures[i]);
177 0 : _idlingValuesPollutants.insert(pollutantIdentifier[i], idlingValuesPollutants[i] * pollutantMultiplyer);
178 : } // end for
179 :
180 0 : _idlingFC = idlingFC * _ratedPower;
181 :
182 0 : } // end of Cep
183 :
184 :
185 0 : PHEMCEP::~PHEMCEP() {
186 : // free power pattern
187 : _powerPatternFC.clear();
188 : _powerPatternPollutants.clear();
189 : _cepCurveFC.clear();
190 : _speedCurveRotational.clear();
191 : _speedPatternRotational.clear();
192 0 : } // end of ~Cep
193 :
194 :
195 : double
196 0 : PHEMCEP::GetEmission(const std::string& pollutant, double power, double speed, bool normalized) const {
197 : std::vector<double> emissionCurve;
198 : std::vector<double> powerPattern;
199 :
200 0 : if (!normalized && fabs(speed) <= ZERO_SPEED_ACCURACY) {
201 0 : if (pollutant == "FC") {
202 0 : return _idlingFC;
203 : } else {
204 0 : return _idlingValuesPollutants.get(pollutant);
205 : }
206 : } // end if
207 :
208 0 : if (pollutant == "FC") {
209 0 : if (normalized) {
210 0 : emissionCurve = _normedCepCurveFC;
211 0 : powerPattern = _normalizedPowerPatternFC;
212 : } else {
213 0 : emissionCurve = _cepCurveFC;
214 0 : powerPattern = _powerPatternFC;
215 : }
216 : } else {
217 0 : if (!_cepCurvePollutants.hasString(pollutant)) {
218 0 : throw InvalidArgument("Emission pollutant " + pollutant + " not found!");
219 : }
220 :
221 0 : if (normalized) {
222 0 : emissionCurve = _normalizedCepCurvePollutants.get(pollutant);
223 0 : powerPattern = _normailzedPowerPatternPollutants;
224 : } else {
225 0 : emissionCurve = _cepCurvePollutants.get(pollutant);
226 0 : powerPattern = _powerPatternPollutants;
227 : }
228 :
229 : } // end if
230 :
231 :
232 :
233 0 : if (emissionCurve.size() == 0) {
234 0 : throw InvalidArgument("Empty emission curve for " + pollutant + " found!");
235 : }
236 :
237 0 : if (emissionCurve.size() == 1) {
238 0 : return emissionCurve[0];
239 : }
240 :
241 : // in case that the demanded power is smaller than the first entry (smallest) in the power pattern the first two entries are extrapolated
242 0 : if (power <= powerPattern.front()) {
243 0 : double calcEmission = PHEMCEP::Interpolate(power, powerPattern[0], powerPattern[1], emissionCurve[0], emissionCurve[1]);
244 :
245 0 : if (calcEmission < 0) {
246 : return 0;
247 : } else {
248 0 : return calcEmission;
249 : }
250 :
251 : } // end if
252 :
253 : // if power bigger than all entries in power pattern the last two values are linearly extrapolated
254 0 : if (power >= powerPattern.back()) {
255 0 : return PHEMCEP::Interpolate(power, powerPattern[powerPattern.size() - 2], powerPattern.back(), emissionCurve[emissionCurve.size() - 2], emissionCurve.back());
256 : } // end if
257 :
258 : // bisection search to find correct position in power pattern
259 : int upperIndex;
260 : int lowerIndex;
261 :
262 0 : PHEMCEP::FindLowerUpperInPattern(lowerIndex, upperIndex, powerPattern, power);
263 :
264 0 : return PHEMCEP::Interpolate(power, powerPattern[lowerIndex], powerPattern[upperIndex], emissionCurve[lowerIndex], emissionCurve[upperIndex]);
265 :
266 : } // end of GetEmission
267 :
268 :
269 : double
270 0 : PHEMCEP::Interpolate(double px, double p1, double p2, double e1, double e2) const {
271 0 : if (p2 == p1) {
272 : return e1;
273 : }
274 0 : return e1 + (px - p1) / (p2 - p1) * (e2 - e1);
275 : } // end of Interpolate
276 :
277 :
278 0 : double PHEMCEP::GetDecelCoast(double speed, double acc, double gradient, double /* vehicleLoading */) const {
279 0 : if (speed < SPEED_DCEL_MIN) {
280 0 : return speed / SPEED_DCEL_MIN * GetDecelCoast(SPEED_DCEL_MIN, acc, gradient, _vehicleLoading); // !!!vehicleLoading
281 : } // end if
282 :
283 0 : double rotCoeff = GetRotationalCoeffecient(speed);
284 :
285 : int upperIndex;
286 : int lowerIndex;
287 :
288 0 : double iGear = GetGearCoeffecient(speed);
289 :
290 0 : double iTot = iGear * _axleRatio;
291 :
292 0 : double n = (30 * speed * iTot) / ((_effictiveWheelDiameter / 2) * M_PI2);
293 0 : double nNorm = (n - _engineIdlingSpeed) / (_engineRatedSpeed - _engineIdlingSpeed);
294 :
295 0 : FindLowerUpperInPattern(lowerIndex, upperIndex, _nNormTable, nNorm);
296 :
297 : double fMot = 0;
298 :
299 0 : if (speed >= 10e-2) {
300 0 : fMot = (-GetDragCoeffecient(nNorm) * _ratedPower * 1000 / speed) / 0.9;
301 : } // end if
302 :
303 0 : double fRoll = (_resistanceF0
304 0 : + _resistanceF1 * speed
305 0 : + pow(_resistanceF2 * speed, 2)
306 0 : + pow(_resistanceF3 * speed, 3)
307 0 : + pow(_resistanceF4 * speed, 4)) * (_massVehicle + _vehicleLoading) * GRAVITY_CONST; // !!!vehicleLoading
308 :
309 0 : double fAir = _cdValue * _crossSectionalArea * 1.2 * 0.5 * pow(speed, 2);
310 :
311 0 : double fGrad = (_massVehicle + _vehicleLoading) * GRAVITY_CONST * gradient / 100; // !!!vehicleLoading
312 :
313 0 : return -(fMot + fRoll + fAir + fGrad) / ((_massVehicle + _vehicleLoading) * rotCoeff); // !!!vehicleLoading
314 : } // end of GetDecelCoast
315 :
316 :
317 : double
318 0 : PHEMCEP::GetRotationalCoeffecient(double speed) const {
319 : int upperIndex;
320 : int lowerIndex;
321 :
322 0 : PHEMCEP::FindLowerUpperInPattern(lowerIndex, upperIndex, _speedPatternRotational, speed);
323 :
324 0 : return PHEMCEP::Interpolate(speed,
325 : _speedPatternRotational[lowerIndex],
326 : _speedPatternRotational[upperIndex],
327 0 : _speedCurveRotational[lowerIndex],
328 0 : _speedCurveRotational[upperIndex]);
329 : } // end of GetRotationalCoeffecient
330 :
331 0 : double PHEMCEP::GetGearCoeffecient(double speed) const {
332 : int upperIndex;
333 : int lowerIndex;
334 :
335 0 : FindLowerUpperInPattern(lowerIndex, upperIndex, _gearTransmissionCurve, speed);
336 :
337 0 : return Interpolate(speed,
338 : _speedPatternRotational[lowerIndex],
339 : _speedPatternRotational[upperIndex],
340 0 : _gearTransmissionCurve[lowerIndex],
341 0 : _gearTransmissionCurve[upperIndex]);
342 : } // end of GetGearCoefficient
343 :
344 0 : double PHEMCEP::GetDragCoeffecient(double nNorm) const {
345 : int upperIndex;
346 : int lowerIndex;
347 :
348 0 : FindLowerUpperInPattern(lowerIndex, upperIndex, _dragNormTable, nNorm);
349 :
350 0 : return Interpolate(nNorm,
351 : _nNormTable[lowerIndex],
352 : _nNormTable[upperIndex],
353 0 : _dragNormTable[lowerIndex],
354 0 : _dragNormTable[upperIndex]);
355 : } // end of GetGearCoefficient
356 :
357 0 : void PHEMCEP::FindLowerUpperInPattern(int& lowerIndex, int& upperIndex, const std::vector<double>& pattern, double value) const {
358 0 : if (value <= pattern.front()) {
359 0 : lowerIndex = 0;
360 0 : upperIndex = 0;
361 0 : return;
362 :
363 : } // end if
364 :
365 0 : if (value >= pattern.back()) {
366 0 : lowerIndex = (int)pattern.size() - 1;
367 0 : upperIndex = (int)pattern.size() - 1;
368 0 : return;
369 : } // end if
370 :
371 : // bisection search to find correct position in power pattern
372 0 : int middleIndex = ((int)pattern.size() - 1) / 2;
373 0 : upperIndex = (int)pattern.size() - 1;
374 0 : lowerIndex = 0;
375 :
376 0 : while (upperIndex - lowerIndex > 1) {
377 0 : if (pattern[middleIndex] == value) {
378 0 : lowerIndex = middleIndex;
379 0 : upperIndex = middleIndex;
380 0 : return;
381 0 : } else if (pattern[middleIndex] < value) {
382 0 : lowerIndex = middleIndex;
383 0 : middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
384 : } else {
385 0 : upperIndex = middleIndex;
386 0 : middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
387 : } // end if
388 : } // end while
389 :
390 0 : if (pattern[lowerIndex] <= value && value < pattern[upperIndex]) {
391 : return;
392 : } else {
393 0 : throw ProcessError("Error during calculation of position in pattern!");
394 : }
395 : } // end of FindLowerUpperInPattern
396 :
397 :
398 : double
399 0 : PHEMCEP::CalcPower(double v, double a, double slope, double /* vehicleLoading */) const {
400 0 : const double rotFactor = GetRotationalCoeffecient(v);
401 0 : double power = (_massVehicle + _vehicleLoading) * GRAVITY_CONST * (_resistanceF0 + _resistanceF1 * v + _resistanceF4 * pow(v, 4)) * v;
402 0 : power += (_crossSectionalArea * _cdValue * AIR_DENSITY_CONST / 2) * pow(v, 3);
403 0 : power += (_massVehicle * rotFactor + _massRot + _vehicleLoading) * a * v;
404 0 : power += (_massVehicle + _vehicleLoading) * slope * 0.01 * v;
405 0 : return power / 950.;
406 : }
407 :
408 :
409 : double
410 0 : PHEMCEP::GetMaxAccel(double v, double a, double gradient, double /* vehicleLoading */) const {
411 : UNUSED_PARAMETER(a);
412 0 : double rotFactor = GetRotationalCoeffecient(v);
413 0 : const double pMaxForAcc = GetPMaxNorm(v) * _ratedPower - PHEMCEP::CalcPower(v, 0, gradient, _vehicleLoading); // !!!vehicleLoading
414 0 : return (pMaxForAcc * 1000) / ((_massVehicle * rotFactor + _massRot + _vehicleLoading) * v); // !!!vehicleLoading
415 : }
416 :
417 :
418 :
419 0 : double PHEMCEP::GetPMaxNorm(double speed) const {
420 : // Linear function between v0 and v1, constant elsewhere
421 0 : if (speed <= _pNormV0) {
422 0 : return _pNormP0;
423 0 : } else if (speed >= _pNormV1) {
424 0 : return _pNormP1;
425 : } else {
426 0 : return PHEMCEP::Interpolate(speed, _pNormV0, _pNormV1, _pNormP0, _pNormP1);
427 : }
428 : } // end of GetPMaxNorm
429 :
430 :
431 : /****************************************************************************/
|
