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