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 <utils/xml/SUMOSAXAttributes.h>
29 : #include <microsim/MSVehicle.h>
30 : #include <microsim/lcmodels/MSAbstractLaneChangeModel.h>
31 : #include "MSCFModel_Rail.h"
32 :
33 : // ===========================================================================
34 : // trainParams method definitions
35 : // ===========================================================================
36 :
37 : double
38 2658921 : MSCFModel_Rail::TrainParams::getResistance(double speed) const {
39 2658921 : if (resCoef_constant != INVALID_DOUBLE) {
40 20704 : return (resCoef_quadratic * speed * speed + resCoef_linear * speed + resCoef_constant); // kN
41 : } else {
42 2638217 : return LinearApproxHelpers::getInterpolatedValue(resistance, speed); // kN
43 : }
44 : }
45 :
46 :
47 : double
48 1714549 : MSCFModel_Rail::TrainParams::getTraction(double speed) const {
49 1714549 : if (maxPower != INVALID_DOUBLE) {
50 6596 : return MIN2(maxPower / speed, maxTraction); // kN
51 : } else {
52 1707953 : return LinearApproxHelpers::getInterpolatedValue(traction, speed); // kN
53 : }
54 : }
55 :
56 :
57 : // ===========================================================================
58 : // RailVehicleVariables method definitions
59 : // ===========================================================================
60 : void
61 10 : MSCFModel_Rail::RailVehicleVariables::saveState(OutputDevice& out, const MSCFModel& /*cfm*/) const {
62 10 : out.openTag(SUMO_TAG_CFM_VARIABLES);
63 10 : out.writeAttr(SUMO_ATTR_ID, "Rail");
64 10 : std::ostringstream internals;
65 10 : internals << odometerAngles.size() << " ";
66 28 : for (auto item : odometerAngles) {
67 18 : internals << item.first << " " << item.second << " ";
68 : }
69 10 : out.writeAttr(SUMO_ATTR_STATE, internals.str());
70 10 : out.closeTag();
71 10 : }
72 :
73 :
74 : void
75 10 : MSCFModel_Rail::RailVehicleVariables::loadState(const SUMOSAXAttributes& attrs) {
76 10 : bool ok = true;
77 10 : const std::string cfmID = attrs.get<std::string>(SUMO_ATTR_ID, nullptr, ok);
78 10 : if (cfmID != "Rail") {
79 0 : throw ProcessError(TLF("incompatible carFollowModel '%' when loading state for Rail", cfmID));
80 : }
81 10 : std::istringstream bis(attrs.getString(SUMO_ATTR_STATE));
82 : int odometerAnglesSize;
83 10 : bis >> odometerAnglesSize;
84 28 : for (int i = 0; i < odometerAnglesSize; i++) {
85 : double o;
86 : double a;
87 : bis >> o;
88 : bis >> a;
89 18 : odometerAngles.push_back(std::make_pair(o, a));
90 : }
91 20 : }
92 :
93 :
94 :
95 : double
96 71494 : MSCFModel_Rail::RailVehicleVariables::getIntegratedRadius(const MSVehicle* veh, double curveIntegration) {
97 71494 : const double odo = veh->getOdometer();
98 : // add new data point
99 71494 : if (odometerAngles.empty() || odometerAngles.back().first != odo) {
100 6034 : odometerAngles.push_back(std::make_pair(odo, veh->getAngle()));
101 : // clean up old data points beyond integration distance
102 11660 : while (odometerAngles.size() > 2) {
103 11584 : double distCleaned = odometerAngles.back().first - odometerAngles[1].first;
104 11584 : if (distCleaned >= curveIntegration) {
105 : odometerAngles.erase(odometerAngles.begin());
106 : } else {
107 : break;
108 : }
109 : }
110 : }
111 71494 : if (odometerAngles.size() > 1) {
112 71052 : const double dist = odometerAngles.back().first - odometerAngles.front().first;
113 71052 : const double angleDiff = GeomHelper::angleDiff(odometerAngles.back().second, odometerAngles.front().second);
114 71052 : if (dist < curveIntegration) {
115 17010 : return veh->getCurveRadius();
116 : }
117 : return angleDiff == 0
118 54042 : ? std::numeric_limits<double>::max()
119 45888 : : dist / fabs(angleDiff);
120 : } else {
121 442 : return veh->getCurveRadius();
122 : }
123 : }
124 :
125 :
126 :
127 : // ===========================================================================
128 : // method definitions
129 : // ===========================================================================
130 :
131 :
132 397 : MSCFModel_Rail::MSCFModel_Rail(const MSVehicleType* vtype) :
133 397 : MSCFModel(vtype) {
134 397 : const std::string trainType = vtype->getParameter().getCFParamString(SUMO_ATTR_TRAIN_TYPE, "NGT400");
135 397 : if (trainType.compare("RB425") == 0) {
136 58 : myTrainParams = initRB425Params();
137 339 : } else if (trainType.compare("RB628") == 0) {
138 19 : myTrainParams = initRB628Params();
139 320 : } else if (trainType.compare("NGT400") == 0) {
140 119 : myTrainParams = initNGT400Params();
141 201 : } else if (trainType.compare("NGT400_16") == 0) {
142 4 : myTrainParams = initNGT400_16Params();
143 197 : } else if (trainType.compare("ICE1") == 0) {
144 9 : myTrainParams = initICE1Params();
145 188 : } else if (trainType.compare("REDosto7") == 0) {
146 89 : myTrainParams = initREDosto7Params();
147 99 : } else if (trainType.compare("Freight") == 0) {
148 74 : myTrainParams = initFreightParams();
149 25 : } else if (trainType.compare("ICE3") == 0) {
150 5 : myTrainParams = initICE3Params();
151 20 : } else if (trainType.compare("MireoPlusB") == 0) {
152 4 : myTrainParams = initMireoPlusB2TParams();
153 16 : } else if (trainType.compare("MireoPlusH") == 0) {
154 4 : myTrainParams = initMireoPlusH2TParams();
155 12 : } else if (trainType.compare("custom") == 0) {
156 12 : myTrainParams = initCustomParams();
157 : } else {
158 0 : WRITE_ERRORF(TL("Unknown train type: %. Exiting!"), trainType);
159 0 : throw ProcessError();
160 : }
161 : // override with user values
162 397 : if (vtype->wasSet(VTYPEPARS_MAXSPEED_SET)) {
163 37 : myTrainParams.vmax = vtype->getMaxSpeed();
164 : }
165 397 : if (vtype->wasSet(VTYPEPARS_LENGTH_SET)) {
166 68 : myTrainParams.length = vtype->getLength();
167 : }
168 397 : myTrainParams.mf = vtype->getParameter().getCFParam(SUMO_ATTR_MASSFACTOR, myTrainParams.mf);
169 397 : myTrainParams.decl = vtype->getParameter().getCFParam(SUMO_ATTR_DECEL, myTrainParams.decl);
170 : setMaxDecel(myTrainParams.decl);
171 397 : setEmergencyDecel(vtype->getParameter().getCFParam(SUMO_ATTR_EMERGENCYDECEL, myTrainParams.decl + 0.3));
172 : // update type parameters so they are shown correctly in the gui (if defaults from trainType are used)
173 397 : const_cast<MSVehicleType*>(vtype)->setMaxSpeed(myTrainParams.vmax);
174 397 : const_cast<MSVehicleType*>(vtype)->setLength(myTrainParams.length);
175 397 : if (!vtype->wasSet(VTYPEPARS_MASS_SET)) {
176 : // tons to kg
177 365 : const_cast<MSVehicleType*>(vtype)->setMass(myTrainParams.weight * 1000);
178 : }
179 :
180 : // init tabular curves
181 397 : myTrainParams.traction = vtype->getParameter().getCFProfile(SUMO_ATTR_TRACTION_TABLE, myTrainParams.traction);
182 397 : myTrainParams.resistance = vtype->getParameter().getCFProfile(SUMO_ATTR_RESISTANCE_TABLE, myTrainParams.resistance);
183 :
184 : // init parametric curves
185 397 : myTrainParams.maxPower = vtype->getParameter().getCFParam(SUMO_ATTR_MAXPOWER, INVALID_DOUBLE);
186 397 : myTrainParams.maxTraction = vtype->getParameter().getCFParam(SUMO_ATTR_MAXTRACTION, INVALID_DOUBLE);
187 397 : myTrainParams.resCoef_constant = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_CONSTANT, INVALID_DOUBLE);
188 397 : myTrainParams.resCoef_linear = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_LINEAR, INVALID_DOUBLE);
189 397 : myTrainParams.resCoef_quadratic = vtype->getParameter().getCFParam(SUMO_ATTR_RESISTANCE_COEFFICIENT_QUADRATIC, INVALID_DOUBLE);
190 : // curve resistance parameters
191 397 : myTrainParams.curveResistance = vtype->getParameter().getCFParam(SUMO_ATTR_CURVE_RESISTANCE, myTrainParams.curveResistance);
192 397 : myTrainParams.curveIntegration = vtype->getParameter().getCFParam(SUMO_ATTR_CURVE_INTEGRATION, myTrainParams.curveIntegration);
193 397 : myTrainParams.roeckl_sharp_radius = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_SHARP_RADIUS, myTrainParams.roeckl_sharp_radius);
194 397 : myTrainParams.roeckl_numerator = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_NUMERATOR, myTrainParams.roeckl_numerator);
195 397 : myTrainParams.roeckl_numerator_sharp = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_NUMERATOR_SHARP, myTrainParams.roeckl_numerator_sharp);
196 397 : myTrainParams.roeckl_offset = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_OFFSET, myTrainParams.roeckl_offset);
197 397 : myTrainParams.roeckl_offset_sharp = vtype->getParameter().getCFParam(SUMO_ATTR_ROECKL_OFFSET_SHARP, myTrainParams.roeckl_offset_sharp);
198 :
199 397 : if (myTrainParams.maxPower != INVALID_DOUBLE && myTrainParams.maxTraction == INVALID_DOUBLE) {
200 0 : throw ProcessError(TLF("Undefined maxPower for vType '%'.", vtype->getID()));
201 397 : } else if (myTrainParams.maxPower == INVALID_DOUBLE && myTrainParams.maxTraction != INVALID_DOUBLE) {
202 0 : throw ProcessError(TLF("Undefined maxTraction for vType '%'.", vtype->getID()));
203 : }
204 409 : if (myTrainParams.maxPower != INVALID_DOUBLE && vtype->getParameter().getCFParamString(SUMO_ATTR_TRACTION_TABLE, "") != "") {
205 0 : WRITE_WARNING(TLF("Ignoring tractionTable because maxPower and maxTraction are set for vType '%'.", vtype->getID()));
206 : }
207 397 : const bool hasSomeResCoef = (myTrainParams.resCoef_constant != INVALID_DOUBLE
208 393 : || myTrainParams.resCoef_linear != INVALID_DOUBLE
209 790 : || myTrainParams.resCoef_quadratic != INVALID_DOUBLE);
210 : const bool hasAllResCoef = (myTrainParams.resCoef_constant != INVALID_DOUBLE
211 4 : && myTrainParams.resCoef_linear != INVALID_DOUBLE
212 401 : && myTrainParams.resCoef_quadratic != INVALID_DOUBLE);
213 397 : if (hasSomeResCoef && !hasAllResCoef) {
214 0 : throw ProcessError(TLF("Some undefined resistance coefficients for vType '%' (requires resCoef_constant, resCoef_linear and resCoef_quadratic)", vtype->getID()));
215 : }
216 409 : if (myTrainParams.resCoef_constant != INVALID_DOUBLE && vtype->getParameter().getCFParamString(SUMO_ATTR_RESISTANCE_TABLE, "") != "") {
217 4 : WRITE_WARNING(TLF("Ignoring resistanceTable because resistance coefficients are set for vType '%'.", vtype->getID()));
218 : }
219 :
220 397 : if (myTrainParams.traction.empty() && myTrainParams.maxPower == INVALID_DOUBLE) {
221 12 : throw ProcessError(TLF("Either tractionTable or maxPower must be defined for vType '%' with Rail model type '%'.", vtype->getID(), trainType));
222 : }
223 393 : if (myTrainParams.resistance.empty() && myTrainParams.resCoef_constant == INVALID_DOUBLE) {
224 0 : throw ProcessError(TLF("Either resistanceTable or resCoef_constant must be defined for vType '%' with Rail model type '%'.", vtype->getID(), trainType));
225 : }
226 401 : }
227 :
228 :
229 784 : MSCFModel_Rail::~MSCFModel_Rail() { }
230 :
231 :
232 111104 : double MSCFModel_Rail::followSpeed(const MSVehicle* const veh, double speed, double gap,
233 : double /* predSpeed */, double /* predMaxDecel*/, const MSVehicle* const /*pred*/, const CalcReason /*usage*/) const {
234 :
235 : // followSpeed module is used for the simulation of moving block operations. The safety gap is chosen similar to the existing german
236 : // system CIR-ELKE (based on LZB). Other implementations of moving block systems may differ, but for now no appropriate parameter
237 : // can be set (would be per lane, not per train) -> hard-coded
238 :
239 : // @note: default train minGap of 5 is already subtracted from gap
240 111104 : if (speed >= 30 / 3.6) {
241 : // safety distance for higher speeds (>= 30 km/h)
242 73824 : gap = MAX2(0.0, gap + veh->getVehicleType().getMinGap() - 50);
243 : }
244 :
245 111104 : const double vsafe = maximumSafeStopSpeed(gap, myDecel, speed, false, TS, false); // absolute breaking distance
246 111104 : const double vmin = minNextSpeed(speed, veh);
247 111104 : const double vmax = maxNextSpeed(speed, veh);
248 :
249 111104 : if (MSGlobals::gSemiImplicitEulerUpdate) {
250 : return MIN2(vsafe, vmax);
251 : } else {
252 : // ballistic
253 : // XXX: the euler variant can break as strong as it wishes immediately! The ballistic cannot, refs. #2575.
254 : return MAX2(MIN2(vsafe, vmax), vmin);
255 : }
256 : }
257 :
258 :
259 : int
260 0 : MSCFModel_Rail::getModelID() const {
261 0 : return SUMO_TAG_CF_RAIL;
262 : }
263 :
264 :
265 : MSCFModel*
266 11 : MSCFModel_Rail::duplicate(const MSVehicleType* vtype) const {
267 11 : return new MSCFModel_Rail(vtype);
268 : }
269 :
270 : double
271 2658921 : MSCFModel_Rail::getRotWeight(const MSVehicle* const veh) const {
272 2658921 : return getWeight(veh) * myTrainParams.mf;
273 : }
274 :
275 : double
276 7976763 : MSCFModel_Rail::getWeight(const MSVehicle* const veh) const {
277 : // kg to tons
278 7976763 : return veh->getVehicleType().getMass() / 1000;
279 : }
280 :
281 : double
282 2658921 : MSCFModel_Rail::getCurveResistance(const MSVehicle* veh) const {
283 2658921 : if (myTrainParams.curveResistance > 0) {
284 : RailVehicleVariables* vars = (RailVehicleVariables*)veh->getCarFollowVariables();
285 : assert(vars != nullptr);
286 71494 : const double r = vars->getIntegratedRadius(veh, myTrainParams.curveIntegration);
287 71494 : if (r == std::numeric_limits<double>::max()) {
288 : return 0;
289 56032 : } else if (r >= myTrainParams.roeckl_sharp_radius) {
290 40316 : return 0.001 * myTrainParams.curveResistance * myTrainParams.roeckl_numerator / (r - myTrainParams.roeckl_offset);
291 15716 : } else if (r > myTrainParams.roeckl_offset_sharp) {
292 13488 : return 0.001 * myTrainParams.curveResistance * myTrainParams.roeckl_numerator_sharp / (r - myTrainParams.roeckl_offset_sharp);
293 : } else {
294 6684 : WRITE_WARNINGF("Cannot compute curve resistance for vehicle '%' with radius % at time %",
295 : veh->getID(), r, time2string(SIMSTEP));
296 2228 : return 0;
297 : }
298 : }
299 : return 0;
300 : }
301 :
302 :
303 2088729 : double MSCFModel_Rail::maxNextSpeed(double speed, const MSVehicle* const veh) const {
304 :
305 2088729 : if (speed >= myTrainParams.vmax) {
306 : return myTrainParams.vmax;
307 : }
308 :
309 : double targetSpeed = myTrainParams.vmax;
310 :
311 1714549 : double res = myTrainParams.getResistance(speed); // kN
312 :
313 1714549 : double slope = veh->getSlope();
314 1714549 : double gr = getWeight(veh) * GRAVITY * sin(DEG2RAD(slope)); //kN
315 1714549 : double cr = getWeight(veh) * getCurveResistance(veh); //kN
316 :
317 1714549 : double totalRes = res + gr + cr; //kN
318 :
319 1714549 : double trac = myTrainParams.getTraction(speed); // kN
320 : double a;
321 1714549 : if (speed < targetSpeed) {
322 1714549 : a = (trac - totalRes) / getRotWeight(veh); //kN/t == N/kg
323 : } else {
324 : a = 0.;
325 0 : if (totalRes > trac) {
326 0 : a = (trac - totalRes) / getRotWeight(veh); //kN/t == N/kg
327 : }
328 : }
329 1714549 : double maxNextSpeed = speed + ACCEL2SPEED(a);
330 :
331 : // std::cout << veh->getID() << " speed: " << (speed*3.6) << std::endl;
332 :
333 1714549 : return MIN2(myTrainParams.vmax, maxNextSpeed);
334 : }
335 :
336 :
337 944372 : double MSCFModel_Rail::minNextSpeed(double speed, const MSVehicle* const veh) const {
338 :
339 944372 : const double slope = veh->getSlope();
340 944372 : const double gr = getWeight(veh) * GRAVITY * sin(DEG2RAD(slope)); //kN
341 944372 : const double cr = getWeight(veh) * getCurveResistance(veh);
342 944372 : const double res = myTrainParams.getResistance(speed); // kN
343 944372 : const double totalRes = res + gr + cr; //kN
344 944372 : const double a = myTrainParams.decl + totalRes / getRotWeight(veh);
345 944372 : const double vMin = speed - ACCEL2SPEED(a);
346 944372 : if (MSGlobals::gSemiImplicitEulerUpdate) {
347 : return MAX2(vMin, 0.);
348 : } else {
349 : // NOTE: ballistic update allows for negative speeds to indicate a stop within the next timestep
350 : return vMin;
351 : }
352 :
353 : }
354 :
355 :
356 : double
357 276132 : MSCFModel_Rail::minNextSpeedEmergency(double speed, const MSVehicle* const veh) const {
358 276132 : return minNextSpeed(speed, veh);
359 : }
360 :
361 :
362 : //void
363 : //MSCFModel_Rail::initVehicleVariables(const MSVehicle *const veh, MSCFModel_Rail::VehicleVariables *pVariables) const {
364 : //
365 : // pVariables->setInitialized();
366 : //
367 : //}
368 :
369 :
370 0 : double MSCFModel_Rail::getSpeedAfterMaxDecel(double /* speed */) const {
371 :
372 : // //TODO: slope not known here
373 : // double gr = 0; //trainParams.weight * GRAVITY * edge.grade
374 : //
375 : // double a = 0;//trainParams.decl - gr/trainParams.rotWeight;
376 : //
377 : // return speed + a * DELTA_T / 1000.;
378 0 : WRITE_ERROR("function call not allowed for rail model. Exiting!");
379 0 : throw ProcessError();
380 : }
381 :
382 :
383 276132 : double MSCFModel_Rail::finalizeSpeed(MSVehicle* const veh, double vPos) const {
384 276132 : return MSCFModel::finalizeSpeed(veh, vPos);
385 : }
386 :
387 :
388 1029710 : double MSCFModel_Rail::freeSpeed(const MSVehicle* const /* veh */, double /* speed */, double dist, double targetSpeed,
389 : const bool onInsertion, const CalcReason /*usage*/) const {
390 :
391 : // MSCFModel_Rail::VehicleVariables *vars = (MSCFModel_Rail::VehicleVariables *) veh->getCarFollowVariables();
392 : // if (vars->isNotYetInitialized()) {
393 : // initVehicleVariables(veh, vars);
394 : // }
395 :
396 : //TODO: signals, coasting, ...
397 :
398 1029710 : if (MSGlobals::gSemiImplicitEulerUpdate) {
399 : // adapt speed to succeeding lane, no reaction time is involved
400 : // when breaking for y steps the following distance g is covered
401 : // (drive with v in the final step)
402 : // g = (y^2 + y) * 0.5 * b + y * v
403 : // y = ((((sqrt((b + 2.0*v)*(b + 2.0*v) + 8.0*b*g)) - b)*0.5 - v)/b)
404 1029709 : const double v = SPEED2DIST(targetSpeed);
405 1029709 : if (dist < v) {
406 : return targetSpeed;
407 : }
408 933759 : const double b = ACCEL2DIST(myDecel);
409 933759 : const double y = MAX2(0.0, ((sqrt((b + 2.0 * v) * (b + 2.0 * v) + 8.0 * b * dist) - b) * 0.5 - v) / b);
410 933759 : const double yFull = floor(y);
411 933759 : const double exactGap = (yFull * yFull + yFull) * 0.5 * b + yFull * v + (y > yFull ? v : 0.0);
412 933759 : const double fullSpeedGain = (yFull + (onInsertion ? 1. : 0.)) * ACCEL2SPEED(myTrainParams.decl);
413 1269716 : return DIST2SPEED(MAX2(0.0, dist - exactGap) / (yFull + 1)) + fullSpeedGain + targetSpeed;
414 : } else {
415 1 : WRITE_ERROR(TL("Anything else than semi implicit euler update is not yet implemented. Exiting!"));
416 1 : throw ProcessError();
417 : }
418 : }
419 :
420 :
421 1094391 : double MSCFModel_Rail::stopSpeed(const MSVehicle* const veh, const double speed, double gap, double decel, const CalcReason /*usage*/) const {
422 1094391 : return MIN2(maximumSafeStopSpeed(gap, decel, speed, false, TS, false), maxNextSpeed(speed, veh));
423 : }
|
