Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2012-2026 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 Laemmel
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 2617707 : MSCFModel_Rail::TrainParams::getResistance(double speed) const {
38 2617707 : if (resCoef_constant != INVALID_DOUBLE) {
39 20704 : return (resCoef_quadratic * speed * speed + resCoef_linear * speed + resCoef_constant); // kN
40 : } else {
41 2597003 : return LinearApproxHelpers::getInterpolatedValue(resistance, speed); // kN
42 : }
43 : }
44 :
45 :
46 : double
47 1687048 : MSCFModel_Rail::TrainParams::getTraction(double speed) const {
48 1687048 : if (maxPower != INVALID_DOUBLE) {
49 6596 : return MIN2(maxPower / speed, maxTraction); // kN
50 : } else {
51 1680452 : return LinearApproxHelpers::getInterpolatedValue(traction, speed); // kN
52 : }
53 : }
54 :
55 : // ===========================================================================
56 : // method definitions
57 : // ===========================================================================
58 :
59 :
60 390 : MSCFModel_Rail::MSCFModel_Rail(const MSVehicleType* vtype) :
61 390 : MSCFModel(vtype) {
62 390 : const std::string trainType = vtype->getParameter().getCFParamString(SUMO_ATTR_TRAIN_TYPE, "NGT400");
63 390 : if (trainType.compare("RB425") == 0) {
64 58 : myTrainParams = initRB425Params();
65 332 : } else if (trainType.compare("RB628") == 0) {
66 19 : myTrainParams = initRB628Params();
67 313 : } else if (trainType.compare("NGT400") == 0) {
68 116 : myTrainParams = initNGT400Params();
69 197 : } else if (trainType.compare("NGT400_16") == 0) {
70 4 : myTrainParams = initNGT400_16Params();
71 193 : } else if (trainType.compare("ICE1") == 0) {
72 9 : myTrainParams = initICE1Params();
73 184 : } else if (trainType.compare("REDosto7") == 0) {
74 89 : myTrainParams = initREDosto7Params();
75 95 : } else if (trainType.compare("Freight") == 0) {
76 70 : 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 0 : WRITE_ERRORF(TL("Unknown train type: %. Exiting!"), trainType);
87 0 : throw ProcessError();
88 : }
89 : // override with user values
90 390 : if (vtype->wasSet(VTYPEPARS_MAXSPEED_SET)) {
91 36 : myTrainParams.vmax = vtype->getMaxSpeed();
92 : }
93 390 : if (vtype->wasSet(VTYPEPARS_LENGTH_SET)) {
94 67 : myTrainParams.length = vtype->getLength();
95 : }
96 390 : myTrainParams.mf = vtype->getParameter().getCFParam(SUMO_ATTR_MASSFACTOR, myTrainParams.mf);
97 390 : myTrainParams.decl = vtype->getParameter().getCFParam(SUMO_ATTR_DECEL, myTrainParams.decl);
98 : setMaxDecel(myTrainParams.decl);
99 390 : setEmergencyDecel(vtype->getParameter().getCFParam(SUMO_ATTR_EMERGENCYDECEL, myTrainParams.decl + 0.3));
100 : // update type parameters so they are shown correctly in the gui (if defaults from trainType are used)
101 390 : const_cast<MSVehicleType*>(vtype)->setMaxSpeed(myTrainParams.vmax);
102 390 : const_cast<MSVehicleType*>(vtype)->setLength(myTrainParams.length);
103 390 : if (!vtype->wasSet(VTYPEPARS_MASS_SET)) {
104 : // tons to kg
105 359 : const_cast<MSVehicleType*>(vtype)->setMass(myTrainParams.weight * 1000);
106 : }
107 :
108 : // init tabular curves
109 390 : myTrainParams.traction = vtype->getParameter().getCFProfile(SUMO_ATTR_TRACTION_TABLE, myTrainParams.traction);
110 390 : myTrainParams.resistance = vtype->getParameter().getCFProfile(SUMO_ATTR_RESISTANCE_TABLE, myTrainParams.resistance);
111 :
112 : // init parametric curves
113 390 : myTrainParams.maxPower = vtype->getParameter().getCFParam(SUMO_ATTR_MAXPOWER, INVALID_DOUBLE);
114 390 : myTrainParams.maxTraction = vtype->getParameter().getCFParam(SUMO_ATTR_MAXTRACTION, INVALID_DOUBLE);
115 390 : myTrainParams.resCoef_constant = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_CONSTANT, INVALID_DOUBLE);
116 390 : myTrainParams.resCoef_linear = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_LINEAR, INVALID_DOUBLE);
117 390 : myTrainParams.resCoef_quadratic = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_QUADRATIC, INVALID_DOUBLE);
118 : // curve resistance parameters
119 390 : myTrainParams.curveResistance = vtype->getParameter().getCFParam(SUMO_ATTR_CURVE_RESISTANCE, myTrainParams.curveResistance);
120 390 : myTrainParams.roeckl_sharp_radius = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_SHARP_RADIUS, myTrainParams.roeckl_sharp_radius);
121 390 : myTrainParams.roeckl_numerator = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_NUMERATOR, myTrainParams.roeckl_numerator);
122 390 : myTrainParams.roeckl_numerator_sharp = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_NUMERATOR_SHARP, myTrainParams.roeckl_numerator_sharp);
123 390 : myTrainParams.roeckl_offset = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_OFFSET, myTrainParams.roeckl_offset);
124 390 : myTrainParams.roeckl_offset_sharp = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_OFFSET_SHARP, myTrainParams.roeckl_offset_sharp);
125 :
126 390 : if (myTrainParams.maxPower != INVALID_DOUBLE && myTrainParams.maxTraction == INVALID_DOUBLE) {
127 0 : throw ProcessError(TLF("Undefined maxPower for vType '%'.", vtype->getID()));
128 390 : } else if (myTrainParams.maxPower == INVALID_DOUBLE && myTrainParams.maxTraction != INVALID_DOUBLE) {
129 0 : throw ProcessError(TLF("Undefined maxTraction for vType '%'.", vtype->getID()));
130 : }
131 402 : if (myTrainParams.maxPower != INVALID_DOUBLE && vtype->getParameter().getCFParamString(SUMO_ATTR_TRACTION_TABLE, "") != "") {
132 0 : WRITE_WARNING(TLF("Ignoring tractionTable because maxPower and maxTraction are set for vType '%'.", vtype->getID()));
133 : }
134 390 : const bool hasSomeResCoef = (myTrainParams.resCoef_constant != INVALID_DOUBLE
135 386 : || myTrainParams.resCoef_linear != INVALID_DOUBLE
136 776 : || myTrainParams.resCoef_quadratic != INVALID_DOUBLE);
137 : const bool hasAllResCoef = (myTrainParams.resCoef_constant != INVALID_DOUBLE
138 4 : && myTrainParams.resCoef_linear != INVALID_DOUBLE
139 394 : && myTrainParams.resCoef_quadratic != INVALID_DOUBLE);
140 390 : if (hasSomeResCoef && !hasAllResCoef) {
141 0 : throw ProcessError(TLF("Some undefined resistance coefficients for vType '%' (requires resCoef_constant, resCoef_linear and resCoef_quadratic)", vtype->getID()));
142 : }
143 402 : if (myTrainParams.resCoef_constant != INVALID_DOUBLE && vtype->getParameter().getCFParamString(SUMO_ATTR_RESISTANCE_TABLE, "") != "") {
144 4 : WRITE_WARNING(TLF("Ignoring resistanceTable because resistance coefficients are set for vType '%'.", vtype->getID()));
145 : }
146 :
147 390 : if (myTrainParams.traction.empty() && myTrainParams.maxPower == INVALID_DOUBLE) {
148 12 : throw ProcessError(TLF("Either tractionTable or maxPower must be defined for vType '%' with Rail model type '%'.", vtype->getID(), trainType));
149 : }
150 386 : if (myTrainParams.resistance.empty() && myTrainParams.resCoef_constant == INVALID_DOUBLE) {
151 0 : throw ProcessError(TLF("Either resistanceTable or resCoef_constant must be defined for vType '%' with Rail model type '%'.", vtype->getID(), trainType));
152 : }
153 394 : }
154 :
155 :
156 770 : MSCFModel_Rail::~MSCFModel_Rail() { }
157 :
158 :
159 108170 : double MSCFModel_Rail::followSpeed(const MSVehicle* const veh, double speed, double gap,
160 : double /* predSpeed */, double /* predMaxDecel*/, const MSVehicle* const /*pred*/, const CalcReason /*usage*/) const {
161 :
162 : // followSpeed module is used for the simulation of moving block operations. The safety gap is chosen similar to the existing german
163 : // system CIR-ELKE (based on LZB). Other implementations of moving block systems may differ, but for now no appropriate parameter
164 : // can be set (would be per lane, not per train) -> hard-coded
165 :
166 : // @note: default train minGap of 5 is already subtracted from gap
167 108170 : if (speed >= 30 / 3.6) {
168 : // safety distance for higher speeds (>= 30 km/h)
169 71486 : gap = MAX2(0.0, gap + veh->getVehicleType().getMinGap() - 50);
170 : }
171 :
172 108170 : const double vsafe = maximumSafeStopSpeed(gap, myDecel, speed, false, TS, false); // absolute breaking distance
173 108170 : const double vmin = minNextSpeed(speed, veh);
174 108170 : const double vmax = maxNextSpeed(speed, veh);
175 :
176 108170 : if (MSGlobals::gSemiImplicitEulerUpdate) {
177 : return MIN2(vsafe, vmax);
178 : } else {
179 : // ballistic
180 : // XXX: the euler variant can break as strong as it wishes immediately! The ballistic cannot, refs. #2575.
181 : return MAX2(MIN2(vsafe, vmax), vmin);
182 : }
183 : }
184 :
185 :
186 : int
187 0 : MSCFModel_Rail::getModelID() const {
188 0 : return SUMO_TAG_CF_RAIL;
189 : }
190 :
191 :
192 : MSCFModel*
193 11 : MSCFModel_Rail::duplicate(const MSVehicleType* vtype) const {
194 11 : return new MSCFModel_Rail(vtype);
195 : }
196 :
197 : double
198 2617707 : MSCFModel_Rail::getRotWeight(const MSVehicle* const veh) const {
199 2617707 : return getWeight(veh) * myTrainParams.mf;
200 : }
201 :
202 : double
203 7853121 : MSCFModel_Rail::getWeight(const MSVehicle* const veh) const {
204 : // kg to tons
205 7853121 : return veh->getVehicleType().getMass() / 1000;
206 : }
207 :
208 : double
209 2617707 : MSCFModel_Rail::getCurveResistance(const MSVehicle* veh) const {
210 2617707 : if (myTrainParams.curveResistance > 0) {
211 30280 : const double r = veh->getCurveRadius();
212 30280 : if (r == std::numeric_limits<double>::max()) {
213 : return 0;
214 23124 : } else if (r >= myTrainParams.roeckl_sharp_radius) {
215 10404 : return 0.001 * myTrainParams.curveResistance * myTrainParams.roeckl_numerator / (r - myTrainParams.roeckl_offset);
216 12720 : } else if (r > myTrainParams.roeckl_offset_sharp) {
217 5404 : return 0.001 * myTrainParams.curveResistance * myTrainParams.roeckl_numerator_sharp / (r - myTrainParams.roeckl_offset_sharp);
218 : } else {
219 21948 : WRITE_WARNINGF("Cannot compute curve resistance for vehicle '%' with radius % at time %",
220 : veh->getID(), r, time2string(SIMSTEP));
221 7316 : return 0;
222 : }
223 : }
224 : return 0;
225 : }
226 :
227 :
228 2055715 : double MSCFModel_Rail::maxNextSpeed(double speed, const MSVehicle* const veh) const {
229 :
230 2055715 : if (speed >= myTrainParams.vmax) {
231 : return myTrainParams.vmax;
232 : }
233 :
234 : double targetSpeed = myTrainParams.vmax;
235 :
236 1687048 : double res = myTrainParams.getResistance(speed); // kN
237 :
238 1687048 : double slope = veh->getSlope();
239 1687048 : double gr = getWeight(veh) * GRAVITY * sin(DEG2RAD(slope)); //kN
240 1687048 : double cr = getWeight(veh) * getCurveResistance(veh); //kN
241 :
242 1687048 : double totalRes = res + gr + cr; //kN
243 :
244 1687048 : double trac = myTrainParams.getTraction(speed); // kN
245 : double a;
246 1687048 : if (speed < targetSpeed) {
247 1687048 : a = (trac - totalRes) / getRotWeight(veh); //kN/t == N/kg
248 : } else {
249 : a = 0.;
250 0 : if (totalRes > trac) {
251 0 : a = (trac - totalRes) / getRotWeight(veh); //kN/t == N/kg
252 : }
253 : }
254 1687048 : double maxNextSpeed = speed + ACCEL2SPEED(a);
255 :
256 : // std::cout << veh->getID() << " speed: " << (speed*3.6) << std::endl;
257 :
258 1687048 : return MIN2(myTrainParams.vmax, maxNextSpeed);
259 : }
260 :
261 :
262 930659 : double MSCFModel_Rail::minNextSpeed(double speed, const MSVehicle* const veh) const {
263 :
264 930659 : const double slope = veh->getSlope();
265 930659 : const double gr = getWeight(veh) * GRAVITY * sin(DEG2RAD(slope)); //kN
266 930659 : const double cr = getWeight(veh) * getCurveResistance(veh);
267 930659 : const double res = myTrainParams.getResistance(speed); // kN
268 930659 : const double totalRes = res + gr + cr; //kN
269 930659 : const double a = myTrainParams.decl + totalRes / getRotWeight(veh);
270 930659 : const double vMin = speed - ACCEL2SPEED(a);
271 930659 : if (MSGlobals::gSemiImplicitEulerUpdate) {
272 : return MAX2(vMin, 0.);
273 : } else {
274 : // NOTE: ballistic update allows for negative speeds to indicate a stop within the next timestep
275 : return vMin;
276 : }
277 :
278 : }
279 :
280 :
281 : double
282 272539 : MSCFModel_Rail::minNextSpeedEmergency(double speed, const MSVehicle* const veh) const {
283 272539 : return minNextSpeed(speed, veh);
284 : }
285 :
286 :
287 : //void
288 : //MSCFModel_Rail::initVehicleVariables(const MSVehicle *const veh, MSCFModel_Rail::VehicleVariables *pVariables) const {
289 : //
290 : // pVariables->setInitialized();
291 : //
292 : //}
293 :
294 :
295 0 : double MSCFModel_Rail::getSpeedAfterMaxDecel(double /* speed */) const {
296 :
297 : // //TODO: slope not known here
298 : // double gr = 0; //trainParams.weight * GRAVITY * edge.grade
299 : //
300 : // double a = 0;//trainParams.decl - gr/trainParams.rotWeight;
301 : //
302 : // return speed + a * DELTA_T / 1000.;
303 0 : WRITE_ERROR("function call not allowed for rail model. Exiting!");
304 0 : throw ProcessError();
305 : }
306 :
307 :
308 272539 : double MSCFModel_Rail::finalizeSpeed(MSVehicle* const veh, double vPos) const {
309 272539 : return MSCFModel::finalizeSpeed(veh, vPos);
310 : }
311 :
312 :
313 1011862 : double MSCFModel_Rail::freeSpeed(const MSVehicle* const /* veh */, double /* speed */, double dist, double targetSpeed,
314 : const bool onInsertion, const CalcReason /*usage*/) const {
315 :
316 : // MSCFModel_Rail::VehicleVariables *vars = (MSCFModel_Rail::VehicleVariables *) veh->getCarFollowVariables();
317 : // if (vars->isNotYetInitialized()) {
318 : // initVehicleVariables(veh, vars);
319 : // }
320 :
321 : //TODO: signals, coasting, ...
322 :
323 1011862 : if (MSGlobals::gSemiImplicitEulerUpdate) {
324 : // adapt speed to succeeding lane, no reaction time is involved
325 : // when breaking for y steps the following distance g is covered
326 : // (drive with v in the final step)
327 : // g = (y^2 + y) * 0.5 * b + y * v
328 : // y = ((((sqrt((b + 2.0*v)*(b + 2.0*v) + 8.0*b*g)) - b)*0.5 - v)/b)
329 1011861 : const double v = SPEED2DIST(targetSpeed);
330 1011861 : if (dist < v) {
331 : return targetSpeed;
332 : }
333 916074 : const double b = ACCEL2DIST(myDecel);
334 916074 : const double y = MAX2(0.0, ((sqrt((b + 2.0 * v) * (b + 2.0 * v) + 8.0 * b * dist) - b) * 0.5 - v) / b);
335 916074 : const double yFull = floor(y);
336 916074 : const double exactGap = (yFull * yFull + yFull) * 0.5 * b + yFull * v + (y > yFull ? v : 0.0);
337 916074 : const double fullSpeedGain = (yFull + (onInsertion ? 1. : 0.)) * ACCEL2SPEED(myTrainParams.decl);
338 1241006 : return DIST2SPEED(MAX2(0.0, dist - exactGap) / (yFull + 1)) + fullSpeedGain + targetSpeed;
339 : } else {
340 1 : WRITE_ERROR(TL("Anything else than semi implicit euler update is not yet implemented. Exiting!"));
341 1 : throw ProcessError();
342 : }
343 : }
344 :
345 :
346 1075314 : double MSCFModel_Rail::stopSpeed(const MSVehicle* const veh, const double speed, double gap, double decel, const CalcReason /*usage*/) const {
347 1075314 : return MIN2(maximumSafeStopSpeed(gap, decel, speed, false, TS, false), maxNextSpeed(speed, veh));
348 : }
|
