Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2016-2024 German Aerospace Center (DLR) and others.
4 : // PHEMlight module
5 : // Copyright (C) 2016-2017 Technische Universitaet Graz, https://www.tugraz.at/
6 : // This program and the accompanying materials are made available under the
7 : // terms of the Eclipse Public License 2.0 which is available at
8 : // https://www.eclipse.org/legal/epl-2.0/
9 : // This Source Code may also be made available under the following Secondary
10 : // Licenses when the conditions for such availability set forth in the Eclipse
11 : // Public License 2.0 are satisfied: GNU General Public License, version 2
12 : // or later which is available at
13 : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
14 : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
15 : /****************************************************************************/
16 : /// @file CEP.cpp
17 : /// @author Martin Dippold
18 : /// @author Michael Behrisch
19 : /// @date July 2016
20 : ///
21 : //
22 : /****************************************************************************/
23 : #include <config.h>
24 :
25 : #include "CEP.h"
26 : #include "Constants.h"
27 : #include "Helpers.h"
28 :
29 :
30 : namespace PHEMlightdll {
31 :
32 109 : CEP::CEP(bool heavyVehicle, double vehicleMass, double vehicleLoading, double vehicleMassRot, double crossArea, double cWValue, double f0, double f1, double f2, double f3, double f4, double axleRatio, std::vector<double>& transmissionGearRatios, double auxPower, double ratedPower, double engineIdlingSpeed, double engineRatedSpeed, double effictiveWheelDiameter, double pNormV0, double pNormP0, double pNormV1, double pNormP1, const std::string& vehicelFuelType, std::vector<std::vector<double> >& matrixFC, std::vector<std::string>& headerLinePollutants, std::vector<std::vector<double> >& matrixPollutants, std::vector<std::vector<double> >& matrixSpeedRotational, std::vector<std::vector<double> >& normedDragTable, double idlingFC, std::vector<double>& idlingPollutants) {
33 : (void)transmissionGearRatios; // just to make the compiler happy about the unused parameter
34 109 : InitializeInstanceFields();
35 109 : _resistanceF0 = f0;
36 109 : _resistanceF1 = f1;
37 109 : _resistanceF2 = f2;
38 109 : _resistanceF3 = f3;
39 109 : _resistanceF4 = f4;
40 109 : _cWValue = cWValue;
41 109 : _crossSectionalArea = crossArea;
42 109 : _massVehicle = vehicleMass;
43 109 : _vehicleLoading = vehicleLoading;
44 109 : _vehicleMassRot = vehicleMassRot;
45 109 : _ratedPower = ratedPower;
46 109 : _engineIdlingSpeed = engineIdlingSpeed;
47 109 : _engineRatedSpeed = engineRatedSpeed;
48 109 : _effectiveWheelDiameter = effictiveWheelDiameter;
49 109 : _heavyVehicle = heavyVehicle;
50 109 : _fuelType = vehicelFuelType;
51 109 : _axleRatio = axleRatio;
52 109 : _auxPower = auxPower;
53 :
54 109 : _pNormV0 = pNormV0 / 3.6;
55 109 : _pNormP0 = pNormP0;
56 109 : _pNormV1 = pNormV1 / 3.6;
57 109 : _pNormP1 = pNormP1;
58 :
59 : std::vector<std::string> pollutantIdentifier;
60 : std::vector<std::vector<double> > pollutantMeasures;
61 : std::vector<std::vector<double> > normalizedPollutantMeasures;
62 :
63 : // init pollutant identifiers
64 763 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
65 654 : pollutantIdentifier.push_back(headerLinePollutants[i]);
66 : }
67 :
68 : // initialize measures
69 763 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
70 654 : pollutantMeasures.push_back(std::vector<double>());
71 654 : normalizedPollutantMeasures.push_back(std::vector<double>());
72 : }
73 :
74 : // looping through matrix and assigning values for speed rotational table
75 109 : _speedCurveRotational = std::vector<double>();
76 109 : _speedPatternRotational = std::vector<double>();
77 109 : _gearTransmissionCurve = std::vector<double>();
78 895 : for (int i = 0; i < (int)matrixSpeedRotational.size(); i++) {
79 786 : if (matrixSpeedRotational[i].size() != 3) {
80 : return;
81 : }
82 :
83 786 : _speedPatternRotational.push_back(matrixSpeedRotational[i][0] / 3.6);
84 786 : _gearTransmissionCurve.push_back(matrixSpeedRotational[i][1]);
85 786 : _speedCurveRotational.push_back(matrixSpeedRotational[i][2]);
86 : }
87 :
88 : // looping through matrix and assigning values for drag table
89 109 : _nNormTable = std::vector<double>();
90 109 : _dragNormTable = std::vector<double>();
91 1744 : for (int i = 0; i < (int)normedDragTable.size(); i++) {
92 1635 : if (normedDragTable[i].size() != 2) {
93 : return;
94 : }
95 :
96 1635 : _nNormTable.push_back(normedDragTable[i][0]);
97 1635 : _dragNormTable.push_back(normedDragTable[i][1]);
98 : }
99 :
100 : // looping through matrix and assigning values for Fuel consumption
101 109 : _cepCurveFC = std::vector<double>();
102 109 : _normedCepCurveFC = std::vector<double>();
103 109 : _powerPatternFC = std::vector<double>();
104 109 : _normalizedPowerPatternFC = std::vector<double>();
105 1526 : for (int i = 0; i < (int)matrixFC.size(); i++) {
106 1417 : if (matrixFC[i].size() != 2) {
107 : return;
108 : }
109 :
110 1417 : _powerPatternFC.push_back(matrixFC[i][0] * _ratedPower);
111 1417 : _normalizedPowerPatternFC.push_back(matrixFC[i][0]);
112 1417 : _cepCurveFC.push_back(matrixFC[i][1] * _ratedPower);
113 1417 : _normedCepCurveFC.push_back(matrixFC[i][1]);
114 :
115 : }
116 :
117 109 : _powerPatternPollutants = std::vector<double>();
118 :
119 : double pollutantMultiplyer = 1;
120 :
121 109 : _drivingPower = _normalizingPower = CalcPower(Constants::NORMALIZING_SPEED, Constants::NORMALIZING_ACCELARATION, 0);
122 :
123 : // looping through matrix and assigning values for pollutants
124 109 : if (heavyVehicle) {
125 0 : _normalizingPower = _ratedPower;
126 0 : _normalizingType = NormalizingType_RatedPower;
127 : pollutantMultiplyer = _ratedPower;
128 : }
129 : else {
130 : _normalizingPower = _drivingPower;
131 109 : _normalizingType = NormalizingType_DrivingPower;
132 : }
133 :
134 109 : _normailzedPowerPatternPollutants = std::vector<double>();
135 :
136 109 : _cepNormalizedCurvePollutants = std::map<std::string, std::vector<double> >();
137 :
138 109 : int headerCount = (int)headerLinePollutants.size();
139 3488 : for (int i = 0; i < (int)matrixPollutants.size(); i++) {
140 27032 : for (int j = 0; j < (int)matrixPollutants[i].size(); j++) {
141 23653 : if ((int)matrixPollutants[i].size() != headerCount + 1) {
142 : return;
143 : }
144 :
145 23653 : if (j == 0) {
146 3379 : _normailzedPowerPatternPollutants.push_back(matrixPollutants[i][j]);
147 3379 : _powerPatternPollutants.push_back(matrixPollutants[i][j] * getNormalizingPower());
148 : }
149 : else {
150 20274 : pollutantMeasures[j - 1].push_back(matrixPollutants[i][j] * pollutantMultiplyer);
151 20274 : normalizedPollutantMeasures[j - 1].push_back(matrixPollutants[i][j]);
152 : }
153 : }
154 : }
155 :
156 109 : _cepCurvePollutants = std::map<std::string, std::vector<double> >();
157 109 : _idlingValuesPollutants = std::map<std::string, double>();
158 :
159 763 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
160 1308 : _cepCurvePollutants.insert(std::make_pair(pollutantIdentifier[i], pollutantMeasures[i]));
161 654 : _cepNormalizedCurvePollutants.insert(std::make_pair(pollutantIdentifier[i], normalizedPollutantMeasures[i]));
162 1308 : _idlingValuesPollutants.insert(std::make_pair(pollutantIdentifier[i], idlingPollutants[i] * pollutantMultiplyer));
163 : }
164 :
165 109 : _idlingValueFC = idlingFC * _ratedPower;
166 109 : }
167 :
168 0 : const bool& CEP::getHeavyVehicle() const {
169 0 : return _heavyVehicle;
170 : }
171 :
172 1863204 : const std::string& CEP::getFuelType() const {
173 1863204 : return _fuelType;
174 : }
175 :
176 0 : const CEP::NormalizingType& CEP::getNormalizingTypeX() const {
177 0 : return _normalizingType;
178 : }
179 :
180 0 : const double& CEP::getRatedPower() const {
181 0 : return _ratedPower;
182 : }
183 :
184 0 : void CEP::setRatedPower(const double& value) {
185 0 : _ratedPower = value;
186 0 : }
187 :
188 3379 : const double& CEP::getNormalizingPower() const {
189 3379 : return _normalizingPower;
190 : }
191 :
192 0 : const double& CEP::getDrivingPower() const {
193 0 : return _drivingPower;
194 : }
195 :
196 0 : void CEP::setDrivingPower(const double& value) {
197 0 : _drivingPower = value;
198 0 : }
199 :
200 1537995 : double CEP::CalcPower(double speed, double acc, double gradient) {
201 : // Declaration
202 : double power = 0;
203 1537995 : double rotFactor = GetRotationalCoeffecient(speed);
204 1537995 : double powerAux = (_auxPower * _ratedPower);
205 :
206 : // Calculate the power
207 1537995 : power += (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST * (_resistanceF0 + _resistanceF1 * speed + _resistanceF4 * std::pow(speed, 4)) * speed;
208 1537995 : power += (_crossSectionalArea * _cWValue * Constants::AIR_DENSITY_CONST / 2) * std::pow(speed, 3);
209 1537995 : power += (_massVehicle * rotFactor + _vehicleMassRot + _vehicleLoading) * acc * speed;
210 1537995 : power += (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST * gradient * 0.01 * speed;
211 1537995 : power /= 1000;
212 1537995 : power /= Constants::getDRIVE_TRAIN_EFFICIENCY();
213 1537995 : power += powerAux;
214 :
215 : // Return result
216 1537995 : return power;
217 : }
218 :
219 0 : double CEP::CalcEngPower(double power) {
220 0 : if (power < _powerPatternFC.front()) {
221 : return _powerPatternFC.front();
222 : }
223 0 : if (power > _powerPatternFC.back()) {
224 0 : return _powerPatternFC.back();
225 : }
226 :
227 : return power;
228 : }
229 :
230 1051968 : double CEP::GetEmission(const std::string& pollutant, double power, double speed, Helpers* VehicleClass) {
231 : //Declaration
232 : std::vector<double> emissionCurve;
233 : std::vector<double> powerPattern;
234 :
235 : // bisection search to find correct position in power pattern
236 : int upperIndex;
237 : int lowerIndex;
238 :
239 1051968 : if (_fuelType != Constants::strBEV) {
240 1051968 : if (std::abs(speed) <= Constants::ZERO_SPEED_ACCURACY) {
241 375712 : if (pollutant == "FC") {
242 93928 : return _idlingValueFC;
243 : }
244 : else {
245 281784 : if (_cepCurvePollutants.find(pollutant) == _cepCurvePollutants.end()) {
246 0 : VehicleClass->setErrMsg(std::string("Emission pollutant ") + pollutant + std::string(" not found!"));
247 0 : return 0;
248 : }
249 :
250 281784 : return _idlingValuesPollutants[pollutant];
251 : }
252 : }
253 : }
254 :
255 676256 : if (pollutant == "FC") {
256 169064 : emissionCurve = _cepCurveFC;
257 169064 : powerPattern = _powerPatternFC;
258 : }
259 : else {
260 507192 : if (_cepCurvePollutants.find(pollutant) == _cepCurvePollutants.end()) {
261 0 : VehicleClass->setErrMsg(std::string("Emission pollutant ") + pollutant + std::string(" not found!"));
262 0 : return 0;
263 : }
264 :
265 507192 : emissionCurve = _cepCurvePollutants[pollutant];
266 507192 : powerPattern = _powerPatternPollutants;
267 : }
268 :
269 676256 : if (emissionCurve.empty()) {
270 0 : VehicleClass->setErrMsg(std::string("Empty emission curve for ") + pollutant + std::string(" found!"));
271 0 : return 0;
272 : }
273 676256 : if (emissionCurve.size() == 1) {
274 0 : return emissionCurve[0];
275 : }
276 :
277 : // in case that the demanded power is smaller than the first entry (smallest) in the power pattern the first is returned (should never happen)
278 676256 : if (power <= powerPattern.front()) {
279 0 : return emissionCurve[0];
280 : }
281 :
282 : // if power bigger than all entries in power pattern return the last (should never happen)
283 676256 : if (power >= powerPattern.back()) {
284 0 : return emissionCurve.back();
285 : }
286 :
287 676256 : FindLowerUpperInPattern(lowerIndex, upperIndex, powerPattern, power);
288 676256 : return Interpolate(power, powerPattern[lowerIndex], powerPattern[upperIndex], emissionCurve[lowerIndex], emissionCurve[upperIndex]);
289 1051968 : }
290 :
291 131496 : double CEP::GetCO2Emission(double _FC, double _CO, double _HC, Helpers* VehicleClass) {
292 : //Declaration
293 : double fCBr;
294 : double fCHC = 0.866;
295 : double fCCO = 0.429;
296 : double fCCO2 = 0.273;
297 :
298 : //C# TO C++ CONVERTER NOTE: The following 'switch' operated on a string variable and was converted to C++ 'if-else' logic:
299 : // switch (_fuelType)
300 : //ORIGINAL LINE: case Constants.strGasoline:
301 131496 : if (_fuelType == Constants::strGasoline) {
302 : fCBr = 0.865;
303 : }
304 : //ORIGINAL LINE: case Constants.strDiesel:
305 5000 : else if (_fuelType == Constants::strDiesel) {
306 : fCBr = 0.863;
307 : }
308 : //ORIGINAL LINE: case Constants.strCNG:
309 0 : else if (_fuelType == Constants::strCNG) {
310 : fCBr = 0.693;
311 : fCHC = 0.803;
312 : }
313 : //ORIGINAL LINE: case Constants.strLPG:
314 0 : else if (_fuelType == Constants::strLPG) {
315 : fCBr = 0.825;
316 : fCHC = 0.825;
317 : }
318 : else {
319 0 : VehicleClass->setErrMsg(std::string("The propolsion type is not known! (") + _fuelType + std::string(")"));
320 0 : return 0;
321 : }
322 :
323 131496 : return (_FC * fCBr - _CO * fCCO - _HC * fCHC) / fCCO2;
324 : }
325 :
326 1317306 : double CEP::GetDecelCoast(double speed, double acc, double gradient) {
327 : //Declaration
328 : int upperIndex;
329 : int lowerIndex;
330 :
331 1317306 : if (speed < Constants::SPEED_DCEL_MIN) {
332 372216 : return speed / Constants::SPEED_DCEL_MIN * GetDecelCoast(Constants::SPEED_DCEL_MIN, acc, gradient);
333 : }
334 :
335 945090 : double rotCoeff = GetRotationalCoeffecient(speed);
336 945090 : FindLowerUpperInPattern(lowerIndex, upperIndex, _speedPatternRotational, speed);
337 945090 : double iGear = Interpolate(speed, _speedPatternRotational[lowerIndex], _speedPatternRotational[upperIndex], _gearTransmissionCurve[lowerIndex], _gearTransmissionCurve[upperIndex]);
338 :
339 945090 : double iTot = iGear * _axleRatio;
340 :
341 945090 : double n = (30 * speed * iTot) / ((_effectiveWheelDiameter / 2) * M_PI);
342 945090 : double nNorm = (n - _engineIdlingSpeed) / (_engineRatedSpeed - _engineIdlingSpeed);
343 :
344 945090 : FindLowerUpperInPattern(lowerIndex, upperIndex, _nNormTable, nNorm);
345 :
346 : double fMot = 0;
347 :
348 945090 : if (speed >= 10e-2) {
349 945090 : fMot = (-Interpolate(nNorm, _nNormTable[lowerIndex], _nNormTable[upperIndex], _dragNormTable[lowerIndex], _dragNormTable[upperIndex]) * _ratedPower * 1000 / speed) / 0.9;
350 : }
351 :
352 945090 : double fRoll = (_resistanceF0 + _resistanceF1 * speed + std::pow(_resistanceF2 * speed, 2) + std::pow(_resistanceF3 * speed, 3) + std::pow(_resistanceF4 * speed, 4)) * (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST;
353 :
354 945090 : double fAir = _cWValue * _crossSectionalArea * 1.2 * 0.5 * std::pow(speed, 2);
355 :
356 945090 : double fGrad = (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST * gradient / 100;
357 :
358 945090 : return -(fMot + fRoll + fAir + fGrad) / ((_massVehicle + _vehicleLoading) * rotCoeff);
359 : }
360 :
361 3100499 : double CEP::GetRotationalCoeffecient(double speed) {
362 : //Declaration
363 : int upperIndex;
364 : int lowerIndex;
365 :
366 3100499 : FindLowerUpperInPattern(lowerIndex, upperIndex, _speedPatternRotational, speed);
367 3100499 : return Interpolate(speed, _speedPatternRotational[lowerIndex], _speedPatternRotational[upperIndex], _speedCurveRotational[lowerIndex], _speedCurveRotational[upperIndex]);
368 : }
369 :
370 5666935 : void CEP::FindLowerUpperInPattern(int& lowerIndex, int& upperIndex, std::vector<double>& pattern, double value) {
371 5666935 : lowerIndex = 0;
372 5666935 : upperIndex = 0;
373 :
374 5666935 : if (value <= pattern.front()) {
375 : lowerIndex = 0;
376 : upperIndex = 0;
377 : return;
378 : }
379 :
380 5341617 : if (value >= pattern.back()) {
381 0 : lowerIndex = (int)pattern.size() - 1;
382 0 : upperIndex = (int)pattern.size() - 1;
383 0 : return;
384 : }
385 :
386 : // bisection search to find correct position in power pattern
387 5341617 : int middleIndex = ((int)pattern.size() - 1) / 2;
388 5341617 : upperIndex = (int)pattern.size() - 1;
389 5341617 : lowerIndex = 0;
390 :
391 22422462 : while (upperIndex - lowerIndex > 1) {
392 17080845 : if (pattern[middleIndex] == value) {
393 0 : lowerIndex = middleIndex;
394 0 : upperIndex = middleIndex;
395 0 : return;
396 : }
397 17080845 : else if (pattern[middleIndex] < value) {
398 8186366 : lowerIndex = middleIndex;
399 8186366 : middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
400 : }
401 : else {
402 8894479 : upperIndex = middleIndex;
403 8894479 : middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
404 : }
405 : }
406 :
407 : if (pattern[lowerIndex] <= value && value < pattern[upperIndex]) {
408 : return;
409 : }
410 : }
411 :
412 6228734 : double CEP::Interpolate(double px, double p1, double p2, double e1, double e2) {
413 6228734 : if (p2 == p1) {
414 : return e1;
415 : }
416 :
417 5903416 : return e1 + (px - p1) / (p2 - p1) * (e2 - e1);
418 : }
419 :
420 617414 : double CEP::GetMaxAccel(double speed, double gradient) {
421 617414 : double rotFactor = GetRotationalCoeffecient(speed);
422 617414 : double pMaxForAcc = GetPMaxNorm(speed) * _ratedPower - CalcPower(speed, 0, gradient);
423 :
424 617414 : return (pMaxForAcc * 1000) / ((_massVehicle * rotFactor + _vehicleMassRot + _vehicleLoading) * speed);
425 : }
426 :
427 617414 : double CEP::GetPMaxNorm(double speed) {
428 : // Linear function between v0 and v1, constant elsewhere
429 617414 : if (speed <= _pNormV0) {
430 14546 : return _pNormP0;
431 : }
432 602868 : else if (speed >= _pNormV1) {
433 41069 : return _pNormP1;
434 : }
435 : else {
436 561799 : return Interpolate(speed, _pNormV0, _pNormV1, _pNormP0, _pNormP1);
437 : }
438 : }
439 :
440 109 : void CEP::InitializeInstanceFields() {
441 109 : _heavyVehicle = false;
442 109 : _normalizingType = static_cast<NormalizingType>(0);
443 109 : _ratedPower = 0;
444 109 : _normalizingPower = 0;
445 109 : _drivingPower = 0;
446 109 : _massVehicle = 0;
447 109 : _vehicleLoading = 0;
448 109 : _vehicleMassRot = 0;
449 109 : _crossSectionalArea = 0;
450 109 : _cWValue = 0;
451 109 : _resistanceF0 = 0;
452 109 : _resistanceF1 = 0;
453 109 : _resistanceF2 = 0;
454 109 : _resistanceF3 = 0;
455 109 : _resistanceF4 = 0;
456 109 : _axleRatio = 0;
457 109 : _auxPower = 0;
458 109 : _pNormV0 = 0;
459 109 : _pNormP0 = 0;
460 109 : _pNormV1 = 0;
461 109 : _pNormP1 = 0;
462 109 : _engineRatedSpeed = 0;
463 109 : _engineIdlingSpeed = 0;
464 109 : _effectiveWheelDiameter = 0;
465 109 : _idlingValueFC = 0;
466 109 : }
467 : }
|
