Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2016-2025 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 121 : 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 121 : InitializeInstanceFields();
35 121 : _resistanceF0 = f0;
36 121 : _resistanceF1 = f1;
37 121 : _resistanceF2 = f2;
38 121 : _resistanceF3 = f3;
39 121 : _resistanceF4 = f4;
40 121 : _cWValue = cWValue;
41 121 : _crossSectionalArea = crossArea;
42 121 : _massVehicle = vehicleMass;
43 121 : _vehicleLoading = vehicleLoading;
44 121 : _vehicleMassRot = vehicleMassRot;
45 121 : _ratedPower = ratedPower;
46 121 : _engineIdlingSpeed = engineIdlingSpeed;
47 121 : _engineRatedSpeed = engineRatedSpeed;
48 121 : _effectiveWheelDiameter = effictiveWheelDiameter;
49 121 : _heavyVehicle = heavyVehicle;
50 121 : _fuelType = vehicelFuelType;
51 121 : _axleRatio = axleRatio;
52 121 : _auxPower = auxPower;
53 :
54 121 : _pNormV0 = pNormV0 / 3.6;
55 121 : _pNormP0 = pNormP0;
56 121 : _pNormV1 = pNormV1 / 3.6;
57 121 : _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 847 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
65 726 : pollutantIdentifier.push_back(headerLinePollutants[i]);
66 : }
67 :
68 : // initialize measures
69 847 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
70 726 : pollutantMeasures.push_back(std::vector<double>());
71 726 : normalizedPollutantMeasures.push_back(std::vector<double>());
72 : }
73 :
74 : // looping through matrix and assigning values for speed rotational table
75 121 : _speedCurveRotational = std::vector<double>();
76 121 : _speedPatternRotational = std::vector<double>();
77 121 : _gearTransmissionCurve = std::vector<double>();
78 991 : for (int i = 0; i < (int)matrixSpeedRotational.size(); i++) {
79 870 : if (matrixSpeedRotational[i].size() != 3) {
80 : return;
81 : }
82 :
83 870 : _speedPatternRotational.push_back(matrixSpeedRotational[i][0] / 3.6);
84 870 : _gearTransmissionCurve.push_back(matrixSpeedRotational[i][1]);
85 870 : _speedCurveRotational.push_back(matrixSpeedRotational[i][2]);
86 : }
87 :
88 : // looping through matrix and assigning values for drag table
89 121 : _nNormTable = std::vector<double>();
90 121 : _dragNormTable = std::vector<double>();
91 1936 : for (int i = 0; i < (int)normedDragTable.size(); i++) {
92 1815 : if (normedDragTable[i].size() != 2) {
93 : return;
94 : }
95 :
96 1815 : _nNormTable.push_back(normedDragTable[i][0]);
97 1815 : _dragNormTable.push_back(normedDragTable[i][1]);
98 : }
99 :
100 : // looping through matrix and assigning values for Fuel consumption
101 121 : _cepCurveFC = std::vector<double>();
102 121 : _normedCepCurveFC = std::vector<double>();
103 121 : _powerPatternFC = std::vector<double>();
104 121 : _normalizedPowerPatternFC = std::vector<double>();
105 1694 : for (int i = 0; i < (int)matrixFC.size(); i++) {
106 1573 : if (matrixFC[i].size() != 2) {
107 : return;
108 : }
109 :
110 1573 : _powerPatternFC.push_back(matrixFC[i][0] * _ratedPower);
111 1573 : _normalizedPowerPatternFC.push_back(matrixFC[i][0]);
112 1573 : _cepCurveFC.push_back(matrixFC[i][1] * _ratedPower);
113 1573 : _normedCepCurveFC.push_back(matrixFC[i][1]);
114 :
115 : }
116 :
117 121 : _powerPatternPollutants = std::vector<double>();
118 :
119 : double pollutantMultiplyer = 1;
120 :
121 121 : _drivingPower = _normalizingPower = CalcPower(Constants::NORMALIZING_SPEED, Constants::NORMALIZING_ACCELARATION, 0);
122 :
123 : // looping through matrix and assigning values for pollutants
124 121 : if (heavyVehicle) {
125 0 : _normalizingPower = _ratedPower;
126 0 : _normalizingType = NormalizingType_RatedPower;
127 : pollutantMultiplyer = _ratedPower;
128 : }
129 : else {
130 : _normalizingPower = _drivingPower;
131 121 : _normalizingType = NormalizingType_DrivingPower;
132 : }
133 :
134 121 : _normailzedPowerPatternPollutants = std::vector<double>();
135 :
136 121 : _cepNormalizedCurvePollutants = std::map<std::string, std::vector<double> >();
137 :
138 121 : int headerCount = (int)headerLinePollutants.size();
139 3872 : for (int i = 0; i < (int)matrixPollutants.size(); i++) {
140 30008 : for (int j = 0; j < (int)matrixPollutants[i].size(); j++) {
141 26257 : if ((int)matrixPollutants[i].size() != headerCount + 1) {
142 : return;
143 : }
144 :
145 26257 : if (j == 0) {
146 3751 : _normailzedPowerPatternPollutants.push_back(matrixPollutants[i][j]);
147 3751 : _powerPatternPollutants.push_back(matrixPollutants[i][j] * getNormalizingPower());
148 : }
149 : else {
150 22506 : pollutantMeasures[j - 1].push_back(matrixPollutants[i][j] * pollutantMultiplyer);
151 22506 : normalizedPollutantMeasures[j - 1].push_back(matrixPollutants[i][j]);
152 : }
153 : }
154 : }
155 :
156 121 : _cepCurvePollutants = std::map<std::string, std::vector<double> >();
157 121 : _idlingValuesPollutants = std::map<std::string, double>();
158 :
159 847 : for (int i = 0; i < (int)headerLinePollutants.size(); i++) {
160 726 : _cepCurvePollutants.insert(std::make_pair(pollutantIdentifier[i], pollutantMeasures[i]));
161 726 : _cepNormalizedCurvePollutants.insert(std::make_pair(pollutantIdentifier[i], normalizedPollutantMeasures[i]));
162 1452 : _idlingValuesPollutants.insert(std::make_pair(pollutantIdentifier[i], idlingPollutants[i] * pollutantMultiplyer));
163 : }
164 :
165 121 : _idlingValueFC = idlingFC * _ratedPower;
166 121 : }
167 :
168 0 : const bool& CEP::getHeavyVehicle() const {
169 0 : return _heavyVehicle;
170 : }
171 :
172 1876420 : const std::string& CEP::getFuelType() const {
173 1876420 : 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 3751 : const double& CEP::getNormalizingPower() const {
189 3751 : 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 1550439 : double CEP::CalcPower(double speed, double acc, double gradient) {
201 : // Declaration
202 : double power = 0;
203 1550439 : double rotFactor = GetRotationalCoeffecient(speed);
204 1550439 : double powerAux = (_auxPower * _ratedPower);
205 :
206 : // Calculate the power
207 1550439 : power += (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST * (_resistanceF0 + _resistanceF1 * speed + _resistanceF4 * std::pow(speed, 4)) * speed;
208 1550439 : power += (_crossSectionalArea * _cWValue * Constants::AIR_DENSITY_CONST / 2) * std::pow(speed, 3);
209 1550439 : power += (_massVehicle * rotFactor + _vehicleMassRot + _vehicleLoading) * acc * speed;
210 1550439 : power += (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST * gradient * 0.01 * speed;
211 1550439 : power /= 1000;
212 1550439 : power /= Constants::getDRIVE_TRAIN_EFFICIENCY();
213 1550439 : power += powerAux;
214 :
215 : // Return result
216 1550439 : 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 1058240 : 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 1058240 : if (_fuelType != Constants::strBEV) {
240 1058240 : if (std::abs(speed) <= Constants::ZERO_SPEED_ACCURACY) {
241 376800 : if (pollutant == "FC") {
242 94200 : return _idlingValueFC;
243 : }
244 : else {
245 282600 : 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 282600 : return _idlingValuesPollutants[pollutant];
251 : }
252 : }
253 : }
254 :
255 681440 : if (pollutant == "FC") {
256 170360 : emissionCurve = _cepCurveFC;
257 170360 : powerPattern = _powerPatternFC;
258 : }
259 : else {
260 511080 : 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 511080 : emissionCurve = _cepCurvePollutants[pollutant];
266 511080 : powerPattern = _powerPatternPollutants;
267 : }
268 :
269 681440 : if (emissionCurve.empty()) {
270 0 : VehicleClass->setErrMsg(std::string("Empty emission curve for ") + pollutant + std::string(" found!"));
271 0 : return 0;
272 : }
273 681440 : 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 681440 : 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 681440 : if (power >= powerPattern.back()) {
284 0 : return emissionCurve.back();
285 : }
286 :
287 681440 : FindLowerUpperInPattern(lowerIndex, upperIndex, powerPattern, power);
288 681440 : return Interpolate(power, powerPattern[lowerIndex], powerPattern[upperIndex], emissionCurve[lowerIndex], emissionCurve[upperIndex]);
289 1058240 : }
290 :
291 132280 : 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 132280 : 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 132280 : return (_FC * fCBr - _CO * fCCO - _HC * fCHC) / fCCO2;
324 : }
325 :
326 1326014 : double CEP::GetDecelCoast(double speed, double acc, double gradient) {
327 : //Declaration
328 : int upperIndex;
329 : int lowerIndex;
330 :
331 1326014 : if (speed < Constants::SPEED_DCEL_MIN) {
332 373196 : return speed / Constants::SPEED_DCEL_MIN * GetDecelCoast(Constants::SPEED_DCEL_MIN, acc, gradient);
333 : }
334 :
335 952818 : double rotCoeff = GetRotationalCoeffecient(speed);
336 952818 : FindLowerUpperInPattern(lowerIndex, upperIndex, _speedPatternRotational, speed);
337 952818 : double iGear = Interpolate(speed, _speedPatternRotational[lowerIndex], _speedPatternRotational[upperIndex], _gearTransmissionCurve[lowerIndex], _gearTransmissionCurve[upperIndex]);
338 :
339 952818 : double iTot = iGear * _axleRatio;
340 :
341 952818 : double n = (30 * speed * iTot) / ((_effectiveWheelDiameter / 2) * M_PI);
342 952818 : double nNorm = (n - _engineIdlingSpeed) / (_engineRatedSpeed - _engineIdlingSpeed);
343 :
344 952818 : FindLowerUpperInPattern(lowerIndex, upperIndex, _nNormTable, nNorm);
345 :
346 : double fMot = 0;
347 :
348 952818 : if (speed >= 10e-2) {
349 952818 : fMot = (-Interpolate(nNorm, _nNormTable[lowerIndex], _nNormTable[upperIndex], _dragNormTable[lowerIndex], _dragNormTable[upperIndex]) * _ratedPower * 1000 / speed) / 0.9;
350 : }
351 :
352 952818 : 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 952818 : double fAir = _cWValue * _crossSectionalArea * 1.2 * 0.5 * std::pow(speed, 2);
355 :
356 952818 : double fGrad = (_massVehicle + _vehicleLoading) * Constants::GRAVITY_CONST * gradient / 100;
357 :
358 952818 : return -(fMot + fRoll + fAir + fGrad) / ((_massVehicle + _vehicleLoading) * rotCoeff);
359 : }
360 :
361 3127615 : double CEP::GetRotationalCoeffecient(double speed) {
362 : //Declaration
363 : int upperIndex;
364 : int lowerIndex;
365 :
366 3127615 : FindLowerUpperInPattern(lowerIndex, upperIndex, _speedPatternRotational, speed);
367 3127615 : return Interpolate(speed, _speedPatternRotational[lowerIndex], _speedPatternRotational[upperIndex], _speedCurveRotational[lowerIndex], _speedCurveRotational[upperIndex]);
368 : }
369 :
370 5714691 : void CEP::FindLowerUpperInPattern(int& lowerIndex, int& upperIndex, std::vector<double>& pattern, double value) {
371 5714691 : lowerIndex = 0;
372 5714691 : upperIndex = 0;
373 :
374 5714691 : if (value <= pattern.front()) {
375 : lowerIndex = 0;
376 : upperIndex = 0;
377 : return;
378 : }
379 :
380 5388589 : 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 5388589 : int middleIndex = ((int)pattern.size() - 1) / 2;
388 5388589 : upperIndex = (int)pattern.size() - 1;
389 5388589 : lowerIndex = 0;
390 :
391 22593782 : while (upperIndex - lowerIndex > 1) {
392 17205193 : if (pattern[middleIndex] == value) {
393 0 : lowerIndex = middleIndex;
394 0 : upperIndex = middleIndex;
395 0 : return;
396 : }
397 17205193 : else if (pattern[middleIndex] < value) {
398 8244426 : lowerIndex = middleIndex;
399 8244426 : middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
400 : }
401 : else {
402 8960767 : upperIndex = middleIndex;
403 8960767 : middleIndex = (upperIndex - lowerIndex) / 2 + lowerIndex;
404 : }
405 : }
406 :
407 : if (pattern[lowerIndex] <= value && value < pattern[upperIndex]) {
408 : return;
409 : }
410 : }
411 :
412 6276742 : double CEP::Interpolate(double px, double p1, double p2, double e1, double e2) {
413 6276742 : if (p2 == p1) {
414 : return e1;
415 : }
416 :
417 5950640 : return e1 + (px - p1) / (p2 - p1) * (e2 - e1);
418 : }
419 :
420 624358 : double CEP::GetMaxAccel(double speed, double gradient) {
421 624358 : double rotFactor = GetRotationalCoeffecient(speed);
422 624358 : double pMaxForAcc = GetPMaxNorm(speed) * _ratedPower - CalcPower(speed, 0, gradient);
423 :
424 624358 : return (pMaxForAcc * 1000) / ((_massVehicle * rotFactor + _vehicleMassRot + _vehicleLoading) * speed);
425 : }
426 :
427 624358 : double CEP::GetPMaxNorm(double speed) {
428 : // Linear function between v0 and v1, constant elsewhere
429 624358 : if (speed <= _pNormV0) {
430 14714 : return _pNormP0;
431 : }
432 609644 : else if (speed >= _pNormV1) {
433 47593 : return _pNormP1;
434 : }
435 : else {
436 562051 : return Interpolate(speed, _pNormV0, _pNormV1, _pNormP0, _pNormP1);
437 : }
438 : }
439 :
440 121 : void CEP::InitializeInstanceFields() {
441 121 : _heavyVehicle = false;
442 121 : _normalizingType = static_cast<NormalizingType>(0);
443 121 : _ratedPower = 0;
444 121 : _normalizingPower = 0;
445 121 : _drivingPower = 0;
446 121 : _massVehicle = 0;
447 121 : _vehicleLoading = 0;
448 121 : _vehicleMassRot = 0;
449 121 : _crossSectionalArea = 0;
450 121 : _cWValue = 0;
451 121 : _resistanceF0 = 0;
452 121 : _resistanceF1 = 0;
453 121 : _resistanceF2 = 0;
454 121 : _resistanceF3 = 0;
455 121 : _resistanceF4 = 0;
456 121 : _axleRatio = 0;
457 121 : _auxPower = 0;
458 121 : _pNormV0 = 0;
459 121 : _pNormP0 = 0;
460 121 : _pNormV1 = 0;
461 121 : _pNormP1 = 0;
462 121 : _engineRatedSpeed = 0;
463 121 : _engineIdlingSpeed = 0;
464 121 : _effectiveWheelDiameter = 0;
465 121 : _idlingValueFC = 0;
466 121 : }
467 : }
|
