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