From d493b4ab32810a518fa018e1f569a1f98bd8ddbc Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Tue, 9 Jun 2026 16:06:06 +0200 Subject: [PATCH 01/13] wip --- mseetc/etcs.py | 43 +++-- simulations/sim_launcher.py | 259 ++++++++++++++++++++++++++++-- tracks/swisstopo/analyzeTracks.py | 4 +- 3 files changed, 275 insertions(+), 31 deletions(-) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index 9d56885..820819d 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -165,10 +165,10 @@ def validateInput(self, target): if not 0 <= target.permittedVelocity < 400 / 3.6: raise ValueError("permittedVelocity must be between 0 and 400 km/h.") - if not 0 < target.EoA < self.track.length: + if not 0 < target.EoA <= self.track.length: raise ValueError("EoA must lie within the track length.") - if not 0 < target.SvL < self.track.length: + if not 0 < target.SvL <= self.track.length + target.overlap: raise ValueError("SvL must lie within the track length.") if not 0 <= target.targetVelocity < target.permittedVelocity: @@ -360,24 +360,19 @@ def computeSBICurve(self, SBI1_curve, SBI2_curve): return SBI_curve - def processCurvesBeforeTarget(self, curves, target): + def trimCurves(self, curves, target): permittedVelocity = target.permittedVelocity - start_position = target.position - self.distancePre speedLimits = computeCeilingSpeedLimits(permittedVelocity) curves["EBI"] = trimCurveToMaxVelocity(curves["EBI"], speedLimits["EBI [m/s]"]) - curves["EBI"] = addStartPointToCurve(curves["EBI"], speedLimits["EBI [m/s]"], start_position) curves["SBI"] = trimCurveToMaxVelocity(curves["SBI"], speedLimits["SBI [m/s]"]) - curves["SBI"] = addStartPointToCurve(curves["SBI"], speedLimits["SBI [m/s]"], start_position) curves["W"] = trimCurveToMaxVelocity(curves["W"], speedLimits["Warning [m/s]"]) - curves["W"] = addStartPointToCurve(curves["W"], speedLimits["Warning [m/s]"], start_position) curves["P"] = trimCurveToMaxVelocity(curves["P"], permittedVelocity) - curves["P"] = addStartPointToCurve(curves["P"], permittedVelocity, start_position) curves["I"] = trimCurveToMaxVelocity(curves["I"], permittedVelocity) @@ -391,6 +386,24 @@ def processCurvesBeforeTarget(self, curves, target): return curves + + def processCurvesBeforeTarget(self, curves, target): + + permittedVelocity = target.permittedVelocity + start_position = target.position - self.distancePre + + speedLimits = computeCeilingSpeedLimits(permittedVelocity) + + curves["EBI"] = addStartPointToCurve(curves["EBI"], speedLimits["EBI [m/s]"], start_position) + + curves["SBI"] = addStartPointToCurve(curves["SBI"], speedLimits["SBI [m/s]"], start_position) + + curves["W"] = addStartPointToCurve(curves["W"], speedLimits["Warning [m/s]"], start_position) + + curves["P"] = addStartPointToCurve(curves["P"], permittedVelocity, start_position) + + return curves + def processCurvcesAfterTarget(self, curves, target): targetVelocity = target.targetVelocity @@ -450,17 +463,19 @@ def computeTarget(self, target): curves["I"] = shiftCurveByTime(curves["P"], T_indication) - curves = self.processCurvesBeforeTarget(curves, target) - - if target.targetVelocity > 0: - - curves = self.processCurvcesAfterTarget(curves, target) + curves = self.trimCurves(curves, target) return curves def plotCurves(self, curves, target): + # curves = self.processCurvesBeforeTarget(curves, target) + # + # if target.targetVelocity > 0: + # + # curves = self.processCurvcesAfterTarget(curves, target) + targetPosition = target.EoA permittedVelocity = target.permittedVelocity targetVelocity = target.targetVelocity @@ -522,7 +537,7 @@ def plotCurves(self, curves, target): position=5000, overlap= 100, permittedVelocity=160/3.6, - targetVelocity=0/3.6 + targetVelocity=80/3.6 ) calculator = EtcsBrakingCurveCalculator(trainBrakingData, track, distancePre=5000, distancePost=1000) diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 515caaa..5ce1ad1 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -1,4 +1,8 @@ +import numpy as np +from matplotlib import pyplot as plt + from mseetc.efficiency import totalLossesFunction +from mseetc.etcs import EtcsBrakingCurveCalculator, BrakingTarget, trimCurveToMaxVelocity from mseetc.ocp import casadiSolver @@ -21,6 +25,185 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 raise ValueError("mode must be one of: 'perfect', 'static', 'dynamic'") +def getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions): + """ + Return stepwise track speed limit at given positions. + Assumes speedLimitPositions are sorted increasingly. + """ + + indices = np.searchsorted(speedLimitPositions, positions, side="right") - 1 + indices = np.clip(indices, 0, len(speedLimits) - 1) + + return speedLimits[indices] + + +def getBrakingTargetsFromSpeedLimits(track): + + speedLimitPositions = track.speedLimits.index.to_numpy(dtype=float) + speedLimits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + + # Add final stop target at end of track + speedLimitPositions = np.append(speedLimitPositions, track.length) + speedLimits = np.append(speedLimits, 0.0) + + v_max = max(speedLimits) + + targets = [] + + for idx in range(1, len(speedLimitPositions)): + + previousVelocity = speedLimits[idx - 1] + targetVelocity = speedLimits[idx] + + # Only speed decreases require a braking curve + if targetVelocity < previousVelocity: + + target = BrakingTarget( + position=speedLimitPositions[idx], + overlap=100, + permittedVelocity=v_max, + targetVelocity=targetVelocity, + ) + + targets.append(target) + + return targets, speedLimitPositions, speedLimits + + +def getEtcsSpeedLimits2(track, trainBrakingData, positionStep=10.0): + + calculator = EtcsBrakingCurveCalculator(trainBrakingData, track) + + targets, speedLimitPositions, speedLimits = getBrakingTargetsFromSpeedLimits(track) + + # Compute P curves for all speed decreases + pCurves = [] + + for target in targets: + + curveSet = calculator.computeTarget(target) + curveSet["P"].loc[target.position, "Velocity [m/s]"] = target.targetVelocity + pCurves.append(curveSet["P"]) + + # Build common position grid + positions = np.arange(0.0, track.length + positionStep, positionStep) + + # Start with the ordinary track speed limit + etcsVelocities = getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions) + + # Apply every ETCS P curve as an additional restriction + for pCurve in pCurves: + + curvePositions = pCurve.index.to_numpy(dtype=float) + curveVelocities = pCurve["Velocity [m/s]"].to_numpy(dtype=float) + + mask = ((curvePositions.min() <= positions) & (positions <= curvePositions.max())) + + interpolatedCurveVelocities = np.interp(positions[mask], curvePositions, curveVelocities) + + etcsVelocities[mask] = np.minimum(etcsVelocities[mask], interpolatedCurveVelocities) + + return positions, etcsVelocities + + + +def getEtcsSpeedLimits(track, trainBrakingData): + + etcsLimitsPositions = [] + etcsLimitsVelocities = [] + + speedLimitPositions = track.speedLimits.index.to_numpy(dtype=float) + speedLimits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + speedLimitPositions = np.append(speedLimitPositions, track.length) + speedLimits = np.append(speedLimits, 0) + + v_max = max(speedLimits) + + calculator = EtcsBrakingCurveCalculator(trainBrakingData, track, distancePre=5000, distancePost=1000) + + rev_idx = len(speedLimitPositions) - 1 + + while rev_idx > 0: + + if speedLimits[rev_idx-1] < speedLimits[rev_idx]: + + # speed limit increase -> no breaking curve needed + if etcsLimitsPositions[0] - 20 < speedLimitPositions[rev_idx] < etcsLimitsPositions[0] + 20: + + etcsLimitsPositions = np.concatenate([np.array([speedLimitPositions[rev_idx]]), etcsLimitsPositions]) + etcsLimitsVelocities = np.concatenate([np.array([speedLimits[rev_idx-1]]), etcsLimitsVelocities]) + + else: + + etcsLimitsPositions = np.concatenate([np.array([speedLimitPositions[rev_idx], speedLimitPositions[rev_idx]+1]), etcsLimitsPositions]) + etcsLimitsVelocities = np.concatenate([np.array([speedLimits[rev_idx-1], speedLimits[rev_idx]]), etcsLimitsVelocities]) + + rev_idx = rev_idx - 1 + continue + + # compute braking curve + target = BrakingTarget( + position=speedLimitPositions[rev_idx], + overlap=100, + permittedVelocity=v_max, + targetVelocity=speedLimits[rev_idx] + ) + + curve_set = calculator.computeTarget(target) + curve = curve_set["P"] + + curvePositions = curve.index.to_numpy(dtype=float) + curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) + + # find intersection with track speed limit profile + found = False + + while not found: + + sectionStart = speedLimitPositions[rev_idx-1] + + # case 1: braking curve intersects speed limit in current section + positionAtVLimit = np.interp(speedLimits[rev_idx-1], curveVelocities[::-1], curvePositions[::-1]) + + if sectionStart < positionAtVLimit: + + curve = trimCurveToMaxVelocity(curve, speedLimits[rev_idx-1]) + + curvePositions = curve.index.to_numpy(dtype=float) + curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) + + etcsLimitsPositions = np.concatenate([np.array([curvePositions[0]-1]), curvePositions, etcsLimitsPositions]) + etcsLimitsVelocities = np.concatenate([np.array([speedLimits[rev_idx-1]]), curveVelocities, etcsLimitsVelocities]) + + rev_idx = rev_idx - 1 + break + + # case 2: speed limit increases at sectionStart and braking curve intersects speed increase + if speedLimits[rev_idx-2] < speedLimits[rev_idx-1]: + + # speed limit increases at sectionStart + vAtSectionStart = np.interp(speedLimitPositions[rev_idx-1], curvePositions, curveVelocities) + + if speedLimits[rev_idx-2] < vAtSectionStart < speedLimits[rev_idx-1]: + + # braking curve intersects speed increase + curve = trimCurveToMaxVelocity(curve, vAtSectionStart) + + curvePositions = curve.index.to_numpy(dtype=float) + curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) + + etcsLimitsPositions = np.concatenate([curvePositions, etcsLimitsPositions]) + etcsLimitsVelocities = np.concatenate([curveVelocities, etcsLimitsVelocities]) + + rev_idx = rev_idx - 1 + break + + rev_idx = rev_idx - 1 + + etcsLimitsPositions = np.concatenate([np.array([0]), etcsLimitsPositions]) + etcsLimitsVelocities = np.concatenate([np.array([speedLimits[0]]), etcsLimitsVelocities]) + return etcsLimitsPositions, etcsLimitsVelocities + if __name__ == '__main__': @@ -29,34 +212,80 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 # Timetable startPosition = 0 # [m] - endPosition = 20000 # [m] - duration = 60*20 # [s] + endPosition = 10000 # [m] + duration = (endPosition-startPosition) / (90/3.6) # [s] train = Train(config={'id':'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') train.forceMinPn = 0 train.withPnBrake = False train.powerLosses = get_power_loss_function(train, "static") - track = Track(config={'id':'CH_ZH_LU'}, pathJSON='../tracks') - track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') + track = Track(config={'id':'00_var_speed_limit_quick_change'}, pathJSON='../tracks') + # track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') track.updateTrainLengthDependentValues(train) track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + trainBrakingData = { + "A_brake_emergency [m/s^2]": { + "velocity [m/s]": [0, 20, 40, 60], + "value [m/s^2]": [-0.9, -0.85, -0.8, -0.75], + }, + "A_brake_service [m/s^2]": { + "velocity [m/s]": [0, 20, 40, 60], + "value [m/s^2]": [-0.5, -0.45, -0.4, -0.35], + }, + "K_dry_rst [-]": 0.8, + "M_NVAVADH [-]": 0, + "K_wet_rst [-]": 0.9, + "T_traction [s]": 1, + "T_be [s]": 4, + "Kt_int [-]": 1.15, + "v_uncertainty [%]": 2.98, + "T_bs [s]": 3, + "T_bs1 [s]": 3, + "T_bs2 [s]": 3, + } + opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} - solver = casadiSolver(train, track, opts) + # solver = casadiSolver(train, track, opts) + # + # df, stats = solver.solve(duration) + # + # # print some info + # if df is not None: + # + # print("") + # print("Objective value = {:.2f} {}".format(stats['Cost'], 'kWh' if solver.opts.energyOptimal else 's')) + # print("") + # print("Maximum acceleration: {:5.2f}, with bound {}".format(df.max()['Acceleration [m/s^2]'], train.accMax if train.accMax is not None else 'None')) + # print("Maximum deceleration: {:5.2f}, with bound {}".format(df.min()['Acceleration [m/s^2]'], train.accMin if train.accMin is not None else 'None')) + # + # else: + # + # print("Solver failed!") - df, stats = solver.solve(duration) - # print some info - if df is not None: + ### Plot Trajectory - print("") - print("Objective value = {:.2f} {}".format(stats['Cost'], 'kWh' if solver.opts.energyOptimal else 's')) - print("") - print("Maximum acceleration: {:5.2f}, with bound {}".format(df.max()['Acceleration [m/s^2]'], train.accMax if train.accMax is not None else 'None')) - print("Maximum deceleration: {:5.2f}, with bound {}".format(df.min()['Acceleration [m/s^2]'], train.accMin if train.accMin is not None else 'None')) + fig, ax = plt.subplots(figsize=(16, 8)) - else: + x = track.speedLimits.index.to_numpy(dtype=float) + v = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + x_plot = np.append(x, track.length) + v_plot = np.append(v, v[-1]) + + etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits2(track, trainBrakingData) + + ax.step(x_plot/1000, v_plot*3.6, where="post", label="Track Speed Limit") + ax.plot(etcsLimitsPositions/1000, etcsLimitsVelocities*3.6, label="ETCS Speed Limit") + # ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, label="non-adjusted speed profile") + ax.set_title("Speed Profile Comparison") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Velocity [km/h]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + # ax.set_xlim(0, df["Position [m]"].max() / 1000) + ax.figure.tight_layout() - print("Solver failed!") \ No newline at end of file + plt.show() \ No newline at end of file diff --git a/tracks/swisstopo/analyzeTracks.py b/tracks/swisstopo/analyzeTracks.py index 250bea5..fd8c036 100644 --- a/tracks/swisstopo/analyzeTracks.py +++ b/tracks/swisstopo/analyzeTracks.py @@ -56,8 +56,8 @@ fig2, ax2 = plt.subplots(figsize=(16, 8)) - ax2.step(SBB_track.speedLimits.index.values / 1000, SBB_track.speedLimits["Speed limit [m/s]"].to_numpy()*3.6, label="SBB") - ax2.step((Topo_track.speedLimits.index.values-shift) / 1000, Topo_track.speedLimits["Speed limit [m/s]"].to_numpy()*3.6, label="Topo") + ax2.step(SBB_track.speedLimits.index.values / 1000, SBB_track.speedLimits["Speed limit [m/s]"].to_numpy()*3.6, where="post", label="SBB") + ax2.step((Topo_track.speedLimits.index.values-shift) / 1000, Topo_track.speedLimits["Speed limit [m/s]"].to_numpy()*3.6, where="post", label="Topo") ax2.set_title("Altitude Comparison") ax2.set_xlabel("Position [km]") ax2.set_ylabel("Velocity [km/h]") From 2a70bb691dc84d219457ab24f54a899a0676978b Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Tue, 9 Jun 2026 16:17:50 +0200 Subject: [PATCH 02/13] compute etcs speed profile for a whole track --- mseetc/etcs.py | 80 +++++++++++++ simulations/sim_launcher.py | 226 ++++-------------------------------- 2 files changed, 103 insertions(+), 203 deletions(-) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index 820819d..a4e20c6 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -8,6 +8,86 @@ from mseetc.track import Track +def getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions): + """ + Return stepwise track speed limit at given positions. + Assumes speedLimitPositions are sorted increasingly. + """ + + indices = np.searchsorted(speedLimitPositions, positions, side="right") - 1 + indices = np.clip(indices, 0, len(speedLimits) - 1) + + return speedLimits[indices] + + +def getBrakingTargetsFromSpeedLimits(track): + + speedLimitPositions = track.speedLimits.index.to_numpy(dtype=float) + speedLimits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + + # Add final stop target at end of track + speedLimitPositions = np.append(speedLimitPositions, track.length) + speedLimits = np.append(speedLimits, 0.0) + + v_max = max(speedLimits) + + targets = [] + + for idx in range(1, len(speedLimitPositions)): + + # Only speed decreases require a braking curve + if speedLimits[idx] < speedLimits[idx-1]: + + targets.append( + BrakingTarget( + position=speedLimitPositions[idx], + overlap=100, + permittedVelocity=v_max, + targetVelocity=speedLimits[idx], + ) + ) + + return targets, speedLimitPositions, speedLimits + + +def getEtcsSpeedLimits(trainBrakingData, track, positionStep=20.0): + + calculator = EtcsBrakingCurveCalculator(trainBrakingData, track) + + targets, speedLimitPositions, speedLimits = getBrakingTargetsFromSpeedLimits(track) + + # Compute P curves for all speed decreases + pCurves = [] + + for target in targets: + + curveSet = calculator.computeTarget(target) + curveSet["P"].loc[target.position, "Velocity [m/s]"] = target.targetVelocity + pCurves.append(curveSet["P"]) + + # Build common position grid + positions = np.arange(0.0, track.length + positionStep, positionStep) + positions = np.union1d(positions, speedLimitPositions) + positions = np.sort(positions) + + # Start with the ordinary track speed limit + etcsVelocities = getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions) + + # Apply every ETCS P curve as an additional restriction + for pCurve in pCurves: + + curvePositions = pCurve.index.to_numpy(dtype=float) + curveVelocities = pCurve["Velocity [m/s]"].to_numpy(dtype=float) + + mask = ((curvePositions.min() <= positions) & (positions <= curvePositions.max())) + + interpolatedCurveVelocities = np.interp(positions[mask], curvePositions, curveVelocities) + + etcsVelocities[mask] = np.minimum(etcsVelocities[mask], interpolatedCurveVelocities) + + return positions, etcsVelocities + + def shiftCurveByTime(dfCurve, timeShift): positionsOriginal = dfCurve.index.to_numpy() diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 5ce1ad1..19702fb 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -2,7 +2,7 @@ from matplotlib import pyplot as plt from mseetc.efficiency import totalLossesFunction -from mseetc.etcs import EtcsBrakingCurveCalculator, BrakingTarget, trimCurveToMaxVelocity +from mseetc.etcs import getEtcsSpeedLimits from mseetc.ocp import casadiSolver @@ -25,186 +25,6 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 raise ValueError("mode must be one of: 'perfect', 'static', 'dynamic'") -def getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions): - """ - Return stepwise track speed limit at given positions. - Assumes speedLimitPositions are sorted increasingly. - """ - - indices = np.searchsorted(speedLimitPositions, positions, side="right") - 1 - indices = np.clip(indices, 0, len(speedLimits) - 1) - - return speedLimits[indices] - - -def getBrakingTargetsFromSpeedLimits(track): - - speedLimitPositions = track.speedLimits.index.to_numpy(dtype=float) - speedLimits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) - - # Add final stop target at end of track - speedLimitPositions = np.append(speedLimitPositions, track.length) - speedLimits = np.append(speedLimits, 0.0) - - v_max = max(speedLimits) - - targets = [] - - for idx in range(1, len(speedLimitPositions)): - - previousVelocity = speedLimits[idx - 1] - targetVelocity = speedLimits[idx] - - # Only speed decreases require a braking curve - if targetVelocity < previousVelocity: - - target = BrakingTarget( - position=speedLimitPositions[idx], - overlap=100, - permittedVelocity=v_max, - targetVelocity=targetVelocity, - ) - - targets.append(target) - - return targets, speedLimitPositions, speedLimits - - -def getEtcsSpeedLimits2(track, trainBrakingData, positionStep=10.0): - - calculator = EtcsBrakingCurveCalculator(trainBrakingData, track) - - targets, speedLimitPositions, speedLimits = getBrakingTargetsFromSpeedLimits(track) - - # Compute P curves for all speed decreases - pCurves = [] - - for target in targets: - - curveSet = calculator.computeTarget(target) - curveSet["P"].loc[target.position, "Velocity [m/s]"] = target.targetVelocity - pCurves.append(curveSet["P"]) - - # Build common position grid - positions = np.arange(0.0, track.length + positionStep, positionStep) - - # Start with the ordinary track speed limit - etcsVelocities = getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions) - - # Apply every ETCS P curve as an additional restriction - for pCurve in pCurves: - - curvePositions = pCurve.index.to_numpy(dtype=float) - curveVelocities = pCurve["Velocity [m/s]"].to_numpy(dtype=float) - - mask = ((curvePositions.min() <= positions) & (positions <= curvePositions.max())) - - interpolatedCurveVelocities = np.interp(positions[mask], curvePositions, curveVelocities) - - etcsVelocities[mask] = np.minimum(etcsVelocities[mask], interpolatedCurveVelocities) - - return positions, etcsVelocities - - - -def getEtcsSpeedLimits(track, trainBrakingData): - - etcsLimitsPositions = [] - etcsLimitsVelocities = [] - - speedLimitPositions = track.speedLimits.index.to_numpy(dtype=float) - speedLimits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) - speedLimitPositions = np.append(speedLimitPositions, track.length) - speedLimits = np.append(speedLimits, 0) - - v_max = max(speedLimits) - - calculator = EtcsBrakingCurveCalculator(trainBrakingData, track, distancePre=5000, distancePost=1000) - - rev_idx = len(speedLimitPositions) - 1 - - while rev_idx > 0: - - if speedLimits[rev_idx-1] < speedLimits[rev_idx]: - - # speed limit increase -> no breaking curve needed - if etcsLimitsPositions[0] - 20 < speedLimitPositions[rev_idx] < etcsLimitsPositions[0] + 20: - - etcsLimitsPositions = np.concatenate([np.array([speedLimitPositions[rev_idx]]), etcsLimitsPositions]) - etcsLimitsVelocities = np.concatenate([np.array([speedLimits[rev_idx-1]]), etcsLimitsVelocities]) - - else: - - etcsLimitsPositions = np.concatenate([np.array([speedLimitPositions[rev_idx], speedLimitPositions[rev_idx]+1]), etcsLimitsPositions]) - etcsLimitsVelocities = np.concatenate([np.array([speedLimits[rev_idx-1], speedLimits[rev_idx]]), etcsLimitsVelocities]) - - rev_idx = rev_idx - 1 - continue - - # compute braking curve - target = BrakingTarget( - position=speedLimitPositions[rev_idx], - overlap=100, - permittedVelocity=v_max, - targetVelocity=speedLimits[rev_idx] - ) - - curve_set = calculator.computeTarget(target) - curve = curve_set["P"] - - curvePositions = curve.index.to_numpy(dtype=float) - curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) - - # find intersection with track speed limit profile - found = False - - while not found: - - sectionStart = speedLimitPositions[rev_idx-1] - - # case 1: braking curve intersects speed limit in current section - positionAtVLimit = np.interp(speedLimits[rev_idx-1], curveVelocities[::-1], curvePositions[::-1]) - - if sectionStart < positionAtVLimit: - - curve = trimCurveToMaxVelocity(curve, speedLimits[rev_idx-1]) - - curvePositions = curve.index.to_numpy(dtype=float) - curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) - - etcsLimitsPositions = np.concatenate([np.array([curvePositions[0]-1]), curvePositions, etcsLimitsPositions]) - etcsLimitsVelocities = np.concatenate([np.array([speedLimits[rev_idx-1]]), curveVelocities, etcsLimitsVelocities]) - - rev_idx = rev_idx - 1 - break - - # case 2: speed limit increases at sectionStart and braking curve intersects speed increase - if speedLimits[rev_idx-2] < speedLimits[rev_idx-1]: - - # speed limit increases at sectionStart - vAtSectionStart = np.interp(speedLimitPositions[rev_idx-1], curvePositions, curveVelocities) - - if speedLimits[rev_idx-2] < vAtSectionStart < speedLimits[rev_idx-1]: - - # braking curve intersects speed increase - curve = trimCurveToMaxVelocity(curve, vAtSectionStart) - - curvePositions = curve.index.to_numpy(dtype=float) - curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) - - etcsLimitsPositions = np.concatenate([curvePositions, etcsLimitsPositions]) - etcsLimitsVelocities = np.concatenate([curveVelocities, etcsLimitsVelocities]) - - rev_idx = rev_idx - 1 - break - - rev_idx = rev_idx - 1 - - etcsLimitsPositions = np.concatenate([np.array([0]), etcsLimitsPositions]) - etcsLimitsVelocities = np.concatenate([np.array([speedLimits[0]]), etcsLimitsVelocities]) - return etcsLimitsPositions, etcsLimitsVelocities - - if __name__ == '__main__': from mseetc.train import Train @@ -212,8 +32,8 @@ def getEtcsSpeedLimits(track, trainBrakingData): # Timetable startPosition = 0 # [m] - endPosition = 10000 # [m] - duration = (endPosition-startPosition) / (90/3.6) # [s] + endPosition = 20000 # [m] + duration = (endPosition-startPosition) / (105/3.6) # [s] train = Train(config={'id':'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') train.forceMinPn = 0 @@ -221,7 +41,7 @@ def getEtcsSpeedLimits(track, trainBrakingData): train.powerLosses = get_power_loss_function(train, "static") track = Track(config={'id':'00_var_speed_limit_quick_change'}, pathJSON='../tracks') - # track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') + track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') track.updateTrainLengthDependentValues(train) track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') @@ -248,22 +68,22 @@ def getEtcsSpeedLimits(track, trainBrakingData): opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} - # solver = casadiSolver(train, track, opts) - # - # df, stats = solver.solve(duration) - # - # # print some info - # if df is not None: - # - # print("") - # print("Objective value = {:.2f} {}".format(stats['Cost'], 'kWh' if solver.opts.energyOptimal else 's')) - # print("") - # print("Maximum acceleration: {:5.2f}, with bound {}".format(df.max()['Acceleration [m/s^2]'], train.accMax if train.accMax is not None else 'None')) - # print("Maximum deceleration: {:5.2f}, with bound {}".format(df.min()['Acceleration [m/s^2]'], train.accMin if train.accMin is not None else 'None')) - # - # else: - # - # print("Solver failed!") + solver = casadiSolver(train, track, opts) + + df, stats = solver.solve(duration) + + # print some info + if df is not None: + + print("") + print("Objective value = {:.2f} {}".format(stats['Cost'], 'kWh' if solver.opts.energyOptimal else 's')) + print("") + print("Maximum acceleration: {:5.2f}, with bound {}".format(df.max()['Acceleration [m/s^2]'], train.accMax if train.accMax is not None else 'None')) + print("Maximum deceleration: {:5.2f}, with bound {}".format(df.min()['Acceleration [m/s^2]'], train.accMin if train.accMin is not None else 'None')) + + else: + + print("Solver failed!") ### Plot Trajectory @@ -275,17 +95,17 @@ def getEtcsSpeedLimits(track, trainBrakingData): x_plot = np.append(x, track.length) v_plot = np.append(v, v[-1]) - etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits2(track, trainBrakingData) + etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits(trainBrakingData, track) ax.step(x_plot/1000, v_plot*3.6, where="post", label="Track Speed Limit") ax.plot(etcsLimitsPositions/1000, etcsLimitsVelocities*3.6, label="ETCS Speed Limit") - # ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, label="non-adjusted speed profile") + ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, label="non-adjusted speed profile") ax.set_title("Speed Profile Comparison") ax.set_xlabel("Position [km]") ax.set_ylabel("Velocity [km/h]") ax.grid(True, which="both", linestyle="--", alpha=0.5) ax.legend(loc="upper right") - # ax.set_xlim(0, df["Position [m]"].max() / 1000) + ax.set_xlim(0, df["Position [m]"].max() / 1000) ax.figure.tight_layout() plt.show() \ No newline at end of file From 5d90f150637585ce1d88092e948d532f54a63b4b Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Tue, 9 Jun 2026 17:51:14 +0200 Subject: [PATCH 03/13] etcs works in mseetc --- mseetc/etcs.py | 8 ++-- mseetc/ocp.py | 11 +++-- mseetc/track.py | 21 ++++++++-- mseetc/utils.py | 18 +++++++- simulations/sim_launcher.py | 84 +++++++++++++++++++++++-------------- 5 files changed, 101 insertions(+), 41 deletions(-) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index a4e20c6..f420bf5 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -5,7 +5,7 @@ import pandas as pd from matplotlib import pyplot as plt -from mseetc.track import Track +# from mseetc.track import Track def getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions): @@ -590,6 +590,8 @@ def plotCurves(self, curves, target): if __name__ == '__main__': + from mseetc.track import Track + trainBrakingData = { "A_brake_emergency [m/s^2]": { "velocity [m/s]": [0, 20, 40, 60], @@ -597,7 +599,7 @@ def plotCurves(self, curves, target): }, "A_brake_service [m/s^2]": { "velocity [m/s]": [0, 20, 40, 60], - "value [m/s^2]": [-0.5, -0.45, -0.4, -0.35], + "value [m/s^2]": [-0.55, -0.5, -0.45, -0.4], }, "K_dry_rst [-]": 0.8, "M_NVAVADH [-]": 0, @@ -617,7 +619,7 @@ def plotCurves(self, curves, target): position=5000, overlap= 100, permittedVelocity=160/3.6, - targetVelocity=80/3.6 + targetVelocity=0/3.6 ) calculator = EtcsBrakingCurveCalculator(trainBrakingData, track, distancePre=5000, distancePost=1000) diff --git a/mseetc/ocp.py b/mseetc/ocp.py index 67ed55b..3a6a427 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -1,6 +1,4 @@ -import numpy as np import pandas as pd -import casadi as ca from mseetc.train import * @@ -27,10 +25,12 @@ def __init__(self, paramsDict): self.integrateLosses = False # integrate losses or take mid-point rule - self.chooseClosestTunnelCrossSection = True # if exact tunnel cross section is not specified in train tunnel resistances, choose the closest value + self.chooseClosestTunnelCrossSection = True # if exact tunnel cross section is not specified in train tunnel resistances, choose the closest value self.pwcLengthDependentTrackAttributes = False + self.withEtcsBrakingCurves = False # use ETCS braking curves + super().__init__(paramsDict) @@ -128,6 +128,11 @@ def __init__(self, train, track, optsDict={}): self.points = computeDiscretizationPoints(track, numIntervals, opts) self.steps = np.diff(self.points.index) + if opts.withEtcsBrakingCurves: + + velocities = np.interp(self.points.index.to_numpy(), track.etcsPositions, track.etcsVelocities) + self.points["Speed limit [m/s]"] = velocities + # real-time parameters self.initialTime = ca.MX.sym('t0') diff --git a/mseetc/track.py b/mseetc/track.py index 086c2d0..8617e36 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -1,13 +1,13 @@ import json -import math import sys from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd +from mseetc.etcs import getEtcsSpeedLimits from mseetc.utils import checkTTOBenchVersion, convertUnit, pickEquallySpacedPoints, plotSpeedLimits, plotGradients, \ - plotCurvatures + plotCurvatures, getRelevantEtcsBrakingPositions plotDebug = False @@ -98,7 +98,14 @@ def computeDiscretizationPoints(track, numIntervals, opts): df1 = track.mergeDataFrames() - pos = pickEquallySpacedPoints(0, track.length, numIntervals, df1.index.to_numpy(dtype=float)) + requiredPoints = df1.index.to_numpy(dtype=float) + + if opts.withEtcsBrakingCurves: + + brakingPoints = getRelevantEtcsBrakingPositions(track) + requiredPoints = np.append(requiredPoints, brakingPoints) + + pos = pickEquallySpacedPoints(0, track.length, numIntervals, requiredPoints) df2 = pd.DataFrame({'position [m]':pos}).set_index('position [m]') df3 = df2.join(df1, how='outer').ffill() @@ -598,6 +605,14 @@ def crop(dfIn): self.crossSections = crop(self.crossSections) + def setEtcsSpeedLimits(self, trainBrakingData): + + etcsPositions, etcsVelocities = getEtcsSpeedLimits(trainBrakingData, self) + + self.etcsPositions = etcsPositions + self.etcsVelocities = etcsVelocities + + def updateTrainLengthDependentValues(self, train): self.updateSpeedLimitsToTrainLength(train.length) diff --git a/mseetc/utils.py b/mseetc/utils.py index d23e0d0..ae7b56e 100644 --- a/mseetc/utils.py +++ b/mseetc/utils.py @@ -541,7 +541,7 @@ def pickEquallySpacedPoints(startPoint, endPoint, numIntervals, requiredPoints): if len(cand_without_required) >= num_of_remaining_points: picked_points = np.random.choice(cand_without_required, size=num_of_remaining_points, replace=False) - return picked_points + return np.sort(np.r_[requiredPoints, picked_points]) m *= 2 @@ -628,6 +628,22 @@ def plotCurvatures(track, pos_adj, c_adj, c_linear): plt.show() +def getRelevantEtcsBrakingPositions(track): + + relevantPositions = [] + + positions = track.etcsPositions + velocities = track.etcsVelocities + + for idx in range(1, len(positions)): + + if velocities[idx - 2] > velocities[idx - 1] == velocities[idx]: + + relevantPositions.append(positions[idx]) + + return relevantPositions + + if __name__ == '__main__': pass diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 19702fb..1537d1d 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -25,34 +25,31 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 raise ValueError("mode must be one of: 'perfect', 'static', 'dynamic'") -if __name__ == '__main__': +def printStats(df, stats, solver): - from mseetc.train import Train - from mseetc.track import Track + if df is not None: - # Timetable - startPosition = 0 # [m] - endPosition = 20000 # [m] - duration = (endPosition-startPosition) / (105/3.6) # [s] + print("") + print("Objective value = {:.2f} {}".format(stats['Cost'], 'kWh' if solver.opts.energyOptimal else 's')) + print("") + print("Maximum acceleration: {:5.2f}, with bound {}".format(df.max()['Acceleration [m/s^2]'], train.accMax if train.accMax is not None else 'None')) + print("Maximum deceleration: {:5.2f}, with bound {}".format(df.min()['Acceleration [m/s^2]'],train.accMin if train.accMin is not None else 'None')) - train = Train(config={'id':'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') - train.forceMinPn = 0 - train.withPnBrake = False - train.powerLosses = get_power_loss_function(train, "static") + else: - track = Track(config={'id':'00_var_speed_limit_quick_change'}, pathJSON='../tracks') - track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') - track.updateTrainLengthDependentValues(train) - track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + print("Solver failed!") - trainBrakingData = { + +def getTrainBrakingData(): + + return { "A_brake_emergency [m/s^2]": { "velocity [m/s]": [0, 20, 40, 60], "value [m/s^2]": [-0.9, -0.85, -0.8, -0.75], }, "A_brake_service [m/s^2]": { "velocity [m/s]": [0, 20, 40, 60], - "value [m/s^2]": [-0.5, -0.45, -0.4, -0.35], + "value [m/s^2]": [-0.55, -0.5, -0.45, -0.4], }, "K_dry_rst [-]": 0.8, "M_NVAVADH [-]": 0, @@ -66,24 +63,42 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 "T_bs2 [s]": 3, } + +if __name__ == '__main__': + + from mseetc.train import Train + from mseetc.track import Track + + # Timetable + startPosition = 0 # [m] + endPosition = 20000 # [m] + duration = (endPosition-startPosition) / (105/3.6) # [s] + + train = Train(config={'id':'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') + train.forceMinPn = 0 + train.withPnBrake = False + train.powerLosses = get_power_loss_function(train, "static") + + track = Track(config={'id':'00_var_speed_limit_quick_change'}, pathJSON='../tracks') + track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') + track.updateTrainLengthDependentValues(train) + track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) - # print some info - if df is not None: + printStats(df, stats, solver) - print("") - print("Objective value = {:.2f} {}".format(stats['Cost'], 'kWh' if solver.opts.energyOptimal else 's')) - print("") - print("Maximum acceleration: {:5.2f}, with bound {}".format(df.max()['Acceleration [m/s^2]'], train.accMax if train.accMax is not None else 'None')) - print("Maximum deceleration: {:5.2f}, with bound {}".format(df.min()['Acceleration [m/s^2]'], train.accMin if train.accMin is not None else 'None')) + trainBrakingData = getTrainBrakingData() + track.setEtcsSpeedLimits(trainBrakingData) + opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} - else: + solverEtcs = casadiSolver(train, track, opts) + dfEtcs, statsEtcs = solverEtcs.solve(duration) - print("Solver failed!") + printStats(dfEtcs, statsEtcs, solverEtcs) ### Plot Trajectory @@ -97,9 +112,12 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits(trainBrakingData, track) - ax.step(x_plot/1000, v_plot*3.6, where="post", label="Track Speed Limit") - ax.plot(etcsLimitsPositions/1000, etcsLimitsVelocities*3.6, label="ETCS Speed Limit") - ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, label="non-adjusted speed profile") + ax.step(x_plot/1000, v_plot*3.6, where="post", color="black", linestyle="-", label="Track Speed Limit") + ax.plot(etcsLimitsPositions/1000, etcsLimitsVelocities*3.6, color="red", linestyle="-", label="ETCS Speed Limit") + + ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, linestyle="--", label="non-adjusted speed profile") + ax.plot(dfEtcs["Position [m]"] / 1000, dfEtcs["Velocity [m/s]"] * 3.6, linestyle="--", label="ETCS-adjusted speed profile") + ax.set_title("Speed Profile Comparison") ax.set_xlabel("Position [km]") ax.set_ylabel("Velocity [km/h]") @@ -108,4 +126,8 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 ax.set_xlim(0, df["Position [m]"].max() / 1000) ax.figure.tight_layout() - plt.show() \ No newline at end of file + plt.show() + + costRatio = (statsEtcs["Cost"] - stats["Cost"]) / stats["Cost"] + + print(f"Cost increase with ETCS: {costRatio:.2%}") \ No newline at end of file From 959895a48f700c1da29f77e06c17f791e9448f28 Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Wed, 10 Jun 2026 11:43:41 +0200 Subject: [PATCH 04/13] wip --- mseetc/etcs.py | 51 +++++++---- mseetc/train.py | 152 +++++++++++++++++++++++++++++-- mseetc/utils.py | 2 +- trains/CH_Stadler_FLIRT_TPF.json | 50 ++++++++++ 4 files changed, 226 insertions(+), 29 deletions(-) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index f420bf5..7f1123b 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -283,18 +283,27 @@ def computeABrakeSafe(self): def computeAGradient(self, currentPosition): - positions = self.track.gradients.index.values - gradients = self.track.gradients["Gradient [permil]"].values + positions = self.track.gradients.index.to_numpy(dtype=float) + gradients = self.track.gradients["Gradient [permil]"].to_numpy(dtype=float) + + if "Gradient linear term [permil/m]" in self.track.gradients.columns: + + gradientsLinearTerm = self.track.gradients["Gradient linear term [permil/m]"].to_numpy(dtype=float) + + else: + + gradientsLinearTerm = np.zeros(len(self.track.gradients)) + idx = bisect_right(positions, currentPosition) - 1 idx = max(0, min(idx, len(positions) - 1)) - if idx == 0: + # If the backward-computed curve extends before the first known gradient point, assume flat track. + if currentPosition < positions[0]: - # If the backward-computed curve extends before the first known gradient point, assume flat track. return 0 - gradient = gradients[idx] + gradient = gradients[idx] + gradientsLinearTerm[idx] * (currentPosition - positions[idx]) return 9.81 * gradient * 0.001 @@ -360,7 +369,7 @@ def computeEBICurve(self, EBD_curve, targetVelocity): T_traction = trainBrakingData["T_traction [s]"] T_be = trainBrakingData["T_be [s]"] Kt_int = trainBrakingData["Kt_int [-]"] - v_uncertainty = trainBrakingData["v_uncertainty [%]"] * 0.01 + v_uncertainty = trainBrakingData["v_uncertainty [%]"] positionsEBD = EBD_curve.index.to_numpy() velocitiesEBD = EBD_curve["Velocity [m/s]"].to_numpy() @@ -529,11 +538,11 @@ def computeTarget(self, target): curves["EBI"] = self.computeEBICurve(curves["EBD"], target.targetVelocity) - curves["SBI2"] = shiftCurveByTime(curves["EBI"], trainBrakingData["T_bs2 [s]"]) + curves["SBI2"] = shiftCurveByTime(curves["EBI"], trainBrakingData["T_bs [s]"]) curves["SBD"] = self.computeBrakingCurve(self.trainBrakingData["A_brake_service [m/s^2]"], target.EoA, target.permittedVelocity, target.targetVelocity) - curves["SBI1"] = shiftCurveByTime(curves["SBD"], trainBrakingData["T_bs1 [s]"]) + curves["SBI1"] = shiftCurveByTime(curves["SBD"], trainBrakingData["T_bs [s]"]) curves["SBI"] = self.computeSBICurve(curves["SBI1"], curves["SBI2"]) @@ -591,29 +600,31 @@ def plotCurves(self, curves, target): if __name__ == '__main__': from mseetc.track import Track + from mseetc.train import Train trainBrakingData = { - "A_brake_emergency [m/s^2]": { + "A_brake_emergency [m/s^2]": { # [train] "velocity [m/s]": [0, 20, 40, 60], "value [m/s^2]": [-0.9, -0.85, -0.8, -0.75], }, - "A_brake_service [m/s^2]": { + "A_brake_service [m/s^2]": { # [train] "velocity [m/s]": [0, 20, 40, 60], "value [m/s^2]": [-0.55, -0.5, -0.45, -0.4], }, - "K_dry_rst [-]": 0.8, - "M_NVAVADH [-]": 0, - "K_wet_rst [-]": 0.9, - "T_traction [s]": 1, - "T_be [s]": 4, - "Kt_int [-]": 1.15, - "v_uncertainty [%]": 2.98, - "T_bs [s]": 3, - "T_bs1 [s]": 3, - "T_bs2 [s]": 3, + "K_dry_rst [-]": 0.8, # [train] + "M_NVAVADH [-]": 0, # [infra] + "K_wet_rst [-]": 0.9, # [train] + "T_traction [s]": 1, # [train] + "T_be [s]": 4, # [train] + "Kt_int [-]": 1.15, # [infra] + "v_uncertainty [%]": 2.98, # [train] + "T_bs [s]": 3 # [train] } + train = Train(config={'id': 'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') + track.updateTrainLengthDependentValues(train) target = BrakingTarget( position=5000, diff --git a/mseetc/train.py b/mseetc/train.py index 9ef6d4c..52e11fc 100644 --- a/mseetc/train.py +++ b/mseetc/train.py @@ -105,7 +105,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> raise ValueError("Both efficiencies need to be specified in json file!") - if 'values' in 'efficiency traction' or 'values' in 'efficiency reg brake': + if 'values' in data['efficiency traction'] or 'values' in data['efficiency reg brake']: raise ValueError("Dynamic efficiency from json file not implemented yet!") @@ -114,23 +114,111 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> if "tunnel resistance" in data: - tunnel_data = data["tunnel resistance"] # additive aerodynamic tunnel drag as dict per tunnel cross section + tunnelData = data["tunnel resistance"] # additive aerodynamic tunnel drag as dict per tunnel cross section - cross_section_unit = tunnel_data["units"]["cross section"] # tunnel cross section [m^2] - resistance_unit = tunnel_data["units"]["coefficient"] # resistance coefficient [N/(m/s)^2] + crossSectionUnit = tunnelData["units"]["cross section"] # tunnel cross section [m^2] + resistanceUnit = tunnelData["units"]["coefficient"] # resistance coefficient [N/(m/s)^2] self.tunnelCoefficients = { - convertUnit(cross_section, cross_section_unit): convertUnit( + convertUnit(cross_section, crossSectionUnit): convertUnit( resistance, - resistance_unit + resistanceUnit ) - for cross_section, resistance in tunnel_data["values"] + for cross_section, resistance in tunnelData["values"] } else: self.tunnelCoefficients = {} + + if "ETCS braking data" in data: + + etcsData = data["ETCS braking data"] + + requiredFields = { + "A_brake_emergency", + "A_brake_service", + "K_dry_rst", + "K_wet_rst", + "T_traction", + "T_be", + "v_uncertainty", + "T_bs" + } + + if set(etcsData) != requiredFields: + + missingFields = requiredFields - set(etcsData) + redundantFields = set(etcsData) - requiredFields + + raise ValueError( + "ETCS braking data must contain exactly the following fields: {}! " + "Missing fields: {}. Redundant fields: {}.".format( + ", ".join(sorted(requiredFields)), + ", ".join(sorted(missingFields)) if missingFields else "none", + ", ".join(sorted(redundantFields)) if redundantFields else "none", + ) + ) + + + def readBrakingCurve(curveName): + + brakingData = etcsData[curveName] + + if set(brakingData) != {"units", "values"}: + + raise ValueError("ETCS braking curve '{}' must contain 'units' and 'values'!".format(curveName)) + + velocityUnit = brakingData["units"]["velocity"] # velocity [m/s] + decelerationUnit = brakingData["units"]["deceleration"] # positive deceleration [m/s^2] + + return { + "velocity [m/s]": [ + convertUnit(velocity, velocityUnit) + for velocity, deceleration in brakingData["values"] + ], + "value [m/s^2]": [ + -convertUnit(deceleration, decelerationUnit) + for velocity, deceleration in brakingData["values"] + ] + } + + self.ABrakeEmergency = readBrakingCurve("A_brake_emergency") # emergency braking acceleration curve [m/s^2] + + self.ABrakeService = readBrakingCurve("A_brake_service") # service braking acceleration curve [m/s^2] + + self.KDryRst = convertUnit(etcsData["K_dry_rst"]["value"], etcsData["K_dry_rst"]["unit"]) # correction factor for dry rail conditions [-] + + self.KWetRst = convertUnit(etcsData["K_wet_rst"]["value"], etcsData["K_wet_rst"]["unit"]) # correction factor for wet rail conditions [-] + + self.TTraction = convertUnit(etcsData["T_traction"]["value"], etcsData["T_traction"]["unit"]) # traction cut-off delay [s] + + self.TBe = convertUnit(etcsData["T_be"]["value"], etcsData["T_be"]["unit"]) # emergency brake build-up time [s] + + self.vUncertainty = convertUnit(etcsData["v_uncertainty"]["value"], etcsData["v_uncertainty"]["unit"]) # train-related speed uncertainty [-] + + self.TBs = convertUnit(etcsData["T_bs"]["value"], etcsData["T_bs"]["unit"]) # service brake build-up time [s] + + else: + + self.ABrakeEmergency = None + + self.ABrakeService = None + + self.KDryRst = None + + self.KWetRst = None + + self.TTraction = None + + self.TBe = None + + self.vUncertainty = None + + self.TBs = None + + self.checkFields() @@ -196,7 +284,6 @@ def checkFields(self): raise ValueError("Rolling resistance coefficient {} must be positive, not {}!".format('r'+ii, coef)) - for crossSection, coef in self.tunnelCoefficients.items(): if crossSection is None or crossSection <= 0 or np.isinf(crossSection): @@ -207,6 +294,55 @@ def checkFields(self): raise ValueError("Tunnel resistance coefficient must be positive, not {}!".format(coef)) + for curveName in ["ABrakeEmergency", "ABrakeService"]: + + curve = getattr(self, curveName) + + if curve is not None: + + velocities = curve["velocity [m/s]"] + accelerations = curve["value [m/s^2]"] + + if len(velocities) == 0 or len(velocities) != len(accelerations): + + raise ValueError("ETCS braking curve {} has invalid length!".format(curveName)) + + if any(velocity < 0 or np.isinf(velocity) for velocity in velocities): + + raise ValueError("ETCS braking curve {} contains invalid velocities!".format(curveName)) + + if any(acceleration >= 0 or np.isinf(acceleration) for acceleration in accelerations): + + raise ValueError("ETCS braking curve {} contains invalid accelerations!".format(curveName)) + + if any(v2 <= v1 for v1, v2 in zip(velocities[:-1], velocities[1:])): + + raise ValueError("ETCS braking curve {} velocities must be strictly increasing!".format(curveName)) + + if self.KDryRst is not None and (self.KDryRst <= 0 or np.isinf(self.KDryRst)): + + raise ValueError("ETCS braking parameter 'K_dry_rst' must be positive, not {}!".format(self.KDryRst)) + + if self.KWetRst is not None and (self.KWetRst <= 0 or np.isinf(self.KWetRst)): + + raise ValueError("ETCS braking parameter 'K_wet_rst' must be positive, not {}!".format(self.KWetRst)) + + if self.TTraction is not None and (self.TTraction <= 0 or np.isinf(self.TTraction)): + + raise ValueError("ETCS braking parameter 'T_traction' must be positive, not {}!".format(self.TTraction)) + + if self.TBe is not None and (self.TBe <= 0 or np.isinf(self.TBe)): + + raise ValueError("ETCS braking parameter 'T_be' must be positive, not {}!".format(self.TBe)) + + if self.TBs is not None and (self.TBs <= 0 or np.isinf(self.TBs)): + + raise ValueError("ETCS braking parameter 'T_bs' must be positive, not {}!".format(self.TBs)) + + if self.vUncertainty is not None and (self.vUncertainty < 0 or self.vUncertainty > 1 or np.isinf(self.vUncertainty)): + + raise ValueError("ETCS braking parameter 'v_uncertainty' must be between 0 and 1, not {}!".format(self.vUncertainty)) + def exportModel(self): "Export train model (ODE and relevant train data)." diff --git a/mseetc/utils.py b/mseetc/utils.py index ae7b56e..b4cb866 100644 --- a/mseetc/utils.py +++ b/mseetc/utils.py @@ -372,7 +372,7 @@ def convertUnit(value, unit): Convert from any known unit to internally used unit. """ - if unit in {'m', 'm/s', 'm^2', 'permil', 'kg', 'W', 'N', 'm/s^2', '-', 'N/(m/s)', 'N/(m/s)^2', 'kg/m'}: + if unit in {'s', 'm', 'm/s', 'm^2', 'permil', 'kg', 'W', 'N', 'm/s^2', '-', 'N/(m/s)', 'N/(m/s)^2', 'kg/m'}: valueOut = value diff --git a/trains/CH_Stadler_FLIRT_TPF.json b/trains/CH_Stadler_FLIRT_TPF.json index 9443e2b..bcd2c76 100644 --- a/trains/CH_Stadler_FLIRT_TPF.json +++ b/trains/CH_Stadler_FLIRT_TPF.json @@ -92,5 +92,55 @@ 0.0005579 ] ] + }, + "ETCS braking data": { + "A_brake_emergency": { + "units": { + "velocity": "m/s", + "deceleration": "m/s^2" + }, + "values": [ + [0, 0.9], + [20, 0.85], + [40, 0.8], + [60, 0.75] + ] + }, + "A_brake_service": { + "units": { + "velocity": "m/s", + "deceleration": "m/s^2" + }, + "values": [ + [0, 0.55], + [20, 0.5], + [40, 0.45], + [60, 0.4] + ] + }, + "K_dry_rst": { + "unit": "-", + "value": 0.8 + }, + "K_wet_rst": { + "unit": "-", + "value": 0.9 + }, + "T_traction": { + "unit": "s", + "value": 1 + }, + "T_be": { + "unit": "s", + "value": 4 + }, + "v_uncertainty": { + "unit": "%", + "value": 2.98 + }, + "T_bs": { + "unit": "s", + "value": 3 + } } } From 6db2411a0fd2ed5b63d67a231d60c981b8e397b2 Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Wed, 10 Jun 2026 13:58:49 +0200 Subject: [PATCH 05/13] new ttobench version with etcs data added --- mseetc/etcs.py | 79 +++++++++++++------------------------ mseetc/track.py | 44 ++++++++++++++++++++- simulations/sim_launcher.py | 31 ++------------- tracks/CH_StGallen_Wil.json | 10 +++++ 4 files changed, 84 insertions(+), 80 deletions(-) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index 7f1123b..194b9f8 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -50,9 +50,9 @@ def getBrakingTargetsFromSpeedLimits(track): return targets, speedLimitPositions, speedLimits -def getEtcsSpeedLimits(trainBrakingData, track, positionStep=20.0): +def getEtcsSpeedLimits(train, track, positionStep=20.0): - calculator = EtcsBrakingCurveCalculator(trainBrakingData, track) + calculator = EtcsBrakingCurveCalculator(train, track) targets, speedLimitPositions, speedLimits = getBrakingTargetsFromSpeedLimits(track) @@ -211,9 +211,9 @@ class EtcsBrakingCurveCalculator: - Curves are computed backwards from the target position. """ - def __init__(self, trainBrakingData, track, distancePre=3000, distancePost=1000): + def __init__(self, train, track, distancePre=3000, distancePost=1000): - self.trainBrakingData = trainBrakingData + self.train = train self.track = track # Compute the braking curve using fixed time steps of length dt. @@ -257,16 +257,12 @@ def validateInput(self, target): def computeABrakeSafe(self): - trainBrakingData = self.trainBrakingData + velocities = self.train.ABrakeEmergency["velocity [m/s]"] + A_emergency_values = self.train.ABrakeEmergency["value [m/s^2]"] - braking = trainBrakingData["A_brake_emergency [m/s^2]"] - - velocities = braking["velocity [m/s]"] - A_emergency_values = braking["value [m/s^2]"] - - K_dry_rst = trainBrakingData["K_dry_rst [-]"] - M_NVAVADH = trainBrakingData["M_NVAVADH [-]"] - K_wet_rst = trainBrakingData["K_wet_rst [-]"] + K_dry_rst = self.train.KDryRst + M_NVAVADH = self.track.MNvavadh + K_wet_rst = self.train.KWetRst K_wet_corr = K_wet_rst + M_NVAVADH * (1 - K_wet_rst) @@ -364,12 +360,10 @@ def computeBrakingCurve(self, brakingProfile, targetPosition, permittedVelocity, def computeEBICurve(self, EBD_curve, targetVelocity): - trainBrakingData = self.trainBrakingData - - T_traction = trainBrakingData["T_traction [s]"] - T_be = trainBrakingData["T_be [s]"] - Kt_int = trainBrakingData["Kt_int [-]"] - v_uncertainty = trainBrakingData["v_uncertainty [%]"] + T_traction = self.train.TTraction + T_be = self.train.TBe + Kt_int = self.track.KtInt + v_uncertainty = self.train.vUncertainty positionsEBD = EBD_curve.index.to_numpy() velocitiesEBD = EBD_curve["Velocity [m/s]"].to_numpy() @@ -493,7 +487,7 @@ def processCurvesBeforeTarget(self, curves, target): return curves - def processCurvcesAfterTarget(self, curves, target): + def processCurvesAfterTarget(self, curves, target): targetVelocity = target.targetVelocity end_position = target.position + self.distancePost @@ -527,10 +521,9 @@ def processCurvcesAfterTarget(self, curves, target): def computeTarget(self, target): self.validateInput(target) - trainBrakingData = self.trainBrakingData ABrakeSafeProfile = self.computeABrakeSafe() - T_indication = max(0.8 * trainBrakingData["T_bs [s]"], 5) + self.T_driver + T_indication = max(0.8 * self.train.TBs, 5) + self.T_driver curves = {} @@ -538,11 +531,11 @@ def computeTarget(self, target): curves["EBI"] = self.computeEBICurve(curves["EBD"], target.targetVelocity) - curves["SBI2"] = shiftCurveByTime(curves["EBI"], trainBrakingData["T_bs [s]"]) + curves["SBI2"] = shiftCurveByTime(curves["EBI"], self.train.TBs) - curves["SBD"] = self.computeBrakingCurve(self.trainBrakingData["A_brake_service [m/s^2]"], target.EoA, target.permittedVelocity, target.targetVelocity) + curves["SBD"] = self.computeBrakingCurve(self.train.ABrakeService, target.EoA, target.permittedVelocity, target.targetVelocity) - curves["SBI1"] = shiftCurveByTime(curves["SBD"], trainBrakingData["T_bs [s]"]) + curves["SBI1"] = shiftCurveByTime(curves["SBD"], self.train.TBs) curves["SBI"] = self.computeSBICurve(curves["SBI1"], curves["SBI2"]) @@ -559,12 +552,6 @@ def computeTarget(self, target): def plotCurves(self, curves, target): - # curves = self.processCurvesBeforeTarget(curves, target) - # - # if target.targetVelocity > 0: - # - # curves = self.processCurvcesAfterTarget(curves, target) - targetPosition = target.EoA permittedVelocity = target.permittedVelocity targetVelocity = target.targetVelocity @@ -602,25 +589,6 @@ def plotCurves(self, curves, target): from mseetc.track import Track from mseetc.train import Train - trainBrakingData = { - "A_brake_emergency [m/s^2]": { # [train] - "velocity [m/s]": [0, 20, 40, 60], - "value [m/s^2]": [-0.9, -0.85, -0.8, -0.75], - }, - "A_brake_service [m/s^2]": { # [train] - "velocity [m/s]": [0, 20, 40, 60], - "value [m/s^2]": [-0.55, -0.5, -0.45, -0.4], - }, - "K_dry_rst [-]": 0.8, # [train] - "M_NVAVADH [-]": 0, # [infra] - "K_wet_rst [-]": 0.9, # [train] - "T_traction [s]": 1, # [train] - "T_be [s]": 4, # [train] - "Kt_int [-]": 1.15, # [infra] - "v_uncertainty [%]": 2.98, # [train] - "T_bs [s]": 3 # [train] - } - train = Train(config={'id': 'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') @@ -633,7 +601,16 @@ def plotCurves(self, curves, target): targetVelocity=0/3.6 ) - calculator = EtcsBrakingCurveCalculator(trainBrakingData, track, distancePre=5000, distancePost=1000) + addConstantVelocitySections = False + + calculator = EtcsBrakingCurveCalculator(train, track, distancePre=5000, distancePost=1000) curve_set = calculator.computeTarget(target) + if addConstantVelocitySections: + + curve_set = calculator.processCurvesBeforeTarget(curve_set, target) + + if target.targetVelocity > 0: + curve_set = calculator.processCurvesAfterTarget(curve_set, target) + calculator.plotCurves(curve_set, target) \ No newline at end of file diff --git a/mseetc/track.py b/mseetc/track.py index 8617e36..b5429a5 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -250,6 +250,38 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'tracks'): data['tunnels']['units']['length'] if 'tunnels' in data else 'm', data['tunnels']['units']['cross section'] if 'tunnels' in data else 'm^2') + if "ETCS braking data" in data: + + etcsData = data["ETCS braking data"] + + requiredFields = { + "M_NVAVADH", + "Kt_int" + } + + if set(etcsData) != requiredFields: + + missingFields = requiredFields - set(etcsData) + redundantFields = set(etcsData) - requiredFields + + raise ValueError( + "ETCS braking data must contain exactly the following fields: {}! " + "Missing fields: {}. Redundant fields: {}.".format( + ", ".join(sorted(requiredFields)), + ", ".join(sorted(missingFields)) if missingFields else "none", + ", ".join(sorted(redundantFields)) if redundantFields else "none", + ) + ) + + self.MNvavadh = convertUnit(etcsData["M_NVAVADH"]["value"], etcsData["M_NVAVADH"]["unit"]) + + self.KtInt = convertUnit(etcsData["Kt_int"]["value"], etcsData["Kt_int"]["unit"]) + + else: + + self.MNvavadh = None + self.KtInt = None + numStops = len(data['stops']['values']) indxDeparture = config['from'] if 'from' in config else 0 @@ -326,6 +358,14 @@ def checkFields(self): raise ValueError("Issue with tunnel cross sections!") + if self.MNvavadh is not None and (self.MNvavadh < 0 or np.isinf(self.MNvavadh)): + + raise ValueError("ETCS braking parameter 'M_NVAVADH' must be positive or zero, not {}!".format(self.MNvavadh)) + + if self.KtInt is not None and (self.KtInt <= 0 or np.isinf(self.KtInt)): + + raise ValueError("ETCS braking parameter 'Kt_int' must be positive, not {}!".format(self.KtInt)) + def importGradientTuples(self, tuples, unit='permil'): @@ -605,9 +645,9 @@ def crop(dfIn): self.crossSections = crop(self.crossSections) - def setEtcsSpeedLimits(self, trainBrakingData): + def setEtcsSpeedLimits(self, train): - etcsPositions, etcsVelocities = getEtcsSpeedLimits(trainBrakingData, self) + etcsPositions, etcsVelocities = getEtcsSpeedLimits(train, self) self.etcsPositions = etcsPositions self.etcsVelocities = etcsVelocities diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 1537d1d..2ff1ae7 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -40,30 +40,6 @@ def printStats(df, stats, solver): print("Solver failed!") -def getTrainBrakingData(): - - return { - "A_brake_emergency [m/s^2]": { - "velocity [m/s]": [0, 20, 40, 60], - "value [m/s^2]": [-0.9, -0.85, -0.8, -0.75], - }, - "A_brake_service [m/s^2]": { - "velocity [m/s]": [0, 20, 40, 60], - "value [m/s^2]": [-0.55, -0.5, -0.45, -0.4], - }, - "K_dry_rst [-]": 0.8, - "M_NVAVADH [-]": 0, - "K_wet_rst [-]": 0.9, - "T_traction [s]": 1, - "T_be [s]": 4, - "Kt_int [-]": 1.15, - "v_uncertainty [%]": 2.98, - "T_bs [s]": 3, - "T_bs1 [s]": 3, - "T_bs2 [s]": 3, - } - - if __name__ == '__main__': from mseetc.train import Train @@ -84,6 +60,7 @@ def getTrainBrakingData(): track.updateTrainLengthDependentValues(train) track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + # non-adjusted speed profile opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} solver = casadiSolver(train, track, opts) @@ -91,8 +68,8 @@ def getTrainBrakingData(): printStats(df, stats, solver) - trainBrakingData = getTrainBrakingData() - track.setEtcsSpeedLimits(trainBrakingData) + # ETCS-adjusted speed profile + track.setEtcsSpeedLimits(train) opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} solverEtcs = casadiSolver(train, track, opts) @@ -110,7 +87,7 @@ def getTrainBrakingData(): x_plot = np.append(x, track.length) v_plot = np.append(v, v[-1]) - etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits(trainBrakingData, track) + etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits(train, track) ax.step(x_plot/1000, v_plot*3.6, where="post", color="black", linestyle="-", label="Track Speed Limit") ax.plot(etcsLimitsPositions/1000, etcsLimitsVelocities*3.6, color="red", linestyle="-", label="ETCS Speed Limit") diff --git a/tracks/CH_StGallen_Wil.json b/tracks/CH_StGallen_Wil.json index c2afc63..fe0915f 100644 --- a/tracks/CH_StGallen_Wil.json +++ b/tracks/CH_StGallen_Wil.json @@ -1894,5 +1894,15 @@ -901.4 ] ] + }, + "ETCS braking data": { + "M_NVAVADH": { + "unit": "-", + "value": 0 + }, + "Kt_int": { + "unit": "-", + "value": 1.15 + } } } \ No newline at end of file From e3cac3acea759f5c37d2bcd6ad79e5fae3efd1a2 Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Wed, 10 Jun 2026 15:24:03 +0200 Subject: [PATCH 06/13] wip --- journeys/CH_StGallen_Wil_Journey_01.json | 43 ++++++++ mseetc/journey.py | 130 +++++++++++++++++++++++ simulations/sim_launcher.py | 5 +- 3 files changed, 177 insertions(+), 1 deletion(-) create mode 100644 journeys/CH_StGallen_Wil_Journey_01.json create mode 100644 mseetc/journey.py diff --git a/journeys/CH_StGallen_Wil_Journey_01.json b/journeys/CH_StGallen_Wil_Journey_01.json new file mode 100644 index 0000000..b1ed319 --- /dev/null +++ b/journeys/CH_StGallen_Wil_Journey_01.json @@ -0,0 +1,43 @@ +{ + "metadata": { + "id": "CH_StGallen_Wil_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License" + }, + "associated track id": { + "id": "CH_StGallen_Wil" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + { + "position": 0, + "lower time constraint": 0, + "upper time constraint": 0, + "lower speed constraint": 0, + "upper speed constraint": 0 + }, + { + "position": 12500, + "lower time constraint": 480, + "upper time constraint": 540, + "lower speed constraint": 10, + "upper speed constraint": null + }, + { + "position": 25000, + "lower time constraint": 1020, + "upper time constraint": 1080, + "lower speed constraint": 0, + "upper speed constraint": 0 + } + ] + } +} \ No newline at end of file diff --git a/mseetc/journey.py b/mseetc/journey.py new file mode 100644 index 0000000..2b3b018 --- /dev/null +++ b/mseetc/journey.py @@ -0,0 +1,130 @@ +import json +import numpy as np +import pandas as pd + +from pathlib import Path + +from mseetc.utils import checkTTOBenchVersion, convertUnit + + +class Journey(): + + def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'journeys') -> None: + """ + Constructor of Journey objects. + """ + + # check config + if not isinstance(config, dict): + + raise ValueError("Journey configuration should be provided as a dictionary!") + + if 'id' not in config: + + raise ValueError("Journey ID must be specified in configuration!") + + # open json file + filename = Path(pathJSON) / (config['id'] + '.json') + + with open(filename) as file: + + data = json.load(file) + + checkTTOBenchVersion(data, ['1.5']) + + # read data + self.id = data['metadata']['id'] + + self.associatedTrackID = data['associated track id']['id'] + + self.timingPoints = self.readTimingPoints(data['timing points']) + + self.checkFields() + + + def readTimingPoints(self, timingPoints): + + units = timingPoints['units'] + + positionUnit = units['position'] + lowerTimeUnit = units['lower time constraint'] + upperTimeUnit = units['upper time constraint'] + lowerSpeedUnit = units['lower speed constraint'] + upperSpeedUnit = units['upper speed constraint'] + + values = { + "Position [m]": [], + "Lower time constraint [s]": [], + "Upper time constraint [s]": [], + "Lower speed constraint [m/s]": [], + "Upper speed constraint [m/s]": [] + } + + for point in timingPoints['values']: + + values["Position [m]"].append(convertUnit(point['position'], positionUnit)) + values["Lower time constraint [s]"].append(self.convertConstraint(point['lower time constraint'], lowerTimeUnit)) + values["Upper time constraint [s]"].append(self.convertConstraint(point['upper time constraint'], upperTimeUnit)) + values["Lower speed constraint [m/s]"].append(self.convertConstraint(point['lower speed constraint'], lowerSpeedUnit)) + values["Upper speed constraint [m/s]"].append(self.convertConstraint(point['upper speed constraint'], upperSpeedUnit)) + + df = pd.DataFrame(values) + df = df.set_index("Position [m]") + + return df + + + def convertConstraint(self, value, unit): + + if value is None: + + return None + + return convertUnit(value, unit) + + + def checkFields(self): + + if len(self.timingPoints) < 2: + + raise ValueError("Journey must contain at least two timing points!") + + positions = self.timingPoints.index.values + + if any(position < 0 or np.isinf(position) for position in positions): + + raise ValueError("Timing point positions must be positive finite numbers!") + + if any(pos2 <= pos1 for pos1, pos2 in zip(positions[:-1], positions[1:])): + + raise ValueError("Timing point positions must be strictly increasing!") + + firstPoint = self.timingPoints.iloc[0] + lastPoint = self.timingPoints.iloc[-1] + + if firstPoint["Lower speed constraint [m/s]"] != 0 or firstPoint["Upper speed constraint [m/s]"] != 0: + + raise ValueError("First timing point must have both speed constraints set to zero!") + + if lastPoint["Lower speed constraint [m/s]"] != 0 or lastPoint["Upper speed constraint [m/s]"] != 0: + + raise ValueError("Last timing point must have both speed constraints set to zero!") + + for ii, point in self.timingPoints.iterrows(): + + tMin = point["Lower time constraint [s]"] + tMax = point["Upper time constraint [s]"] + vMin = point["Lower speed constraint [m/s]"] + vMax = point["Upper speed constraint [m/s]"] + + if tMin is not None and tMax is not None and tMin > tMax: + + raise ValueError("Lower time constraint must be smaller than or equal to upper time constraint!") + + if vMin is not None and vMax is not None and vMin > vMax: + + raise ValueError("Lower speed constraint must be smaller than or equal to upper speed constraint!") + + if any(value is not None and (value < 0 or np.isinf(value)) for value in [tMin, tMax, vMin, vMax]): + + raise ValueError("Timing point constraints must be positive finite numbers or None!") \ No newline at end of file diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 2ff1ae7..b8afef2 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -3,6 +3,7 @@ from mseetc.efficiency import totalLossesFunction from mseetc.etcs import getEtcsSpeedLimits +from mseetc.journey import Journey from mseetc.ocp import casadiSolver @@ -55,9 +56,11 @@ def printStats(df, stats, solver): train.withPnBrake = False train.powerLosses = get_power_loss_function(train, "static") - track = Track(config={'id':'00_var_speed_limit_quick_change'}, pathJSON='../tracks') + # track = Track(config={'id':'00_var_speed_limit_quick_change'}, pathJSON='../tracks') track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') track.updateTrainLengthDependentValues(train) + + journey = Journey(config={'id':'CH_StGallen_Wil_Journey_01'}, pathJSON='../journeys') track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') # non-adjusted speed profile From fd28fcfcd295558b3527916321ba26e619d591bf Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Thu, 11 Jun 2026 12:12:13 +0200 Subject: [PATCH 07/13] timingPoints added --- journeys/CH_StGallen_Wil_Journey_02.json | 50 ++++++++++ mseetc/journey.py | 51 +++++++--- mseetc/ocp.py | 78 +++++++++++---- mseetc/track.py | 12 +-- mseetc/utils.py | 6 ++ simulations/sim_launcher.py | 12 +-- simulations/sim_launcher_timingPoints.py | 118 +++++++++++++++++++++++ 7 files changed, 283 insertions(+), 44 deletions(-) create mode 100644 journeys/CH_StGallen_Wil_Journey_02.json create mode 100644 simulations/sim_launcher_timingPoints.py diff --git a/journeys/CH_StGallen_Wil_Journey_02.json b/journeys/CH_StGallen_Wil_Journey_02.json new file mode 100644 index 0000000..a955c77 --- /dev/null +++ b/journeys/CH_StGallen_Wil_Journey_02.json @@ -0,0 +1,50 @@ +{ + "metadata": { + "id": "CH_StGallen_Wil_Journey_02", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License" + }, + "associated track id": { + "id": "CH_StGallen_Wil" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + { + "position": 1000, + "lower time constraint": 100, + "upper time constraint": 200, + "lower speed constraint": 0, + "upper speed constraint": 0 + }, + { + "position": 12500, + "lower time constraint": 480, + "upper time constraint": 540, + "lower speed constraint": 110, + "upper speed constraint": 120 + }, + { + "position": 18000, + "lower time constraint": 800, + "upper time constraint": null, + "lower speed constraint": null, + "upper speed constraint": 70 + }, + { + "position": 25000, + "lower time constraint": 900, + "upper time constraint": 1100, + "lower speed constraint": 0, + "upper speed constraint": 0 + } + ] + } +} \ No newline at end of file diff --git a/mseetc/journey.py b/mseetc/journey.py index 2b3b018..09cd631 100644 --- a/mseetc/journey.py +++ b/mseetc/journey.py @@ -41,17 +41,15 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'journeys') - self.checkFields() + self.computeStartAndEndPoints() + + self.computeInitialAndTerminalStates() + def readTimingPoints(self, timingPoints): units = timingPoints['units'] - positionUnit = units['position'] - lowerTimeUnit = units['lower time constraint'] - upperTimeUnit = units['upper time constraint'] - lowerSpeedUnit = units['lower speed constraint'] - upperSpeedUnit = units['upper speed constraint'] - values = { "Position [m]": [], "Lower time constraint [s]": [], @@ -62,11 +60,11 @@ def readTimingPoints(self, timingPoints): for point in timingPoints['values']: - values["Position [m]"].append(convertUnit(point['position'], positionUnit)) - values["Lower time constraint [s]"].append(self.convertConstraint(point['lower time constraint'], lowerTimeUnit)) - values["Upper time constraint [s]"].append(self.convertConstraint(point['upper time constraint'], upperTimeUnit)) - values["Lower speed constraint [m/s]"].append(self.convertConstraint(point['lower speed constraint'], lowerSpeedUnit)) - values["Upper speed constraint [m/s]"].append(self.convertConstraint(point['upper speed constraint'], upperSpeedUnit)) + values["Position [m]"].append(convertUnit(point['position'], units['position'])) + values["Lower time constraint [s]"].append(self.convertConstraint(point['lower time constraint'], units['lower time constraint'])) + values["Upper time constraint [s]"].append(self.convertConstraint(point['upper time constraint'], units['upper time constraint'])) + values["Lower speed constraint [m/s]"].append(self.convertConstraint(point['lower speed constraint'], units['lower speed constraint'])) + values["Upper speed constraint [m/s]"].append(self.convertConstraint(point['upper speed constraint'], units['upper speed constraint'])) df = pd.DataFrame(values) df = df.set_index("Position [m]") @@ -127,4 +125,33 @@ def checkFields(self): if any(value is not None and (value < 0 or np.isinf(value)) for value in [tMin, tMax, vMin, vMax]): - raise ValueError("Timing point constraints must be positive finite numbers or None!") \ No newline at end of file + raise ValueError("Timing point constraints must be positive finite numbers or None!") + + + def computeStartAndEndPoints(self): + + self.positionStart = self.timingPoints.index.values[0] + + self.positionEnd = self.timingPoints.index.values[-1] + + + def computeInitialAndTerminalStates(self): + + self.initialTime = self.timingPoints.iloc[0]["Lower time constraint [s]"] + + self.terminalTime = self.timingPoints.iloc[-1]["Upper time constraint [s]"] + + self.initialVelocity = self.timingPoints.iloc[0]["Lower speed constraint [m/s]"] + + self.terminalVelocity = self.timingPoints.iloc[-1]["Upper speed constraint [m/s]"] + + + def getTimingPoint(self, position): + + idx = np.where(np.isclose(self.timingPoints.index.values, position))[0] + + if len(idx) == 0: + return None + + return self.timingPoints.iloc[idx[0]] + diff --git a/mseetc/ocp.py b/mseetc/ocp.py index 3a6a427..f8b0d06 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -4,7 +4,7 @@ from mseetc.track import computeDiscretizationPoints -from mseetc.utils import Options, var, postProcessDataFrame, splitLosses, computeTunnelFactor +from mseetc.utils import Options, var, postProcessDataFrame, splitLosses, computeTunnelFactor, isSet class OptionsCasadiSolver(Options): @@ -81,11 +81,12 @@ def checkValues(self): class casadiSolver(): "NLP solver object using casadi and ipopt." - def __init__(self, train, track, optsDict={}): + def __init__(self, train, track, journey, optsDict={}): # input checking track.checkFields() train.checkFields() + journey.checkFields() opts = OptionsCasadiSolver(optsDict) @@ -125,7 +126,8 @@ def __init__(self, train, track, optsDict={}): # track parameters - self.points = computeDiscretizationPoints(track, numIntervals, opts) + timingPointPositions = journey.timingPoints.index.to_numpy(dtype=float)[1:-1] - journey.positionStart + self.points = computeDiscretizationPoints(track, numIntervals, opts, timingPointPositions) self.steps = np.diff(self.points.index) if opts.withEtcsBrakingCurves: @@ -269,28 +271,60 @@ def __init__(self, train, track, optsDict={}): z += [time[i]] z += [velSq[i]] - if i == 0: + position = self.points.index.values[i] + timingPoint = journey.getTimingPoint(position + journey.positionStart) - # initial state constraints - lbz += [self.initialTime, self.initialVelocitySquared] - ubz += [self.initialTime, self.initialVelocitySquared] + speedLimit = self.points.iloc[i]['Speed limit [m/s]'] + speedLimit = min(speedLimit, velocityMax) - elif i == numIntervals: + if i > 0: + + speedLimit = min(speedLimit, self.points.iloc[i-1]['Speed limit [m/s]']) # do not accelerate before speed limit increase - # terminal state constraints - lbz += [self.initialTime, self.terminalVelocitySquared] - ubz += [self.terminalTime, self.terminalVelocitySquared] + timeLower = self.initialTime + timeUpper = self.terminalTime - else: + velocityLower = velocityMin ** 2 + velocityUpper = speedLimit ** 2 - # state constraints - speedLimit = self.points.iloc[i]['Speed limit [m/s]'] - speedLimit = min(speedLimit, velocityMax) + if timingPoint is not None: - speedLimit = min(speedLimit, self.points.iloc[i-1]['Speed limit [m/s]']) # do not accelerate before speed limit increase + timeMinPoint = timingPoint["Lower time constraint [s]"] + timeMaxPoint = timingPoint["Upper time constraint [s]"] + velocityMinPoint = timingPoint["Lower speed constraint [m/s]"] + velocityMaxPoint = timingPoint["Upper speed constraint [m/s]"] + + if isSet(timeMinPoint): + + timeLower = timeMinPoint + + if isSet(timeMaxPoint): + + timeUpper = timeMaxPoint + + if isSet(velocityMinPoint): + + velocityLower = max(velocityLower, velocityMinPoint ** 2) + + if isSet(velocityMaxPoint): + + velocityUpper = min(velocityUpper, velocityMaxPoint ** 2) + + if i == 0: + + timeLower = self.initialTime + velocityLower = self.initialVelocitySquared + velocityUpper = self.initialVelocitySquared + + elif i == numIntervals: + + timeUpper = self.terminalTime + velocityLower = self.terminalVelocitySquared + velocityUpper = self.terminalVelocitySquared + + lbz += [timeLower, velocityLower] + ubz += [timeUpper, velocityUpper] - lbz += [self.initialTime, velocityMin**2] - ubz += [self.terminalTime, speedLimit**2] # scaling of objective function (fixes convergence issues when using powerLosses) @@ -328,7 +362,13 @@ def __init__(self, train, track, optsDict={}): self.ubg = ca.vcat(ubg) - def solve(self, terminalTime, initialTime=0, terminalVelocity=1, initialVelocity=1): + def solve(self, journey): + + initialTime = journey.initialTime + terminalTime = journey.terminalTime + initialVelocity = journey.initialVelocity + terminalVelocity = journey.terminalVelocity + # check boundary conditions on time diff --git a/mseetc/track.py b/mseetc/track.py index b5429a5..e193f1a 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -91,7 +91,7 @@ def computeAltitude(gradients, length, altitudeStart=0): return df -def computeDiscretizationPoints(track, numIntervals, opts): +def computeDiscretizationPoints(track, numIntervals, opts, timingPointPositions): """ Compute the space discretization points based on track characteristics and horizon length. """ @@ -99,6 +99,7 @@ def computeDiscretizationPoints(track, numIntervals, opts): df1 = track.mergeDataFrames() requiredPoints = df1.index.to_numpy(dtype=float) + requiredPoints = np.append(requiredPoints, timingPointPositions) if opts.withEtcsBrakingCurves: @@ -616,16 +617,13 @@ def updateLimits(self, positionStart=None, positionEnd=None, unit='m'): Truncate track to given positions. """ - positionStart = 0 if positionStart is None else positionStart - positionEnd = self.length if positionEnd is None else positionEnd + positionStart = 0 if positionStart is None else convertUnit(positionStart, unit) + positionEnd = self.length if positionEnd is None else convertUnit(positionEnd, unit) if (not 0 <= positionStart < self.length) or (not 0 < positionEnd <= self.length) : raise ValueError("Given positions must be between limits of track!") - positionStart = convertUnit(positionStart, unit) - positionEnd = convertUnit(positionEnd, unit) - newPos = pd.DataFrame({'Position [m]':[positionStart]}).set_index('Position [m]') def crop(dfIn): @@ -637,7 +635,7 @@ def crop(dfIn): return dfOut - self.length -= positionStart + (self.length - positionEnd) + self.length = positionEnd - positionStart self.speedLimits = crop(self.speedLimits) self.gradients = crop(self.gradients) diff --git a/mseetc/utils.py b/mseetc/utils.py index b4cb866..c39fa51 100644 --- a/mseetc/utils.py +++ b/mseetc/utils.py @@ -7,6 +7,7 @@ from types import MethodType from distutils.spawn import find_executable +import pandas as pd from matplotlib import pyplot as plt @@ -644,6 +645,11 @@ def getRelevantEtcsBrakingPositions(track): return relevantPositions +def isSet( value): + + return value is not None and not pd.isna(value) + + if __name__ == '__main__': pass diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index b8afef2..9b235b2 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -26,7 +26,7 @@ def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000 raise ValueError("mode must be one of: 'perfect', 'static', 'dynamic'") -def printStats(df, stats, solver): +def printStats(df, stats, solver, train): if df is not None: @@ -61,15 +61,15 @@ def printStats(df, stats, solver): track.updateTrainLengthDependentValues(train) journey = Journey(config={'id':'CH_StGallen_Wil_Journey_01'}, pathJSON='../journeys') - track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') # non-adjusted speed profile opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} - solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve(terminalTime=journey.terminalTime, initialTime=journey.initialTime) - printStats(df, stats, solver) + printStats(df, stats, solver, train) # ETCS-adjusted speed profile track.setEtcsSpeedLimits(train) @@ -78,7 +78,7 @@ def printStats(df, stats, solver): solverEtcs = casadiSolver(train, track, opts) dfEtcs, statsEtcs = solverEtcs.solve(duration) - printStats(dfEtcs, statsEtcs, solverEtcs) + printStats(dfEtcs, statsEtcs, solverEtcs, train) ### Plot Trajectory diff --git a/simulations/sim_launcher_timingPoints.py b/simulations/sim_launcher_timingPoints.py new file mode 100644 index 0000000..53f83bd --- /dev/null +++ b/simulations/sim_launcher_timingPoints.py @@ -0,0 +1,118 @@ +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.patches import Rectangle + +from mseetc.ocp import casadiSolver +from mseetc.utils import isSet +from simulations.sim_launcher import get_power_loss_function, printStats + + +def plotTpConstraintsTime(ax, journey): + + posShift = journey.positionEnd - journey.positionStart + timeShift = journey.terminalTime - journey.initialTime + + for pos, tP in journey.timingPoints.iterrows(): + + lowerTime = tP["Lower time constraint [s]"] + upperTime = tP["Upper time constraint [s]"] + + if isSet(lowerTime): + + ax.add_patch(Rectangle((pos/1000, lowerTime-timeShift), posShift/1000, timeShift, facecolor='lightgrey', edgecolor='none')) + + if isSet(upperTime): + + ax.add_patch(Rectangle(((pos-posShift) / 1000, upperTime), posShift / 1000, timeShift, facecolor='lightgrey', edgecolor='none')) + + +def plotTpConstraintsVelocity(ax, journey): + + maxVel = 400 # [km/h] + + for pos, tP in journey.timingPoints.iterrows(): + + lowerVel = tP["Lower speed constraint [m/s]"] + upperVel = tP["Upper speed constraint [m/s]"] + + if isSet(lowerVel): + + ax.vlines(pos/1000, 0, lowerVel*3.6, color='red', linewidth=1) + + if isSet(upperVel): + + ax.vlines(pos/1000, upperVel*3.6, maxVel, color='red', linewidth=1) + + +def plotPositionTimeProfile(df, journey): + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.plot((df["Position [m]"]+journey.positionStart)/ 1000, df.index.to_numpy(), linestyle="--") + plotTpConstraintsTime(ax, journey) + + ax.set_title("Position - Time - Profile") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Time [s]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.set_xlim(journey.positionStart / 1000 - 1, journey.positionEnd / 1000 + 1) + ax.set_ylim(journey.initialTime - 100, journey.terminalTime + 100) + ax.figure.tight_layout() + + plt.show() + + +def plotPositionVelocityProfile(df, journey): + + x = track.speedLimits.index.to_numpy(dtype=float) + v = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + x_plot = np.append(x, track.length) + v_plot = np.append(v, v[-1]) + + vMax = max(v.max(), df["Velocity [m/s]"].max()) * 3.6 + 20 + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.plot((df["Position [m]"]+journey.positionStart) / 1000, df["Velocity [m/s]"] * 3.6, linestyle="--") + ax.step((x_plot+journey.positionStart) / 1000, v_plot * 3.6, where="post", color="black", linestyle="-", label="Track Speed Limit") + plotTpConstraintsVelocity(ax, journey) + + ax.set_title("Position - Velocity - Profile") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Velocity [km/h]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.set_xlim(journey.positionStart / 1000 - 1, journey.positionEnd / 1000 + 1) + ax.set_ylim(0, vMax) + ax.legend(loc="upper right") + ax.figure.tight_layout() + + plt.show() + + +if __name__ == '__main__': + + from mseetc.train import Train + from mseetc.track import Track + from mseetc.journey import Journey + + + train = Train(config={'id':'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') + train.forceMinPn = 0 + train.withPnBrake = False + train.powerLosses = get_power_loss_function(train, "static") + + track = Track(config={'id':'CH_StGallen_Wil'}, pathJSON='../tracks') + track.updateTrainLengthDependentValues(train) + + journey = Journey(config={'id':'CH_StGallen_Wil_Journey_02'}, pathJSON='../journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + + opts = {'numIntervals':800, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} + + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve(journey) + + printStats(df, stats, solver, train) + + plotPositionTimeProfile(df, journey) + plotPositionVelocityProfile(df, journey) \ No newline at end of file From 529a1c4dd03fdce0c1dfc86719db58857256a2cd Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Thu, 11 Jun 2026 13:45:28 +0200 Subject: [PATCH 08/13] good code --- .../00_var_speed_limit_100_Journey_01.json | 36 +++++++++++++++++++ journeys/CH_StGallen_Wil_Journey_01.json | 13 ++----- mseetc/ocp.py | 8 +++-- mseetc/train.py | 6 ++-- simulations/sim_launcher.py | 11 +++--- 5 files changed, 52 insertions(+), 22 deletions(-) create mode 100644 journeys/00_var_speed_limit_100_Journey_01.json diff --git a/journeys/00_var_speed_limit_100_Journey_01.json b/journeys/00_var_speed_limit_100_Journey_01.json new file mode 100644 index 0000000..8708964 --- /dev/null +++ b/journeys/00_var_speed_limit_100_Journey_01.json @@ -0,0 +1,36 @@ +{ + "metadata": { + "id": "00_var_speed_limit_100_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License" + }, + "associated track id": { + "id": "00_var_speed_limit_100" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + { + "position": 0, + "lower time constraint": 0, + "upper time constraint": 0, + "lower speed constraint": 0, + "upper speed constraint": 0 + }, + { + "position": 20000, + "lower time constraint": 685, + "upper time constraint": 685, + "lower speed constraint": 0, + "upper speed constraint": 0 + } + ] + } +} \ No newline at end of file diff --git a/journeys/CH_StGallen_Wil_Journey_01.json b/journeys/CH_StGallen_Wil_Journey_01.json index b1ed319..b8bef56 100644 --- a/journeys/CH_StGallen_Wil_Journey_01.json +++ b/journeys/CH_StGallen_Wil_Journey_01.json @@ -25,16 +25,9 @@ "upper speed constraint": 0 }, { - "position": 12500, - "lower time constraint": 480, - "upper time constraint": 540, - "lower speed constraint": 10, - "upper speed constraint": null - }, - { - "position": 25000, - "lower time constraint": 1020, - "upper time constraint": 1080, + "position": 20000, + "lower time constraint": 685, + "upper time constraint": 685, "lower speed constraint": 0, "upper speed constraint": 0 } diff --git a/mseetc/ocp.py b/mseetc/ocp.py index f8b0d06..27fc011 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -475,6 +475,7 @@ def solve(self, journey): from mseetc.train import Train from mseetc.track import Track + from mseetc.journey import Journey # Example on how to solve an OCP @@ -482,11 +483,14 @@ def solve(self, journey): track = Track(config={'id':'00_var_speed_limit_100'}) + journey = Journey(config={'id':'00_var_speed_limit_100_Journey_01'}, pathJSON='../journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + opts = {'numIntervals':200, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} - solver = casadiSolver(train, track, opts) + solver = casadiSolver(train, track, journey, opts) - df, stats = solver.solve(1541) + df, stats = solver.solve(journey) # print some info if df is not None: diff --git a/mseetc/train.py b/mseetc/train.py index 52e11fc..6a77652 100644 --- a/mseetc/train.py +++ b/mseetc/train.py @@ -66,7 +66,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> raise ValueError("Redundant fields in train configuration: {}!".format(', '.join(set(config) - usedFields))) # read data - self.length = convertUnit(data['length']['value'], data['length']['unit']) # train length [m] + self.length = convertUnit(data['length']['value'], data['length']['unit']) if "length" in data else None # train length [m] self.mass = convertUnit(data['mass']['value'], data['mass']['unit']) # train mass [kg] @@ -96,7 +96,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> self.r1 = convertUnit(data['rolling resistance r1']['value'], data['rolling resistance r1']['unit']) # linear term [N/(m/s)] - self.r2 = convertUnit(data['rolling resistance r2']['value'], data['rolling resistance r2']['unit']) # quadratic term [N/(m/s)^2] + self.r2 = convertUnit(data['rolling resistance r2']['value'], data['rolling resistance r2']['unit']) # quadratic term [N/(m/s)^2 # TODO: unify with case of dynamic efficiency if 'efficiency traction' in data or 'efficiency reg brake' in data: @@ -224,7 +224,7 @@ def readBrakingCurve(curveName): def checkFields(self): - if self.length is None or self.length < 0 or np.isinf(self.length): + if self.length is not None and (self.length < 0 or np.isinf(self.length)): raise ValueError("Train length must be a positive number, not {}!".format(self.length)) diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 9b235b2..d95a613 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -45,11 +45,8 @@ def printStats(df, stats, solver, train): from mseetc.train import Train from mseetc.track import Track + from mseetc.journey import Journey - # Timetable - startPosition = 0 # [m] - endPosition = 20000 # [m] - duration = (endPosition-startPosition) / (105/3.6) # [s] train = Train(config={'id':'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') train.forceMinPn = 0 @@ -67,7 +64,7 @@ def printStats(df, stats, solver, train): opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} solver = casadiSolver(train, track, journey, opts) - df, stats = solver.solve(terminalTime=journey.terminalTime, initialTime=journey.initialTime) + df, stats = solver.solve(journey) printStats(df, stats, solver, train) @@ -75,8 +72,8 @@ def printStats(df, stats, solver, train): track.setEtcsSpeedLimits(train) opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} - solverEtcs = casadiSolver(train, track, opts) - dfEtcs, statsEtcs = solverEtcs.solve(duration) + solverEtcs = casadiSolver(train, track, journey, opts) + dfEtcs, statsEtcs = solverEtcs.solve(journey) printStats(dfEtcs, statsEtcs, solverEtcs, train) From 661d794a37928d1bc617d3c6715892609dc4b7cc Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Thu, 11 Jun 2026 15:12:06 +0200 Subject: [PATCH 09/13] big fix missing position unit in track data --- mseetc/track.py | 43 ++++++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/mseetc/track.py b/mseetc/track.py index e193f1a..de99362 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -237,17 +237,22 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'tracks'): self.altitude = convertUnit(data['altitude']['value'], data['altitude']['unit']) if 'altitude' in data else 0 self.title = data['metadata']['id'] - self.importSpeedLimitTuples(data['speed limits']['values'], data['speed limits']['units']['velocity']) + self.importSpeedLimitTuples(data['speed limits']['values'], + data['speed limits']['units']['position'], + data['speed limits']['units']['velocity']) self.importGradientTuples(data['gradients']['values'] if 'gradients' in data else [(0.0, 0.0)], + data['gradients']['units']['position'] if 'gradients' in data else 'm', data['gradients']['units']['slope'] if 'gradients' in data else 'permil') self.importCurvatureTuples(data['curvatures']['values'] if 'curvatures' in data else [(0.0, "infinity", "infinity")], + data['curvatures']['units']['position'] if 'curvatures' in data else 'm', data['curvatures']['units']['radius at start'] if 'curvatures' in data else "m", data['curvatures']['units']['radius at end'] if 'curvatures' in data else "m", config['clothoidSamplingInterval'] if 'clothoidSamplingInterval' in config else None) self.importTunnelTuples(data['tunnels']['values'] if 'tunnels' in data else [(0.0, 0.0, "infinity")], + data['tunnels']['units']['position'] if 'tunnels' in data else 'm', data['tunnels']['units']['length'] if 'tunnels' in data else 'm', data['tunnels']['units']['cross section'] if 'tunnels' in data else 'm^2') @@ -368,49 +373,61 @@ def checkFields(self): raise ValueError("ETCS braking parameter 'Kt_int' must be positive, not {}!".format(self.KtInt)) - def importGradientTuples(self, tuples, unit='permil'): + def importGradientTuples(self, tuples, unitPosition='m', unitGradient='permil'): if not self.lengthOk(): raise ValueError("Cannot import gradients without a valid track length!") - if unit not in {'permil'}: + if unitPosition not in {'m', 'km'}: + + raise ValueError("Specified gradient position unit not supported!") + + if unitGradient not in {'permil'}: raise ValueError("Specified gradient unit not supported!") + tuples = [(convertUnit(p, unitPosition), g) for p, g in tuples] self.gradients = importTuples(tuples, 'Position [m]', 'Gradient [permil]') checkDataFrame(self.gradients, self.length) - def importSpeedLimitTuples(self, tuples, unit='km/h'): + def importSpeedLimitTuples(self, tuples, unitPosition='m', unitVelocity='km/h'): if not self.lengthOk(): raise ValueError("Cannot import speed limits without a valid track length!") - if unit not in {'km/h', 'm/s'}: + if unitPosition not in {'m', 'km'}: + + raise ValueError("Specified speed limit position unit not supported!") + + if unitVelocity not in {'km/h', 'm/s'}: raise ValueError("Specified speed unit not supported!") - tuples = [(p, convertUnit(v, unit)) for p,v in tuples] + tuples = [(convertUnit(p, unitPosition), convertUnit(v, unitVelocity)) for p, v in tuples] self.speedLimits = importTuples(tuples, 'Position [m]', 'Speed limit [m/s]') checkDataFrame(self.speedLimits, self.length) - - def importCurvatureTuples(self, tuples, unitRadiusStart='m', unitRadiusEnd='m', clothoidSamplingInterval=None): + def importCurvatureTuples(self, tuples, unitPosition='m', unitRadiusStart='m', unitRadiusEnd='m', clothoidSamplingInterval=None): if not self.lengthOk(): raise ValueError("Cannot import curvature without a valid track length!") + if unitPosition not in {'m', 'km'}: + + raise ValueError("Specified curvature position unit not supported!") + if unitRadiusStart not in {'m', 'km'} or unitRadiusEnd not in {'m', 'km'}: raise ValueError("Specified curvature radius unit not supported!") # if radius is 'infinity', the casting to float produces the float inf - tuples = [(p, convertUnit(float(radiusStart), unitRadiusStart), convertUnit(float(radiusEnd), unitRadiusEnd)) for p, radiusStart, radiusEnd in tuples] + tuples = [(convertUnit(p, unitPosition), convertUnit(float(radiusStart), unitRadiusStart), convertUnit(float(radiusEnd), unitRadiusEnd)) for p, radiusStart, radiusEnd in tuples] tuples = self.sampleClothoid(tuples, clothoidSamplingInterval) @@ -501,17 +518,21 @@ def sampleClothoid(self, tuples, ds=None): return result - def importTunnelTuples(self, tuples, unitLength='m', unitCrossSection='m^2'): + def importTunnelTuples(self, tuples, unitPosition='m', unitLength='m', unitCrossSection='m^2'): if not self.lengthOk(): raise ValueError("Cannot import tunnels without a valid track length!") + if unitPosition not in {'m', 'km'}: + + raise ValueError("Specified tunnel position unit not supported!") + if unitLength not in {'m', 'km'} or unitCrossSection not in {'m^2'}: raise ValueError("Specified tunnel units not supported!") - tuples = [(p, convertUnit(l, unitLength), convertUnit(c, unitCrossSection)) for p,l,c in tuples] + tuples = [(convertUnit(p, unitPosition), convertUnit(l, unitLength), convertUnit(c, unitCrossSection)) for p,l,c in tuples] self.crossSections = importTuples(tuples, 'Position [m]', ['Length [m]', 'CrossSection [m^2]']) From b886e6d653b5d5769eeb71bb245a33f6d8afc5e9 Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Fri, 12 Jun 2026 16:20:02 +0200 Subject: [PATCH 10/13] good code visualization and equal shooting nodes --- journeys/CH_StGallen_Wil_Journey_01.json | 4 +- mseetc/etcs.py | 4 +- mseetc/ocp.py | 2 +- mseetc/track.py | 19 +++-- mseetc/utils.py | 91 +++++++++++++++--------- simulations/sim_launcher.py | 35 +++++++-- 6 files changed, 108 insertions(+), 47 deletions(-) diff --git a/journeys/CH_StGallen_Wil_Journey_01.json b/journeys/CH_StGallen_Wil_Journey_01.json index b8bef56..006137d 100644 --- a/journeys/CH_StGallen_Wil_Journey_01.json +++ b/journeys/CH_StGallen_Wil_Journey_01.json @@ -26,8 +26,8 @@ }, { "position": 20000, - "lower time constraint": 685, - "upper time constraint": 685, + "lower time constraint": 700, + "upper time constraint": 700, "lower speed constraint": 0, "upper speed constraint": 0 } diff --git a/mseetc/etcs.py b/mseetc/etcs.py index 194b9f8..13fe289 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -597,8 +597,8 @@ def plotCurves(self, curves, target): target = BrakingTarget( position=5000, overlap= 100, - permittedVelocity=160/3.6, - targetVelocity=0/3.6 + permittedVelocity=130/3.6, + targetVelocity=100/3.6 ) addConstantVelocitySections = False diff --git a/mseetc/ocp.py b/mseetc/ocp.py index 27fc011..71aa6ef 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -133,7 +133,7 @@ def __init__(self, train, track, journey, optsDict={}): if opts.withEtcsBrakingCurves: velocities = np.interp(self.points.index.to_numpy(), track.etcsPositions, track.etcsVelocities) - self.points["Speed limit [m/s]"] = velocities + self.points["Speed limit [m/s]"] = np.maximum(velocities, velocityMin) # real-time parameters diff --git a/mseetc/track.py b/mseetc/track.py index de99362..a8094b1 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -6,8 +6,8 @@ import pandas as pd from mseetc.etcs import getEtcsSpeedLimits -from mseetc.utils import checkTTOBenchVersion, convertUnit, pickEquallySpacedPoints, plotSpeedLimits, plotGradients, \ - plotCurvatures, getRelevantEtcsBrakingPositions +from mseetc.utils import checkTTOBenchVersion, convertUnit, plotSpeedLimits, plotGradients, \ + plotCurvatures, introduceSufficientShootingNodesForETCSBrakingCurves, getEquallyDistributedPoints plotDebug = False @@ -98,15 +98,22 @@ def computeDiscretizationPoints(track, numIntervals, opts, timingPointPositions) df1 = track.mergeDataFrames() - requiredPoints = df1.index.to_numpy(dtype=float) - requiredPoints = np.append(requiredPoints, timingPointPositions) + requiredPoints = np.concatenate([df1.index.to_numpy(dtype=float), timingPointPositions]) if opts.withEtcsBrakingCurves: - brakingPoints = getRelevantEtcsBrakingPositions(track) + brakingPoints = introduceSufficientShootingNodesForETCSBrakingCurves(track, requiredPoints) requiredPoints = np.append(requiredPoints, brakingPoints) - pos = pickEquallySpacedPoints(0, track.length, numIntervals, requiredPoints) + requiredPoints = np.unique(np.append(requiredPoints, [0, track.length])) + numMissingPoints = numIntervals + 1 - len(requiredPoints) + + if numMissingPoints < 0: + + raise ValueError(f"Increase num of intervals by {numMissingPoints}!") + + newPositions = getEquallyDistributedPoints(0.0, track.length, numMissingPoints, requiredPoints) + pos = np.sort(np.unique(np.concatenate([requiredPoints, newPositions]))) df2 = pd.DataFrame({'position [m]':pos}).set_index('position [m]') df3 = df2.join(df1, how='outer').ffill() diff --git a/mseetc/utils.py b/mseetc/utils.py index c39fa51..23da746 100644 --- a/mseetc/utils.py +++ b/mseetc/utils.py @@ -521,32 +521,6 @@ def computeTunnelFactor(cross_section, train, opts): return tunnelCoefficient/total_mass -def pickEquallySpacedPoints(startPoint, endPoint, numIntervals, requiredPoints): - - np.random.seed(42) - - if len(requiredPoints) > numIntervals + 1: - - raise ValueError(f"Too many required points ({len(requiredPoints)}) for N.") - - num_of_remaining_points = numIntervals + 1 - len(requiredPoints) - m = 1 # number of points to oversample - - while True: - - cand = np.linspace(startPoint, endPoint, num_of_remaining_points + m + 2)[1:-1] # oversample to avoid overlaps - cand = np.round(cand, 0) - out = np.unique(np.r_[requiredPoints, cand]) - cand_without_required = out[~np.isin(out, requiredPoints)] - - if len(cand_without_required) >= num_of_remaining_points: - - picked_points = np.random.choice(cand_without_required, size=num_of_remaining_points, replace=False) - return np.sort(np.r_[requiredPoints, picked_points]) - - m *= 2 - - def plotSpeedLimits(track, pos_adj, v_adj): v_limits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) @@ -629,20 +603,73 @@ def plotCurvatures(track, pos_adj, c_adj, c_linear): plt.show() -def getRelevantEtcsBrakingPositions(track): +def getEquallyDistributedPoints(startPoint, endPoint, numMissingPoints, existingPointsInInterval): + + points = np.sort(np.unique(existingPointsInInterval)) + + assert startPoint < endPoint + assert np.isclose(startPoint, points[0]) + assert np.isclose(endPoint, points[-1]) + assert len(points) == len(existingPointsInInterval) + assert numMissingPoints >= 0 + + additionalPoints = [] + + for _ in range(numMissingPoints): + + gaps = np.diff(points) + gapIdx = np.argmax(gaps) - relevantPositions = [] + leftPoint = points[gapIdx] + rightPoint = points[gapIdx + 1] + + newPoint = np.round((leftPoint + rightPoint) / 2, 0) + + if not leftPoint < newPoint < rightPoint or np.any(np.isclose(newPoint, points)): + + raise ValueError("Could not insert another unique rounded point.") + + points = np.sort(np.append(points, newPoint)) + additionalPoints.append(newPoint) + + return additionalPoints + + +def introduceSufficientShootingNodesForETCSBrakingCurves(track, existingPoints): + + minPointsPer100Meters = 4 + + additionalPoints = [] + + startPoints, endPoints = [], [] positions = track.etcsPositions velocities = track.etcsVelocities - for idx in range(1, len(positions)): + for idx in range(1, len(positions)-1): + + if np.isclose(velocities[idx-1], velocities[idx]) and velocities[idx] > velocities[idx+1]: + + startPoints.append(positions[idx]) + + if np.isclose(velocities[idx], velocities[idx+1]) and velocities[idx-1] > velocities[idx]: + + endPoints.append(positions[idx]) + + for startPoint, endPoint in zip(startPoints, endPoints): + + numRequiredPoints = int(np.ceil(minPointsPer100Meters / 100 * (endPoint - startPoint))) + mask = (existingPoints >= startPoint) & (existingPoints <= endPoint) + existingPointsInInterval = existingPoints[mask] + existingPointsInInterval = np.unique(np.append(existingPointsInInterval, [startPoint, endPoint])) + additionalPoints.extend([startPoint, endPoint]) + numMissingPoints = max(numRequiredPoints - len(existingPointsInInterval), 0) - if velocities[idx - 2] > velocities[idx - 1] == velocities[idx]: + if numMissingPoints > 0: - relevantPositions.append(positions[idx]) + additionalPoints.extend(getEquallyDistributedPoints(startPoint, endPoint, numMissingPoints, existingPointsInInterval)) - return relevantPositions + return np.asarray(additionalPoints, dtype=float) def isSet( value): diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index d95a613..52b04f8 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -3,7 +3,6 @@ from mseetc.efficiency import totalLossesFunction from mseetc.etcs import getEtcsSpeedLimits -from mseetc.journey import Journey from mseetc.ocp import casadiSolver @@ -41,6 +40,22 @@ def printStats(df, stats, solver, train): print("Solver failed!") +def plotShootingNodes(ax, positions): + + for pos in positions: + + ax.vlines(pos / 1000, 0, 400, color='red', linewidth=0.1) + + +def plotOCPSpeedLimit(ax, solver): + + positions = solver.points.index.to_numpy() + speedLimits = solver.points["Speed limit [m/s]"].to_numpy() + + + ax.step(positions/1000, speedLimits*3.6, where="post", color="green", linestyle="-", linewidth=1, label="OCP Speed Limit") + + if __name__ == '__main__': from mseetc.train import Train @@ -61,7 +76,7 @@ def printStats(df, stats, solver, train): track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') # non-adjusted speed profile - opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} + opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':5}, 'energyOptimal':True} solver = casadiSolver(train, track, journey, opts) df, stats = solver.solve(journey) @@ -70,7 +85,7 @@ def printStats(df, stats, solver, train): # ETCS-adjusted speed profile track.setEtcsSpeedLimits(train) - opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} + opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':5}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} solverEtcs = casadiSolver(train, track, journey, opts) dfEtcs, statsEtcs = solverEtcs.solve(journey) @@ -80,7 +95,10 @@ def printStats(df, stats, solver, train): ### Plot Trajectory - fig, ax = plt.subplots(figsize=(16, 8)) + withShootingNodes = True + withOCPSpeedLimit = True + + fig, ax = plt.subplots(figsize=(24, 12)) x = track.speedLimits.index.to_numpy(dtype=float) v = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) @@ -95,12 +113,21 @@ def printStats(df, stats, solver, train): ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, linestyle="--", label="non-adjusted speed profile") ax.plot(dfEtcs["Position [m]"] / 1000, dfEtcs["Velocity [m/s]"] * 3.6, linestyle="--", label="ETCS-adjusted speed profile") + if withShootingNodes: + + plotShootingNodes(ax, dfEtcs["Position [m]"].to_numpy()) + + if withOCPSpeedLimit: + + plotOCPSpeedLimit(ax, solverEtcs) + ax.set_title("Speed Profile Comparison") ax.set_xlabel("Position [km]") ax.set_ylabel("Velocity [km/h]") ax.grid(True, which="both", linestyle="--", alpha=0.5) ax.legend(loc="upper right") ax.set_xlim(0, df["Position [m]"].max() / 1000) + ax.set_ylim(0, v_plot.max() * 3.6 * 1.3) ax.figure.tight_layout() plt.show() From 62f5e50d1ac62cfd51f7cb1523d350d36bb027f0 Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Tue, 16 Jun 2026 11:53:18 +0200 Subject: [PATCH 11/13] wip --- mseetc/etcs.py | 84 ++++++---- mseetc/track.py | 3 + .../fixtures/tracks/test_ETCS_flat_track.json | 41 +++++ .../tracks/test_ETCS_non_flat_track.json | 61 ++++++++ .../fixtures/tracks}/test_flat_no_tunnel.json | 0 .../tracks}/test_flat_with_tunnel.json | 0 .../fixtures/tracks}/test_one_hill.json | 0 .../tracks}/test_one_speed_increase.json | 0 .../fixtures/tracks}/test_two_radii.json | 0 .../fixtures/trains/CH_Stadler_FLIRT_TPF.json | 146 ++++++++++++++++++ .../testCurvatureResistance.py | 0 tests/units/etcs/testEtcs.py | 119 ++++++++++++++ .../units/integrators/testIntegrators.py | 0 .../testCurvature.py | 0 .../testGradient.py | 0 .../testSpeedLimit.py | 0 .../tunnelResistance/testTunnelResistance.py | 0 trains/CH_Stadler_FLIRT_TPF.json | 4 +- 18 files changed, 425 insertions(+), 33 deletions(-) create mode 100644 tests/fixtures/tracks/test_ETCS_flat_track.json create mode 100644 tests/fixtures/tracks/test_ETCS_non_flat_track.json rename {tracks => tests/fixtures/tracks}/test_flat_no_tunnel.json (100%) rename {tracks => tests/fixtures/tracks}/test_flat_with_tunnel.json (100%) rename {tracks => tests/fixtures/tracks}/test_one_hill.json (100%) rename {tracks => tests/fixtures/tracks}/test_one_speed_increase.json (100%) rename {tracks => tests/fixtures/tracks}/test_two_radii.json (100%) create mode 100644 tests/fixtures/trains/CH_Stadler_FLIRT_TPF.json rename unitTests/curvatureResistance/curvatureResistance.py => tests/units/curvatureResistance/testCurvatureResistance.py (100%) create mode 100644 tests/units/etcs/testEtcs.py rename unitTests/integrators/integrators.py => tests/units/integrators/testIntegrators.py (100%) rename unitTests/trainLengthDependentTrackAttributes/curvature.py => tests/units/trainLengthDependentTrackAttributes/testCurvature.py (100%) rename unitTests/trainLengthDependentTrackAttributes/gradient.py => tests/units/trainLengthDependentTrackAttributes/testGradient.py (100%) rename unitTests/trainLengthDependentTrackAttributes/speedLimit.py => tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py (100%) rename unitTests/tunnelResistance/tunnelResistance.py => tests/units/tunnelResistance/testTunnelResistance.py (100%) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index 13fe289..8257331 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -5,8 +5,6 @@ import pandas as pd from matplotlib import pyplot as plt -# from mseetc.track import Track - def getTrackVelocityAtPositions(speedLimitPositions, speedLimits, positions): """ @@ -52,16 +50,16 @@ def getBrakingTargetsFromSpeedLimits(track): def getEtcsSpeedLimits(train, track, positionStep=20.0): - calculator = EtcsBrakingCurveCalculator(train, track) - targets, speedLimitPositions, speedLimits = getBrakingTargetsFromSpeedLimits(track) + calculator = EtcsBrakingCurveCalculator(train, track) + # Compute P curves for all speed decreases pCurves = [] for target in targets: - curveSet = calculator.computeTarget(target) + curveSet, _ = calculator.computeTarget(target) curveSet["P"].loc[target.position, "Velocity [m/s]"] = target.targetVelocity pCurves.append(curveSet["P"]) @@ -161,9 +159,12 @@ def addEndPointToCurve(curve, velocity, end_position): def trimCurveToMaxVelocity(curve, maxVelocity): velocities = curve["Velocity [m/s]"].to_numpy(dtype=float) - keep_mask = velocities <= maxVelocity - return curve[keep_mask].copy() + keepMask = velocities <= maxVelocity + firstKeptIdx = np.where(keepMask)[0][0] + keepMask[max(firstKeptIdx - 1, 0)] = True + + return curve[keepMask].copy() def trimCurveFromMinVelocity(curve, minVelocity): @@ -279,27 +280,22 @@ def computeABrakeSafe(self): def computeAGradient(self, currentPosition): - positions = self.track.gradients.index.to_numpy(dtype=float) - gradients = self.track.gradients["Gradient [permil]"].to_numpy(dtype=float) - - if "Gradient linear term [permil/m]" in self.track.gradients.columns: - - gradientsLinearTerm = self.track.gradients["Gradient linear term [permil/m]"].to_numpy(dtype=float) - - else: - - gradientsLinearTerm = np.zeros(len(self.track.gradients)) - - - idx = bisect_right(positions, currentPosition) - 1 - idx = max(0, min(idx, len(positions) - 1)) + positions = self.track.gradientsTrainLengthIndependent.index.to_numpy(dtype=float) + gradients = self.track.gradientsTrainLengthIndependent["Gradient [permil]"].to_numpy(dtype=float) # If the backward-computed curve extends before the first known gradient point, assume flat track. if currentPosition < positions[0]: return 0 - gradient = gradients[idx] + gradientsLinearTerm[idx] * (currentPosition - positions[idx]) + idxFront = bisect_right(positions, currentPosition) - 1 + idxRear = bisect_right(positions, currentPosition - self.train.length) - 1 + + idxFront = max(0, min(idxFront, len(positions) - 1)) + idxRear = max(0, min(idxRear, len(positions) - 1)) + + gradient = np.min(gradients[idxRear:idxFront + 1]) + return 9.81 * gradient * 0.001 @@ -376,7 +372,7 @@ def computeEBICurve(self, EBD_curve, targetVelocity): for pos, vel in zip(positionsEBD, velocitiesEBD): - A_est1 = 0.1 # todo + A_est1 = 0.0 # todo A_est2 = min(A_est1, 0.4) V_est = (vel - A_est1 * T_traction - A_est2 * T_berem) / (1 + v_uncertainty) @@ -547,7 +543,9 @@ def computeTarget(self, target): curves = self.trimCurves(curves, target) - return curves + interventionPoints = self.computeInterventionPoints(curves, target) + + return curves, interventionPoints def plotCurves(self, curves, target): @@ -584,6 +582,26 @@ def plotCurves(self, curves, target): plt.show() + def computeInterventionPoints(self, curves, target): + + interventionPoints = {} + + for name, curve in curves.items(): + + curvePositions = curve.index.to_numpy(dtype=float) + curveVelocities = curve["Velocity [m/s]"].to_numpy(dtype=float) + interventionPoints[name] = target.position - np.interp(target.permittedVelocity, curveVelocities[::-1], curvePositions[::-1]) + + return interventionPoints + + + def printInterventionPoints(self, interventionPoints): + + for name, value in interventionPoints.items(): + + print(name, " point: ", round(value, 2)) + + if __name__ == '__main__': from mseetc.track import Track @@ -591,20 +609,23 @@ def plotCurves(self, curves, target): train = Train(config={'id': 'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') - track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') - track.updateTrainLengthDependentValues(train) + # track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') + track = Track(config={'id': 'test_ETCS_non_flat_track'}, pathJSON='../tests/fixtures//tracks') + # track.updateTrainLengthDependentValues(train) target = BrakingTarget( position=5000, overlap= 100, - permittedVelocity=130/3.6, - targetVelocity=100/3.6 + permittedVelocity=140/3.6, + targetVelocity=00/3.6 ) - addConstantVelocitySections = False + addConstantVelocitySections = True calculator = EtcsBrakingCurveCalculator(train, track, distancePre=5000, distancePost=1000) - curve_set = calculator.computeTarget(target) + curve_set, interventionPoints = calculator.computeTarget(target) + + calculator.printInterventionPoints(interventionPoints) if addConstantVelocitySections: @@ -613,4 +634,5 @@ def plotCurves(self, curves, target): if target.targetVelocity > 0: curve_set = calculator.processCurvesAfterTarget(curve_set, target) - calculator.plotCurves(curve_set, target) \ No newline at end of file + calculator.plotCurves(curve_set, target) + diff --git a/mseetc/track.py b/mseetc/track.py index a8094b1..d850089 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -399,6 +399,8 @@ def importGradientTuples(self, tuples, unitPosition='m', unitGradient='permil'): checkDataFrame(self.gradients, self.length) + self.gradientsTrainLengthIndependent = self.gradients.copy() + def importSpeedLimitTuples(self, tuples, unitPosition='m', unitVelocity='km/h'): @@ -669,6 +671,7 @@ def crop(dfIn): self.gradients = crop(self.gradients) self.curvatures = crop(self.curvatures) self.crossSections = crop(self.crossSections) + self.gradientsTrainLengthIndependent = crop(self.gradientsTrainLengthIndependent) def setEtcsSpeedLimits(self, train): diff --git a/tests/fixtures/tracks/test_ETCS_flat_track.json b/tests/fixtures/tracks/test_ETCS_flat_track.json new file mode 100644 index 0000000..896e0ba --- /dev/null +++ b/tests/fixtures/tracks/test_ETCS_flat_track.json @@ -0,0 +1,41 @@ +{ + "metadata": { + "id": "test_ETCS_flat_track", + "created by": "Roland Staerk", + "library version": "TTOBench v1.4", + "license": "BSD 2-Clause License" + }, + "altitude": { + "unit": "m", + "value": 0 + }, + "stops": { + "unit": "m", + "values": [ + 0.0, + 30000.0 + ] + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": [ + [ + 0.0, + 160 + ] + ] + }, + "ETCS braking data": { + "M_NVAVADH": { + "unit": "-", + "value": 0 + }, + "Kt_int": { + "unit": "-", + "value": 1 + } + } +} \ No newline at end of file diff --git a/tests/fixtures/tracks/test_ETCS_non_flat_track.json b/tests/fixtures/tracks/test_ETCS_non_flat_track.json new file mode 100644 index 0000000..ee5c051 --- /dev/null +++ b/tests/fixtures/tracks/test_ETCS_non_flat_track.json @@ -0,0 +1,61 @@ +{ + "metadata": { + "id": "test_ETCS_flat_track", + "created by": "Roland Staerk", + "library version": "TTOBench v1.4", + "license": "BSD 2-Clause License" + }, + "altitude": { + "unit": "m", + "value": 0 + }, + "stops": { + "unit": "m", + "values": [ + 0.0, + 30000.0 + ] + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": [ + [ + 0.0, + 160 + ] + ] + }, + "gradients": { + "units": { + "position": "m", + "slope": "permil" + }, + "values": [ + [ + 0.0, + 0.0 + ], + [ + 3500.0, + 20.0 + ], + [ + 4200.0, + -20.0 + ] + ] + }, + "ETCS braking data": { + "M_NVAVADH": { + "unit": "-", + "value": 0 + }, + "Kt_int": { + "unit": "-", + "value": 1 + } + } +} \ No newline at end of file diff --git a/tracks/test_flat_no_tunnel.json b/tests/fixtures/tracks/test_flat_no_tunnel.json similarity index 100% rename from tracks/test_flat_no_tunnel.json rename to tests/fixtures/tracks/test_flat_no_tunnel.json diff --git a/tracks/test_flat_with_tunnel.json b/tests/fixtures/tracks/test_flat_with_tunnel.json similarity index 100% rename from tracks/test_flat_with_tunnel.json rename to tests/fixtures/tracks/test_flat_with_tunnel.json diff --git a/tracks/test_one_hill.json b/tests/fixtures/tracks/test_one_hill.json similarity index 100% rename from tracks/test_one_hill.json rename to tests/fixtures/tracks/test_one_hill.json diff --git a/tracks/test_one_speed_increase.json b/tests/fixtures/tracks/test_one_speed_increase.json similarity index 100% rename from tracks/test_one_speed_increase.json rename to tests/fixtures/tracks/test_one_speed_increase.json diff --git a/tracks/test_two_radii.json b/tests/fixtures/tracks/test_two_radii.json similarity index 100% rename from tracks/test_two_radii.json rename to tests/fixtures/tracks/test_two_radii.json diff --git a/tests/fixtures/trains/CH_Stadler_FLIRT_TPF.json b/tests/fixtures/trains/CH_Stadler_FLIRT_TPF.json new file mode 100644 index 0000000..2ee668f --- /dev/null +++ b/tests/fixtures/trains/CH_Stadler_FLIRT_TPF.json @@ -0,0 +1,146 @@ +{ + "metadata": { + "id": "CH_Stadler_FLIRT_TPF", + "created by": "Dimitris Kouzoupis", + "library version": "TTOBench v1.4", + "license": "BSD 2-Clause License" + }, + "num seats": { + "unit": "-", + "value":167 + }, + "num coaches": { + "unit": "-", + "value": 4 + }, + "length": { + "unit": "m", + "value": 58.6 + }, + "mass": { + "unit": "kg", + "value": 122000.0 + }, + "rho": { + "unit": "%", + "value": 10 + }, + "max traction power": { + "unit": "kW", + "value": 2600.0 + }, + "max reg braking power": { + "unit": "kW", + "value": 2600.0 + }, + "max traction force": { + "unit": "kN", + "value": 200.0 + }, + "max reg braking force": { + "unit": "kN", + "value": 200.0 + }, + "max pn braking force": { + "unit": "kN", + "value": 1220.0 + }, + "max acceleration": { + "unit": "m/s^2", + "value": 1.1 + }, + "max deceleration": { + "unit": "m/s^2", + "value": 1.1 + }, + "max speed": { + "unit": "km/h", + "value": 160 + }, + "rolling resistance r0": { + "unit": "kN", + "value": 2.37888 + }, + "rolling resistance r1": { + "unit": "kN/(km/h)", + "value": 0.01698165 + }, + "rolling resistance r2": { + "unit": "kN/(km/h)^2", + "value": 0.00093264 + }, + "efficiency traction": { + "unit": "%", + "value": 90.0 + }, + "efficiency reg brake": { + "unit": "%", + "value": 90.0 + }, + "tunnel resistance": { + "units": { + "cross section": "m^2", + "coefficient": "kN/(km/h)^2" + }, + "values": [ + [ + 24, + 0.0011698 + ], + [ + 40, + 0.0005579 + ] + ] + }, + "ETCS braking data": { + "A_brake_emergency": { + "units": { + "velocity": "m/s", + "deceleration": "m/s^2" + }, + "values": [ + [0, 0.9], + [20, 0.85], + [40, 0.8], + [60, 0.75] + ] + }, + "A_brake_service": { + "units": { + "velocity": "m/s", + "deceleration": "m/s^2" + }, + "values": [ + [0, 0.55], + [20, 0.5], + [40, 0.45], + [60, 0.4] + ] + }, + "K_dry_rst": { + "unit": "-", + "value": 0.75 + }, + "K_wet_rst": { + "unit": "-", + "value": 0.9 + }, + "T_traction": { + "unit": "s", + "value": 1 + }, + "T_be": { + "unit": "s", + "value": 3 + }, + "v_uncertainty": { + "unit": "%", + "value": 2.98 + }, + "T_bs": { + "unit": "s", + "value": 3 + } + } +} diff --git a/unitTests/curvatureResistance/curvatureResistance.py b/tests/units/curvatureResistance/testCurvatureResistance.py similarity index 100% rename from unitTests/curvatureResistance/curvatureResistance.py rename to tests/units/curvatureResistance/testCurvatureResistance.py diff --git a/tests/units/etcs/testEtcs.py b/tests/units/etcs/testEtcs.py new file mode 100644 index 0000000..1173bfc --- /dev/null +++ b/tests/units/etcs/testEtcs.py @@ -0,0 +1,119 @@ +import unittest +from pathlib import Path + +from mseetc.etcs import BrakingTarget, EtcsBrakingCurveCalculator +from mseetc.track import Track +from mseetc.train import Train + + +class TestETCS(unittest.TestCase): + + def testInterventionPointsForFlatTrack(self): + ''' + Verify that the ETCS braking curve calculator computes the expected + intervention points for a flat test track. + + The test checks the main supervision limits I, P, W, SBI, and EBI + against known reference positions. A small absolute tolerance is used + because the values may vary slightly due to numerical computation. + ''' + + Path(__file__).resolve() + + train = Train(config={'id': 'CH_Stadler_FLIRT_TPF'}, pathJSON='tests/fixtures/trains') + + track = Track(config={'id': 'test_ETCS_flat_track'}, pathJSON='tests/fixtures/tracks') + track.updateTrainLengthDependentValues(train) + + target = BrakingTarget( + position=5000, + overlap=100, + permittedVelocity=140/3.6, + targetVelocity=0 + ) + + + calculator = EtcsBrakingCurveCalculator(train, track, distancePre=5000, distancePost=1000) + _, interventionPoints = calculator.computeTarget(target) + + tol = 1.6 + + expectedInterventionPoints = { + "I": 2098.10, + "P": 1748.10, + "W": 1670.33, + "SBI": 1592.55, + "EBI": 1397.00, + } + + for pointName, expectedValue in expectedInterventionPoints.items(): + + actualValue = interventionPoints[pointName] + + self.assertAlmostEqual( + actualValue, + expectedValue, + delta=tol, + msg=( + f"Intervention point {pointName} is not within the expected tolerance. " + f"Expected: {expectedValue:.2f} m, " + f"actual: {actualValue:.2f} m, " + f"allowed tolerance: ±{tol:.2f} m." + ) + ) + + + def testInterventionPointsForNonFlatTrack(self): + ''' + Verify that the ETCS braking curve calculator computes the expected + intervention points for a non-flat test track. + + The test checks the main supervision limits I, P, W, SBI, and EBI + against known reference positions. A small absolute tolerance is used + because the values may vary slightly due to numerical computation. + ''' + + Path(__file__).resolve() + + train = Train(config={'id': 'CH_Stadler_FLIRT_TPF'}, pathJSON='tests/fixtures/trains') + + track = Track(config={'id': 'test_ETCS_non_flat_track'}, pathJSON='tests/fixtures/tracks') + track.updateTrainLengthDependentValues(train) + + target = BrakingTarget( + position=5000, + overlap=100, + permittedVelocity=140/3.6, + targetVelocity=0 + ) + + + calculator = EtcsBrakingCurveCalculator(train, track, distancePre=5000, distancePost=1000) + _, interventionPoints = calculator.computeTarget(target) + + tol = 33 + + expectedInterventionPoints = { + "I": 2139.76, + "P": 1789.76, + "W": 1711.98, + "SBI": 1634.21, + "EBI": 1491.69, + } + + for pointName, expectedValue in expectedInterventionPoints.items(): + + actualValue = interventionPoints[pointName] + + self.assertAlmostEqual( + actualValue, + expectedValue, + delta=tol, + msg=( + f"Intervention point {pointName} is not within the expected tolerance. " + f"Expected: {expectedValue:.2f} m, " + f"actual: {actualValue:.2f} m, " + f"allowed tolerance: ±{tol:.2f} m." + ) + ) + diff --git a/unitTests/integrators/integrators.py b/tests/units/integrators/testIntegrators.py similarity index 100% rename from unitTests/integrators/integrators.py rename to tests/units/integrators/testIntegrators.py diff --git a/unitTests/trainLengthDependentTrackAttributes/curvature.py b/tests/units/trainLengthDependentTrackAttributes/testCurvature.py similarity index 100% rename from unitTests/trainLengthDependentTrackAttributes/curvature.py rename to tests/units/trainLengthDependentTrackAttributes/testCurvature.py diff --git a/unitTests/trainLengthDependentTrackAttributes/gradient.py b/tests/units/trainLengthDependentTrackAttributes/testGradient.py similarity index 100% rename from unitTests/trainLengthDependentTrackAttributes/gradient.py rename to tests/units/trainLengthDependentTrackAttributes/testGradient.py diff --git a/unitTests/trainLengthDependentTrackAttributes/speedLimit.py b/tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py similarity index 100% rename from unitTests/trainLengthDependentTrackAttributes/speedLimit.py rename to tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py diff --git a/unitTests/tunnelResistance/tunnelResistance.py b/tests/units/tunnelResistance/testTunnelResistance.py similarity index 100% rename from unitTests/tunnelResistance/tunnelResistance.py rename to tests/units/tunnelResistance/testTunnelResistance.py diff --git a/trains/CH_Stadler_FLIRT_TPF.json b/trains/CH_Stadler_FLIRT_TPF.json index bcd2c76..2ee668f 100644 --- a/trains/CH_Stadler_FLIRT_TPF.json +++ b/trains/CH_Stadler_FLIRT_TPF.json @@ -120,7 +120,7 @@ }, "K_dry_rst": { "unit": "-", - "value": 0.8 + "value": 0.75 }, "K_wet_rst": { "unit": "-", @@ -132,7 +132,7 @@ }, "T_be": { "unit": "s", - "value": 4 + "value": 3 }, "v_uncertainty": { "unit": "%", From 97accd1a8048e1d29d476d2bbb82cd4964049aad Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Tue, 16 Jun 2026 14:30:02 +0200 Subject: [PATCH 12/13] wip 2 --- mseetc/etcs.py | 5 ++--- mseetc/ocp.py | 8 ++++++-- simulations/sim_launcher.py | 6 ++---- tests/units/etcs/testEtcs.py | 26 +++++++++++++++++++++++--- 4 files changed, 33 insertions(+), 12 deletions(-) diff --git a/mseetc/etcs.py b/mseetc/etcs.py index 8257331..704a9bf 100644 --- a/mseetc/etcs.py +++ b/mseetc/etcs.py @@ -609,9 +609,8 @@ def printInterventionPoints(self, interventionPoints): train = Train(config={'id': 'CH_Stadler_FLIRT_TPF'}, pathJSON='../trains') - # track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') - track = Track(config={'id': 'test_ETCS_non_flat_track'}, pathJSON='../tests/fixtures//tracks') - # track.updateTrainLengthDependentValues(train) + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') + track.updateTrainLengthDependentValues(train) target = BrakingTarget( position=5000, diff --git a/mseetc/ocp.py b/mseetc/ocp.py index 71aa6ef..adbc36c 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -126,14 +126,18 @@ def __init__(self, train, track, journey, optsDict={}): # track parameters + if opts.withEtcsBrakingCurves: + + track.setEtcsSpeedLimits(train) + timingPointPositions = journey.timingPoints.index.to_numpy(dtype=float)[1:-1] - journey.positionStart self.points = computeDiscretizationPoints(track, numIntervals, opts, timingPointPositions) self.steps = np.diff(self.points.index) if opts.withEtcsBrakingCurves: - velocities = np.interp(self.points.index.to_numpy(), track.etcsPositions, track.etcsVelocities) - self.points["Speed limit [m/s]"] = np.maximum(velocities, velocityMin) + etcsVelocities = np.interp(self.points.index.to_numpy(), track.etcsPositions, track.etcsVelocities) + self.points["Speed limit [m/s]"] = np.maximum(etcsVelocities, velocityMin) # real-time parameters diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 52b04f8..7c28ab0 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -83,8 +83,8 @@ def plotOCPSpeedLimit(ax, solver): printStats(df, stats, solver, train) + # ETCS-adjusted speed profile - track.setEtcsSpeedLimits(train) opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':5}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} solverEtcs = casadiSolver(train, track, journey, opts) @@ -105,10 +105,8 @@ def plotOCPSpeedLimit(ax, solver): x_plot = np.append(x, track.length) v_plot = np.append(v, v[-1]) - etcsLimitsPositions, etcsLimitsVelocities = getEtcsSpeedLimits(train, track) - ax.step(x_plot/1000, v_plot*3.6, where="post", color="black", linestyle="-", label="Track Speed Limit") - ax.plot(etcsLimitsPositions/1000, etcsLimitsVelocities*3.6, color="red", linestyle="-", label="ETCS Speed Limit") + ax.plot(track.etcsPositions/1000, track.etcsVelocities*3.6, color="red", linestyle="-", label="ETCS Speed Limit") ax.plot(df["Position [m]"] / 1000, df["Velocity [m/s]"] * 3.6, linestyle="--", label="non-adjusted speed profile") ax.plot(dfEtcs["Position [m]"] / 1000, dfEtcs["Velocity [m/s]"] * 3.6, linestyle="--", label="ETCS-adjusted speed profile") diff --git a/tests/units/etcs/testEtcs.py b/tests/units/etcs/testEtcs.py index 1173bfc..d910a28 100644 --- a/tests/units/etcs/testEtcs.py +++ b/tests/units/etcs/testEtcs.py @@ -91,14 +91,13 @@ def testInterventionPointsForNonFlatTrack(self): calculator = EtcsBrakingCurveCalculator(train, track, distancePre=5000, distancePost=1000) _, interventionPoints = calculator.computeTarget(target) - tol = 33 + tol = 1 expectedInterventionPoints = { "I": 2139.76, "P": 1789.76, "W": 1711.98, - "SBI": 1634.21, - "EBI": 1491.69, + "SBI": 1634.21 } for pointName, expectedValue in expectedInterventionPoints.items(): @@ -117,3 +116,24 @@ def testInterventionPointsForNonFlatTrack(self): ) ) + tol = 8.5 + + expectedInterventionPoints = { + "EBI": 1491.69 + } + + for pointName, expectedValue in expectedInterventionPoints.items(): + actualValue = interventionPoints[pointName] + + self.assertAlmostEqual( + actualValue, + expectedValue, + delta=tol, + msg=( + f"Intervention point {pointName} is not within the expected tolerance. " + f"Expected: {expectedValue:.2f} m, " + f"actual: {actualValue:.2f} m, " + f"allowed tolerance: ±{tol:.2f} m." + ) + ) + From d71706b83e3ba19b9b4b55cc4b07a7984b518dcc Mon Sep 17 00:00:00 2001 From: Roli15 <81290196+Roli15@users.noreply.github.com> Date: Tue, 16 Jun 2026 16:37:26 +0200 Subject: [PATCH 13/13] good code --- .../00_var_speed_limit_100_Journey_01.json | 34 +- journeys/CH_StGallen_Wil_Journey_01.json | 34 +- journeys/CH_StGallen_Wil_Journey_02.json | 62 +- mseetc/journey.py | 25 +- mseetc/ocp.py | 11 +- mseetc/train.py | 2 +- simulations/sim_launcher.py | 4 +- simulations/sim_launcher_timingPoints.py | 2 +- .../journeys/CH_StGallen_Wil_Journey_01.json | 34 + .../journeys/CH_StGallen_Wil_Journey_02.json | 34 + .../test_flat_no_tunnel_Journey_01.json | 34 + .../test_flat_with_tunnel_Journey_01.json | 34 + .../journeys/test_one_hill_Journey_01.json | 34 + .../test_one_speed_increase_Journey_01.json | 34 + .../journeys/test_two_radii_Journey_01.json | 34 + tests/fixtures/tracks/CH_StGallen_Wil.json | 1908 +++++++++++++++++ tests/units/integrators/testIntegrators.py | 33 +- .../testCurvature.py | 49 +- .../testGradient.py | 51 +- .../testSpeedLimit.py | 20 +- .../tunnelResistance/testTunnelResistance.py | 27 +- 21 files changed, 2323 insertions(+), 177 deletions(-) create mode 100644 tests/fixtures/journeys/CH_StGallen_Wil_Journey_01.json create mode 100644 tests/fixtures/journeys/CH_StGallen_Wil_Journey_02.json create mode 100644 tests/fixtures/journeys/test_flat_no_tunnel_Journey_01.json create mode 100644 tests/fixtures/journeys/test_flat_with_tunnel_Journey_01.json create mode 100644 tests/fixtures/journeys/test_one_hill_Journey_01.json create mode 100644 tests/fixtures/journeys/test_one_speed_increase_Journey_01.json create mode 100644 tests/fixtures/journeys/test_two_radii_Journey_01.json create mode 100644 tests/fixtures/tracks/CH_StGallen_Wil.json diff --git a/journeys/00_var_speed_limit_100_Journey_01.json b/journeys/00_var_speed_limit_100_Journey_01.json index 8708964..aa4b1dd 100644 --- a/journeys/00_var_speed_limit_100_Journey_01.json +++ b/journeys/00_var_speed_limit_100_Journey_01.json @@ -3,10 +3,8 @@ "id": "00_var_speed_limit_100_Journey_01", "created by": "Roland Staerk", "library version": "TTOBench v1.5", - "license": "BSD 2-Clause License" - }, - "associated track id": { - "id": "00_var_speed_limit_100" + "license": "BSD 2-Clause License", + "associated track": "00_var_speed_limit_100" }, "timing points": { "units": { @@ -17,20 +15,20 @@ "upper speed constraint": "km/h" }, "values": [ - { - "position": 0, - "lower time constraint": 0, - "upper time constraint": 0, - "lower speed constraint": 0, - "upper speed constraint": 0 - }, - { - "position": 20000, - "lower time constraint": 685, - "upper time constraint": 685, - "lower speed constraint": 0, - "upper speed constraint": 0 - } + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 20000, + 685, + 685, + 0, + 0 + ] ] } } \ No newline at end of file diff --git a/journeys/CH_StGallen_Wil_Journey_01.json b/journeys/CH_StGallen_Wil_Journey_01.json index 006137d..3655d54 100644 --- a/journeys/CH_StGallen_Wil_Journey_01.json +++ b/journeys/CH_StGallen_Wil_Journey_01.json @@ -3,10 +3,8 @@ "id": "CH_StGallen_Wil_Journey_01", "created by": "Roland Staerk", "library version": "TTOBench v1.5", - "license": "BSD 2-Clause License" - }, - "associated track id": { - "id": "CH_StGallen_Wil" + "license": "BSD 2-Clause License", + "associated track": "CH_StGallen_Wil" }, "timing points": { "units": { @@ -17,20 +15,20 @@ "upper speed constraint": "km/h" }, "values": [ - { - "position": 0, - "lower time constraint": 0, - "upper time constraint": 0, - "lower speed constraint": 0, - "upper speed constraint": 0 - }, - { - "position": 20000, - "lower time constraint": 700, - "upper time constraint": 700, - "lower speed constraint": 0, - "upper speed constraint": 0 - } + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 20000, + 700, + 700, + 0, + 0 + ] ] } } \ No newline at end of file diff --git a/journeys/CH_StGallen_Wil_Journey_02.json b/journeys/CH_StGallen_Wil_Journey_02.json index a955c77..f382d7d 100644 --- a/journeys/CH_StGallen_Wil_Journey_02.json +++ b/journeys/CH_StGallen_Wil_Journey_02.json @@ -3,10 +3,8 @@ "id": "CH_StGallen_Wil_Journey_02", "created by": "Roland Staerk", "library version": "TTOBench v1.5", - "license": "BSD 2-Clause License" - }, - "associated track id": { - "id": "CH_StGallen_Wil" + "license": "BSD 2-Clause License", + "associated track": "CH_StGallen_Wil" }, "timing points": { "units": { @@ -17,34 +15,34 @@ "upper speed constraint": "km/h" }, "values": [ - { - "position": 1000, - "lower time constraint": 100, - "upper time constraint": 200, - "lower speed constraint": 0, - "upper speed constraint": 0 - }, - { - "position": 12500, - "lower time constraint": 480, - "upper time constraint": 540, - "lower speed constraint": 110, - "upper speed constraint": 120 - }, - { - "position": 18000, - "lower time constraint": 800, - "upper time constraint": null, - "lower speed constraint": null, - "upper speed constraint": 70 - }, - { - "position": 25000, - "lower time constraint": 900, - "upper time constraint": 1100, - "lower speed constraint": 0, - "upper speed constraint": 0 - } + [ + 1000, + 100, + 200, + 0, + 0 + ], + [ + 12500, + 480, + 540, + 110, + 120 + ], + [ + 18000, + 800, + null, + null, + 70 + ], + [ + 25000, + 900, + 1100, + 0, + 0 + ] ] } } \ No newline at end of file diff --git a/mseetc/journey.py b/mseetc/journey.py index 09cd631..7faf5b1 100644 --- a/mseetc/journey.py +++ b/mseetc/journey.py @@ -35,7 +35,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'journeys') - # read data self.id = data['metadata']['id'] - self.associatedTrackID = data['associated track id']['id'] + self.associatedTrackID = data['metadata']['associated track'] self.timingPoints = self.readTimingPoints(data['timing points']) @@ -45,7 +45,6 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'journeys') - self.computeInitialAndTerminalStates() - def readTimingPoints(self, timingPoints): units = timingPoints['units'] @@ -59,12 +58,13 @@ def readTimingPoints(self, timingPoints): } for point in timingPoints['values']: + pos, tMin, tMax, vMin, vMax = point - values["Position [m]"].append(convertUnit(point['position'], units['position'])) - values["Lower time constraint [s]"].append(self.convertConstraint(point['lower time constraint'], units['lower time constraint'])) - values["Upper time constraint [s]"].append(self.convertConstraint(point['upper time constraint'], units['upper time constraint'])) - values["Lower speed constraint [m/s]"].append(self.convertConstraint(point['lower speed constraint'], units['lower speed constraint'])) - values["Upper speed constraint [m/s]"].append(self.convertConstraint(point['upper speed constraint'], units['upper speed constraint'])) + values["Position [m]"].append(convertUnit(pos, units['position'])) + values["Lower time constraint [s]"].append(self.convertConstraint(tMin, units['lower time constraint'])) + values["Upper time constraint [s]"].append(self.convertConstraint(tMax, units['upper time constraint'])) + values["Lower speed constraint [m/s]"].append(self.convertConstraint(vMin, units['lower speed constraint'])) + values["Upper speed constraint [m/s]"].append(self.convertConstraint(vMax, units['upper speed constraint'])) df = pd.DataFrame(values) df = df.set_index("Position [m]") @@ -81,6 +81,11 @@ def convertConstraint(self, value, unit): return convertUnit(value, unit) + def isSet(self, value): + + return value is not None and not pd.isna(value) + + def checkFields(self): if len(self.timingPoints) < 2: @@ -115,15 +120,15 @@ def checkFields(self): vMin = point["Lower speed constraint [m/s]"] vMax = point["Upper speed constraint [m/s]"] - if tMin is not None and tMax is not None and tMin > tMax: + if self.isSet(tMin) and self.isSet(tMax) and tMin > tMax: raise ValueError("Lower time constraint must be smaller than or equal to upper time constraint!") - if vMin is not None and vMax is not None and vMin > vMax: + if self.isSet(vMin) and self.isSet(vMax) and vMin > vMax: raise ValueError("Lower speed constraint must be smaller than or equal to upper speed constraint!") - if any(value is not None and (value < 0 or np.isinf(value)) for value in [tMin, tMax, vMin, vMax]): + if any(self.isSet(value) and (value < 0 or np.isinf(value)) for value in [tMin, tMax, vMin, vMax]): raise ValueError("Timing point constraints must be positive finite numbers or None!") diff --git a/mseetc/ocp.py b/mseetc/ocp.py index adbc36c..9b1b253 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -356,6 +356,7 @@ def __init__(self, train, track, journey, optsDict={}): self.withRgBrake = withRgBrake self.withPnBrake = withPnBrake self.train = train + self.journey = journey self.energyOptimal = opts.energyOptimal self.scalingFactorObjective = scalingFactorObjective self.opts = opts @@ -366,12 +367,12 @@ def __init__(self, train, track, journey, optsDict={}): self.ubg = ca.vcat(ubg) - def solve(self, journey): + def solve(self): - initialTime = journey.initialTime - terminalTime = journey.terminalTime - initialVelocity = journey.initialVelocity - terminalVelocity = journey.terminalVelocity + initialTime = self.journey.initialTime + terminalTime = self.journey.terminalTime + initialVelocity = self.journey.initialVelocity + terminalVelocity = self.journey.terminalVelocity # check boundary conditions on time diff --git a/mseetc/train.py b/mseetc/train.py index 6a77652..b4eb971 100644 --- a/mseetc/train.py +++ b/mseetc/train.py @@ -70,7 +70,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> self.mass = convertUnit(data['mass']['value'], data['mass']['unit']) # train mass [kg] - self.rho = convertUnit(data['rho']['value'], data['rho']['unit']) # rotating-mass factor [-] + self.rho = convertUnit(data['rho']['value'], data['rho']['unit']) if 'max traction force' in data else 1 # rotating-mass factor [-] if self.rho < 1: diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py index 7c28ab0..bcbc44e 100644 --- a/simulations/sim_launcher.py +++ b/simulations/sim_launcher.py @@ -79,7 +79,7 @@ def plotOCPSpeedLimit(ax, solver): opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':5}, 'energyOptimal':True} solver = casadiSolver(train, track, journey, opts) - df, stats = solver.solve(journey) + df, stats = solver.solve() printStats(df, stats, solver, train) @@ -88,7 +88,7 @@ def plotOCPSpeedLimit(ax, solver): opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':5}, 'energyOptimal':True, 'withEtcsBrakingCurves': True} solverEtcs = casadiSolver(train, track, journey, opts) - dfEtcs, statsEtcs = solverEtcs.solve(journey) + dfEtcs, statsEtcs = solverEtcs.solve() printStats(dfEtcs, statsEtcs, solverEtcs, train) diff --git a/simulations/sim_launcher_timingPoints.py b/simulations/sim_launcher_timingPoints.py index 53f83bd..8cbe842 100644 --- a/simulations/sim_launcher_timingPoints.py +++ b/simulations/sim_launcher_timingPoints.py @@ -110,7 +110,7 @@ def plotPositionVelocityProfile(df, journey): opts = {'numIntervals':800, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':1}, 'energyOptimal':True} solver = casadiSolver(train, track, journey, opts) - df, stats = solver.solve(journey) + df, stats = solver.solve() printStats(df, stats, solver, train) diff --git a/tests/fixtures/journeys/CH_StGallen_Wil_Journey_01.json b/tests/fixtures/journeys/CH_StGallen_Wil_Journey_01.json new file mode 100644 index 0000000..a1d94d4 --- /dev/null +++ b/tests/fixtures/journeys/CH_StGallen_Wil_Journey_01.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "CH_StGallen_Wil_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "CH_StGallen_Wil" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 5000, + 300, + 300, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/journeys/CH_StGallen_Wil_Journey_02.json b/tests/fixtures/journeys/CH_StGallen_Wil_Journey_02.json new file mode 100644 index 0000000..9d07f49 --- /dev/null +++ b/tests/fixtures/journeys/CH_StGallen_Wil_Journey_02.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "CH_StGallen_Wil_Journey_02", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "CH_StGallen_Wil" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 29500, + 1327.5, + 1327.5, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/journeys/test_flat_no_tunnel_Journey_01.json b/tests/fixtures/journeys/test_flat_no_tunnel_Journey_01.json new file mode 100644 index 0000000..6d016c1 --- /dev/null +++ b/tests/fixtures/journeys/test_flat_no_tunnel_Journey_01.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "test_flat_no_tunnel_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "test_flat_no_tunnel" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 28000, + 695, + 695, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/journeys/test_flat_with_tunnel_Journey_01.json b/tests/fixtures/journeys/test_flat_with_tunnel_Journey_01.json new file mode 100644 index 0000000..0396a3e --- /dev/null +++ b/tests/fixtures/journeys/test_flat_with_tunnel_Journey_01.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "test_flat_with_tunnel_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "test_flat_with_tunnel" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 28000, + 695, + 695, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/journeys/test_one_hill_Journey_01.json b/tests/fixtures/journeys/test_one_hill_Journey_01.json new file mode 100644 index 0000000..55f973c --- /dev/null +++ b/tests/fixtures/journeys/test_one_hill_Journey_01.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "test_one_hill_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "test_one_hill" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 5000, + 300, + 300, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/journeys/test_one_speed_increase_Journey_01.json b/tests/fixtures/journeys/test_one_speed_increase_Journey_01.json new file mode 100644 index 0000000..178df59 --- /dev/null +++ b/tests/fixtures/journeys/test_one_speed_increase_Journey_01.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "test_one_speed_increase_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "test_one_speed_increase" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 12000, + 375, + 375, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/journeys/test_two_radii_Journey_01.json b/tests/fixtures/journeys/test_two_radii_Journey_01.json new file mode 100644 index 0000000..973d5ae --- /dev/null +++ b/tests/fixtures/journeys/test_two_radii_Journey_01.json @@ -0,0 +1,34 @@ +{ + "metadata": { + "id": "test_two_radii_Journey_01", + "created by": "Roland Staerk", + "library version": "TTOBench v1.5", + "license": "BSD 2-Clause License", + "associated track": "test_two_radii" + }, + "timing points": { + "units": { + "position": "m", + "lower time constraint": "s", + "upper time constraint": "s", + "lower speed constraint": "km/h", + "upper speed constraint": "km/h" + }, + "values": [ + [ + 0, + 0, + 0, + 0, + 0 + ], + [ + 5000, + 300, + 300, + 0, + 0 + ] + ] + } +} \ No newline at end of file diff --git a/tests/fixtures/tracks/CH_StGallen_Wil.json b/tests/fixtures/tracks/CH_StGallen_Wil.json new file mode 100644 index 0000000..fe0915f --- /dev/null +++ b/tests/fixtures/tracks/CH_StGallen_Wil.json @@ -0,0 +1,1908 @@ +{ + "metadata": { + "id": "CH_StGallen_Wil", + "created by": "Dimitris Kouzoupis, Juxhino Kavaja", + "library version": "TTOBench v1.2", + "license": "BSD 2-Clause License" + }, + "altitude": { + "unit": "m", + "value": 675.2 + }, + "stops": { + "unit": "m", + "values": [ + 0.0, + 29556.1 + ] + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": [ + [ + 0.0, + 90 + ], + [ + 49.6, + 100 + ], + [ + 125.6, + 110 + ], + [ + 504.3, + 125 + ], + [ + 2140.4, + 110 + ], + [ + 3475.8, + 120 + ], + [ + 4824.1, + 125 + ], + [ + 13497.5, + 115 + ], + [ + 14872.2, + 105 + ], + [ + 18679.4, + 90 + ], + [ + 20761.3, + 105 + ], + [ + 28178.7, + 90 + ], + [ + 28456.1, + 80 + ] + ] + }, + "gradients": { + "units": { + "position": "m", + "slope": "permil" + }, + "values": [ + [ + 0.0, + 11.9 + ], + [ + 145.1, + 2.0 + ], + [ + 239.5, + -5.4 + ], + [ + 296.7, + -7.7 + ], + [ + 385.6, + -9.9 + ], + [ + 465.4, + -11.2 + ], + [ + 775.2, + -11.4 + ], + [ + 909.0, + -10.0 + ], + [ + 1047.5, + -6.2 + ], + [ + 1129.2, + -9.8 + ], + [ + 1235.4, + -7.8 + ], + [ + 1260.3, + -9.8 + ], + [ + 1475.1, + -10.1 + ], + [ + 1917.3, + -10.2 + ], + [ + 2335.6, + -9.8 + ], + [ + 2502.4, + -12.4 + ], + [ + 2523.9, + -9.7 + ], + [ + 2821.4, + -0.7 + ], + [ + 3037.6, + -5.8 + ], + [ + 3205.4, + -4.4 + ], + [ + 3361.1, + 0.0 + ], + [ + 3700.5, + 7.2 + ], + [ + 3875.7, + 5.6 + ], + [ + 3965.9, + 8.2 + ], + [ + 4212.0, + 0.8 + ], + [ + 4481.8, + -5.6 + ], + [ + 4520.9, + 0.0 + ], + [ + 4788.0, + 5.2 + ], + [ + 4871.8, + 5.0 + ], + [ + 4933.5, + 5.8 + ], + [ + 5057.8, + 0.0 + ], + [ + 5204.6, + 0.1 + ], + [ + 5457.0, + -5.3 + ], + [ + 5638.1, + -6.9 + ], + [ + 5778.3, + -1.4 + ], + [ + 5932.2, + -4.4 + ], + [ + 6146.4, + -5.4 + ], + [ + 6837.9, + -4.8 + ], + [ + 7474.0, + -9.6 + ], + [ + 7629.9, + -12.5 + ], + [ + 7791.9, + -4.2 + ], + [ + 7905.0, + -9.6 + ], + [ + 7954.8, + -10.9 + ], + [ + 8052.6, + -9.4 + ], + [ + 8152.9, + -3.9 + ], + [ + 8221.4, + -0.6 + ], + [ + 8297.8, + -1.8 + ], + [ + 8360.2, + -0.5 + ], + [ + 8508.1, + -0.7 + ], + [ + 8533.9, + 0.0 + ], + [ + 8658.3, + -0.4 + ], + [ + 8737.2, + 0.9 + ], + [ + 8759.6, + -0.2 + ], + [ + 8856.0, + -0.6 + ], + [ + 9037.7, + -5.3 + ], + [ + 9130.1, + -0.4 + ], + [ + 9280.5, + -3.7 + ], + [ + 9528.3, + -3.0 + ], + [ + 9749.0, + -3.3 + ], + [ + 10483.2, + -9.4 + ], + [ + 11085.7, + -6.4 + ], + [ + 11386.9, + -2.6 + ], + [ + 11540.4, + -3.2 + ], + [ + 11889.0, + -5.0 + ], + [ + 12385.5, + -10.0 + ], + [ + 12669.8, + -8.0 + ], + [ + 13191.6, + -15.4 + ], + [ + 13313.9, + -8.8 + ], + [ + 13521.9, + 0.0 + ], + [ + 13691.5, + -0.2 + ], + [ + 13765.4, + 0.2 + ], + [ + 13900.2, + 0.1 + ], + [ + 13966.9, + -0.3 + ], + [ + 14066.3, + 0.2 + ], + [ + 14213.9, + 3.0 + ], + [ + 14332.7, + 2.7 + ], + [ + 14437.1, + 3.6 + ], + [ + 14523.5, + 0.2 + ], + [ + 14603.7, + -1.5 + ], + [ + 14642.0, + -1.2 + ], + [ + 14720.6, + -9.2 + ], + [ + 14809.4, + -6.6 + ], + [ + 14842.6, + -5.6 + ], + [ + 14935.4, + -13.7 + ], + [ + 14992.1, + -9.7 + ], + [ + 15193.8, + -9.2 + ], + [ + 15608.1, + -10.2 + ], + [ + 16230.1, + -8.7 + ], + [ + 16391.4, + -10.2 + ], + [ + 16592.0, + -10.7 + ], + [ + 16799.9, + -10.0 + ], + [ + 18406.3, + -9.8 + ], + [ + 18546.4, + -8.0 + ], + [ + 18661.2, + -14.0 + ], + [ + 18743.1, + -8.6 + ], + [ + 18795.1, + -10.6 + ], + [ + 18835.5, + -8.9 + ], + [ + 18877.0, + -10.0 + ], + [ + 18965.9, + -11.2 + ], + [ + 19271.9, + -9.5 + ], + [ + 19521.0, + 1.0 + ], + [ + 19638.1, + -0.5 + ], + [ + 19715.6, + 3.5 + ], + [ + 19798.2, + -1.4 + ], + [ + 19932.4, + 3.3 + ], + [ + 19986.9, + -3.0 + ], + [ + 20124.2, + -2.9 + ], + [ + 20245.1, + -11.5 + ], + [ + 20295.5, + -10.6 + ], + [ + 20554.6, + -9.4 + ], + [ + 20579.0, + -10.2 + ], + [ + 20645.3, + -8.7 + ], + [ + 20714.6, + -9.9 + ], + [ + 20929.5, + -1.4 + ], + [ + 21033.8, + -2.3 + ], + [ + 21086.5, + -1.2 + ], + [ + 21272.6, + -1.3 + ], + [ + 21827.0, + -1.7 + ], + [ + 21960.7, + -4.2 + ], + [ + 22485.2, + -5.9 + ], + [ + 22634.7, + -1.2 + ], + [ + 22840.1, + -2.4 + ], + [ + 23491.7, + -8.0 + ], + [ + 24109.4, + -8.8 + ], + [ + 24276.9, + -7.7 + ], + [ + 24681.4, + -2.4 + ], + [ + 24923.5, + -1.4 + ], + [ + 25012.2, + -2.9 + ], + [ + 25030.1, + -1.7 + ], + [ + 25115.8, + -1.6 + ], + [ + 25302.6, + 0.0 + ], + [ + 25625.1, + 1.8 + ], + [ + 25777.4, + 0.0 + ], + [ + 25956.4, + 10.0 + ], + [ + 26498.7, + 9.4 + ], + [ + 26814.7, + 6.2 + ], + [ + 26906.5, + 15.9 + ], + [ + 26944.2, + 11.2 + ], + [ + 27114.0, + 9.8 + ], + [ + 27775.1, + 9.4 + ], + [ + 28102.0, + 11.0 + ], + [ + 28228.9, + 9.9 + ], + [ + 28315.2, + 10.0 + ], + [ + 28600.4, + 11.1 + ], + [ + 28682.5, + 15.0 + ], + [ + 28789.9, + 10.7 + ], + [ + 28952.0, + 1.9 + ], + [ + 29073.3, + -0.6 + ], + [ + 29219.7, + 0.0 + ], + [ + 29318.1, + -1.1 + ], + [ + 29381.8, + 0.0 + ], + [ + 29480.5, + 5.1 + ], + [ + 29535.6, + -4.6 + ] + ] + }, + "curvatures": { + "units": { + "position": "m", + "radius at start": "m", + "radius at end": "m" + }, + "values": [ + [ + 0.0, + 502.0, + 502.0 + ], + [ + 49.6, + 502.0, + 3570.0 + ], + [ + 125.6, + 3570.0, + 3570.0 + ], + [ + 172.5, + 3570.0, + 1250.0 + ], + [ + 198.5, + 1250.0, + 1250.0 + ], + [ + 232.1, + 1250.0, + "infinity" + ], + [ + 287.1, + "infinity", + "infinity" + ], + [ + 330.2, + -5700.0, + -5700.0 + ], + [ + 393.1, + "infinity", + "infinity" + ], + [ + 445.4, + "infinity", + 1567.0 + ], + [ + 594.4, + 1567.0, + 1567.0 + ], + [ + 948.8, + 1567.0, + "infinity" + ], + [ + 1018.8, + "infinity", + -850.0 + ], + [ + 1106.1, + -850.0, + -850.0 + ], + [ + 1234.7, + -850.0, + -2600.0 + ], + [ + 1314.7, + -2600.0, + -2600.0 + ], + [ + 1425.1, + -2600.0, + "infinity" + ], + [ + 1505.1, + "infinity", + "infinity" + ], + [ + 2140.4, + "infinity", + 508.0 + ], + [ + 2235.8, + 508.0, + 508.0 + ], + [ + 2326.0, + 508.0, + "infinity" + ], + [ + 2422.0, + "infinity", + "infinity" + ], + [ + 2552.6, + "infinity", + -752.0 + ], + [ + 2627.6, + -752.0, + -752.0 + ], + [ + 2782.3, + -752.0, + -605.0 + ], + [ + 2812.3, + -605.0, + -605.0 + ], + [ + 2953.0, + -605.0, + "infinity" + ], + [ + 3043.0, + "infinity", + "infinity" + ], + [ + 3103.0, + "infinity", + 520.2 + ], + [ + 3195.0, + 520.2, + 520.2 + ], + [ + 3380.8, + 520.2, + "infinity" + ], + [ + 3475.8, + "infinity", + "infinity" + ], + [ + 3557.0, + "infinity", + 2896.2 + ], + [ + 3597.0, + 2896.2, + 2896.2 + ], + [ + 3705.7, + 2896.2, + 6246.2 + ], + [ + 3755.7, + 6246.2, + 6246.2 + ], + [ + 3790.6, + 6246.2, + 1020.0 + ], + [ + 3884.6, + 1020.0, + 1020.0 + ], + [ + 3982.6, + 1020.0, + "infinity" + ], + [ + 4058.6, + "infinity", + "infinity" + ], + [ + 4499.2, + "infinity", + 657.2 + ], + [ + 4609.2, + 657.2, + 657.2 + ], + [ + 4717.1, + 657.2, + "infinity" + ], + [ + 4824.1, + "infinity", + "infinity" + ], + [ + 4926.7, + "infinity", + 5000.0 + ], + [ + 4943.7, + 5000.0, + 5000.0 + ], + [ + 4977.4, + 5000.0, + "infinity" + ], + [ + 4994.4, + "infinity", + -5000.0 + ], + [ + 5011.4, + -5000.0, + -5000.0 + ], + [ + 5045.0, + -5000.0, + "infinity" + ], + [ + 5062.0, + "infinity", + "infinity" + ], + [ + 5486.3, + "infinity", + -6000.0 + ], + [ + 5514.3, + -6000.0, + -6000.0 + ], + [ + 5537.9, + -6000.0, + "infinity" + ], + [ + 5565.9, + "infinity", + 6000.0 + ], + [ + 5593.9, + 6000.0, + 6000.0 + ], + [ + 5617.5, + 6000.0, + "infinity" + ], + [ + 5645.5, + "infinity", + "infinity" + ], + [ + 7385.4, + "infinity", + -3003.8 + ], + [ + 7425.5, + -3003.8, + -3003.8 + ], + [ + 7496.1, + -3003.8, + "infinity" + ], + [ + 7536.2, + "infinity", + "infinity" + ], + [ + 7793.9, + "infinity", + -962.0 + ], + [ + 7894.7, + -962.0, + -962.0 + ], + [ + 8105.6, + -962.0, + "infinity" + ], + [ + 8194.0, + "infinity", + "infinity" + ], + [ + 8464.2, + "infinity", + 2800.0 + ], + [ + 8507.1, + 2800.0, + 2800.0 + ], + [ + 8582.9, + 2800.0, + "infinity" + ], + [ + 8625.8, + "infinity", + -2780.0 + ], + [ + 8657.8, + -2780.0, + -2780.0 + ], + [ + 8737.2, + -2780.0, + "infinity" + ], + [ + 8784.2, + "infinity", + "infinity" + ], + [ + 8974.8, + "infinity", + -1350.0 + ], + [ + 9034.8, + -1350.0, + -1350.0 + ], + [ + 9079.3, + -1350.0, + -6130.8 + ], + [ + 9135.4, + -6130.8, + "infinity" + ], + [ + 9151.3, + "infinity", + "infinity" + ], + [ + 9212.0, + "infinity", + -803.6 + ], + [ + 9310.0, + -803.6, + -803.6 + ], + [ + 9341.8, + -803.6, + "infinity" + ], + [ + 9439.8, + "infinity", + "infinity" + ], + [ + 9629.7, + "infinity", + -778.0 + ], + [ + 9734.7, + -778.0, + -778.0 + ], + [ + 9943.2, + -778.0, + "infinity" + ], + [ + 10048.2, + "infinity", + "infinity" + ], + [ + 10336.8, + "infinity", + 700.0 + ], + [ + 10476.8, + 700.0, + 700.0 + ], + [ + 10986.5, + 700.0, + "infinity" + ], + [ + 11126.5, + "infinity", + "infinity" + ], + [ + 11501.9, + "infinity", + -6800.0 + ], + [ + 11531.9, + -6800.0, + -6800.0 + ], + [ + 11618.6, + -6800.0, + "infinity" + ], + [ + 11648.6, + "infinity", + "infinity" + ], + [ + 11714.3, + "infinity", + -6600.0 + ], + [ + 11744.3, + -6600.0, + -6600.0 + ], + [ + 11838.2, + -6600.0, + "infinity" + ], + [ + 11868.2, + "infinity", + "infinity" + ], + [ + 12241.3, + "infinity", + 736.3 + ], + [ + 12350.3, + 736.3, + 736.3 + ], + [ + 12410.9, + 736.3, + "infinity" + ], + [ + 12519.9, + "infinity", + "infinity" + ], + [ + 13645.7, + "infinity", + 585.0 + ], + [ + 13751.7, + 585.0, + 585.0 + ], + [ + 13863.8, + 585.0, + "infinity" + ], + [ + 13969.8, + "infinity", + "infinity" + ], + [ + 14000.9, + "infinity", + -5000.0 + ], + [ + 14020.9, + -5000.0, + -5000.0 + ], + [ + 14111.9, + -5000.0, + "infinity" + ], + [ + 14141.9, + "infinity", + 5000.0 + ], + [ + 14171.9, + 5000.0, + 5000.0 + ], + [ + 14231.1, + 5000.0, + "infinity" + ], + [ + 14251.1, + "infinity", + "infinity" + ], + [ + 14351.7, + "infinity", + -3000.0 + ], + [ + 14381.7, + -3000.0, + -3000.0 + ], + [ + 14456.5, + -3000.0, + "infinity" + ], + [ + 14486.5, + "infinity", + 3000.0 + ], + [ + 14516.5, + 3000.0, + 3000.0 + ], + [ + 14591.6, + 3000.0, + "infinity" + ], + [ + 14621.6, + "infinity", + "infinity" + ], + [ + 14796.0, + "infinity", + -610.4 + ], + [ + 14906.0, + -610.4, + -610.4 + ], + [ + 14935.8, + -634.8, + -634.8 + ], + [ + 15252.8, + -634.8, + -583.8 + ], + [ + 15283.0, + -583.8, + -583.8 + ], + [ + 15424.4, + -583.8, + "infinity" + ], + [ + 15506.7, + "infinity", + "infinity" + ], + [ + 15625.2, + "infinity", + 488.2 + ], + [ + 15724.8, + 488.2, + 488.2 + ], + [ + 15812.8, + 488.2, + 498.1 + ], + [ + 15837.6, + 498.1, + 498.1 + ], + [ + 16256.9, + 498.1, + "infinity" + ], + [ + 16346.6, + "infinity", + "infinity" + ], + [ + 16521.0, + "infinity", + -588.8 + ], + [ + 16603.3, + -588.8, + -588.8 + ], + [ + 16799.4, + -588.8, + -610.8 + ], + [ + 16824.6, + -610.8, + -610.8 + ], + [ + 17018.6, + -610.8, + "infinity" + ], + [ + 17100.8, + "infinity", + "infinity" + ], + [ + 17718.9, + "infinity", + -533.8 + ], + [ + 17811.2, + -533.8, + -533.8 + ], + [ + 17860.7, + -533.8, + "infinity" + ], + [ + 17957.8, + "infinity", + "infinity" + ], + [ + 18042.0, + "infinity", + 526.2 + ], + [ + 18133.6, + 526.2, + 526.2 + ], + [ + 18200.7, + 526.2, + "infinity" + ], + [ + 18298.3, + "infinity", + "infinity" + ], + [ + 18438.1, + "infinity", + -591.8 + ], + [ + 18534.4, + -591.8, + -591.8 + ], + [ + 18658.1, + -591.8, + "infinity" + ], + [ + 18738.3, + "infinity", + "infinity" + ], + [ + 18787.4, + "infinity", + 347.8 + ], + [ + 18864.4, + 347.8, + 347.8 + ], + [ + 19201.7, + 350.0, + 350.0 + ], + [ + 19223.1, + 343.8, + 343.8 + ], + [ + 19355.8, + 343.8, + 400.0 + ], + [ + 19385.8, + 400.0, + 400.0 + ], + [ + 19477.6, + 400.0, + "infinity" + ], + [ + 19549.6, + "infinity", + "infinity" + ], + [ + 19655.4, + "infinity", + -403.1 + ], + [ + 19752.1, + -403.1, + -403.1 + ], + [ + 19845.4, + -388.6, + -388.6 + ], + [ + 19870.3, + -399.3, + -399.3 + ], + [ + 19974.2, + -399.3, + -1900.0 + ], + [ + 20056.0, + -1900.0, + -1900.0 + ], + [ + 20093.6, + -1900.0, + -700.0 + ], + [ + 20138.6, + -700.0, + -700.0 + ], + [ + 20183.2, + -700.0, + -393.8 + ], + [ + 20228.2, + -393.8, + -393.8 + ], + [ + 20473.0, + -393.8, + "infinity" + ], + [ + 20563.4, + "infinity", + "infinity" + ], + [ + 20618.0, + "infinity", + 650.0 + ], + [ + 20717.8, + 650.0, + 650.0 + ], + [ + 20915.1, + 650.0, + "infinity" + ], + [ + 21014.8, + "infinity", + "infinity" + ], + [ + 21375.7, + "infinity", + 864.2 + ], + [ + 21460.5, + 864.2, + 864.2 + ], + [ + 21617.7, + 864.2, + "infinity" + ], + [ + 21702.6, + "infinity", + "infinity" + ], + [ + 21801.7, + "infinity", + -1103.8 + ], + [ + 21881.8, + -1103.8, + -1103.8 + ], + [ + 22036.8, + -1103.8, + "infinity" + ], + [ + 22116.9, + "infinity", + "infinity" + ], + [ + 22496.4, + "infinity", + -645.8 + ], + [ + 22586.7, + -645.8, + -645.8 + ], + [ + 22814.7, + -645.8, + "infinity" + ], + [ + 22904.9, + "infinity", + "infinity" + ], + [ + 23114.4, + "infinity", + -538.8 + ], + [ + 23224.8, + -538.8, + -538.8 + ], + [ + 23403.5, + -538.8, + "infinity" + ], + [ + 23513.9, + "infinity", + "infinity" + ], + [ + 23932.4, + "infinity", + 576.2 + ], + [ + 24029.0, + 576.2, + 576.2 + ], + [ + 24120.7, + 576.2, + "infinity" + ], + [ + 24217.4, + "infinity", + "infinity" + ], + [ + 24287.7, + "infinity", + -585.8 + ], + [ + 24384.0, + -585.8, + -585.8 + ], + [ + 24606.4, + -585.8, + "infinity" + ], + [ + 24702.7, + "infinity", + "infinity" + ], + [ + 24776.0, + "infinity", + 8400.0 + ], + [ + 24816.0, + 8400.0, + 8400.0 + ], + [ + 24879.3, + 8400.0, + "infinity" + ], + [ + 24919.3, + "infinity", + "infinity" + ], + [ + 25179.4, + "infinity", + 525.0 + ], + [ + 25269.4, + 525.0, + 525.0 + ], + [ + 25703.3, + 525.0, + "infinity" + ], + [ + 25793.3, + "infinity", + "infinity" + ], + [ + 26242.1, + "infinity", + -491.0 + ], + [ + 26336.1, + -491.0, + -491.0 + ], + [ + 26588.8, + -491.0, + "infinity" + ], + [ + 26688.8, + "infinity", + "infinity" + ], + [ + 26752.4, + "infinity", + 457.2 + ], + [ + 26852.4, + 457.2, + 457.2 + ], + [ + 27208.3, + 457.2, + 461.0 + ], + [ + 27232.3, + 461.0, + 461.0 + ], + [ + 27546.7, + 461.0, + 453.0 + ], + [ + 27570.7, + 453.0, + 453.0 + ], + [ + 27731.2, + 453.0, + "infinity" + ], + [ + 27827.2, + "infinity", + "infinity" + ], + [ + 28456.1, + "infinity", + -372.1 + ], + [ + 28532.1, + -372.1, + -372.1 + ], + [ + 28638.4, + -372.1, + -397.8 + ], + [ + 28678.4, + -397.8, + -397.8 + ], + [ + 28748.8, + -397.8, + -340.1 + ], + [ + 28788.8, + -340.1, + -340.1 + ], + [ + 28879.4, + -340.1, + -520.0 + ], + [ + 28919.4, + -520.0, + -520.0 + ], + [ + 28948.8, + -520.0, + -351.0 + ], + [ + 28988.8, + -351.0, + -351.0 + ], + [ + 29090.0, + -351.0, + "infinity" + ], + [ + 29160.0, + "infinity", + "infinity" + ], + [ + 29311.1, + "infinity", + -1600.0 + ], + [ + 29331.1, + -1600.0, + -1600.0 + ], + [ + 29383.3, + -1600.0, + "infinity" + ], + [ + 29403.3, + "infinity", + "infinity" + ], + [ + 29457.2, + "infinity", + -490.0 + ], + [ + 29507.2, + -490.0, + -490.0 + ], + [ + 29531.0, + -490.0, + -901.4 + ] + ] + }, + "ETCS braking data": { + "M_NVAVADH": { + "unit": "-", + "value": 0 + }, + "Kt_int": { + "unit": "-", + "value": 1.15 + } + } +} \ No newline at end of file diff --git a/tests/units/integrators/testIntegrators.py b/tests/units/integrators/testIntegrators.py index 062aa64..2d71108 100644 --- a/tests/units/integrators/testIntegrators.py +++ b/tests/units/integrators/testIntegrators.py @@ -1,5 +1,6 @@ import unittest +from mseetc.journey import Journey from mseetc.ocp import casadiSolver from mseetc.track import Track from mseetc.train import Train @@ -17,47 +18,45 @@ def testAllIntegratorTypesWork(self): differ by a small relative tolerance. ''' - startPosition = 0 # [m] - endPosition = 5000 # [m] - duration = 5000 / (60 / 3.6) # [s] - tol = 0.1 numIntervals = 200 - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = 600 - track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') - track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tests/fixtures/tracks') track.updateTrainLengthDependentValues(train) + journey = Journey(config={'id': 'CH_StGallen_Wil_Journey_01'}, pathJSON='tests/fixtures/journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + opts = {'numIntervals': numIntervals, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}} - solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve() energy_RK_Approx = stats['Cost'] opts = {'numIntervals': numIntervals, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 0}} - solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve() energy_RK = stats['Cost'] opts = {'numIntervals': numIntervals, 'integrationMethod': 'IRK', 'integrationOptions': {'numApproxSteps': 1}} - solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve() energy_IRK_Approx = stats['Cost'] opts = {'numIntervals': numIntervals, 'integrationMethod': 'IRK', 'integrationOptions': {'numApproxSteps': 0}} - solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve() energy_IRK = stats['Cost'] opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES'} - solver = casadiSolver(train, track, opts) - df, stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + df, stats = solver.solve() energy_CVODES = stats['Cost'] diff --git a/tests/units/trainLengthDependentTrackAttributes/testCurvature.py b/tests/units/trainLengthDependentTrackAttributes/testCurvature.py index 1b3bba7..1f6c74c 100644 --- a/tests/units/trainLengthDependentTrackAttributes/testCurvature.py +++ b/tests/units/trainLengthDependentTrackAttributes/testCurvature.py @@ -4,6 +4,7 @@ import pandas as pd from matplotlib import pyplot as plt +from mseetc.journey import Journey from mseetc.ocp import casadiSolver, OptionsCasadiSolver from mseetc.track import Track from mseetc.train import Train, TrainIntegrator @@ -55,7 +56,7 @@ def test_cvodes_pwc_midpoint_curvature_converges_to_pwl_curvature(self): CVODES is used as the integrator. ''' - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') optsDict = {'integrationMethod': 'CVODES'} opts = OptionsCasadiSolver(optsDict) @@ -191,7 +192,7 @@ def test_rk_pwc_midpoint_curvature_converges_to_pwl_curvature(self): RK substeps are increased until convergence ''' - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') optsDict = {'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps': 0, 'numSteps': 1}} opts = OptionsCasadiSolver(optsDict) @@ -342,7 +343,7 @@ def test_rk_time_approx_pwc_midpoint_curvature_converges_to_pwl_curvature(self): Time approx steps are increased until convergence ''' - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') trainModel = train.exportModel() # Scenario @@ -488,7 +489,7 @@ def test_train_length_dependent_curvature_preserves_target_heading(self): headingTolerance = 1e-6 trainLength = 800 # [m] - track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tests/fixtures/tracks') # track needs to be straight at least train length meters before the end of the track track.curvatures = track.curvatures[track.curvatures.index < track.length - trainLength] @@ -496,7 +497,7 @@ def test_train_length_dependent_curvature_preserves_target_heading(self): df_heading_indep = computeHeadingFromCurvature(track.curvatures, track.length) - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = trainLength track.updateTrainLengthDependentValues(train) @@ -546,31 +547,29 @@ def test_pwc_curvature_approximation_matches_pwl_curvature_energy(self): linear curvatures or equivalent piecewise constant curvatures. ''' - startPosition = 0 # [m] - endPosition = 5000 # [m] - duration = 5000/(60/3.6) # [s] - energyRelativeTolerance = 1e-4 numIntervals = 100 - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = 600 - track = Track(config={'id': 'test_two_radii'}, pathJSON='tracks') - track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track = Track(config={'id': 'test_two_radii'}, pathJSON='tests/fixtures/tracks') track.updateTrainLengthDependentValues(train) + journey = Journey(config={'id': 'CH_StGallen_Wil_Journey_01'}, pathJSON='tests/fixtures/journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + # PWL Curvatures opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES'} - solver = casadiSolver(train, track, opts) - pwl_df, pwl_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + pwl_df, pwl_stats = solver.solve() energyConsumptionWithLinearTerms = pwl_stats['Cost'] # PWC Curvatures opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES', 'pwcLengthDependentTrackAttributes': True} - solver = casadiSolver(train, track, opts) - pwc_df, pwc_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + pwc_df, pwc_stats = solver.solve() energyConsumptionWithPwcTerms= pwc_stats['Cost'] @@ -616,27 +615,29 @@ def test_short_train_length_has_negligible_effect_on_curvature_energy_consumptio relativeTolerance = 1e-2 trainLength = 10 # [m] - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = trainLength - track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tests/fixtures/tracks') track.gradients = track.gradients.iloc[[0]] track.gradients["Gradient [permil]"].iloc[0] = 0 + journey = Journey(config={'id': 'CH_StGallen_Wil_Journey_02'}, pathJSON='tests/fixtures/journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + # track needs to be straight at least train length meters before the end of the track track.curvatures = track.curvatures[track.curvatures.index < track.length - trainLength] track.curvatures.loc[track.length - trainLength] = {"Curvature [1/m]": 0.0} - duration = track.length / (80/3.6) - # train-length-independent opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':0}, 'energyOptimal':True} - solver = casadiSolver(train, track, opts) - indep_df, indep_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + indep_df, indep_stats = solver.solve() + # train-length-dependent track.updateTrainLengthDependentValues(train) - solver = casadiSolver(train, track, opts) - dep_df, dep_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + dep_df, dep_stats = solver.solve() energyConsumptionIndependentOfTrainLength = indep_stats['Cost'] energyConsumptionDependentOfTrainLength = dep_stats['Cost'] diff --git a/tests/units/trainLengthDependentTrackAttributes/testGradient.py b/tests/units/trainLengthDependentTrackAttributes/testGradient.py index 6fd5018..d23a3fb 100644 --- a/tests/units/trainLengthDependentTrackAttributes/testGradient.py +++ b/tests/units/trainLengthDependentTrackAttributes/testGradient.py @@ -3,6 +3,7 @@ import numpy as np from matplotlib import pyplot as plt +from mseetc.journey import Journey from mseetc.ocp import casadiSolver, OptionsCasadiSolver from mseetc.track import Track, computeAltitude from mseetc.train import Train, TrainIntegrator @@ -26,7 +27,7 @@ def test_cvodes_pwc_midpoint_gradient_converges_to_pwl_gradient(self): CVODES is used as the integrator. ''' - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') optsDict = {'integrationMethod': 'CVODES'} opts = OptionsCasadiSolver(optsDict) @@ -164,7 +165,7 @@ def test_rk_pwc_midpoint_gradient_converges_to_pwl_gradient(self): RK substeps are increased until convergence ''' - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') optsDict = {'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps': 0, 'numSteps': 1}} opts = OptionsCasadiSolver(optsDict) @@ -317,7 +318,7 @@ def test_rk_time_approx_pwc_midpoint_gradient_converges_to_pwl_gradient(self): Time approx steps are increased until convergence ''' - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') trainModel = train.exportModel() # Scenario @@ -465,7 +466,7 @@ def test_train_length_dependent_gradient_preserves_target_altitude(self): altitudeTolerance = 1e-6 trainLength = 800 # [m] - track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tests/fixtures/tracks') # track needs to be flat at least train length meters before the end of the track track.gradients = track.gradients[track.gradients.index < track.length - trainLength] @@ -473,7 +474,7 @@ def test_train_length_dependent_gradient_preserves_target_altitude(self): df_alt = computeAltitude(track.gradients, track.length) - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = trainLength track.updateTrainLengthDependentValues(train) @@ -538,32 +539,29 @@ def test_pwc_gradient_approximation_matches_pwl_gradient_energy(self): Altitude should be 0 m at the end. ''' - startPosition = 0 # [m] - endPosition = 5000 # [m] - duration = 5000/(60/3.6) # [s] - altitudeTolerance = 1e-4 energyRelativeTolerance = 1e-3 - numIntervals = 50 - - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + numIntervals = 25 + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = 600 - track = Track(config={'id': 'test_one_hill'}, pathJSON='tracks') - track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track = Track(config={'id': 'test_one_hill'}, pathJSON='tests/fixtures/tracks') track.updateTrainLengthDependentValues(train) + journey = Journey(config={'id': 'test_one_hill_Journey_01'}, pathJSON='tests/fixtures/journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + # PWL Gradients opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES'} - solver = casadiSolver(train, track, opts) - pwl_df, pwl_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + pwl_df, pwl_stats = solver.solve() energyConsumptionWithLinearTerms = pwl_stats['Cost'] # PWC Gradients opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES', 'pwcLengthDependentTrackAttributes': True} - solver = casadiSolver(train, track, opts) - pwc_df, pwc_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + pwc_df, pwc_stats = solver.solve() energyConsumptionWithPwcTerms= pwc_stats['Cost'] @@ -624,27 +622,28 @@ def test_short_train_length_has_negligible_effect_on_energy_consumption(self): relativeTolerance = 0.01 trainLength = 10 # [m] - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') train.length = trainLength - track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tests/fixtures/tracks') track.curvatures = track.curvatures.iloc[[0]] track.curvatures["Curvature [1/m]"].iloc[0] = 0 + journey = Journey(config={'id': 'CH_StGallen_Wil_Journey_02'}, pathJSON='tests/fixtures/journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') + # track needs to be flat at least train length meters before the end of the track track.gradients = track.gradients[track.gradients.index < track.length - trainLength] track.gradients.loc[track.length - trainLength] = {"Gradient [permil]": 0.0} - duration = track.length / (80/3.6) - # train-length-independent opts = {'numIntervals':600, 'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps':0}, 'energyOptimal':True} - solver = casadiSolver(train, track, opts) - indep_df, indep_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + indep_df, indep_stats = solver.solve() track.updateTrainLengthDependentValues(train) - solver = casadiSolver(train, track, opts) - dep_df, dep_stats = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + dep_df, dep_stats = solver.solve() energyConsumptionIndependentOfTrainLength = indep_stats['Cost'] energyConsumptionDependentOfTrainLength = dep_stats['Cost'] diff --git a/tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py b/tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py index 8b99ffe..6ed3941 100644 --- a/tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py +++ b/tests/units/trainLengthDependentTrackAttributes/testSpeedLimit.py @@ -1,5 +1,6 @@ import unittest +from mseetc.journey import Journey from mseetc.ocp import casadiSolver from mseetc.track import Track from mseetc.train import Train @@ -16,29 +17,28 @@ def test_train_length_dependent_speed_limit_delays_acceleration(self): 22 m/s after the whole train has passed the speed-increase position. ''' - startPosition = 0 # [m] - endPosition = 12000 # [m] - duration = endPosition/(115/3.6) # [s] speedIncreasePosition = 1000 # [m] speedLimitBeforeIncrease = 22 # [m/s] speedTolerance = 0.001 # [m/s] - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') - track = Track(config={'id': 'test_one_speed_increase'}, pathJSON='tracks') - track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track = Track(config={'id': 'test_one_speed_increase'}, pathJSON='tests/fixtures/tracks') + + journey = Journey(config={'id': 'test_one_speed_increase_Journey_01'}, pathJSON='tests/fixtures/journeys') + track.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} - solver = casadiSolver(train, track, opts) - dfWithoutTrainLength, statsWithoutTrainLength = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + dfWithoutTrainLength, statsWithoutTrainLength = solver.solve() speedWithoutTrainLengthAfterIncrease = dfWithoutTrainLength[dfWithoutTrainLength['Position [m]'] > speedIncreasePosition].iloc[0]['Velocity [m/s]'] track.updateTrainLengthDependentValues(train) opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} - solver = casadiSolver(train, track, opts) - dfWithTrainLength, statsWithTrainLength = solver.solve(duration) + solver = casadiSolver(train, track, journey, opts) + dfWithTrainLength, statsWithTrainLength = solver.solve() speedWithTrainLengthAfterIncrease = dfWithTrainLength[dfWithTrainLength['Position [m]'] > speedIncreasePosition].iloc[0]['Velocity [m/s]'] speedWithTrainLengthAfterTrainPassedIncrease = dfWithTrainLength[dfWithTrainLength['Position [m]'] > (speedIncreasePosition+train.length)].iloc[0]['Velocity [m/s]'] diff --git a/tests/units/tunnelResistance/testTunnelResistance.py b/tests/units/tunnelResistance/testTunnelResistance.py index e3f600f..fffe500 100644 --- a/tests/units/tunnelResistance/testTunnelResistance.py +++ b/tests/units/tunnelResistance/testTunnelResistance.py @@ -1,5 +1,6 @@ import unittest +from mseetc.journey import Journey from mseetc.ocp import casadiSolver from mseetc.track import Track from mseetc.train import Train @@ -12,29 +13,29 @@ def test_tunnel_resistance_increases_energy_consumption(self): 26 km long small tunnel with cross section of 24 m^2 on a track of 28 km results in significant higher energy consumption. ''' - startPosition = 0 # [m] - endPosition = 28000 # [m] - duration = 28000/(145/3.6) # [s] - minEnergyRatio = 1.5 - train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='tests/fixtures/trains') + + trackWithoutTunnel = Track(config={'id': 'test_flat_no_tunnel'}, pathJSON='tests/fixtures/tracks') - trackWithoutTunnel = Track(config={'id': 'test_flat_no_tunnel'}, pathJSON='tracks') - trackWithoutTunnel.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + journey = Journey(config={'id': 'test_flat_no_tunnel_Journey_01'}, pathJSON='tests/fixtures/journeys') + trackWithoutTunnel.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} - solver = casadiSolver(train, trackWithoutTunnel, opts) - dfWithoutTunnel, statsWithoutTunnel = solver.solve(duration) + solver = casadiSolver(train, trackWithoutTunnel, journey, opts) + dfWithoutTunnel, statsWithoutTunnel = solver.solve() energyConsumptionWithoutTunnel = statsWithoutTunnel['Cost'] - trackWithTunnel = Track(config={'id': 'test_flat_with_tunnel'}, pathJSON='tracks') - trackWithTunnel.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + trackWithTunnel = Track(config={'id': 'test_flat_with_tunnel'}, pathJSON='tests/fixtures/tracks') + + journey = Journey(config={'id': 'test_flat_with_tunnel_Journey_01'}, pathJSON='tests/fixtures/journeys') + trackWithTunnel.updateLimits(positionStart=journey.positionStart, positionEnd=journey.positionEnd, unit='m') opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} - solver = casadiSolver(train, trackWithTunnel, opts) - dfWithTunnel, statsWithTunnel = solver.solve(duration) + solver = casadiSolver(train, trackWithTunnel, journey, opts) + dfWithTunnel, statsWithTunnel = solver.solve() energyConsumptionWithTunnel = statsWithTunnel['Cost']