Line data    Source code

       1             : /****************************************************************************/
       2             : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
       3             : // Copyright (C) 2012-2024 German Aerospace Center (DLR) and others.
       4             : // This program and the accompanying materials are made available under the
       5             : // terms of the Eclipse Public License 2.0 which is available at
       6             : // https://www.eclipse.org/legal/epl-2.0/
       7             : // This Source Code may also be made available under the following Secondary
       8             : // Licenses when the conditions for such availability set forth in the Eclipse
       9             : // Public License 2.0 are satisfied: GNU General Public License, version 2
      10             : // or later which is available at
      11             : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
      12             : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
      13             : /****************************************************************************/
      14             : /// @file    MSCFModel_Rail.cpp
      15             : /// @author  Gregor L\"ammel
      16             : /// @author  Leander Flamm
      17             : /// @date    Tue, 08 Feb 2017
      18             : ///
      19             : // <description missing>
      20             : /****************************************************************************/
      21             : #include <config.h>
      22             : 
      23             : #include <iostream>
      24             : #include <utils/common/MsgHandler.h>
      25             : #include <utils/common/StringUtils.h>
      26             : #include <utils/common/StringTokenizer.h>
      27             : #include <utils/geom/GeomHelper.h>
      28             : #include <microsim/MSVehicle.h>
      29             : #include <microsim/lcmodels/MSAbstractLaneChangeModel.h>
      30             : #include "MSCFModel_Rail.h"
      31             : 
      32             : // ===========================================================================
      33             : // trainParams method definitions
      34             : // ===========================================================================
      35             : 
      36             : double
      37     2406959 : MSCFModel_Rail::TrainParams::getResistance(double speed) const {
      38     2406959 :     if (resCoef_constant != INVALID_DOUBLE) {
      39       20704 :         return (resCoef_quadratic * speed * speed + resCoef_linear * speed + resCoef_constant); // kN
      40             :     } else {
      41     2386255 :         return getInterpolatedValueFromLookUpMap(speed, &resistance); // kN
      42             :     }
      43             : }
      44             : 
      45             : 
      46             : double
      47     1527965 : MSCFModel_Rail::TrainParams::getTraction(double speed) const {
      48     1527965 :     if (maxPower != INVALID_DOUBLE) {
      49        6596 :         return MIN2(maxPower / speed, maxTraction); // kN
      50             :     } else {
      51     1521369 :         return getInterpolatedValueFromLookUpMap(speed, &traction); // kN
      52             :     }
      53             : }
      54             : 
      55             : // ===========================================================================
      56             : // method definitions
      57             : // ===========================================================================
      58             : 
      59             : 
      60         319 : MSCFModel_Rail::MSCFModel_Rail(const MSVehicleType* vtype) :
      61         319 :     MSCFModel(vtype) {
      62         319 :     const std::string trainType = vtype->getParameter().getCFParamString(SUMO_ATTR_TRAIN_TYPE, "NGT400");
      63         319 :     if (trainType.compare("RB425") == 0) {
      64          38 :         myTrainParams = initRB425Params();
      65         281 :     } else if (trainType.compare("RB628") == 0) {
      66          19 :         myTrainParams = initRB628Params();
      67         262 :     } else if (trainType.compare("NGT400") == 0) {
      68         116 :         myTrainParams = initNGT400Params();
      69         146 :     } else if (trainType.compare("NGT400_16") == 0) {
      70           4 :         myTrainParams = initNGT400_16Params();
      71         142 :     } else if (trainType.compare("ICE1") == 0) {
      72           9 :         myTrainParams = initICE1Params();
      73         133 :     } else if (trainType.compare("REDosto7") == 0) {
      74          72 :         myTrainParams = initREDosto7Params();
      75          61 :     } else if (trainType.compare("Freight") == 0) {
      76          36 :         myTrainParams = initFreightParams();
      77          25 :     } else if (trainType.compare("ICE3") == 0) {
      78           5 :         myTrainParams = initICE3Params();
      79          20 :     } else if (trainType.compare("MireoPlusB") == 0) {
      80           4 :         myTrainParams = initMireoPlusB2TParams();
      81          16 :     } else if (trainType.compare("MireoPlusH") == 0) {
      82           4 :         myTrainParams = initMireoPlusH2TParams();
      83          12 :     } else if (trainType.compare("custom") == 0) {
      84          12 :         myTrainParams = initCustomParams();
      85             :     } else {
      86           4 :         WRITE_ERRORF(TL("Unknown train type: %. Exiting!"), trainType);
      87           0 :         throw ProcessError();
      88             :     }
      89             :     // override with user values
      90         319 :     if (vtype->wasSet(VTYPEPARS_MAXSPEED_SET)) {
      91          25 :         myTrainParams.vmax = vtype->getMaxSpeed();
      92             :     }
      93         319 :     if (vtype->wasSet(VTYPEPARS_LENGTH_SET)) {
      94          38 :         myTrainParams.length = vtype->getLength();
      95             :     }
      96         319 :     if (vtype->wasSet(VTYPEPARS_MASS_SET)) {
      97             :         // kg to tons
      98           4 :         myTrainParams.weight = vtype->getMass() / 1000;
      99             :     }
     100         319 :     myTrainParams.mf = vtype->getParameter().getCFParam(SUMO_ATTR_MASSFACTOR, myTrainParams.mf);
     101         319 :     myTrainParams.decl = vtype->getParameter().getCFParam(SUMO_ATTR_DECEL, myTrainParams.decl);
     102             :     setMaxDecel(myTrainParams.decl);
     103         319 :     setEmergencyDecel(vtype->getParameter().getCFParam(SUMO_ATTR_EMERGENCYDECEL, myTrainParams.decl + 0.3));
     104             :     // update type parameters so they are shown correctly in the gui (if defaults from trainType are used)
     105         319 :     const_cast<MSVehicleType*>(vtype)->setMaxSpeed(myTrainParams.vmax);
     106         319 :     const_cast<MSVehicleType*>(vtype)->setLength(myTrainParams.length);
     107             : 
     108             :     // init tabular curves
     109         319 :     std::vector<double> speedTable = getValueTable(vtype, SUMO_ATTR_SPEED_TABLE);
     110         319 :     std::vector<double> tractionTable = getValueTable(vtype, SUMO_ATTR_TRACTION_TABLE);
     111         319 :     std::vector<double> resistanceTable = getValueTable(vtype, SUMO_ATTR_RESISTANCE_TABLE);
     112         319 :     if (speedTable.size() > 0 || tractionTable.size() > 0 || resistanceTable.size() > 0) {
     113           8 :         if (speedTable.size() == 1) {
     114           0 :             throw ProcessError(TLF("Invalid size of speedTable for vType '%' (at least 2 values are required).", vtype->getID()));
     115           8 :         } else if (speedTable.size() != tractionTable.size()) {
     116           0 :             throw ProcessError(TLF("Mismatching size of speedTable and tractionTable for vType '%'.", vtype->getID()));
     117           8 :         } else if (speedTable.size() != resistanceTable.size()) {
     118           0 :             throw ProcessError(TLF("Mismatching size of speedTable and resistanceTable for vType '%'.", vtype->getID()));
     119             :         }
     120           8 :         myTrainParams.traction.clear();
     121           8 :         myTrainParams.resistance.clear();
     122         112 :         for (int i = 0; i < (int)speedTable.size(); i++) {
     123         104 :             myTrainParams.traction[speedTable[i]] = tractionTable[i];
     124         104 :             myTrainParams.resistance[speedTable[i]] = resistanceTable[i];
     125             :         }
     126             :     }
     127             : 
     128             :     // init parametric curves
     129         319 :     myTrainParams.maxPower = vtype->getParameter().getCFParam(SUMO_ATTR_MAXPOWER, INVALID_DOUBLE);
     130         319 :     myTrainParams.maxTraction = vtype->getParameter().getCFParam(SUMO_ATTR_MAXTRACTION, INVALID_DOUBLE);
     131         319 :     myTrainParams.resCoef_constant = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_CONSTANT, INVALID_DOUBLE);
     132         319 :     myTrainParams.resCoef_linear = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_LINEAR, INVALID_DOUBLE);
     133         319 :     myTrainParams.resCoef_quadratic = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_QUADRATIC, INVALID_DOUBLE);
     134             : 
     135         319 :     if (myTrainParams.maxPower != INVALID_DOUBLE && myTrainParams.maxTraction == INVALID_DOUBLE) {
     136           0 :         throw ProcessError(TLF("Undefined maxPower for vType '%'.", vtype->getID()));
     137         319 :     } else if (myTrainParams.maxPower == INVALID_DOUBLE && myTrainParams.maxTraction != INVALID_DOUBLE) {
     138           0 :         throw ProcessError(TLF("Undefined maxTraction for vType '%'.", vtype->getID()));
     139             :     }
     140         319 :     if (myTrainParams.maxPower != INVALID_DOUBLE && tractionTable.size() > 0) {
     141           0 :         WRITE_WARNING(TLF("Ignoring tractionTable because maxPower and maxTraction are set for vType '%'.", vtype->getID()));
     142             :     }
     143         319 :     const bool hasSomeResCoef = (myTrainParams.resCoef_constant != INVALID_DOUBLE
     144         315 :             || myTrainParams.resCoef_linear != INVALID_DOUBLE
     145         634 :             || myTrainParams.resCoef_quadratic != INVALID_DOUBLE);
     146             :     const bool hasAllResCoef = (myTrainParams.resCoef_constant != INVALID_DOUBLE
     147           4 :             && myTrainParams.resCoef_linear != INVALID_DOUBLE
     148         323 :             && myTrainParams.resCoef_quadratic != INVALID_DOUBLE);
     149         319 :     if (hasSomeResCoef && !hasAllResCoef) {
     150           0 :         throw ProcessError(TLF("Some undefined resistance coefficients for vType '%' (requires resCoef_constant, resCoef_linear and resCoef_quadratic)", vtype->getID()));
     151             :     }
     152         319 :     if (myTrainParams.resCoef_constant != INVALID_DOUBLE && resistanceTable.size() > 0) {
     153           4 :         WRITE_WARNING(TLF("Ignoring resistanceTable because resistance coefficents are set for vType '%'.", vtype->getID()));
     154             :     }
     155             : 
     156         319 :     if (myTrainParams.traction.empty() && myTrainParams.maxPower == INVALID_DOUBLE) {
     157          16 :         throw ProcessError(TLF("Either tractionTable or maxPower must be defined for vType '%' with Rail model type '%'.", vtype->getID(), trainType));
     158             :     }
     159         315 :     if (myTrainParams.resistance.empty() && myTrainParams.resCoef_constant == INVALID_DOUBLE) {
     160           0 :         throw ProcessError(TLF("Either resistanceTable or resCoef_constant must be defined for vType '%' with Rail model type '%'.", vtype->getID(), trainType));
     161             :     }
     162         323 : }
     163             : 
     164             : 
     165         628 : MSCFModel_Rail::~MSCFModel_Rail() { }
     166             : 
     167             : 
     168             : std::vector<double>
     169         957 : MSCFModel_Rail::getValueTable(const MSVehicleType* vtype, SumoXMLAttr attr) {
     170             :     std::vector<double> result;
     171        1914 :     const std::string values = vtype->getParameter().getCFParamString(attr, "");
     172         957 :     if (!values.empty()) {
     173         360 :         for (std::string value : StringTokenizer(values).getVector()) {
     174         312 :             result.push_back(StringUtils::toDouble(value));
     175          24 :         }
     176             :     }
     177         957 :     return result;
     178             : }
     179             : 
     180             : 
     181      143097 : double MSCFModel_Rail::followSpeed(const MSVehicle* const veh, double speed, double gap,
     182             :                                    double /* predSpeed */, double /* predMaxDecel*/, const MSVehicle* const /*pred*/, const CalcReason /*usage*/) const {
     183             : 
     184             :     // followSpeed module is used for the simulation of moving block operations. The safety gap is chosen similar to the existing german
     185             :     // system CIR-ELKE (based on LZB). Other implementations of moving block systems may differ, but for now no appropriate parameter
     186             :     // can be set (would be per lane, not per train) -> hard-coded
     187             : 
     188             :     // @note: default train minGap of 5 is already subtracted from gap
     189      143097 :     if (speed >= 30 / 3.6) {
     190             :         // safety distance for higher speeds (>= 30 km/h)
     191       96479 :         gap = MAX2(0.0, gap + veh->getVehicleType().getMinGap() - 50);
     192             :     }
     193             : 
     194      143097 :     const double vsafe = maximumSafeStopSpeed(gap, myDecel, speed, false, TS, false); // absolute breaking distance
     195      143097 :     const double vmin = minNextSpeed(speed, veh);
     196      143097 :     const double vmax = maxNextSpeed(speed, veh);
     197             : 
     198      143097 :     if (MSGlobals::gSemiImplicitEulerUpdate) {
     199             :         return MIN2(vsafe, vmax);
     200             :     } else {
     201             :         // ballistic
     202             :         // XXX: the euler variant can break as strong as it wishes immediately! The ballistic cannot, refs. #2575.
     203             :         return MAX2(MIN2(vsafe, vmax), vmin);
     204             :     }
     205             : }
     206             : 
     207             : int
     208           0 : MSCFModel_Rail::getModelID() const {
     209           0 :     return SUMO_TAG_CF_RAIL;
     210             : }
     211             : 
     212             : MSCFModel*
     213           3 : MSCFModel_Rail::duplicate(const MSVehicleType* vtype) const {
     214           3 :     return new MSCFModel_Rail(vtype);
     215             : }
     216             : 
     217     1656554 : double MSCFModel_Rail::maxNextSpeed(double speed, const MSVehicle* const veh) const {
     218             : 
     219     1656554 :     if (speed >= myTrainParams.vmax) {
     220             :         return myTrainParams.vmax;
     221             :     }
     222             : 
     223             :     double targetSpeed = myTrainParams.vmax;
     224             : 
     225     1527965 :     double res = myTrainParams.getResistance(speed); // kN
     226             : 
     227     1527965 :     double slope = veh->getSlope();
     228     1527965 :     double gr = myTrainParams.weight * GRAVITY * sin(DEG2RAD(slope)); //kN
     229             : 
     230     1527965 :     double totalRes = res + gr; //kN
     231             : 
     232     1527965 :     double trac = myTrainParams.getTraction(speed); // kN
     233             : 
     234             :     double a;
     235     1527965 :     if (speed < targetSpeed) {
     236     1527965 :         a = (trac - totalRes) / myTrainParams.getRotWeight(); //kN/t == N/kg
     237             :     } else {
     238             :         a = 0.;
     239           0 :         if (totalRes > trac) {
     240           0 :             a = (trac - totalRes) / myTrainParams.getRotWeight(); //kN/t == N/kg
     241             :         }
     242             :     }
     243     1527965 :     double maxNextSpeed = speed + ACCEL2SPEED(a);
     244             : 
     245             : //    std::cout << veh->getID() << " speed: " << (speed*3.6) << std::endl;
     246             : 
     247     1527965 :     return maxNextSpeed;
     248             : }
     249             : 
     250      878994 : double MSCFModel_Rail::minNextSpeed(double speed, const MSVehicle* const veh) const {
     251             : 
     252      878994 :     const double slope = veh->getSlope();
     253      878994 :     const double gr = myTrainParams.weight * GRAVITY * sin(DEG2RAD(slope)); //kN
     254      878994 :     const double res = myTrainParams.getResistance(speed); // kN
     255      878994 :     const double totalRes = res + gr; //kN
     256      878994 :     const double a = myTrainParams.decl + totalRes / myTrainParams.getRotWeight();
     257      878994 :     const double vMin = speed - ACCEL2SPEED(a);
     258      878994 :     if (MSGlobals::gSemiImplicitEulerUpdate) {
     259             :         return MAX2(vMin, 0.);
     260             :     } else {
     261             :         // NOTE: ballistic update allows for negative speeds to indicate a stop within the next timestep
     262             :         return vMin;
     263             :     }
     264             : 
     265             : }
     266             : 
     267             : 
     268             : double
     269      243675 : MSCFModel_Rail::minNextSpeedEmergency(double speed, const MSVehicle* const veh) const {
     270      243675 :     return minNextSpeed(speed, veh);
     271             : }
     272             : 
     273             : 
     274     3907624 : double MSCFModel_Rail::getInterpolatedValueFromLookUpMap(double speed, const LookUpMap* lookUpMap) {
     275             :     std::map<double, double>::const_iterator low, prev;
     276             :     low = lookUpMap->lower_bound(speed);
     277             : 
     278     3907624 :     if (low == lookUpMap->end()) { //speed > max speed
     279        9160 :         return (lookUpMap->rbegin())->second;
     280             :     }
     281             : 
     282     3898464 :     if (low == lookUpMap->begin()) {
     283      602267 :         return low->second;
     284             :     }
     285             : 
     286             :     prev = low;
     287             :     --prev;
     288             : 
     289     3296197 :     double range = low->first - prev->first;
     290     3296197 :     double dist = speed - prev->first;
     291             :     assert(range > 0);
     292             :     assert(dist > 0);
     293             : 
     294     3296197 :     double weight = dist / range;
     295             : 
     296     3296197 :     double res = (1 - weight) * prev->second + weight * low->second;
     297             : 
     298     3296197 :     return res;
     299             : 
     300             : }
     301             : 
     302             : 
     303             : 
     304             : //void
     305             : //MSCFModel_Rail::initVehicleVariables(const MSVehicle *const veh, MSCFModel_Rail::VehicleVariables *pVariables) const {
     306             : //
     307             : //    pVariables->setInitialized();
     308             : //
     309             : //}
     310             : 
     311           0 : double MSCFModel_Rail::getSpeedAfterMaxDecel(double /* speed */) const {
     312             : 
     313             : //    //TODO: slope not known here
     314             : //    double gr = 0; //trainParams.weight * GRAVITY * edge.grade
     315             : //
     316             : //    double a = 0;//trainParams.decl - gr/trainParams.rotWeight;
     317             : //
     318             : //    return speed + a * DELTA_T / 1000.;
     319           0 :     WRITE_ERROR("function call not allowed for rail model. Exiting!");
     320           0 :     throw ProcessError();
     321             : }
     322             : 
     323         354 : MSCFModel::VehicleVariables* MSCFModel_Rail::createVehicleVariables() const {
     324         354 :     VehicleVariables* ret = new VehicleVariables();
     325         354 :     return ret;
     326             : }
     327             : 
     328             : 
     329      243675 : double MSCFModel_Rail::finalizeSpeed(MSVehicle* const veh, double vPos) const {
     330      243675 :     return MSCFModel::finalizeSpeed(veh, vPos);
     331             : }
     332             : 
     333      914082 : double MSCFModel_Rail::freeSpeed(const MSVehicle* const /* veh */, double /* speed */, double dist, double targetSpeed,
     334             :                                  const bool onInsertion, const CalcReason /*usage*/) const {
     335             : 
     336             : //    MSCFModel_Rail::VehicleVariables *vars = (MSCFModel_Rail::VehicleVariables *) veh->getCarFollowVariables();
     337             : //    if (vars->isNotYetInitialized()) {
     338             : //        initVehicleVariables(veh, vars);
     339             : //    }
     340             : 
     341             :     //TODO: signals, coasting, ...
     342             : 
     343      914082 :     if (MSGlobals::gSemiImplicitEulerUpdate) {
     344             :         // adapt speed to succeeding lane, no reaction time is involved
     345             :         // when breaking for y steps the following distance g is covered
     346             :         // (drive with v in the final step)
     347             :         // g = (y^2 + y) * 0.5 * b + y * v
     348             :         // y = ((((sqrt((b + 2.0*v)*(b + 2.0*v) + 8.0*b*g)) - b)*0.5 - v)/b)
     349      914081 :         const double v = SPEED2DIST(targetSpeed);
     350      914081 :         if (dist < v) {
     351             :             return targetSpeed;
     352             :         }
     353      820839 :         const double b = ACCEL2DIST(myDecel);
     354      820839 :         const double y = MAX2(0.0, ((sqrt((b + 2.0 * v) * (b + 2.0 * v) + 8.0 * b * dist) - b) * 0.5 - v) / b);
     355      820839 :         const double yFull = floor(y);
     356      820839 :         const double exactGap = (yFull * yFull + yFull) * 0.5 * b + yFull * v + (y > yFull ? v : 0.0);
     357      820839 :         const double fullSpeedGain = (yFull + (onInsertion ? 1. : 0.)) * ACCEL2SPEED(myTrainParams.decl);
     358     1093593 :         return DIST2SPEED(MAX2(0.0, dist - exactGap) / (yFull + 1)) + fullSpeedGain + targetSpeed;
     359             :     } else {
     360           1 :         WRITE_ERROR(TL("Anything else than semi implicit euler update is not yet implemented. Exiting!"));
     361           1 :         throw ProcessError();
     362             :     }
     363             : }
     364             : 
     365      972176 : double MSCFModel_Rail::stopSpeed(const MSVehicle* const veh, const double speed, double gap, double decel, const CalcReason /*usage*/) const {
     366      972176 :     return MIN2(maximumSafeStopSpeed(gap, decel, speed, false, TS, false), maxNextSpeed(speed, veh));
     367             : }
     368             : 
     369             : 
     370             : void
     371         614 : MSCFModel_Rail::convertMap(LookUpMap& map, double keyFactor, double valueFactor) {
     372             :     LookUpMap map2;
     373       18586 :     for (auto item : map) {
     374       17972 :         map2[item.first * keyFactor] = item.second * valueFactor;
     375             :     }
     376             :     map.clear();
     377         614 :     map.insert(map2.begin(), map2.end());
     378         614 : }