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