diff --git a/mseetc/etcs.py b/mseetc/etcs.py new file mode 100644 index 0000000..9d56885 --- /dev/null +++ b/mseetc/etcs.py @@ -0,0 +1,531 @@ +from bisect import bisect_right +from dataclasses import dataclass + +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt + +from mseetc.track import Track + + +def shiftCurveByTime(dfCurve, timeShift): + + positionsOriginal = dfCurve.index.to_numpy() + velocitiesOriginal = dfCurve["Velocity [m/s]"].to_numpy() + + velocitiesShifted = velocitiesOriginal.copy() + positionsShifted = positionsOriginal - velocitiesShifted * timeShift + + dfCurveShifted = pd.DataFrame( + {"Velocity [m/s]": velocitiesShifted}, + index=positionsShifted, + ) + + dfCurveShifted.index.name = "Position [m]" + + return dfCurveShifted + + +def computeCeilingSpeedLimits(V_permitted_mps): + + dV_ebi_min = 7.5/3.6 + dV_ebi_max = 15.0/3.6 + V_ebi_min = 110.0/3.6 + V_ebi_max = 210.0/3.6 + + C_ebi = (dV_ebi_max - dV_ebi_min) / (V_ebi_max - V_ebi_min) + + if V_permitted_mps <= V_ebi_min: + + dV_ebi = dV_ebi_min + + else: + dV_ebi = min(dV_ebi_min + C_ebi * (V_permitted_mps - V_ebi_min),dV_ebi_max,) + + dV_warning = 0.5 * dV_ebi + dV_sbi = 0.75 * dV_ebi + + return { + "Warning [m/s]": V_permitted_mps + dV_warning, + "SBI [m/s]": V_permitted_mps + dV_sbi, + "EBI [m/s]": V_permitted_mps + dV_ebi, + } + + +def addStartPointToCurve(curve, velocity, start_position): + + delta_s = curve.index.to_numpy(dtype=float)[1] - curve.index.to_numpy(dtype=float)[0] + + start_point = pd.DataFrame( + {"Velocity [m/s]": [velocity, velocity]}, + index=[start_position, curve.index.to_numpy(dtype=float)[0] - delta_s], + ) + start_point.index.name = "Position [m]" + + return pd.concat([start_point, curve]) + + +def addEndPointToCurve(curve, velocity, end_position): + + delta_s = curve.index.to_numpy(dtype=float)[-1] - curve.index.to_numpy(dtype=float)[-2] + + end_point = pd.DataFrame( + {"Velocity [m/s]": [velocity, velocity]}, + index=[curve.index.to_numpy(dtype=float)[-1] + delta_s, end_position], + ) + end_point.index.name = "Position [m]" + + return pd.concat([curve, end_point]) + + +def trimCurveToMaxVelocity(curve, maxVelocity): + + velocities = curve["Velocity [m/s]"].to_numpy(dtype=float) + keep_mask = velocities <= maxVelocity + + return curve[keep_mask].copy() + + +def trimCurveFromMinVelocity(curve, minVelocity): + + velocities = curve["Velocity [m/s]"].to_numpy(dtype=float) + keep_mask = velocities >= minVelocity + + return curve[keep_mask].copy() + + +def trimCurveFromMinPosition(curve, minPosition): + + positions = curve.index.to_numpy(dtype=float) + keep_mask = positions >= minPosition + + return curve[keep_mask].copy() + + +@dataclass(frozen=True) +class BrakingTarget: + position: float # [m] + overlap: float # [m] + permittedVelocity: float # [m/s] + targetVelocity: float # [m/s] + + # EoA: End of authority + @property + def EoA(self): + return self.position + + # SvL: Supervised location + @property + def SvL(self): + return self.position + self.overlap + + +class EtcsBrakingCurveCalculator: + """ + Conventions + ----------- + - Position increases in the train running direction. + - Braking accelerations are stored as negative values. + - Gradient is positive for uphill and negative for downhill. + - A_gradient = g * gradient / 1000. + - Curves are computed backwards from the target position. + """ + + def __init__(self, trainBrakingData, track, distancePre=3000, distancePost=1000): + + self.trainBrakingData = trainBrakingData + self.track = track + + # Compute the braking curve using fixed time steps of length dt. + self.dt = 0.1 # [s] + + # speed decrease is plotted from BrakingTarget.position - distancePre until BrakingTarget.position + distancePost + self.distancePre = distancePre # [m] + self.distancePost = distancePost # [m] + + # ETCS Constants + self.T_warning = 2.0 # [s] + self.T_driver = 4.0 # [s] + + self.curveStyles = { + "EBD": {"color": "blue", "linestyle": "-", "linewidth": 2.0}, + "EBI": {"color": "blue", "linestyle": ":", "linewidth": 1.5}, + "SBD": {"color": "green", "linestyle": "-", "linewidth": 2.0}, + "SBI1": {"color": "green", "linestyle": "--", "linewidth": 1.5}, + "SBI2": {"color": "blue", "linestyle": "--", "linewidth": 1.5}, + "SBI": {"color": "red", "linestyle": "-", "linewidth": 2.0}, + "W": {"color": "orange", "linestyle": "-", "linewidth": 2.0}, + "P": {"color": "grey", "linestyle": "-", "linewidth": 2.0}, + "I": {"color": "gold", "linestyle": "-", "linewidth": 2.0}, + } + + + 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: + raise ValueError("EoA must lie within the track length.") + + if not 0 < target.SvL < self.track.length: + raise ValueError("SvL must lie within the track length.") + + if not 0 <= target.targetVelocity < target.permittedVelocity: + raise ValueError("targetVelocity must be lower than permittedVelocity.") + + + def computeABrakeSafe(self): + + trainBrakingData = self.trainBrakingData + + 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_wet_corr = K_wet_rst + M_NVAVADH * (1 - K_wet_rst) + + A_brake_safe_values = [ + A_emergency * K_dry_rst * K_wet_corr + for A_emergency in A_emergency_values + ] + + return { + "velocity [m/s]": velocities, + "value [m/s^2]": A_brake_safe_values, + } + + + def computeAGradient(self, currentPosition): + + positions = self.track.gradients.index.values + gradients = self.track.gradients["Gradient [permil]"].values + + 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. + return 0 + + gradient = gradients[idx] + return 9.81 * gradient * 0.001 + + + def computeBrakingCurve(self, brakingProfile, targetPosition, permittedVelocity, targetVelocity): + """ + Compute a braking curve backwards from a target position and target velocity. + + The braking profile is defined by velocity thresholds and corresponding braking decelerations. + The curve is integrated backwards in fixed time steps until the maximum relevant velocity is reached. + """ + + maxVelocity = permittedVelocity * 1.4 + + positions = [targetPosition] + velocities = [targetVelocity] + + thresholdVelocities = list(brakingProfile["velocity [m/s]"]) + brakingValues = list(brakingProfile["value [m/s^2]"]) + + openEndedVelocityThreshold = 200 # [m/s], practical upper bound for open-ended last interval of the braking profile, basically inf velocity for a train + thresholdVelocities.append(openEndedVelocityThreshold ) + + for idx in range(len(thresholdVelocities) - 1): + + upperThreshold = thresholdVelocities[idx + 1] + + # The curve starts at targetVelocity, so lower velocity ranges are not relevant. + if targetVelocity > upperThreshold: + + continue + + A_brake = brakingValues[idx] + + while velocities[-1] < upperThreshold and velocities[-1] < maxVelocity: + + A_gradient = self.computeAGradient(positions[-1]) + + v_old = velocities[-1] + v_new = v_old - (A_brake - A_gradient) * self.dt + x_new = positions[-1] - 0.5 * (v_new + v_old) * self.dt + + positions.append(x_new) + velocities.append(v_new) + + if velocities[-1] >= maxVelocity: + + break + + curve = pd.DataFrame( + {"Velocity [m/s]": velocities[::-1]}, + index=positions[::-1], + ) + + curve.index.name = "Position [m]" + + return curve + + + 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 [%]"] * 0.01 + + positionsEBD = EBD_curve.index.to_numpy() + velocitiesEBD = EBD_curve["Velocity [m/s]"].to_numpy() + + t_be = T_be * Kt_int + T_berem = max(t_be - T_traction, 0) + + positionsEBI = [] + velocitiesEBI = [] + + for pos, vel in zip(positionsEBD, velocitiesEBD): + + A_est1 = 0.1 # todo + A_est2 = min(A_est1, 0.4) + + V_est = (vel - A_est1 * T_traction - A_est2 * T_berem) / (1 + v_uncertainty) + V_est = max(V_est, 0.0) + + if V_est < targetVelocity: + + # Stop once the estimated velocity drops below the target velocity. + # Add the target velocity as the final point, using a small position offset as a simplified position. + # This is only used to terminate the plotted EBI curve at targetVelocity. + velocitiesEBI.append(targetVelocity) + positionsEBI.append(positionsEBI[-1] + 1) + + break + + velocitiesEBI.append(V_est) + + V_delta_0 = V_est * v_uncertainty + V_delta1 = A_est1 * T_traction + V_delta2 = A_est2 * T_berem + D_bec = T_traction * (V_est + V_delta_0 + 0.5 * V_delta1) + T_berem * (V_est + V_delta_0 + V_delta1 + 0.5 * V_delta2) + positionsEBI.append(pos - D_bec) + + EBI_curve = pd.DataFrame( + {"Velocity [m/s]": velocitiesEBI}, + index=positionsEBI, + ) + + EBI_curve.index.name = "Position [m]" + + return EBI_curve + + + def computeSBICurve(self, SBI1_curve, SBI2_curve): + + positionsSBI1 = SBI1_curve.index.to_numpy() + velocitiesSBI1 = SBI1_curve["Velocity [m/s]"].to_numpy() + + positionsSBI2 = SBI2_curve.index.to_numpy() + velocitiesSBI2 = SBI2_curve["Velocity [m/s]"].to_numpy() + + # Use only the overlapping position range + minPosition = max(positionsSBI1.min(), positionsSBI2.min()) + maxPosition = min(positionsSBI1.max(), positionsSBI2.max()) + + step = 10.0 # [m] + positionsSBI = np.arange(minPosition, maxPosition + step, step) + positionsSBI[-1] = maxPosition + + velocitiesSBI1_interpol = np.interp(positionsSBI, positionsSBI1, velocitiesSBI1) + velocitiesSBI2_interpol = np.interp(positionsSBI, positionsSBI2, velocitiesSBI2) + + # At each position, take the lower speed of SBI1 and SBI2. + # This gives the more restrictive plotted SBI curve. + velocitiesSBI = np.minimum(velocitiesSBI1_interpol, velocitiesSBI2_interpol) + + SBI_curve = pd.DataFrame( + {"Velocity [m/s]": velocitiesSBI}, + index=positionsSBI, + ) + + SBI_curve.index.name = "Position [m]" + + return SBI_curve + + + def processCurvesBeforeTarget(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) + + curves["EBD"] = trimCurveFromMinPosition(curves["EBD"], curves["EBI"].index.to_numpy(dtype=float)[2]) + + curves["SBI2"] = trimCurveFromMinPosition(curves["SBI2"], curves["EBI"].index.to_numpy(dtype=float)[2]) + + curves["SBD"] = trimCurveFromMinPosition(curves["SBD"], curves["SBI"].index.to_numpy(dtype=float)[2]) + + curves["SBI1"] = trimCurveFromMinPosition(curves["SBI1"], curves["SBI"].index.to_numpy(dtype=float)[2]) + + return curves + + def processCurvcesAfterTarget(self, curves, target): + + targetVelocity = target.targetVelocity + end_position = target.position + self.distancePost + + speedLimits = computeCeilingSpeedLimits(targetVelocity) + + curves["EBI"] = trimCurveFromMinVelocity(curves["EBI"], speedLimits["EBI [m/s]"]) + curves["EBI"] = addEndPointToCurve(curves["EBI"], speedLimits["EBI [m/s]"], end_position) + + curves["SBI"] = trimCurveFromMinVelocity(curves["SBI"], speedLimits["SBI [m/s]"]) + curves["SBI"] = addEndPointToCurve(curves["SBI"], speedLimits["SBI [m/s]"], end_position) + + curves["W"] = trimCurveFromMinVelocity(curves["W"], speedLimits["Warning [m/s]"]) + curves["W"] = addEndPointToCurve(curves["W"], speedLimits["Warning [m/s]"], end_position) + + curves["P"] = addEndPointToCurve(curves["P"], targetVelocity, end_position) + + curves["I"] = addEndPointToCurve(curves["I"], targetVelocity, curves["P"].index.to_numpy(dtype=float)[-4]) + + curves["EBD"] = trimCurveFromMinVelocity(curves["EBD"], speedLimits["EBI [m/s]"]) + + curves["SBI2"] = trimCurveFromMinVelocity(curves["SBI2"], speedLimits["SBI [m/s]"]) + + curves["SBD"] = trimCurveFromMinVelocity(curves["SBD"], speedLimits["EBI [m/s]"]) + + curves["SBI1"] = trimCurveFromMinVelocity(curves["SBI1"], speedLimits["SBI [m/s]"]) + + return curves + + + 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 + + curves = {} + + curves["EBD"] = self.computeBrakingCurve(ABrakeSafeProfile, target.SvL, target.permittedVelocity, target.targetVelocity) + + curves["EBI"] = self.computeEBICurve(curves["EBD"], target.targetVelocity) + + curves["SBI2"] = shiftCurveByTime(curves["EBI"], trainBrakingData["T_bs2 [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["SBI"] = self.computeSBICurve(curves["SBI1"], curves["SBI2"]) + + curves["W"] = shiftCurveByTime(curves["SBI"], self.T_warning) + + curves["P"] = shiftCurveByTime(curves["SBI"], self.T_driver) + + curves["I"] = shiftCurveByTime(curves["P"], T_indication) + + curves = self.processCurvesBeforeTarget(curves, target) + + if target.targetVelocity > 0: + + curves = self.processCurvcesAfterTarget(curves, target) + + return curves + + + def plotCurves(self, curves, target): + + targetPosition = target.EoA + permittedVelocity = target.permittedVelocity + targetVelocity = target.targetVelocity + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.step( + np.array( + [targetPosition - self.distancePre, targetPosition, targetPosition + self.distancePost]) / 1000, + np.array([permittedVelocity, permittedVelocity, targetVelocity]) * 3.6, + label="Speed limit", color="black", linewidth=2.0 + ) + + for name, curve in curves.items(): + + style = {} + + if self.curveStyles is not None and name in self.curveStyles: + style = self.curveStyles[name] + + ax.plot(curve.index.values / 1000, curve["Velocity [m/s]"] * 3.6, label=name, **style) + + ax.set_title("ETCS Braking Curves") + 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.figure.tight_layout() + + plt.show() + + +if __name__ == '__main__': + + 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, + } + + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='../tracks') + + target = BrakingTarget( + position=5000, + overlap= 100, + permittedVelocity=160/3.6, + targetVelocity=0/3.6 + ) + + calculator = EtcsBrakingCurveCalculator(trainBrakingData, track, distancePre=5000, distancePost=1000) + curve_set = calculator.computeTarget(target) + + calculator.plotCurves(curve_set, target) \ No newline at end of file diff --git a/mseetc/ocp.py b/mseetc/ocp.py index a9e321d..67ed55b 100644 --- a/mseetc/ocp.py +++ b/mseetc/ocp.py @@ -6,7 +6,7 @@ from mseetc.track import computeDiscretizationPoints -from mseetc.utils import Options, var, postProcessDataFrame, splitLosses +from mseetc.utils import Options, var, postProcessDataFrame, splitLosses, computeTunnelFactor class OptionsCasadiSolver(Options): @@ -27,6 +27,10 @@ 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.pwcLengthDependentTrackAttributes = False + super().__init__(paramsDict) @@ -121,7 +125,7 @@ def __init__(self, train, track, optsDict={}): # track parameters - self.points = computeDiscretizationPoints(track, numIntervals) + self.points = computeDiscretizationPoints(track, numIntervals, opts) self.steps = np.diff(self.points.index) # real-time parameters @@ -193,16 +197,28 @@ def __init__(self, train, track, optsDict={}): # gradient and curvature of current index grad = self.points.iloc[i]['Gradient [permil]']/1e3 + gradLinearTerm = self.points.iloc[i]["Gradient linear term [permil/m]"]/1e3 curv = self.points.iloc[i]['Curvature [1/m]'] + curvLinearTerm = self.points.iloc[i]["Curvature linear term [1/m^2]"] + crossSection = self.points.iloc[i]['CrossSection [m^2]'] + tunnelFactor = computeTunnelFactor(crossSection, train, opts) # acceleration constraints - g += [trainModel.accelerationFun(ca.vertcat(time[i], velSq[i]), ca.vcat(u), grad, curv)] + g += [trainModel.accelerationFun(ca.vertcat(time[i], velSq[i], 0), ca.vcat(u), grad, gradLinearTerm, curv, curvLinearTerm, tunnelFactor)] lbg += [accMin] ubg += [accMax] # coupling constraints - out = trainIntegrator.solve(time=time[i], velocitySquared=velSq[i], ds=self.steps[i], - traction=Fel[i], pnBrake=Fpb[i], gradient=grad, curvature=curv) + out = trainIntegrator.solve(time=time[i], + velocitySquared=velSq[i], + ds=self.steps[i], + traction=Fel[i], + pnBrake=Fpb[i], + gradient=grad, + gradientLinearTerm=gradLinearTerm, + curvature=curv, + curvatureLinearTerm=curvLinearTerm, + tunnelFactor=tunnelFactor) xNxt1 = ca.vertcat(time[i+1], velSq[i+1]) xNxt2 = ca.vertcat(out['time'], out['velSquared']) @@ -322,7 +338,8 @@ def solve(self, terminalTime, initialTime=0, terminalVelocity=1, initialVelocity # initial guess # NOTE: good idea vel0 to be compatible with f0 (power-wise) to avoid nans at first iteration - vel0 = (60/3.6)**2 + vel_avg = (self.points.index[-1]-self.points.index[0])/(terminalTime-initialTime) * 0.95 + vel0 = vel_avg*vel_avg dt = (terminalTime - initialTime)/self.numIntervals t0 = initialTime @@ -437,4 +454,4 @@ def solve(self, terminalTime, initialTime=0, terminalVelocity=1, initialVelocity else: - print("Solver failed!") \ No newline at end of file + print("Solver failed!") diff --git a/mseetc/track.py b/mseetc/track.py index a214aa5..086c2d0 100644 --- a/mseetc/track.py +++ b/mseetc/track.py @@ -6,7 +6,10 @@ import numpy as np import pandas as pd -from mseetc.utils import checkTTOBenchVersion, convertUnit +from mseetc.utils import checkTTOBenchVersion, convertUnit, pickEquallySpacedPoints, plotSpeedLimits, plotGradients, \ + plotCurvatures + +plotDebug = False def importTuples(tuples, xLabel, yLabels): """ @@ -88,14 +91,14 @@ def computeAltitude(gradients, length, altitudeStart=0): return df -def computeDiscretizationPoints(track, numIntervals): +def computeDiscretizationPoints(track, numIntervals, opts): """ Compute the space discretization points based on track characteristics and horizon length. """ df1 = track.mergeDataFrames() - pos = np.linspace(0, track.length, numIntervals + 1 - (len(df1) - 1)) + pos = pickEquallySpacedPoints(0, track.length, numIntervals, df1.index.to_numpy(dtype=float)) df2 = pd.DataFrame({'position [m]':pos}).set_index('position [m]') df3 = df2.join(df1, how='outer').ffill() @@ -104,9 +107,96 @@ def computeDiscretizationPoints(track, numIntervals): raise ValueError("Wrong number of computed discretization intervals!") + adaptConstantTrackAttributesToNewShootingNodes(df3, numIntervals) + + if opts.pwcLengthDependentTrackAttributes: + + makePwcLengthDependentTrackAttibutes(df3) + return df3 +def adaptConstantTrackAttributesToNewShootingNodes(df, numIntervals): + ''' + Adapt gradient and curvature values to the current shooting nodes. + + For piecewise linear profiles, the constant term at each node is updated using the linear term from the previous interval. + + Example: + Gradient at 0 m is 10 permil and the linear term is 0.01 permil/m. + At a new shooting node at 100 m, the gradient becomes: + + 10 + 100 * 0.01 = 11 permil + + Without a linear term, the value would remain 10 permil. + ''' + + if "Gradient linear term [permil/m]" in df.columns: + + positions = df.index.to_numpy(dtype=float) + grads = [df["Gradient [permil]"].iloc[0]] + + for idx in range(1, numIntervals + 1): + + if np.isclose(df["Gradient [permil]"].iloc[idx - 1], df["Gradient [permil]"].iloc[idx]): + + grads.append(grads[-1] + (positions[idx] - positions[idx - 1]) * df["Gradient linear term [permil/m]"].iloc[idx - 1]) + + else: + + grads.append(df["Gradient [permil]"].iloc[idx]) + + df["Gradient [permil]"] = grads + + else: + + df["Gradient linear term [permil/m]"] = np.zeros(len(df)) + + + if "Curvature linear term [1/m^2]" in df.columns: + + positions = df.index.to_numpy(dtype=float) + curvs = [df["Curvature [1/m]"].iloc[0]] + + for idx in range(1, numIntervals + 1): + + if np.isclose(df["Curvature [1/m]"].iloc[idx - 1], df["Curvature [1/m]"].iloc[idx]): + + curvs.append(curvs[-1] + (positions[idx] - positions[idx - 1]) * df["Curvature linear term [1/m^2]"].iloc[idx - 1]) + + else: + + curvs.append(df["Curvature [1/m]"].iloc[idx]) + + df["Curvature [1/m]"] = curvs + + else: + + df["Curvature linear term [1/m^2]"] = np.zeros(len(df)) + + +def makePwcLengthDependentTrackAttibutes(df): + + positions = df.index.to_numpy(dtype=float) + + g_pwl = df["Gradient [permil]"].to_numpy(dtype=float) + g_linear = df["Gradient linear term [permil/m]"].to_numpy(dtype=float) + + c_pwl = df["Curvature [1/m]"].to_numpy(dtype=float) + c_linear = df["Curvature linear term [1/m^2]"].to_numpy(dtype=float) + + ds = positions[1:] - positions[:-1] + + g_pwc = g_pwl[:-1] + 0.5 * g_linear[:-1] * ds + c_pwc = c_pwl[:-1] + 0.5 * c_linear[:-1] * ds + + df["Gradient [permil]"] = np.r_[g_pwc, g_pwc[-1]] + df["Curvature [1/m]"] = np.r_[c_pwc, c_pwc[-1]] + + df["Gradient linear term [permil/m]"] = np.zeros(len(df)) + df["Curvature linear term [1/m^2]"] = np.zeros(len(df)) + + class Track(): CURVATURE_THRESHOLD = 1/150 # absolute value of maximum allowed cruvature [1/m] @@ -132,7 +222,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'tracks'): data = json.load(file) - checkTTOBenchVersion(data, ['1.1', '1.2', '1.3']) + checkTTOBenchVersion(data, ['1.1', '1.2', '1.3', '1.4']) # read data self.length = convertUnit(data['stops']['values'][-1], data['stops']['unit']) @@ -149,6 +239,11 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'tracks'): 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']['length'] if 'tunnels' in data else 'm', + data['tunnels']['units']['cross section'] if 'tunnels' in data else 'm^2') + + numStops = len(data['stops']['values']) indxDeparture = config['from'] if 'from' in config else 0 indxDestination = config['to'] if 'to' in config else numStops-1 @@ -193,6 +288,11 @@ def curvaturesOk(self): return True if self.curvatures.shape[0] > 0 and checkDataFrame(self.curvatures, self.length) else False + def crossSectionsOk(self): + + return True if self.crossSections.shape[0] > 0 and checkDataFrame(self.crossSections, self.length) else False + + def checkFields(self): if not self.lengthOk(): @@ -215,6 +315,10 @@ def checkFields(self): raise ValueError("Issue with track curvatures!") + if not self.crossSectionsOk(): + + raise ValueError("Issue with tunnel cross sections!") + def importGradientTuples(self, tuples, unit='permil'): @@ -263,6 +367,7 @@ def importCurvatureTuples(self, tuples, unitRadiusStart='m', unitRadiusEnd='m', tuples = self.sampleClothoid(tuples, clothoidSamplingInterval) self.curvatures = importTuples(tuples, 'Position [m]', ['Curvature [1/m]']) + self.curvatures["Curvature [1/m]"] = self.curvatures["Curvature [1/m]"].abs() checkDataFrame(self.curvatures, self.length) @@ -348,6 +453,46 @@ def sampleClothoid(self, tuples, ds=None): return result + def importTunnelTuples(self, tuples, unitLength='m', unitCrossSection='m^2'): + + if not self.lengthOk(): + + raise ValueError("Cannot import tunnels without a valid track length!") + + 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] + self.crossSections = importTuples(tuples, 'Position [m]', ['Length [m]', 'CrossSection [m^2]']) + + + # get end of tunnel positions and assign them a cross section of inf + positions = self.crossSections.index.astype(float) + tunnelLengths = self.crossSections["Length [m]"].astype(float) + + endOfTunnelPositions = positions + tunnelLengths + + # a tunnel may change its cross section, therefore some end of tunnel positions need to be removed + endOfTunnelPositions = [ + e for e in endOfTunnelPositions + if not any(abs((p - e)) < 0.1 for p in positions) + ] + + nonTunnelSections_df = pd.DataFrame({"Length [m]": 0.0, "CrossSection [m^2]": float("inf")}, index=endOfTunnelPositions) + nonTunnelSections_df.index.name = self.crossSections.index.name + self.crossSections = pd.concat([self.crossSections, nonTunnelSections_df]).sort_index() + + if positions[0] != 0.0: + first_row = {"Position [m]": 0.0, "Length [m]": 0.0, "CrossSection [m^2]": float("inf")} + self.crossSections.loc[0] = first_row + self.crossSections = self.crossSections.sort_index() + + self.crossSections.drop(columns=["Length [m]"], inplace=True) + + checkDataFrame(self.crossSections, self.length) + + def reverse(self): # switch to opposite trip @@ -368,6 +513,7 @@ def flipData(df): self.gradients = -flipData(self.gradients) self.speedLimits = flipData(self.speedLimits) self.curvatures = -flipData(self.curvatures) + self.crossSections = flipData(self.crossSections) self.title = self.title + ' (reversed)' @@ -380,7 +526,8 @@ def mergeDataFrames(self): """ joinedGradientsAndSpeedLimits = self.gradients.join(self.speedLimits, how='outer').fillna(method='ffill') - return self.curvatures.join(joinedGradientsAndSpeedLimits, how='outer').fillna(method='ffill') + joined = self.curvatures.join(joinedGradientsAndSpeedLimits, how='outer').fillna(method='ffill') + return self.crossSections.join(joined, how='outer').fillna(method='ffill') def print(self): """ @@ -448,6 +595,151 @@ def crop(dfIn): self.speedLimits = crop(self.speedLimits) self.gradients = crop(self.gradients) self.curvatures = crop(self.curvatures) + self.crossSections = crop(self.crossSections) + + + def updateTrainLengthDependentValues(self, train): + + self.updateSpeedLimitsToTrainLength(train.length) + self.updateGradientsToTrainLength(train.length) + self.updateCurvaturesToTrainLength(train.length) + + + def updateSpeedLimitsToTrainLength(self, trainLength): + + v = self.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + pos = self.speedLimits.index.to_numpy(dtype=float) + + if len(pos) > 1: + + pos_adj = [] + v_adj = [] + + for i in range(len(pos)): + new_pos = pos[i] + + # Delay speed increases by train length + if i > 0 and v[i] > v[i - 1]: + new_pos += trainLength + + # Skip points outside the track + if new_pos >= self.length: + continue + + # Remove previous points that are now after this point + while pos_adj and pos_adj[-1] > new_pos: + pos_adj.pop() + v_adj.pop() + + pos_adj.append(new_pos) + v_adj.append(v[i]) + + if plotDebug: + + plotSpeedLimits(self, np.asarray(pos_adj, dtype=float), np.asarray(v_adj, dtype=float)) + + self.speedLimits = pd.DataFrame( + {"Speed limit [m/s]": v_adj}, + index=pos_adj, + ) + + + def updateTrackAttributeToTrainLength(self, df, trainLength): + """ + Convert pointwise track attribute changes into train-length-dependent piecewise linear values. + + A step at position s does not affect the whole train at once. Instead, its + effect is spread over one train length, from s to s + trainLength. + """ + + if "Gradient [permil]" in df.columns: + + valueColumn = "Gradient [permil]" + linearTermColumn = "Gradient linear term [permil/m]" + plotFunction = plotGradients + + elif "Curvature [1/m]" in df.columns: + + valueColumn = "Curvature [1/m]" + linearTermColumn = "Curvature linear term [1/m^2]" + plotFunction = plotCurvatures + + else: + + raise ValueError("Unknown track attribute DataFrame!") + + values = df[valueColumn].to_numpy(dtype=float) + positions = df.index.to_numpy(dtype=float) + + if len(positions) <= 1: + + return df + + # Assume the train starts before the track with a flat and straight track. + values = np.r_[0.0, values] + positions = np.r_[-trainLength, positions] + + # Each original step is spread linearly over one train length assuming uniform mass. + # The first point has no previous value, so its slope contribution is zero. + jumpSlopes = np.r_[0.0, (values[1:] - values[:-1]) / trainLength] + + # New breakpoints occur both when the front of the train reaches a step and when the rear of the train has passed it. + adjustedPositions = np.sort(np.unique(np.r_[positions, positions + trainLength])) + adjustedPositions = adjustedPositions[adjustedPositions < self.length] + + adjustedValues = [0.0] + linearTerms = [0.0] + + epsilon = 1e-3 + + for idx in range(1, len(adjustedPositions)): + + currentPosition = adjustedPositions[idx] + previousPosition = adjustedPositions[idx - 1] + + intervalLength = currentPosition - previousPosition + + # Continue the previous linear value to the current position. + currentValue = adjustedValues[-1] + intervalLength * linearTerms[-1] + + # Active jumps are those currently within one train length behind the train front. + activeJumpMask = ( + (currentPosition - trainLength + epsilon < positions) + & (positions < currentPosition + epsilon) + ) + + currentLinearTerm = np.sum(jumpSlopes[activeJumpMask]) + + adjustedValues.append(currentValue) + linearTerms.append(currentLinearTerm) + + # Remove artificial point at -trainLength. + adjustedPositions = adjustedPositions[1:] + adjustedValues = adjustedValues[1:] + linearTerms = linearTerms[1:] + + if plotDebug and plotFunction is not None: + + plotFunction(self, np.asarray(adjustedPositions, dtype=float), np.asarray(adjustedValues, dtype=float), np.asarray(linearTerms, dtype=float) +) + + return pd.DataFrame( + { + valueColumn: adjustedValues, + linearTermColumn: linearTerms + }, + index=adjustedPositions + ) + + + def updateGradientsToTrainLength(self, trainLength): + + self.gradients = self.updateTrackAttributeToTrainLength(self.gradients, trainLength) + + + def updateCurvaturesToTrainLength(self, trainLength): + + self.curvatures = self.updateTrackAttributeToTrainLength(self.curvatures, trainLength) if __name__ == '__main__': diff --git a/mseetc/train.py b/mseetc/train.py index b9b0e36..9ef6d4c 100644 --- a/mseetc/train.py +++ b/mseetc/train.py @@ -31,7 +31,7 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> data = json.load(file) - checkTTOBenchVersion(data, ['1.1', '1.2', '1.3']) + checkTTOBenchVersion(data, ['1.1', '1.2', '1.3', '1.4']) # overwrite json data with config values if applicable @@ -66,6 +66,8 @@ 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.mass = convertUnit(data['mass']['value'], data['mass']['unit']) # train mass [kg] self.rho = convertUnit(data['rho']['value'], data['rho']['unit']) # rotating-mass factor [-] @@ -110,11 +112,34 @@ def __init__(self, config, pathJSON=Path(__file__).parent.parent / 'trains') -> self.etaTraction = convertUnit(data['efficiency traction']['value'], data['efficiency traction']['unit']) self.etaRgBrake = convertUnit(data['efficiency reg brake']['value'], data['efficiency reg brake']['unit']) + if "tunnel resistance" in data: + + tunnel_data = 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] + + self.tunnelCoefficients = { + convertUnit(cross_section, cross_section_unit): convertUnit( + resistance, + resistance_unit + ) + for cross_section, resistance in tunnel_data["values"] + } + + else: + + self.tunnelCoefficients = {} + self.checkFields() def checkFields(self): + if self.length is None or self.length < 0 or np.isinf(self.length): + + raise ValueError("Train length must be a positive number, not {}!".format(self.length)) + if self.mass is None or self.mass < 0 or np.isinf(self.mass): raise ValueError("Train mass must be a positive number, not {}!".format(self.mass)) @@ -172,6 +197,17 @@ 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): + + raise ValueError("Tunnel cross section must be positive, not {}!".format(crossSection)) + + if coef is None or coef <= 0 or np.isinf(coef): + + raise ValueError("Tunnel resistance coefficient must be positive, not {}!".format(coef)) + + def exportModel(self): "Export train model (ODE and relevant train data)." @@ -226,45 +262,54 @@ def __init__(self, sr0, sr1, sr2, rho=1, g=9.81, withPnBrake=True) -> None: # states - time = ca.MX.sym('time') - velocitySquared = ca.MX.sym('velocitySquared') + time = ca.MX.sym('time') # [s] + velocitySquared = ca.MX.sym('velocitySquared') # [m^2/s^2] + position = ca.MX.sym('position') # [m] - x = ca.vertcat(time, velocitySquared) + x = ca.vertcat(time, velocitySquared, position) # controls - traction = ca.MX.sym('traction') - pnBrake = ca.MX.sym('pnBrake') + traction = ca.MX.sym('traction') # [N/kg] + pnBrake = ca.MX.sym('pnBrake') # [N/kg] u = ca.vertcat(traction, pnBrake if withPnBrake else []) # parameters - gradient = ca.MX.sym('gradient') - curvature = ca.MX.sym('curvature') + gradient = ca.MX.sym('gradient') # [-] + gradientLinearTerm = ca.MX.sym('gradientLinearTerm') # [1/m] + curvature = ca.MX.sym('curvature') # [1/m] + curvatureLinearTerm = ca.MX.sym('curvatureLinearTerm') # [1/m^2] + tunnelFactor = ca.MX.sym('tunnelFactor') # [1/m] ds = ca.MX.sym('ds') - p = ca.vertcat(gradient, curvature, ds) + p = ca.vertcat(gradient, gradientLinearTerm, curvature, curvatureLinearTerm, tunnelFactor, ds) # ODE - rollingResistance = sr0 + sr1*ca.sqrt(velocitySquared) + sr2*velocitySquared - curvatureResistance = ca.if_else(ca.fabs(curvature)<=1/300, g*0.5*ca.fabs(curvature)/(1-30*ca.fabs(curvature)), - g*0.65*ca.fabs(curvature)/(1-55*ca.fabs(curvature))) - acceleration = traction + (pnBrake if withPnBrake else 0) - rollingResistance - g*gradient*(1/rho) - curvatureResistance*(1/rho) + rollingResistance = sr0 + sr1*ca.sqrt(velocitySquared) + sr2*velocitySquared # [N/kg] + gradientResistance = g*(1/rho)*(gradient+gradientLinearTerm*position) # [N/kg] + curvatureResistance = (1/rho)*(5.07*(curvature+curvatureLinearTerm*position)) # [N/kg] + tunnelResistance = tunnelFactor*velocitySquared # [N/kg] + + acceleration = traction + (pnBrake if withPnBrake else 0) - rollingResistance - gradientResistance - curvatureResistance - tunnelResistance # [m/s^2] + timeODE = 1/ca.sqrt(velocitySquared) velocityODE = 2*acceleration + positionODE = 1 timeODE *= ds velocityODE *= ds + positionODE *= ds - fExplicit = ca.vertcat(timeODE, velocityODE) + fExplicit = ca.vertcat(timeODE, velocityODE, positionODE) # model self.ode = fExplicit self.acceleration = acceleration - self.accelerationFun = ca.Function('a', [x, u, gradient, curvature], [acceleration]) + self.accelerationFun = ca.Function('a', [x, u, gradient, gradientLinearTerm, curvature, curvatureLinearTerm, tunnelFactor], [acceleration]) self.rollingResistance = rollingResistance self.parameters = p self.controls = u @@ -295,8 +340,12 @@ def __init__(self, model, solver, optsDict={}) -> None: opts = OptionsRK(optsDict) - ode = model.ode if opts.numApproxSteps == 0 else model.ode[1] - states = model.states if opts.numApproxSteps == 0 else model.states[1] + if opts.numApproxSteps == 0: + ode = model.ode + states = model.states + else: + ode = ca.vertcat(model.ode[1], model.ode[2]) + states = ca.vertcat(model.states[1], model.states[2]) self.eval = ca.simpleRK(ca.Function('ode', [states, params], [ode]), opts.numSteps, opts.order) @@ -304,20 +353,22 @@ def __init__(self, model, solver, optsDict={}) -> None: opts = OptionsIRK(optsDict) - ode = model.ode if opts.numApproxSteps == 0 else model.ode[1] - states = model.states if opts.numApproxSteps == 0 else model.states[1] + if opts.numApproxSteps == 0: + ode = model.ode + states = model.states + else: + ode = ca.vertcat(model.ode[1], model.ode[2]) + states = ca.vertcat(model.states[1], model.states[2]) - self.eval = ca.simpleIRK(ca.Function('ode', [states, params], [ode]), opts.numSteps, opts.order, opts.collMethod, 'fast_newton', {'max_iter':opts.maxIter, 'jit':opts.jit, 'error_on_fail':False}) + self.eval = ca.simpleIRK(ca.Function('ode', [states, params], [ode]), opts.numSteps, opts.order, opts.collMethod, 'fast_newton', {'max_iter': opts.maxIter, 'jit': opts.jit, 'error_on_fail': False}) elif solver == 'CVODES': opts = OptionsCVODES(optsDict) opts.numApproxSteps = 0 - states = model.states - t0, tf = 0, 1 - cvodesFun = ca.integrator('integrator', 'cvodes', {'x':model.states, 'p':params, 'ode':model.ode}, t0, tf, {'abstol':opts.absTol, 'reltol':opts.relTol}) + cvodesFun = ca.integrator('integrator', 'cvodes', {'x': model.states, 'p': params, 'ode': model.ode}, t0, tf, {'abstol': opts.absTol, 'reltol': opts.relTol}) self.eval = lambda x, p, dummy: cvodesFun(x0=x, p=p)['xf'] @@ -327,24 +378,24 @@ def __init__(self, model, solver, optsDict={}) -> None: evalPoints = [0] + [i/ns for i in range(1, ns+1)] - b0 = model.states[1] + z0 = ca.vertcat(model.states[1], model.states[2]) p0 = ca.vertcat(model.controls, model.parameters) - bf = self.eval(b0, p0, ca.hcat(evalPoints)) + zf = self.eval(z0, p0, ca.hcat(evalPoints)) # zf[0, idx]: velSq, zf[1, idx]: pos tApprox = model.states[0] for idx in range(ns): - vCurr = ca.sqrt(bf[idx]) - vNext = ca.sqrt(bf[idx+1]) - tApprox += 2*model.parameters[2]*(evalPoints[idx+1]-evalPoints[idx])/(vCurr + vNext) + vCurr = ca.sqrt(zf[0, idx]) + vNext = ca.sqrt(zf[0, idx+1]) + tApprox += 2*model.parameters[-1]*(evalPoints[idx+1]-evalPoints[idx])/(vCurr + vNext) - eval = ca.vertcat(tApprox, bf[-1]) + eval = ca.vertcat(tApprox, zf[0, ns], zf[1, ns]) self.eval = ca.Function('xNxt', [model.states, ca.vertcat(model.controls, model.parameters), ca.MX.sym('ds')], [eval]) - def solve(self, time, velocitySquared, ds, traction=0, pnBrake=0, gradient=0, curvature=0): + def solve(self, time, velocitySquared, ds, position=0, traction=0, pnBrake=0, gradient=0, gradientLinearTerm=0, curvature=0, curvatureLinearTerm=0, tunnelFactor=0): withPnBrake = self.model.withPnBrake @@ -352,14 +403,15 @@ def solve(self, time, velocitySquared, ds, traction=0, pnBrake=0, gradient=0, cu raise ValueError("Cannot define value for pneumatic braking when this brake is deactivated!") - x0 = ca.vertcat(time, velocitySquared) + x0 = ca.vertcat(time, velocitySquared, position) u0 = ca.vertcat(traction, pnBrake if withPnBrake else []) - p0 = ca.vertcat(gradient, curvature, ds) + p0 = ca.vertcat(gradient, gradientLinearTerm, curvature, curvatureLinearTerm, tunnelFactor, ds) x1 = self.eval(x0, ca.vertcat(u0, p0), 1) out = {} out['time'] = x1[0] out['velSquared'] = x1[1] + out['position'] = x1[2] return out @@ -371,16 +423,18 @@ def initLosses(self, lossesTrFun, lossesRgbFun, totalMass, solver='CVODES'): mdl = self.model vel = ca.MX.sym('v') + pos = ca.MX.sym('pos') velDot = ca.substitute(mdl.acceleration, mdl.states[1], vel**2) + velDot = ca.substitute(velDot, mdl.states[2], pos) energyTrDot = lossesTrFun(mdl.controls[0]*totalMass, vel)/totalMass # tractive energy energyBrDot = lossesRgbFun(mdl.controls[0]*totalMass, vel)/totalMass # braking energy dt = ca.MX.sym('dt') - x = ca.vertcat(vel, ca.MX.sym('eTr'), ca.MX.sym('eBr')) - p = ca.vertcat(mdl.controls, mdl.parameters[0], mdl.parameters[1], dt) - xdot = ca.vertcat(velDot, energyTrDot, energyBrDot) + x = ca.vertcat(vel, pos, ca.MX.sym('eTr'), ca.MX.sym('eBr')) + p = ca.vertcat(mdl.controls, mdl.parameters[0], mdl.parameters[1], mdl.parameters[2], mdl.parameters[3], mdl.parameters[4], dt) + xdot = ca.vertcat(velDot, vel, energyTrDot, energyBrDot) if solver == 'RK': @@ -402,13 +456,13 @@ def initLosses(self, lossesTrFun, lossesRgbFun, totalMass, solver='CVODES'): raise ValueError("Unknown solver!") - def calcLosses(self, velocity, dt, traction=0, pnBrake=0, gradient=0, curvature=0): + def calcLosses(self, velocity, dt, position=0, traction=0, pnBrake=0, gradient=0, gradientLinearTerm=0, curvature=0, curvatureLinearTerm=0, tunnelFactor=0): mdl = self.model - out = self.lossesIntegrator(ca.vertcat(velocity, 0, 0), ca.vertcat(traction, pnBrake if mdl.withPnBrake else [], gradient, curvature), dt) + out = self.lossesIntegrator(ca.vertcat(velocity, position, 0, 0), ca.vertcat(traction, pnBrake if mdl.withPnBrake else [], gradient, gradientLinearTerm, curvature, curvatureLinearTerm, tunnelFactor), dt) - lossesTr, lossesRgb = out[1], out[2] + lossesTr, lossesRgb = out[2], out[3] return lossesTr, lossesRgb @@ -418,11 +472,12 @@ def initRollingResistance(self, solver='CVODES'): mdl = self.model bDot = mdl.ode[1] - eDot = mdl.rollingResistance*mdl.parameters[2] + sDot = mdl.ode[2] + eDot = mdl.rollingResistance*mdl.parameters[-1] - x = ca.vertcat(mdl.states[1], ca.MX.sym('e')) + x = ca.vertcat(mdl.states[1], mdl.states[2], ca.MX.sym('e')) p = ca.vertcat(mdl.controls, mdl.parameters) - xdot = ca.vertcat(bDot, eDot) + xdot = ca.vertcat(bDot, sDot, eDot) if solver == 'RK': @@ -442,14 +497,14 @@ def initRollingResistance(self, solver='CVODES'): raise ValueError("Unknown solver!") - def calcRollingResistance(self, velocity, ds, traction=0, pnBrake=0, gradient=0, curvature=0): + def calcRollingResistance(self, velocity, ds, position=0, traction=0, pnBrake=0, gradient=0, gradientLinearTerm=0, curvature=0, curvatureLinearTerm=0, tunnelFactor=0): mdl = self.model - out = self.rollingResistanceIntegrator(ca.vertcat(velocity**2, 0), ca.vertcat(traction, pnBrake if mdl.withPnBrake else [], gradient, - curvature, ds), 1) + out = self.rollingResistanceIntegrator(ca.vertcat(velocity**2, position, 0), ca.vertcat(traction, pnBrake if mdl.withPnBrake else [], + gradient, gradientLinearTerm, curvature, curvatureLinearTerm, tunnelFactor, ds), 1) - losses = out[1] + losses = out[2] return losses, ca.sqrt(out[0]) @@ -548,6 +603,6 @@ def checkValues(self): trainSpecs = Train(config={'id':'NL_intercity_VIRM6'}) integrator = TrainIntegrator(trainSpecs.exportModel(), 'RK', optsDict={'numApproxSteps':2}) - solution = integrator.solve(t0, v0**2, ds, f0, gradient=gd, curvature=cr) + solution = integrator.solve(t0, v0**2, ds, traction=f0, gradient=gd, curvature=cr) print(solution) diff --git a/mseetc/utils.py b/mseetc/utils.py index 978adc1..d23e0d0 100644 --- a/mseetc/utils.py +++ b/mseetc/utils.py @@ -7,6 +7,9 @@ from types import MethodType from distutils.spawn import find_executable +from matplotlib import pyplot as plt + + def var(tag, dim=None): "Wrapper to create symbolic variables in casadi." @@ -369,7 +372,7 @@ def convertUnit(value, unit): Convert from any known unit to internally used unit. """ - if unit in {'m', 'm/s', 'permil', 'kg', 'W', 'N', 'm/s^2', '-', 'N/(m/s)', 'N/(m/s)^2', 'kg/m'}: + if unit in {'m', 'm/s', 'm^2', 'permil', 'kg', 'W', 'N', 'm/s^2', '-', 'N/(m/s)', 'N/(m/s)^2', 'kg/m'}: valueOut = value @@ -479,6 +482,152 @@ def latexify(): return latexFound +def computeTunnelFactor(cross_section, train, opts): + + if cross_section == float("inf"): + + return 0 # no tunnel + + total_mass = train.mass * train.rho + tunnelCoefficients = train.tunnelCoefficients + + if not tunnelCoefficients: + raise ValueError("Tunnel cross section was specified, but train has no tunnel resistance data.") + + availableCrossSections = list(tunnelCoefficients.keys()) + + if opts.chooseClosestTunnelCrossSection: + + closestCrossSection = min(availableCrossSections, key=lambda cs: abs(cs - cross_section)) + tunnelCoefficient = tunnelCoefficients[closestCrossSection] + + else: + + if cross_section not in tunnelCoefficients: + + raise ValueError( + "Tunnel resistance coefficient for cross section {} m^2 is not available. " + "Available cross sections are: {}. " + "Set chooseClosestTunnelCrossSection=True to use the closest available value.".format( + cross_section, + availableCrossSections + ) + ) + + tunnelCoefficient = tunnelCoefficients[cross_section] + + + 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 picked_points + + m *= 2 + + +def plotSpeedLimits(track, pos_adj, v_adj): + + v_limits = track.speedLimits["Speed limit [m/s]"].to_numpy(dtype=float) + pos_v_limits = track.speedLimits.index.to_numpy(dtype=float) + + v_limits = np.append(v_limits, v_limits[-1]) + pos_v_limits = np.append(pos_v_limits, track.length) + + v_adj = np.append(v_adj, v_adj[-1]) + pos_adj = np.append(pos_adj, track.length) + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.step(pos_v_limits / 1000, v_limits * 3.6, where='post', label="piecewise constant speed limit") + ax.step(pos_adj/1000, v_adj * 3.6, where='post', linestyle="--", label="train length adjusted speed limit") + + ax.set_title("Speed Limits") + 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, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + +def plotGradients(track, pos_adj, g_adj, g_linear): + + x = track.gradients.index.to_numpy(dtype=float) + g = track.gradients["Gradient [permil]"].to_numpy(dtype=float) + + x = np.append(x, track.length) + g = np.append(g, g[-1]) + + pos_adj = np.append(pos_adj, track.length) + g_adj = np.append(g_adj, g_adj[-1]) + g_linear = np.append(g_linear, g_linear[-1]) + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.step(x/1000, g, where='post', label="piecewise constant gradients") + ax.plot(pos_adj / 1000, g_adj, linestyle="--", label="train length averaged gradients") + + g_computed = np.r_[g_adj[0], g_adj[:-1] + g_linear[:-1] * (pos_adj[1:] - pos_adj[:-1])] + ax.scatter(pos_adj / 1000, g_computed, marker="o", label="computed gradients") + + ax.set_title("Gradients") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Gradient [‰]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + +def plotCurvatures(track, pos_adj, c_adj, c_linear): + + x = track.curvatures.index.to_numpy(dtype=float) + c = track.curvatures["Curvature [1/m]"].to_numpy(dtype=float) + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.step(x/1000, c, where='post', label="piecewise constant curvatures") + ax.plot(pos_adj / 1000, c_adj, linestyle="--", label="train length averaged curvatures") + + c_computed = np.r_[c_adj[0], c_adj[:-1] + c_linear[:-1] * (pos_adj[1:] - pos_adj[:-1])] + ax.scatter(pos_adj / 1000, c_computed, marker="o", label="computed curvatures") + + ax.set_title("Curvatures") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Curvature [1/m]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + if __name__ == '__main__': pass diff --git a/simulations/sim_launcher.py b/simulations/sim_launcher.py new file mode 100644 index 0000000..515caaa --- /dev/null +++ b/simulations/sim_launcher.py @@ -0,0 +1,62 @@ +from mseetc.efficiency import totalLossesFunction +from mseetc.ocp import casadiSolver + + +def get_power_loss_function(train, mode="perfect",* ,auxiliaries: float = 27_000, eta_gear: float = 0.96): + + if mode == "perfect": + + return lambda f, v: 0 + + elif mode == "static": + + return lambda f, v: (f>0)*f*v*(1-train.etaTraction)/train.etaTraction - (f<0)*f*v*(1-train.etaRgBrake) + + elif mode == "dynamic": + + return totalLossesFunction(train, auxiliaries=auxiliaries, etaGear=eta_gear) + + else: + + raise ValueError("mode must be one of: 'perfect', 'static', 'dynamic'") + + + +if __name__ == '__main__': + + from mseetc.train import Train + from mseetc.track import Track + + # Timetable + startPosition = 0 # [m] + endPosition = 20000 # [m] + duration = 60*20 # [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.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: + + 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!") \ No newline at end of file diff --git a/tracks/CH_StGallen_Wil.json b/tracks/CH_StGallen_Wil.json index 5d4199c..c2afc63 100644 --- a/tracks/CH_StGallen_Wil.json +++ b/tracks/CH_StGallen_Wil.json @@ -1,8 +1,8 @@ { "metadata": { "id": "CH_StGallen_Wil", - "created by": "Dimitris Kouzoupis", - "library version": "TTOBench v1.1", + "created by": "Dimitris Kouzoupis, Juxhino Kavaja", + "library version": "TTOBench v1.2", "license": "BSD 2-Clause License" }, "altitude": { @@ -695,5 +695,1204 @@ -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 + ] + ] } } \ No newline at end of file diff --git a/tracks/swisstopo/analyzeTracks.py b/tracks/swisstopo/analyzeTracks.py new file mode 100644 index 0000000..250bea5 --- /dev/null +++ b/tracks/swisstopo/analyzeTracks.py @@ -0,0 +1,117 @@ +import numpy as np +from matplotlib import pyplot as plt + +from mseetc.ocp import casadiSolver +from simulations.sim_launcher import get_power_loss_function +from mseetc.track import Track +from mseetc.train import Train + + +if __name__ == '__main__': + + + SBB_track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='') + + SBB_positions = SBB_track.gradients.index.values + SBB_gradients = SBB_track.gradients["Gradient [permil]"].to_numpy() + + initial_altitude = SBB_track.altitude + delta_s = np.diff(SBB_positions) + delta_h = SBB_gradients[:-1] / 1000 * delta_s + + SBB_altitude = np.insert(initial_altitude + np.cumsum(delta_h),0, initial_altitude) + + Topo_track = Track(config={'id': 'CH_StGallen_Wil_Swisstopo'}, pathJSON='') + + Topo_positions = Topo_track.gradients.index.values + Topo_gradients = Topo_track.gradients["Gradient [permil]"].to_numpy() + + initial_altitude = Topo_track.altitude + delta_s = np.diff(Topo_positions) + delta_h = Topo_gradients[:-1] / 1000 * delta_s + + Topo_altitude = np.insert(initial_altitude + np.cumsum(delta_h),0, initial_altitude) + + shift = 800 # 770 + + + ### Plot Altitude Comparison + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.plot(SBB_positions / 1000, SBB_altitude, label="SBB") + ax.plot((Topo_positions - shift) / 1000, Topo_altitude, label="Topo") + ax.set_title("Altitude Comparison") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Altitude [m]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, SBB_track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + + ### Plot Speed Limit Comparison + + 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.set_title("Altitude Comparison") + ax2.set_xlabel("Position [km]") + ax2.set_ylabel("Velocity [km/h]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + ax2.set_xlim(0, SBB_track.length / 1000) + ax2.figure.tight_layout() + + plt.show() + + + ### Compute Energy Comparison + + # Timetable + startPosition = 0 # [m] + endPosition = 23000 # [m] + duration = 23000/(80/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") + opts = {'numIntervals': 1000, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 2}, 'energyOptimal': True} + + SBB_track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + SBB_track.updateTrainLengthDependentValues(train) + solver = casadiSolver(train, SBB_track, opts) + dfSBB, statsSBB = solver.solve(duration) + + Topo_track.updateLimits(positionStart=startPosition + shift, positionEnd=endPosition + shift, unit='m') + Topo_track.updateTrainLengthDependentValues(train) + solver = casadiSolver(train, Topo_track, opts) + dfTopo, statsTopo = solver.solve(duration) + + print(f"Cost SBB: {statsSBB['Cost']:.2f}") + print(f"Cost Topo: {statsTopo['Cost']:.2f}") + + print(f"{abs(statsSBB['Cost'] - statsTopo['Cost']) / statsSBB['Cost'] * 100:.2f}%") + + + ### Plot Trajectory + + fig3, ax3 = plt.subplots(figsize=(16, 8)) + + ax3.plot(dfSBB["Position [m]"] / 1000, dfSBB["Velocity [m/s]"] * 3.6, label="SBB") + ax3.plot(dfTopo["Position [m]"] / 1000, dfTopo["Velocity [m/s]"] * 3.6, label="Topo") + ax3.set_title("Speed Profile Comparison") + ax3.set_xlabel("Position [km]") + ax3.set_ylabel("Velocity [km/h]") + ax3.grid(True, which="both", linestyle="--", alpha=0.5) + ax3.legend(loc="upper right") + ax3.set_xlim(0, dfSBB["Position [m]"].max() / 1000) + ax3.figure.tight_layout() + + plt.show() + + diff --git a/tracks/swisstopo/csvImporter.py b/tracks/swisstopo/csvImporter.py new file mode 100644 index 0000000..f448c28 --- /dev/null +++ b/tracks/swisstopo/csvImporter.py @@ -0,0 +1,108 @@ +import json +from pathlib import Path + +import numpy as np +import pandas as pd + +if __name__ == '__main__': + + + ### Read CSV + + filePath = r"C:\Users\rolan\Documents\ms-eetc-innocheque\tracks\swisstopo\Track_StGallen_Wil.csv" + + df = pd.read_csv(filePath,na_values=["", "null", ""]) + + print(df.head()) + print(df.dtypes) + + + ### Speed limits + + speedProfile = df.loc[ + df["V_max"].notna(), + ["Total_Distance", "V_max"] + ].copy() + + + ### Gradients + + totalDistance = df["Total_Distance"] + altitude = df["Altitude"] + + window_size = 7 + altitude = altitude.rolling(window=window_size, center=True, min_periods=1).mean() + + spacing = 50 # [m] + + positions = np.arange(0, totalDistance.max(), spacing) + altitude_interp = np.interp(positions, totalDistance, altitude) + + gradientPerMille = np.insert(1000 * np.diff(altitude_interp) / np.diff(positions),0,0 ) + gradientPerMille = np.round(gradientPerMille, 1) + + + ### Parse to Json + + track_id = "CH_StGallen_Wil_Swisstopo" + author = "Roland Staerk" + + name = "CH_StGallen_Wil_Swisstopo" + + output_dir = Path(r"C:\Users\rolan\Documents\ms-eetc-innocheque\tracks\swisstopo") + output_path = output_dir / f"{name}.json" + + stops = [ + 0.0, + float(totalDistance.iloc[-1]) + ] + + speed_limits = [ + [float(pos), float(vmax)] + for pos, vmax in zip( + speedProfile["Total_Distance"], + speedProfile["V_max"] + ) + ] + + gradients = [ + [float(pos), float(grad)] + for pos, grad in zip(positions, gradientPerMille) + ] + + track_data = { + "metadata": { + "id": track_id, + "created by": author, + "library version": "TTOBench v1.4", + "license": "BSD 2-Clause License" + }, + "altitude": { + "unit": "m", + "value": float(altitude_interp[0]) + }, + "stops": { + "unit": "m", + "values": stops + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": speed_limits + }, + "gradients": { + "units": { + "position": "m", + "slope": "permil" + }, + "values": gradients + } + } + + + ### Save Json + + with open(output_path, "w", encoding="utf-8") as f: + json.dump(track_data, f, indent=4) \ No newline at end of file diff --git a/tracks/test_flat_no_tunnel.json b/tracks/test_flat_no_tunnel.json new file mode 100644 index 0000000..1e47af0 --- /dev/null +++ b/tracks/test_flat_no_tunnel.json @@ -0,0 +1,31 @@ +{ + "metadata": { + "id": "test_flat_no_tunnel", + "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 + ] + ] + } +} \ No newline at end of file diff --git a/tracks/test_flat_with_tunnel.json b/tracks/test_flat_with_tunnel.json new file mode 100644 index 0000000..1a40288 --- /dev/null +++ b/tracks/test_flat_with_tunnel.json @@ -0,0 +1,45 @@ +{ + "metadata": { + "id": "test_flat_with_tunnel", + "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 + ] + ] + }, + "tunnels": { + "units": { + "position": "m", + "length": "m", + "cross section": "m^2" + }, + "values": [ + [ + 1000.0, + 26000.0, + 24.0 + ] + ] + } +} \ No newline at end of file diff --git a/tracks/test_one_hill.json b/tracks/test_one_hill.json new file mode 100644 index 0000000..42f751a --- /dev/null +++ b/tracks/test_one_hill.json @@ -0,0 +1,59 @@ +{ + "metadata": { + "id": "test_one_hill", + "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, + 12000.0 + ] + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": [ + [ + 0.0, + 79.2 + ] + ] + }, + "gradients": { + "units": { + "position": "m", + "slope": "permil" + }, + "values": [ + [ + 0.0, + 0.0 + ], + [ + 1000.0, + 20.0 + ], + [ + 2000.0, + 0.0 + ], + [ + 3000.0, + -20.0 + ], + [ + 4000.0, + 0.0 + ] + ] + } +} \ No newline at end of file diff --git a/tracks/test_one_speed_increase.json b/tracks/test_one_speed_increase.json new file mode 100644 index 0000000..dbc4a5e --- /dev/null +++ b/tracks/test_one_speed_increase.json @@ -0,0 +1,35 @@ +{ + "metadata": { + "id": "test_one_speed_increase", + "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, + 12000.0 + ] + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": [ + [ + 0.0, + 79.2 + ], + [ + 1000.0, + 144 + ] + ] + } +} \ No newline at end of file diff --git a/tracks/test_two_radii.json b/tracks/test_two_radii.json new file mode 100644 index 0000000..27fd8b7 --- /dev/null +++ b/tracks/test_two_radii.json @@ -0,0 +1,65 @@ +{ + "metadata": { + "id": "test_two_radii", + "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, + 12000.0 + ] + }, + "speed limits": { + "units": { + "position": "m", + "velocity": "km/h" + }, + "values": [ + [ + 0.0, + 79.2 + ] + ] + }, + "curvatures": { + "units": { + "position": "m", + "radius at start": "m", + "radius at end": "m" + }, + "values": [ + [ + 0.0, + "infinity", + "infinity" + ], + [ + 1000.0, + 300.0, + 300.0 + ], + [ + 2000.0, + "infinity", + "infinity" + ], + [ + 3000.0, + -300.0, + -300.0 + ], + [ + 4000.0, + "infinity", + "infinity" + ] + ] + } +} \ No newline at end of file diff --git a/trains/CH_Stadler_FLIRT_TPF.json b/trains/CH_Stadler_FLIRT_TPF.json new file mode 100644 index 0000000..9443e2b --- /dev/null +++ b/trains/CH_Stadler_FLIRT_TPF.json @@ -0,0 +1,96 @@ +{ + "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 + ] + ] + } +} diff --git a/unitTests/curvatureResistance/curvatureResistance.py b/unitTests/curvatureResistance/curvatureResistance.py index cf6bbac..e7e8ec0 100644 --- a/unitTests/curvatureResistance/curvatureResistance.py +++ b/unitTests/curvatureResistance/curvatureResistance.py @@ -1,6 +1,5 @@ import copy import numpy as np -import sys import unittest from mseetc.efficiency import totalLossesFunction diff --git a/unitTests/integrators/integrators.py b/unitTests/integrators/integrators.py new file mode 100644 index 0000000..062aa64 --- /dev/null +++ b/unitTests/integrators/integrators.py @@ -0,0 +1,99 @@ +import unittest + +from mseetc.ocp import casadiSolver +from mseetc.track import Track +from mseetc.train import Train + + +class TestGradient(unittest.TestCase): + + def testAllIntegratorTypesWork(self): + ''' + Verify that all supported integration methods produce consistent results + for the same train, track, and optimization setup. + + The test compares RK, IRK, and CVODES, including the approximate time + integration option for RK and IRK. The resulting energy costs should only + 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.length = 600 + + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track.updateTrainLengthDependentValues(train) + + opts = {'numIntervals': numIntervals, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}} + solver = casadiSolver(train, track, opts) + df, stats = solver.solve(duration) + + energy_RK_Approx = stats['Cost'] + + opts = {'numIntervals': numIntervals, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 0}} + solver = casadiSolver(train, track, opts) + df, stats = solver.solve(duration) + + energy_RK = stats['Cost'] + + opts = {'numIntervals': numIntervals, 'integrationMethod': 'IRK', 'integrationOptions': {'numApproxSteps': 1}} + solver = casadiSolver(train, track, opts) + df, stats = solver.solve(duration) + + energy_IRK_Approx = stats['Cost'] + + opts = {'numIntervals': numIntervals, 'integrationMethod': 'IRK', 'integrationOptions': {'numApproxSteps': 0}} + solver = casadiSolver(train, track, opts) + df, stats = solver.solve(duration) + + energy_IRK = stats['Cost'] + + opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES'} + solver = casadiSolver(train, track, opts) + df, stats = solver.solve(duration) + + energy_CVODES = stats['Cost'] + + relDiff_RKApprox_IRKApprox = abs(energy_RK_Approx - energy_IRK_Approx) / energy_RK_Approx + relDiff_RKApprox_CVODES = abs(energy_RK_Approx - energy_CVODES) / energy_RK_Approx + relDiff_RK_IRK = abs(energy_RK - energy_IRK) / energy_RK + + self.assertLess( + relDiff_RKApprox_IRKApprox, + tol, + msg=( + "RK and IRK with numApproxSteps=1 should give similar costs. " + f"RK approx: {energy_RK_Approx:.6f}, " + f"IRK approx: {energy_IRK_Approx:.6f}, " + f"relative difference: {relDiff_RKApprox_IRKApprox:.6f}." + ) + ) + + self.assertLess( + relDiff_RKApprox_CVODES, + tol, + msg=( + "RK with numApproxSteps=1 and CVODES should give similar costs. " + f"RK approx: {energy_RK_Approx:.6f}, " + f"CVODES: {energy_CVODES:.6f}, " + f"relative difference: {relDiff_RKApprox_CVODES:.6f}." + ) + ) + + self.assertLess( + relDiff_RK_IRK, + tol, + msg=( + "RK and IRK with numApproxSteps=0 should give similar costs. " + f"RK: {energy_RK:.6f}, " + f"IRK: {energy_IRK:.6f}, " + f"relative difference: {relDiff_RK_IRK:.6f}." + ) + ) \ No newline at end of file diff --git a/unitTests/trainLengthDependentTrackAttributes/curvature.py b/unitTests/trainLengthDependentTrackAttributes/curvature.py new file mode 100644 index 0000000..1b3bba7 --- /dev/null +++ b/unitTests/trainLengthDependentTrackAttributes/curvature.py @@ -0,0 +1,655 @@ +import unittest + +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt + +from mseetc.ocp import casadiSolver, OptionsCasadiSolver +from mseetc.track import Track +from mseetc.train import Train, TrainIntegrator + + +plotDebug = False + + +def computeHeadingFromCurvature(curvatures, trackLength): + + positions = curvatures.index.to_numpy(dtype=float) + + if positions[-1] < trackLength: + positions = np.r_[positions, trackLength] + + curvature = curvatures["Curvature [1/m]"].to_numpy(dtype=float) + + if "Curvature linear term [1/m^2]" in curvatures.columns: + + curvatureLinear = curvatures["Curvature linear term [1/m^2]"].to_numpy(dtype=float) + + else: + + curvatureLinear = np.zeros(len(curvature)) + + ds = positions[1:] - positions[:-1] + + heading = np.zeros(len(positions)) + heading[1:] = np.cumsum(curvature * ds + 0.5 * curvatureLinear * ds**2) + + return pd.DataFrame( + {"Heading [rad]": heading}, + index=positions + ) + + +class TestCurvature(unittest.TestCase): + + def test_cvodes_pwc_midpoint_curvature_converges_to_pwl_curvature(self): + ''' + Track with a linearly increasing curvature over 1000 m. + + The result obtained using the piecewise linear curvature model is compared + against a piecewise constant midpoint approximation of the curvature with increasing numbers of intervals. + + The piecewise constant approximation should converge to the piecewise linear + result for both duration and final velocity. + + CVODES is used as the integrator. + ''' + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + + optsDict = {'integrationMethod': 'CVODES'} + opts = OptionsCasadiSolver(optsDict) + + trainModel = train.exportModel() + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + # Scenario + time0 = 0 + velSq0 = 1 + ds = 1000 + traction = 0.8 * (train.forceMax / train.mass) + + initialCurvature = 0 + finalCurvature = 0.005 + + maxIntervals = 50 + relativeTolerance = 1e-4 + + # PWL curvature reference + out = trainIntegrator.solve( + time=time0, + velocitySquared=velSq0, + ds=ds, + traction=traction, + pnBrake=0, + gradient=0, + gradientLinearTerm=0, + curvature=initialCurvature, + curvatureLinearTerm=(finalCurvature - initialCurvature) / ds, + tunnelFactor=0 + ) + + pwlDuration = float(out['time']) + pwlVelocity = np.sqrt(float(out['velSquared'])) + + # PWC curvature using midpoint rule + times = [] + velocities = [] + intervalCounts = [] + + for numIntervals in range(1, maxIntervals + 1): + + time = time0 + velSq = velSq0 + + for idx in range(numIntervals): + + curvature = (initialCurvature + (idx + 0.5) * (finalCurvature - initialCurvature) / numIntervals) + + out = trainIntegrator.solve( + time=time, + velocitySquared=velSq, + ds=ds / numIntervals, + traction=traction, + pnBrake=0, + gradient=0, + gradientLinearTerm=0, + curvature=curvature, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + time = out['time'] + velSq = out['velSquared'] + + intervalCounts.append(numIntervals) + times.append(float(time)) + velocities.append(np.sqrt(float(velSq))) + + finalPwcDuration = times[-1] + finalPwcVelocity = velocities[-1] + + relativeDurationError = abs(finalPwcDuration - pwlDuration) / pwlDuration + relativeVelocityError = abs(finalPwcVelocity - pwlVelocity) / pwlVelocity + + self.assertLess( + relativeDurationError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently to the PWL duration. " + f"PWL duration: {pwlDuration:.6f}, " + f"PWC duration: {finalPwcDuration:.6f}, " + f"relative error: {relativeDurationError:.6e}." + ) + ) + + self.assertLess( + relativeVelocityError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently to the PWL final velocity. " + f"PWL velocity: {pwlVelocity:.6f}, " + f"PWC velocity: {finalPwcVelocity:.6f}, " + f"relative error: {relativeVelocityError:.6e}." + ) + ) + + if plotDebug: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10)) + + fig.suptitle("CVODES: PWC midpoint curvature approximation compared to PWL curvature", fontsize=14) + + ax1.axhline(pwlDuration, label="pwl") + ax1.plot(intervalCounts, times, marker="o", color="orange", label="pwc midpoint") + ax1.set_xlabel("Number of intervals of pwc curvature approximation") + ax1.set_ylabel("Duration [s]") + ax1.grid(True, which="both", linestyle="--", alpha=0.5) + ax1.legend(loc="upper right") + + ax2.axhline(pwlVelocity, label="pwl") + ax2.plot(intervalCounts, velocities, marker="o", color="orange", label="pwc midpoint") + ax2.set_xlabel("Number of intervals of pwc curvature approximation") + ax2.set_ylabel("Velocity [m/s]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + + fig.tight_layout() + plt.show() + + + def test_rk_pwc_midpoint_curvature_converges_to_pwl_curvature(self): + ''' + Track with a linearly increasing curvature over 1000 m. + + The result obtained using the piecewise linear curvature model is compared + against a piecewise constant midpoint approximation of the curvature with increasing numbers of intervals. + + The piecewise constant approximation should converge to the piecewise linear + result for both duration and final velocity. + + RK without time approximation is used as the integrator. + RK substeps are increased until convergence + ''' + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + + optsDict = {'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps': 0, 'numSteps': 1}} + opts = OptionsCasadiSolver(optsDict) + + trainModel = train.exportModel() + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + # Scenario + time0 = 0 + velSq0 = 1 + ds = 1000 + traction = 0.8 * (train.forceMax / train.mass) + + initialCurvature = 0 + finalCurvature = 0.005 + + maxIntervals = 50 + relativeTolerance = 1e-5 + + # PWL gradient reference + pwlTimes = [] + pwlVelocities = [] + pwlIntervalCounts = [] + + for numStep in range(1, maxIntervals + 1): + + optsDict = {'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 0, 'numSteps': numStep}} + opts = OptionsCasadiSolver(optsDict) + pwlTrainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + out = pwlTrainIntegrator.solve( + time=time0, + velocitySquared=velSq0, + ds=ds, + traction=traction, + pnBrake=0, + gradient=0, + gradientLinearTerm=0, + curvature=initialCurvature, + curvatureLinearTerm=(finalCurvature - initialCurvature) / ds, + tunnelFactor=0 + ) + + pwlTimes.append(float(out['time'])) + pwlVelocities.append(np.sqrt(float(out['velSquared']))) + pwlIntervalCounts.append(numStep) + + finalPwlDuration = pwlTimes[-1] + finalPwlVelocity = pwlVelocities[-1] + + # PWC curvature using midpoint rule + pwcTimes = [] + pwcVelocities = [] + pwcIntervalCounts = [] + + for numIntervals in range(1, maxIntervals + 1): + + time = time0 + velSq = velSq0 + + for idx in range(numIntervals): + + curvature = (initialCurvature + (idx + 0.5) * (finalCurvature - initialCurvature) / numIntervals) + + out = trainIntegrator.solve( + time=time, + velocitySquared=velSq, + ds=ds / numIntervals, + traction=traction, + pnBrake=0, + gradient=0, + gradientLinearTerm=0, + curvature=curvature, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + time = out['time'] + velSq = out['velSquared'] + + pwcIntervalCounts.append(numIntervals) + pwcTimes.append(float(time)) + pwcVelocities.append(np.sqrt(float(velSq))) + + finalPwcDuration = pwcTimes[-1] + finalPwcVelocity = pwcVelocities[-1] + + relativeDurationError = abs(finalPwcDuration - finalPwlDuration) / finalPwlDuration + relativeVelocityError = abs(finalPwcVelocity - finalPwlVelocity) / finalPwlVelocity + + self.assertLess( + relativeDurationError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently to the PWL duration. " + f"PWL duration: {finalPwlDuration:.6f}, " + f"PWC duration: {finalPwcDuration:.6f}, " + f"relative error: {relativeDurationError:.6e}." + ) + ) + + self.assertLess( + relativeVelocityError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently to the PWL final velocity. " + f"PWL velocity: {finalPwlVelocity:.6f}, " + f"PWC velocity: {finalPwcVelocity:.6f}, " + f"relative error: {relativeVelocityError:.6e}." + ) + ) + + if plotDebug: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10)) + + fig.suptitle("Explicit RK: PWL curvature with RK substeps vs PWC midpoint approximation", fontsize=14) + + ax1.plot(pwlIntervalCounts, pwlTimes, marker="o", label="PWL: RK substeps") + ax1.plot(pwcIntervalCounts, pwcTimes, marker="o", color="orange", label="PWC: intervals") + ax1.set_xlabel("Refinement level: PWC intervals and RK substeps") + ax1.set_ylabel("Duration [s]") + ax1.grid(True, which="both", linestyle="--", alpha=0.5) + ax1.legend(loc="upper right") + + ax2.plot(pwlIntervalCounts, pwlVelocities, marker="o", label="PWL: RK substeps") + ax2.plot(pwcIntervalCounts, pwcVelocities, marker="o", color="orange", label="PWC: intervals") + ax2.set_xlabel("Refinement level: PWC intervals and RK substeps") + ax2.set_ylabel("Velocity [m/s]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + + fig.tight_layout() + plt.show() + + + def test_rk_time_approx_pwc_midpoint_curvature_converges_to_pwl_curvature(self): + ''' + Track with a linearly increasing curvature over 1000 m. + + The result obtained using the piecewise linear curvature model is compared + against a piecewise constant midpoint approximation of the curvature. + + The piecewise constant approximation should converge to the piecewise linear + result for both duration and final velocity. + + RK with time approximation is used as the integrator. + RK uses 50 substeps. + Time approx steps are increased until convergence + ''' + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + trainModel = train.exportModel() + + # Scenario + time0 = 0 + velSq0 = 1 + ds = 1000 + traction = 0.8 * (train.forceMax / train.mass) + + initialCurvature = 0 + finalCurvature = 0.005 + + numIntervals = 50 + timeApproxSteps = 30 + relativeTolerance = 1e-3 + + # PWL gradient reference + pwlTimes = [] + pwlVelocities = [] + timeApproxStepCounts = [] + + for timeSteps in range(1, timeApproxSteps + 1): + + optsDict = {'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': timeSteps, 'numSteps': 50}} + opts = OptionsCasadiSolver(optsDict) + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + out = trainIntegrator.solve( + time=time0, + velocitySquared=velSq0, + ds=ds, + traction=traction, + pnBrake=0, + gradient=0, + gradientLinearTerm=0, + curvature=initialCurvature, + curvatureLinearTerm=(finalCurvature - initialCurvature) / ds, + tunnelFactor=0 + ) + + pwlTimes.append(float(out['time'])) + pwlVelocities.append(np.sqrt(float(out['velSquared']))) + timeApproxStepCounts.append(timeSteps) + + finalPwlDuration = pwlTimes[-1] + finalPwlVelocity = pwlVelocities[-1] + + # PWC curvature using midpoint rule + pwcTimes = [] + pwcVelocities = [] + + for timeSteps in range(1, timeApproxSteps + 1): + + optsDict = {'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': timeSteps, 'numSteps': 50}} + opts = OptionsCasadiSolver(optsDict) + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + time = time0 + velSq = velSq0 + + for idx in range(numIntervals): + + curvature = (initialCurvature + (idx + 0.5) * (finalCurvature - initialCurvature) / numIntervals) + + out = trainIntegrator.solve( + time=time, + velocitySquared=velSq, + ds=ds / numIntervals, + traction=traction, + pnBrake=0, + gradient=0, + gradientLinearTerm=0, + curvature=curvature, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + time = out['time'] + velSq = out['velSquared'] + + pwcTimes.append(float(time)) + pwcVelocities.append(np.sqrt(float(velSq))) + + finalPwcDuration = pwcTimes[-1] + finalPwcVelocity = pwcVelocities[-1] + + relativeDurationError = abs(finalPwcDuration - finalPwlDuration) / finalPwlDuration + relativeVelocityError = abs(finalPwcVelocity - finalPwlVelocity) / finalPwlVelocity + + self.assertLess( + relativeDurationError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently to the PWL duration. " + f"PWL duration: {finalPwlDuration:.6f}, " + f"PWC duration: {finalPwcDuration:.6f}, " + f"relative error: {relativeDurationError:.6e}." + ) + ) + + self.assertLess( + relativeVelocityError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently to the PWL final velocity. " + f"PWL velocity: {finalPwlVelocity:.6f}, " + f"PWC velocity: {finalPwcVelocity:.6f}, " + f"relative error: {relativeVelocityError:.6e}." + ) + ) + + if plotDebug: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10)) + + fig.suptitle("RK with Time Approx: PWC midpoint curvature approximation compared to PWL curvature", fontsize=14) + + ax1.plot(timeApproxStepCounts , pwlTimes, marker="o", label="pwl") + ax1.plot(timeApproxStepCounts, pwcTimes, marker="o", color="orange", label="pwc midpoint") + ax1.set_xlabel("Number of time approx steps") + ax1.set_ylabel("Duration [s]") + ax1.grid(True, which="both", linestyle="--", alpha=0.5) + ax1.legend(loc="upper right") + + ax2.plot(timeApproxStepCounts , pwlVelocities, marker="o", label="pwl") + ax2.plot(timeApproxStepCounts, pwcVelocities, marker="o", color="orange", label="pwc midpoint") + ax2.set_xlabel("Number of time approx steps") + ax2.set_ylabel("Velocity [m/s]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + ax2.ticklabel_format(axis="y", style="plain", useOffset=False) + + fig.tight_layout() + plt.show() + + + def test_train_length_dependent_curvature_preserves_target_heading(self): + ''' + Compare the final heading of the original length-independent curvature profile + with the train-length-dependent piecewise linear curvature profile. + + Both profiles should start from the same heading and end at the same target heading. + ''' + + headingTolerance = 1e-6 + + trainLength = 800 # [m] + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='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] + track.curvatures.loc[track.length - trainLength] = {"Curvature [1/m]": 0.0} + + df_heading_indep = computeHeadingFromCurvature(track.curvatures, track.length) + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train.length = trainLength + track.updateTrainLengthDependentValues(train) + + df_heading_dep = computeHeadingFromCurvature(track.curvatures, track.length) + + self.assertAlmostEqual( + df_heading_indep["Heading [rad]"].iloc[0], + df_heading_dep["Heading [rad]"].iloc[0], + delta=headingTolerance, + msg="Length-independent and train-length-dependent curvature profiles should start with the same heading." + ) + + self.assertAlmostEqual( + df_heading_indep["Heading [rad]"].iloc[-1], + df_heading_dep["Heading [rad]"].iloc[-1], + delta=headingTolerance, + msg=( + "Length-independent and train-length-dependent curvature profiles " + "should end with the same heading. " + f"Length-independent final heading: {df_heading_indep['Heading [rad]'].iloc[-1]:.12f} rad, " + f"train-length-dependent final heading: {df_heading_dep['Heading [rad]'].iloc[-1]:.12f} rad." + ) + ) + + if plotDebug: + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.plot(df_heading_indep.index.values / 1000, df_heading_indep["Heading [rad]"].to_numpy(), label="length-independent heading") + ax.plot(df_heading_dep.index.values / 1000, df_heading_dep["Heading [rad]"].to_numpy(), label="length-dependent heading") + ax.set_title("Heading Comparison") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Heading [rad]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + + def test_pwc_curvature_approximation_matches_pwl_curvature_energy(self): + ''' + Track with right turn from 1000 m to 2000 m and a left turn from 3000 m to 4000 m. + + Energy consumption should be roughly equal if computed using piecewise + 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.length = 600 + + track = Track(config={'id': 'test_two_radii'}, pathJSON='tracks') + track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track.updateTrainLengthDependentValues(train) + + # PWL Curvatures + opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES'} + solver = casadiSolver(train, track, opts) + pwl_df, pwl_stats = solver.solve(duration) + + 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) + + energyConsumptionWithPwcTerms= pwc_stats['Cost'] + + relativeEnergyDifference = (abs(energyConsumptionWithLinearTerms - energyConsumptionWithPwcTerms) / energyConsumptionWithLinearTerms) + + self.assertLess( + relativeEnergyDifference, + energyRelativeTolerance, + msg=( + "Energy consumption differs too much between piecewise linear and " + f"piecewise constant curvatures. Relative difference: " + f"{relativeEnergyDifference:.6f}, tolerance: {energyRelativeTolerance:.6f}. " + f"PWL cost: {energyConsumptionWithLinearTerms:.6f}, " + f"PWC cost: {energyConsumptionWithPwcTerms:.6f}." + ) + ) + + if plotDebug: + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.plot(pwl_df["Position [m]"] / 1000, pwl_df["Curvature [1/m]"], label="pwl curvatures") + ax.step(pwc_df["Position [m]"] / 1000, pwc_df["Curvature [1/m]"], "--", where="post", label="pwc curvatures") + ax.set_title("Curvatures") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Curvature [1/m]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + + def test_short_train_length_has_negligible_effect_on_curvature_energy_consumption(self): + ''' + Use an artificially short train on a real track profile. + + For a very small train length, the train-length-dependent curvature profile + should be almost identical to the original curvature profile. Therefore, the + energy consumption should remain within a small relative tolerance. + ''' + + relativeTolerance = 1e-2 + trainLength = 10 # [m] + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train.length = trainLength + + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track.gradients = track.gradients.iloc[[0]] + track.gradients["Gradient [permil]"].iloc[0] = 0 + + # 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) + + track.updateTrainLengthDependentValues(train) + solver = casadiSolver(train, track, opts) + dep_df, dep_stats = solver.solve(duration) + + energyConsumptionIndependentOfTrainLength = indep_stats['Cost'] + energyConsumptionDependentOfTrainLength = dep_stats['Cost'] + + relativeDifference = (abs(energyConsumptionDependentOfTrainLength - energyConsumptionIndependentOfTrainLength) / energyConsumptionIndependentOfTrainLength) + + self.assertLess( + relativeDifference, + relativeTolerance, + msg=( + "Energy consumption with and without train-length-dependent curvature should be roughly equal. " + f"Independent: {energyConsumptionIndependentOfTrainLength:.6f}, " + f"dependent: {energyConsumptionDependentOfTrainLength:.6f}, " + f"relative difference: {relativeDifference:.6e}." + ) + ) \ No newline at end of file diff --git a/unitTests/trainLengthDependentTrackAttributes/gradient.py b/unitTests/trainLengthDependentTrackAttributes/gradient.py new file mode 100644 index 0000000..6fd5018 --- /dev/null +++ b/unitTests/trainLengthDependentTrackAttributes/gradient.py @@ -0,0 +1,663 @@ +import unittest + +import numpy as np +from matplotlib import pyplot as plt + +from mseetc.ocp import casadiSolver, OptionsCasadiSolver +from mseetc.track import Track, computeAltitude +from mseetc.train import Train, TrainIntegrator + + +plotDebug = False + + +class TestGradient(unittest.TestCase): + + def test_cvodes_pwc_midpoint_gradient_converges_to_pwl_gradient(self): + ''' + Track with a linearly increasing gradient over 1000 m. + + The result obtained using the piecewise linear gradient model is compared + against a piecewise constant midpoint approximation of the gradient with increasing numbers of intervals. + + The piecewise constant approximation should converge to the piecewise linear + result for both duration and final velocity. + + CVODES is used as the integrator. + ''' + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + + optsDict = {'integrationMethod': 'CVODES'} + opts = OptionsCasadiSolver(optsDict) + + trainModel = train.exportModel() + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + # Scenario + time0 = 0 + velSq0 = 1 + ds = 1000 + traction = 0.8 * (train.forceMax / train.mass) + + initialGradient = 0 + finalGradient = 0.07 + + maxIntervals = 50 + relativeTolerance = 1e-3 + + # PWL gradient reference + out = trainIntegrator.solve( + time=time0, + velocitySquared=velSq0, + ds=ds, + traction=traction, + pnBrake=0, + gradient=initialGradient, + gradientLinearTerm=(finalGradient - initialGradient) / ds, + curvature=0, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + pwlDuration = float(out['time']) + pwlVelocity = np.sqrt(float(out['velSquared'])) + + # PWC gradients using midpoint rule + times = [] + velocities = [] + intervalCounts = [] + + for numIntervals in range(1, maxIntervals + 1): + + time = time0 + velSq = velSq0 + + for idx in range(numIntervals): + + gradient = (initialGradient + (idx + 0.5) * (finalGradient - initialGradient) / numIntervals) + + out = trainIntegrator.solve( + time=time, + velocitySquared=velSq, + ds=ds / numIntervals, + traction=traction, + pnBrake=0, + gradient=gradient, + gradientLinearTerm=0, + curvature=0, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + time = out['time'] + velSq = out['velSquared'] + + intervalCounts.append(numIntervals) + times.append(float(time)) + velocities.append(np.sqrt(float(velSq))) + + finalPwcDuration = times[-1] + finalPwcVelocity = velocities[-1] + + relativeDurationError = abs(finalPwcDuration - pwlDuration) / pwlDuration + relativeVelocityError = abs(finalPwcVelocity - pwlVelocity) / pwlVelocity + + self.assertLess( + relativeDurationError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently " + "to the PWL duration. " + f"PWL duration: {pwlDuration:.6f}, " + f"PWC duration: {finalPwcDuration:.6f}, " + f"relative error: {relativeDurationError:.6e}." + ) + ) + + self.assertLess( + relativeVelocityError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently " + "to the PWL final velocity. " + f"PWL velocity: {pwlVelocity:.6f}, " + f"PWC velocity: {finalPwcVelocity:.6f}, " + f"relative error: {relativeVelocityError:.6e}." + ) + ) + + if plotDebug: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10)) + + fig.suptitle("CVODES: PWC midpoint gradient approximation compared to PWL gradient", fontsize=14) + + ax1.axhline(pwlDuration, label="pwl") + ax1.plot(intervalCounts, times, marker="o", color="orange", label="pwc midpoint") + ax1.set_xlabel("Number of intervals of pwc gradient approximation") + ax1.set_ylabel("Duration [s]") + ax1.grid(True, which="both", linestyle="--", alpha=0.5) + ax1.legend(loc="upper right") + + ax2.axhline(pwlVelocity, label="pwl") + ax2.plot(intervalCounts, velocities, marker="o", color="orange", label="pwc midpoint") + ax2.set_xlabel("Number of intervals of pwc gradient approximation") + ax2.set_ylabel("Velocity [m/s]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + + fig.tight_layout() + plt.show() + + + def test_rk_pwc_midpoint_gradient_converges_to_pwl_gradient(self): + ''' + Track with a linearly increasing gradient over 1000 m. + + The result obtained using the piecewise linear gradient model is compared + against a piecewise constant midpoint approximation of the gradient with increasing numbers of intervals. + + The piecewise constant approximation should converge to the piecewise linear + result for both duration and final velocity. + + RK without time approximation is used as the integrator. + RK substeps are increased until convergence + ''' + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + + optsDict = {'integrationMethod':'RK', 'integrationOptions':{'numApproxSteps': 0, 'numSteps': 1}} + opts = OptionsCasadiSolver(optsDict) + + trainModel = train.exportModel() + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + # Scenario + time0 = 0 + velSq0 = 1 + ds = 1000 + traction = 0.8 * (train.forceMax / train.mass) + + initialGradient = 0 + finalGradient = 0.07 + + maxIntervals = 50 + relativeTolerance = 1e-3 + + # PWL gradient reference + pwlTimes = [] + pwlVelocities = [] + pwlIntervalCounts = [] + + for numStep in range(1, maxIntervals + 1): + + optsDict = {'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 0, 'numSteps': numStep}} + opts = OptionsCasadiSolver(optsDict) + pwlTrainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + out = pwlTrainIntegrator.solve( + time=time0, + velocitySquared=velSq0, + ds=ds, + traction=traction, + pnBrake=0, + gradient=initialGradient, + gradientLinearTerm=(finalGradient - initialGradient) / ds, + curvature=0, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + pwlTimes.append(float(out['time'])) + pwlVelocities.append(np.sqrt(float(out['velSquared']))) + pwlIntervalCounts.append(numStep) + + finalPwlDuration = pwlTimes[-1] + finalPwlVelocity = pwlVelocities[-1] + + # PWC gradients using midpoint rule + pwcTimes = [] + pwcVelocities = [] + pwcIntervalCounts = [] + + for numIntervals in range(1, maxIntervals + 1): + + time = time0 + velSq = velSq0 + + for idx in range(numIntervals): + + gradient = (initialGradient + (idx + 0.5) * (finalGradient - initialGradient) / numIntervals) + + out = trainIntegrator.solve( + time=time, + velocitySquared=velSq, + ds=ds / numIntervals, + traction=traction, + pnBrake=0, + gradient=gradient, + gradientLinearTerm=0, + curvature=0, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + time = out['time'] + velSq = out['velSquared'] + + pwcIntervalCounts.append(numIntervals) + pwcTimes.append(float(time)) + pwcVelocities.append(np.sqrt(float(velSq))) + + finalPwcDuration = pwcTimes[-1] + finalPwcVelocity = pwcVelocities[-1] + + relativeDurationError = abs(finalPwcDuration - finalPwlDuration) / finalPwlDuration + relativeVelocityError = abs(finalPwcVelocity - finalPwlVelocity) / finalPwlVelocity + + self.assertLess( + relativeDurationError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently " + "to the PWL duration. " + f"PWL duration: {finalPwlDuration:.6f}, " + f"PWC duration: {finalPwcDuration:.6f}, " + f"relative error: {relativeDurationError:.6e}." + ) + ) + + self.assertLess( + relativeVelocityError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently " + "to the PWL final velocity. " + f"PWL velocity: {finalPwlVelocity:.6f}, " + f"PWC velocity: {finalPwcVelocity:.6f}, " + f"relative error: {relativeVelocityError:.6e}." + ) + ) + + if plotDebug: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10)) + + fig.suptitle("Explicit RK: PWL gradient with RK substeps vs PWC midpoint approximation", fontsize=14) + + ax1.plot(pwlIntervalCounts, pwlTimes, marker="o", label="PWL: RK substeps") + ax1.plot(pwcIntervalCounts, pwcTimes, marker="o", color="orange", label="PWC: intervals") + ax1.set_xlabel("Refinement level: PWC intervals and RK substeps") + ax1.set_ylabel("Duration [s]") + ax1.grid(True, which="both", linestyle="--", alpha=0.5) + ax1.legend(loc="upper right") + + ax2.plot(pwlIntervalCounts, pwlVelocities, marker="o", label="PWL: RK substeps") + ax2.plot(pwcIntervalCounts, pwcVelocities, marker="o", color="orange", label="PWC: intervals") + ax2.set_xlabel("Refinement level: PWC intervals and RK substeps") + ax2.set_ylabel("Velocity [m/s]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + + fig.tight_layout() + plt.show() + + + def test_rk_time_approx_pwc_midpoint_gradient_converges_to_pwl_gradient(self): + ''' + Track with a linearly increasing gradient over 1000 m. + + The result obtained using the piecewise linear gradient model is compared + against a piecewise constant midpoint approximation of the gradient. + + The piecewise constant approximation should converge to the piecewise linear + result for both duration and final velocity. + + RK with time approximation is used as the integrator. + RK uses 50 substeps. + Time approx steps are increased until convergence + ''' + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + trainModel = train.exportModel() + + # Scenario + time0 = 0 + velSq0 = 1 + ds = 1000 + traction = 0.8 * (train.forceMax / train.mass) + + initialGradient = 0 + finalGradient = 0.07 + + numIntervals = 50 + timeApproxSteps = 30 + relativeTolerance = 1e-3 + + # PWL gradient reference + pwlTimes = [] + pwlVelocities = [] + timeApproxStepCounts = [] + + for timeSteps in range(1, timeApproxSteps + 1): + + optsDict = {'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': timeSteps, 'numSteps': 50}} + opts = OptionsCasadiSolver(optsDict) + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + out = trainIntegrator.solve( + time=time0, + velocitySquared=velSq0, + ds=ds, + traction=traction, + pnBrake=0, + gradient=initialGradient, + gradientLinearTerm=(finalGradient - initialGradient) / ds, + curvature=0, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + pwlTimes.append(float(out['time'])) + pwlVelocities.append(np.sqrt(float(out['velSquared']))) + timeApproxStepCounts.append(timeSteps) + + finalPwlDuration = pwlTimes[-1] + finalPwlVelocity = pwlVelocities[-1] + + # PWC gradients using midpoint rule + pwcTimes = [] + pwcVelocities = [] + + for timeSteps in range(1, timeApproxSteps + 1): + + optsDict = {'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': timeSteps, 'numSteps': 50}} + opts = OptionsCasadiSolver(optsDict) + trainIntegrator = TrainIntegrator(trainModel, opts.integrationMethod, opts.integrationOptions.toDict()) + + time = time0 + velSq = velSq0 + + for idx in range(numIntervals): + + gradient = (initialGradient + (idx + 0.5) * (finalGradient - initialGradient) / numIntervals) + + out = trainIntegrator.solve( + time=time, + velocitySquared=velSq, + ds=ds / numIntervals, + traction=traction, + pnBrake=0, + gradient=gradient, + gradientLinearTerm=0, + curvature=0, + curvatureLinearTerm=0, + tunnelFactor=0 + ) + + time = out['time'] + velSq = out['velSquared'] + + pwcTimes.append(float(time)) + pwcVelocities.append(np.sqrt(float(velSq))) + + finalPwcDuration = pwcTimes[-1] + finalPwcVelocity = pwcVelocities[-1] + + relativeDurationError = abs(finalPwcDuration - finalPwlDuration) / finalPwlDuration + relativeVelocityError = abs(finalPwcVelocity - finalPwlVelocity) / finalPwlVelocity + + self.assertLess( + relativeDurationError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently " + "to the PWL duration. " + f"PWL duration: {finalPwlDuration:.6f}, " + f"PWC duration: {finalPwcDuration:.6f}, " + f"relative error: {relativeDurationError:.6e}." + ) + ) + + self.assertLess( + relativeVelocityError, + relativeTolerance, + msg=( + "PWC midpoint approximation did not converge sufficiently " + "to the PWL final velocity. " + f"PWL velocity: {finalPwlVelocity:.6f}, " + f"PWC velocity: {finalPwcVelocity:.6f}, " + f"relative error: {relativeVelocityError:.6e}." + ) + ) + + if plotDebug: + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10)) + + fig.suptitle("RK with Time Approx: PWC midpoint gradient approximation compared to PWL gradient", fontsize=14) + + ax1.plot(timeApproxStepCounts , pwlTimes, marker="o", label="pwl") + ax1.plot(timeApproxStepCounts, pwcTimes, marker="o", color="orange", label="pwc midpoint") + ax1.set_xlabel("Number of time approx steps") + ax1.set_ylabel("Duration [s]") + ax1.grid(True, which="both", linestyle="--", alpha=0.5) + ax1.legend(loc="upper right") + + ax2.plot(timeApproxStepCounts , pwlVelocities, marker="o", label="pwl") + ax2.plot(timeApproxStepCounts, pwcVelocities, marker="o", color="orange", label="pwc midpoint") + ax2.set_xlabel("Number of time approx steps") + ax2.set_ylabel("Velocity [m/s]") + ax2.grid(True, which="both", linestyle="--", alpha=0.5) + ax2.legend(loc="upper right") + ax2.ticklabel_format(axis="y", style="plain", useOffset=False) + + fig.tight_layout() + plt.show() + + + def test_train_length_dependent_gradient_preserves_target_altitude(self): + ''' + Compare the final altitude of the original length-independent gradient profile + with the train-length-dependent piecewise linear profile. + + Both profiles should start from the same altitude and end at the same target altitude. + ''' + + altitudeTolerance = 1e-6 + + trainLength = 800 # [m] + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='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] + track.gradients.loc[track.length - trainLength] = {"Gradient [permil]": 0.0} + + df_alt = computeAltitude(track.gradients, track.length) + + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train.length = trainLength + track.updateTrainLengthDependentValues(train) + + positions = track.gradients.index.to_numpy(dtype=float) + positions = np.r_[positions, track.length] + + gradient = track.gradients["Gradient [permil]"].to_numpy(dtype=float) + gradientLinear = track.gradients["Gradient linear term [permil/m]"].to_numpy(dtype=float) + + ds = positions[1:] - positions[:-1] + + pwlAltitude = np.zeros(len(positions)) + pwlAltitude[1:] = np.cumsum((gradient * ds + 0.5 * gradientLinear * ds ** 2) / 1000) + + self.assertAlmostEqual( + df_alt["Altitude [m]"].iloc[0], + pwlAltitude[0], + delta=altitudeTolerance, + msg="Length-independent and train-length-dependent altitude profiles should start at the same altitude." + ) + + lengthIndependentFinalAltitude = df_alt["Altitude [m]"].iloc[-1] + lengthDependentFinalAltitude = pwlAltitude[-1] + + self.assertAlmostEqual( + lengthIndependentFinalAltitude, + lengthDependentFinalAltitude, + delta=altitudeTolerance, + msg=( + "Length-independent and train-length-dependent altitude profiles " + "should end at the same altitude. " + f"Length-independent final altitude: {lengthIndependentFinalAltitude:.8f} m, " + f"train-length-dependent final altitude: {lengthDependentFinalAltitude:.8f} m." + ) + ) + + if plotDebug: + + fig, ax = plt.subplots(figsize=(16, 8)) + + ax.plot(df_alt.index.values / 1000, df_alt["Altitude [m]"].to_numpy(), label="length-independent altitude") + ax.plot(positions / 1000, pwlAltitude, label="length-dependent altitude") + ax.set_title("Altitude Comparison") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Altitude [m]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + + def test_pwc_gradient_approximation_matches_pwl_gradient_energy(self): + ''' + Track with 20 permil increase from 1000 m to 2000 m and 20 permil + decrease from 3000 m to 4000 m. + + Energy consumption should be roughly equal if computed using piecewise + linear gradients or equivalent piecewise constant gradients due to a high number of shooting intervals. + + 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') + train.length = 600 + + track = Track(config={'id': 'test_one_hill'}, pathJSON='tracks') + track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + track.updateTrainLengthDependentValues(train) + + # PWL Gradients + opts = {'numIntervals': numIntervals, 'integrationMethod': 'CVODES'} + solver = casadiSolver(train, track, opts) + pwl_df, pwl_stats = solver.solve(duration) + + 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) + + energyConsumptionWithPwcTerms= pwc_stats['Cost'] + + relativeEnergyDifference = (abs(energyConsumptionWithLinearTerms - energyConsumptionWithPwcTerms) / energyConsumptionWithLinearTerms) + + pwc_df_grads = pwc_df.set_index("Position [m]")[["Gradient [permil]"]] + pwc_df_alt = computeAltitude(pwc_df_grads, track.length) + finalAltitude = pwc_df_alt.iloc[-1]["Altitude [m]"] + + self.assertLess( + relativeEnergyDifference, + energyRelativeTolerance, + msg=( + "Energy consumption differs too much between piecewise linear and " + f"piecewise constant gradients. Relative difference: " + f"{relativeEnergyDifference:.6f}, tolerance: {energyRelativeTolerance:.6f}. " + f"PWL cost: {energyConsumptionWithLinearTerms:.6f}, " + f"PWC cost: {energyConsumptionWithPwcTerms:.6f}." + ) + ) + + self.assertLess( + abs(finalAltitude), + altitudeTolerance, + msg=( + "Final altitude should be close to 0 m. " + f"Final altitude: {finalAltitude:.8f} m, " + f"tolerance: {altitudeTolerance:.8f} m." + ) + ) + + if plotDebug: + + fig, ax = plt.subplots(figsize=(16, 8)) + + df_grads_1 = pwl_df.set_index("Position [m]")[["Gradient [permil]"]] + ax.plot(df_grads_1.index.values / 1000, df_grads_1["Gradient [permil]"],label="pwl gradients") + ax.step(pwc_df_grads.index.values / 1000, pwc_df_grads["Gradient [permil]"], '--', where='post', label="pwc gradients") + ax.set_title("Gradients") + ax.set_xlabel("Position [km]") + ax.set_ylabel("Gradient [‰]") + ax.grid(True, which="both", linestyle="--", alpha=0.5) + ax.legend(loc="upper right") + ax.set_xlim(0, track.length / 1000) + ax.figure.tight_layout() + + plt.show() + + + def test_short_train_length_has_negligible_effect_on_energy_consumption(self): + ''' + Use an artificially short train on a real track profile. + + For a very small train length, the train-length-dependent gradient profile + should be almost identical to the original gradient profile. Therefore, the + energy consumption should remain within a small relative tolerance. + ''' + + relativeTolerance = 0.01 + trainLength = 10 # [m] + train = Train(config={'id': 'CH_Stadler_Flirt_TPF'}, pathJSON='trains') + train.length = trainLength + + track = Track(config={'id': 'CH_StGallen_Wil'}, pathJSON='tracks') + track.curvatures = track.curvatures.iloc[[0]] + track.curvatures["Curvature [1/m]"].iloc[0] = 0 + + # 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) + + track.updateTrainLengthDependentValues(train) + solver = casadiSolver(train, track, opts) + dep_df, dep_stats = solver.solve(duration) + + energyConsumptionIndependentOfTrainLength = indep_stats['Cost'] + energyConsumptionDependentOfTrainLength = dep_stats['Cost'] + + relativeDifference = (abs(energyConsumptionDependentOfTrainLength - energyConsumptionIndependentOfTrainLength) / energyConsumptionIndependentOfTrainLength) + + self.assertLess( + relativeDifference, + relativeTolerance, + msg=( + "Energy consumption with and without train-length-dependent gradients should be roughly equal. " + f"Independent: {energyConsumptionIndependentOfTrainLength:.6f}, " + f"dependent: {energyConsumptionDependentOfTrainLength:.6f}, " + f"relative difference: {relativeDifference:.6e}." + ) + ) \ No newline at end of file diff --git a/unitTests/trainLengthDependentTrackAttributes/speedLimit.py b/unitTests/trainLengthDependentTrackAttributes/speedLimit.py new file mode 100644 index 0000000..8b99ffe --- /dev/null +++ b/unitTests/trainLengthDependentTrackAttributes/speedLimit.py @@ -0,0 +1,62 @@ +import unittest + +from mseetc.ocp import casadiSolver +from mseetc.track import Track +from mseetc.train import Train + + +class TestSpeedLimit(unittest.TestCase): + + def test_train_length_dependent_speed_limit_delays_acceleration(self): + ''' + Speed limit increases from 22 m/s to 40 m/s at position 1000 m. + + Without train-length-dependent values, the optimizer may accelerate too early. + With train-length-dependent values, the front of the train may only exceed + 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') + + track = Track(config={'id': 'test_one_speed_increase'}, pathJSON='tracks') + track.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + + opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} + solver = casadiSolver(train, track, opts) + dfWithoutTrainLength, statsWithoutTrainLength = solver.solve(duration) + + 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) + + speedWithTrainLengthAfterIncrease = dfWithTrainLength[dfWithTrainLength['Position [m]'] > speedIncreasePosition].iloc[0]['Velocity [m/s]'] + speedWithTrainLengthAfterTrainPassedIncrease = dfWithTrainLength[dfWithTrainLength['Position [m]'] > (speedIncreasePosition+train.length)].iloc[0]['Velocity [m/s]'] + + self.assertGreater( + speedWithoutTrainLengthAfterIncrease, + speedLimitBeforeIncrease, + msg="Without train-length-dependent values, the train should accelerate immediately after the speed limit increase." + ) + + self.assertLessEqual( + speedWithTrainLengthAfterIncrease, + speedLimitBeforeIncrease + speedTolerance, + msg="With train-length-dependent values, the train should not accelerate before the whole train has passed the speed limit increase." + ) + + self.assertGreater( + speedWithTrainLengthAfterTrainPassedIncrease, + speedLimitBeforeIncrease, + msg="With train-length-dependent values, the train should accelerate after the whole train has passed the speed limit increase." + ) \ No newline at end of file diff --git a/unitTests/tunnelResistance/tunnelResistance.py b/unitTests/tunnelResistance/tunnelResistance.py new file mode 100644 index 0000000..e3f600f --- /dev/null +++ b/unitTests/tunnelResistance/tunnelResistance.py @@ -0,0 +1,56 @@ +import unittest + +from mseetc.ocp import casadiSolver +from mseetc.track import Track +from mseetc.train import Train + + +class TestTunnelResistance(unittest.TestCase): + + 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') + + trackWithoutTunnel = Track(config={'id': 'test_flat_no_tunnel'}, pathJSON='tracks') + trackWithoutTunnel.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + + opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} + solver = casadiSolver(train, trackWithoutTunnel, opts) + dfWithoutTunnel, statsWithoutTunnel = solver.solve(duration) + + energyConsumptionWithoutTunnel = statsWithoutTunnel['Cost'] + + trackWithTunnel = Track(config={'id': 'test_flat_with_tunnel'}, pathJSON='tracks') + trackWithTunnel.updateLimits(positionStart=startPosition, positionEnd=endPosition, unit='m') + + opts = {'numIntervals': 300, 'integrationMethod': 'RK', 'integrationOptions': {'numApproxSteps': 1}, 'energyOptimal': True} + solver = casadiSolver(train, trackWithTunnel, opts) + dfWithTunnel, statsWithTunnel = solver.solve(duration) + + energyConsumptionWithTunnel = statsWithTunnel['Cost'] + + self.assertGreater( + energyConsumptionWithTunnel, + energyConsumptionWithoutTunnel, + msg="Energy consumption with tunnel should be higher than without tunnel." + ) + + self.assertGreater( + energyConsumptionWithTunnel / energyConsumptionWithoutTunnel, + minEnergyRatio, + msg=( + "Energy consumption with tunnel should be significantly higher. " + f"Expected ratio > {minEnergyRatio}, got " + f"{energyConsumptionWithTunnel / energyConsumptionWithoutTunnel:.3f}." + ) + + ) \ No newline at end of file