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