1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2012-2019 German Aerospace Center (DLR) and others.
4 // This program and the accompanying materials
5 // are made available under the terms of the Eclipse Public License v2.0
6 // which accompanies this distribution, and is available at
7 // http://www.eclipse.org/legal/epl-v20.html
8 // SPDX-License-Identifier: EPL-2.0
9 /****************************************************************************/
10 /// @file    MSCFModel_Rail.cpp
11 /// @author  Gregor L\"ammel
12 /// @author  Leander Flamm
13 /// @date    Tue, 08 Feb 2017
14 /// @version $Id$
15 ///
16 // <description missing>
17 /****************************************************************************/
18 // ===========================================================================
19 // included modules
20 // ===========================================================================
21 #include <config.h>
22 
23 #include <iostream>
24 #include <utils/common/MsgHandler.h>
25 #include <utils/geom/GeomHelper.h>
26 #include <microsim/MSVehicle.h>
27 #include <microsim/lcmodels/MSAbstractLaneChangeModel.h>
28 #include "MSCFModel_Rail.h"
29 
30 #define G  9.80665
31 
MSCFModel_Rail(const MSVehicleType * vtype)32 MSCFModel_Rail::MSCFModel_Rail(const MSVehicleType* vtype) :
33     MSCFModel(vtype) {
34     const std::string trainType = vtype->getParameter().getCFParamString(SUMO_ATTR_TRAIN_TYPE, "NGT400");
35     if (trainType.compare("RB425") == 0) {
36         myTrainParams = initRB425Params();
37     } else if (trainType.compare("RB628") == 0) {
38         myTrainParams = initRB628Params();
39     } else if (trainType.compare("NGT400") == 0) {
40         myTrainParams = initNGT400Params();
41     } else if (trainType.compare("NGT400_16") == 0) {
42         myTrainParams = initNGT400_16Params();
43     } else if (trainType.compare("ICE1") == 0) {
44         myTrainParams = initICE1Params();
45     } else if (trainType.compare("REDosto7") == 0) {
46         myTrainParams = initREDosto7Params();
47     } else if (trainType.compare("Freight") == 0) {
48         myTrainParams = initFreightParams();
49     } else if (trainType.compare("ICE3") == 0) {
50         myTrainParams = initICE3Params();
51     } else {
52         WRITE_ERROR("Unknown train type: " + trainType + ". Exiting!");
53         throw ProcessError();
54     }
55     // override with user values
56     if (vtype->wasSet(VTYPEPARS_MAXSPEED_SET)) {
57         myTrainParams.vmax = vtype->getMaxSpeed();
58     }
59     if (vtype->wasSet(VTYPEPARS_LENGTH_SET)) {
60         myTrainParams.length = vtype->getLength();
61     }
62     myTrainParams.decl = vtype->getParameter().getCFParam(SUMO_ATTR_DECEL, myTrainParams.decl);
63     setMaxDecel(myTrainParams.decl);
64     setEmergencyDecel(vtype->getParameter().getCFParam(SUMO_ATTR_EMERGENCYDECEL, myTrainParams.decl + 0.3));
65     // update type parameters so they are shown correctly in the gui (if defaults from trainType are used)
66     const_cast<MSVehicleType*>(vtype)->setMaxSpeed(myTrainParams.vmax);
67     const_cast<MSVehicleType*>(vtype)->setLength(myTrainParams.length);
68 
69 }
70 
~MSCFModel_Rail()71 MSCFModel_Rail::~MSCFModel_Rail() { }
72 
followSpeed(const MSVehicle * const veh,double speed,double gap,double,double,const MSVehicle * const) const73 double MSCFModel_Rail::followSpeed(const MSVehicle* const veh, double speed, double gap,
74                                    double /* predSpeed */, double /* predMaxDecel*/, const MSVehicle* const /*pred*/) const {
75 
76     // followSpeed module is used for the simulation of moving block operations. The safety gap is chosen similar to the existing german
77     // system CIR-ELKE (based on LZB). Other implementations of moving block systems may differ, but for now no appropriate parameter
78     // can be set (would be per lane, not per train) -> hard-coded
79 
80     double safetyGap = 5.0; // default value for low speeds (< 30 km/h)
81     if (speed >= 30 / 3.6) {
82         safetyGap = 50.0; // safety distance for higher speeds (>= 30 km/h)
83     }
84 
85     const double vsafe = maximumSafeStopSpeed(gap - safetyGap, speed, false, TS); // absolute breaking distance
86     const double vmin = minNextSpeed(speed, veh);
87     const double vmax = maxNextSpeed(speed, veh);
88 
89     if (MSGlobals::gSemiImplicitEulerUpdate) {
90         return MIN2(vsafe, vmax);
91     } else {
92         // ballistic
93         // XXX: the euler variant can break as strong as it wishes immediately! The ballistic cannot, refs. #2575.
94         return MAX2(MIN2(vsafe, vmax), vmin);
95     }
96 }
97 
98 int
getModelID() const99 MSCFModel_Rail::getModelID() const {
100     return SUMO_TAG_CF_RAIL;
101 }
102 
103 MSCFModel*
duplicate(const MSVehicleType * vtype) const104 MSCFModel_Rail::duplicate(const MSVehicleType* vtype) const {
105     return new MSCFModel_Rail(vtype);
106 }
107 
maxNextSpeed(double speed,const MSVehicle * const veh) const108 double MSCFModel_Rail::maxNextSpeed(double speed, const MSVehicle* const veh) const {
109 
110     if (speed >= myTrainParams.vmax) {
111         return myTrainParams.vmax;
112     }
113 
114     double targetSpeed = myTrainParams.vmax;
115 
116     double res = getInterpolatedValueFromLookUpMap(speed, &(myTrainParams.resistance)); // kN
117 
118     double slope = veh->getSlope();
119     double gr = myTrainParams.weight * G * sin(DEG2RAD(slope)); //kN
120 
121     double totalRes = res + gr; //kN
122 
123     double trac = getInterpolatedValueFromLookUpMap(speed, &(myTrainParams.traction)); // kN
124 
125     double a;
126     if (speed < targetSpeed) {
127         a = (trac - totalRes) / myTrainParams.rotWeight; //kN/t == N/kg
128     } else {
129         a = 0.;
130         if (totalRes > trac) {
131             a = (trac - totalRes) / myTrainParams.rotWeight; //kN/t == N/kg
132         }
133     }
134 
135     double maxNextSpeed = speed + a * DELTA_T / 1000.;
136 
137 //    std::cout << veh->getID() << " speed: " << (speed*3.6) << std::endl;
138 
139     return maxNextSpeed;
140 }
141 
minNextSpeed(double speed,const MSVehicle * const veh) const142 double MSCFModel_Rail::minNextSpeed(double speed, const MSVehicle* const veh) const {
143 
144     const double slope = veh->getSlope();
145     const double gr = myTrainParams.weight * G * sin(DEG2RAD(slope)); //kN
146     const double res = getInterpolatedValueFromLookUpMap(speed, &(myTrainParams.resistance)); // kN
147     const double totalRes = res + gr; //kN
148     const double a = myTrainParams.decl + totalRes / myTrainParams.rotWeight;
149     const double vMin = speed - a * DELTA_T / 1000.;
150     if (MSGlobals::gSemiImplicitEulerUpdate) {
151         return MAX2(vMin, 0.);
152     } else {
153         // NOTE: ballistic update allows for negative speeds to indicate a stop within the next timestep
154         return vMin;
155     }
156 
157 }
158 
159 
160 double
minNextSpeedEmergency(double speed,const MSVehicle * const veh) const161 MSCFModel_Rail::minNextSpeedEmergency(double speed, const MSVehicle* const veh) const {
162     return minNextSpeed(speed, veh);
163 }
164 
165 
getInterpolatedValueFromLookUpMap(double speed,const LookUpMap * lookUpMap) const166 double MSCFModel_Rail::getInterpolatedValueFromLookUpMap(double speed, const LookUpMap* lookUpMap) const {
167     speed = speed * 3.6; // lookup values in km/h
168     std::map<double, double>::const_iterator low, prev;
169     low = lookUpMap->lower_bound(speed);
170 
171     if (low == lookUpMap->end()) { //speed > max speed
172         return (lookUpMap->rbegin())->second;
173     }
174 
175     if (low == lookUpMap->begin()) {
176         return low->second;
177     }
178 
179     prev = low;
180     --prev;
181 
182     double range = low->first - prev->first;
183     double dist = speed - prev->first;
184     assert(range > 0);
185     assert(dist > 0);
186 
187     double weight = dist / range;
188 
189     double res = (1 - weight) * prev->second + weight * low->second;
190 
191     return res;
192 
193 }
194 
195 
196 
197 //void
198 //MSCFModel_Rail::initVehicleVariables(const MSVehicle *const veh, MSCFModel_Rail::VehicleVariables *pVariables) const {
199 //
200 //    pVariables->setInitialized();
201 //
202 //}
203 
getSpeedAfterMaxDecel(double) const204 double MSCFModel_Rail::getSpeedAfterMaxDecel(double /* speed */) const {
205 
206 //    //TODO: slope not known here
207 //    double gr = 0; //trainParams.weight * 9.81 * edge.grade
208 //
209 //    double a = 0;//trainParams.decl - gr/trainParams.rotWeight;
210 //
211 //    return speed + a * DELTA_T / 1000.;
212     WRITE_ERROR("function call not allowd for rail model. Exiting!");
213     throw ProcessError();
214 }
215 
createVehicleVariables() const216 MSCFModel::VehicleVariables* MSCFModel_Rail::createVehicleVariables() const {
217     VehicleVariables* ret = new VehicleVariables();
218     return ret;
219 }
220 
221 
finalizeSpeed(MSVehicle * const veh,double vPos) const222 double MSCFModel_Rail::finalizeSpeed(MSVehicle* const veh, double vPos) const {
223     return MSCFModel::finalizeSpeed(veh, vPos);
224 }
225 
freeSpeed(const MSVehicle * const,double,double dist,double targetSpeed,const bool onInsertion) const226 double MSCFModel_Rail::freeSpeed(const MSVehicle* const /* veh */, double /* speed */, double dist, double targetSpeed,
227                                  const bool onInsertion) const {
228 
229 //    MSCFModel_Rail::VehicleVariables *vars = (MSCFModel_Rail::VehicleVariables *) veh->getCarFollowVariables();
230 //    if (vars->isNotYetInitialized()) {
231 //        initVehicleVariables(veh, vars);
232 //    }
233 
234     //TODO: signals, coasting, ...
235 
236     if (MSGlobals::gSemiImplicitEulerUpdate) {
237         // adapt speed to succeeding lane, no reaction time is involved
238         // when breaking for y steps the following distance g is covered
239         // (drive with v in the final step)
240         // g = (y^2 + y) * 0.5 * b + y * v
241         // y = ((((sqrt((b + 2.0*v)*(b + 2.0*v) + 8.0*b*g)) - b)*0.5 - v)/b)
242         const double v = SPEED2DIST(targetSpeed);
243         if (dist < v) {
244             return targetSpeed;
245         }
246         const double b = ACCEL2DIST(myDecel);
247         const double y = MAX2(0.0, ((sqrt((b + 2.0 * v) * (b + 2.0 * v) + 8.0 * b * dist) - b) * 0.5 - v) / b);
248         const double yFull = floor(y);
249         const double exactGap = (yFull * yFull + yFull) * 0.5 * b + yFull * v + (y > yFull ? v : 0.0);
250         const double fullSpeedGain = (yFull + (onInsertion ? 1. : 0.)) * ACCEL2SPEED(myTrainParams.decl);
251         return DIST2SPEED(MAX2(0.0, dist - exactGap) / (yFull + 1)) + fullSpeedGain + targetSpeed;
252     } else {
253         WRITE_ERROR("Anything else than semi implicit euler update is not yet implemented. Exiting!");
254         throw ProcessError();
255     }
256 }
257 
stopSpeed(const MSVehicle * const veh,const double speed,double gap) const258 double MSCFModel_Rail::stopSpeed(const MSVehicle* const veh, const double speed, double gap) const {
259     return MIN2(maximumSafeStopSpeed(gap, speed, false, TS), maxNextSpeed(speed, veh));
260 }
261