From 61f6cc817335fc6ab1991d5af6faffa3143da93c Mon Sep 17 00:00:00 2001 From: Erashin Date: Mon, 13 Jan 2025 10:18:55 +0100 Subject: [PATCH] core: implement etcs braking simulator Signed-off-by: Erashin --- core/envelope-sim/build.gradle | 1 + .../sncf/osrd/envelope/EnvelopePhysics.java | 3 +- .../osrd/envelope_sim/EnvelopeSimPath.java | 25 ++ .../sncf/osrd/envelope_sim/PhysicsPath.java | 5 +- .../envelope_sim/PhysicsRollingStock.java | 18 ++ .../envelope_sim/TrainPhysicsIntegrator.java | 92 ++++-- .../overlays/EnvelopeDeceleration.java | 15 +- .../sncf/osrd/envelope_sim/etcs/Constants.kt | 6 + .../envelope_sim/etcs/ETCSBrakingCurves.kt | 263 ++++++++++++++++++ .../envelope_sim/etcs/ETCSBrakingSimulator.kt | 51 ++-- .../pipelines/MaxSpeedEnvelope.kt | 7 +- .../TrainPhysicsIntegratorTest.java | 8 +- .../fr/sncf/osrd/envelope_sim/FlatPath.java | 5 + .../osrd/envelope_sim/SimpleRollingStock.java | 6 + core/osrd-geom/build.gradle | 1 + core/osrd-railjson/build.gradle | 5 +- .../rollingstock/RJSEtcsBrakeParams.java | 60 +++- .../osrd/standalone_sim/StandaloneSim.java | 2 +- .../java/fr/sncf/osrd/train/RollingStock.java | 13 +- .../standalone_sim/StandaloneSimulation.kt | 1 - 20 files changed, 517 insertions(+), 70 deletions(-) create mode 100644 core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/Constants.kt create mode 100644 core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingCurves.kt diff --git a/core/envelope-sim/build.gradle b/core/envelope-sim/build.gradle index 4aafebb47cf..5a7422c00d4 100644 --- a/core/envelope-sim/build.gradle +++ b/core/envelope-sim/build.gradle @@ -20,6 +20,7 @@ dependencies { implementation project(':osrd-reporting') implementation project(":kt-osrd-sim-infra") + api project(":osrd-railjson") api project(":kt-osrd-utils") implementation libs.guava implementation libs.slf4j diff --git a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope/EnvelopePhysics.java b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope/EnvelopePhysics.java index 2247e381111..16b67b9ac84 100644 --- a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope/EnvelopePhysics.java +++ b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope/EnvelopePhysics.java @@ -1,5 +1,6 @@ package fr.sncf.osrd.envelope; +import static fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator.GRAVITY_ACCELERATION; import static fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator.areAccelerationsEqual; import static java.lang.Math.max; @@ -217,7 +218,7 @@ public static double getPartMechanicalEnergyConsumed( var meanGrade = 0.001 * path.getAverageGrade(part.getBeginPos(), part.getEndPos()); var altitudeDelta = meanGrade * part.getTotalDistance(); - var workGravity = -mass * 9.81 * altitudeDelta; + var workGravity = -mass * GRAVITY_ACCELERATION * altitudeDelta; var kineticEnergyDelta = 0.5 * inertia * (part.getEndSpeed() * part.getEndSpeed() - part.getBeginSpeed() * part.getBeginSpeed()); diff --git a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/EnvelopeSimPath.java b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/EnvelopeSimPath.java index a98b95ff519..ab2a9613b8b 100644 --- a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/EnvelopeSimPath.java +++ b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/EnvelopeSimPath.java @@ -106,6 +106,31 @@ public double getAverageGrade(double begin, double end) { return (getCumGrade(end) - getCumGrade(begin)) / (end - begin); } + @Override + public double getMinGrade(double begin, double end) { + // TODO: Optimise method by adding in a cache. + int indexBegin = getIndexBeforePos(begin); + int indexEnd = getIndexBeforePos(end); + var lowestGradient = gradeValues[indexBegin]; + for (int i = indexBegin; i < indexEnd; i++) { + var grad = gradeValues[i]; + if (grad < lowestGradient) lowestGradient = grad; + } + return lowestGradient; + } + + /** For a given position, return the index of the position just before in gradePositions */ + public int getIndexBeforePos(double position) { + // TODO: Optimise method by using binary search. + if (position <= gradePositions[0]) return 0; + if (position >= gradePositions[gradePositions.length - 1]) return gradePositions.length - 1; + for (int i = 0; i < gradePositions.length; i++) { + var pos = gradePositions[i]; + if (pos > position) return i - 1; + } + return gradePositions.length - 1; + } + private RangeMap getModeAndProfileMap( String powerClass, Range range, boolean ignoreElectricalProfiles) { if (ignoreElectricalProfiles) powerClass = null; diff --git a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsPath.java b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsPath.java index 811b4987a07..8acd441e78c 100644 --- a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsPath.java +++ b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsPath.java @@ -4,6 +4,9 @@ public interface PhysicsPath { /** The length of the path, in meters */ double getLength(); - /** The average slope on a given range, in meters per kilometers */ + /** The average slope on a given range, in m/km */ double getAverageGrade(double begin, double end); + + /** The lowest slope on a given range, in m/km */ + double getMinGrade(double begin, double end); } diff --git a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsRollingStock.java b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsRollingStock.java index 1a49e7afe9d..401b6d93f70 100644 --- a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsRollingStock.java +++ b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/PhysicsRollingStock.java @@ -1,5 +1,11 @@ package fr.sncf.osrd.envelope_sim; +import static fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator.GRAVITY_ACCELERATION; +import static fr.sncf.osrd.envelope_sim.etcs.ConstantsKt.mRotatingMax; +import static fr.sncf.osrd.envelope_sim.etcs.ConstantsKt.mRotatingMin; + +import fr.sncf.osrd.railjson.schema.rollingstock.RJSEtcsBrakeParams; + public interface PhysicsRollingStock { /** The mass of the train, in kilograms */ double getMass(); @@ -19,6 +25,8 @@ public interface PhysicsRollingStock { /** The first derivative of the resistance to movement at a given speed, in kg/s */ double getRollingResistanceDeriv(double speed); + RJSEtcsBrakeParams getRJSEtcsBrakeParams(); + /** Get the effort the train can apply at a given speed, in newtons */ static double getMaxEffort(double speed, TractiveEffortPoint[] tractiveEffortCurve) { int index = 0; @@ -50,6 +58,16 @@ static double getMaxEffort(double speed, TractiveEffortPoint[] tractiveEffortCur return previousPoint.maxEffort() + coeff * (Math.abs(speed) - previousPoint.speed()); } + /** + * The gradient acceleration of the rolling stock taking its rotating mass into account, in m/s². + * Grade is in m/km. + * mRotating (Max or Min) is in %, as seen in ERA braking curves simulation tool v5.1. + */ + static double getGradientAcceleration(double grade) { + var mRotating = grade >= 0 ? mRotatingMax : mRotatingMin; + return -GRAVITY_ACCELERATION * grade / (1000.0 + 10.0 * mRotating); + } + /** The maximum constant deceleration, in m/s^2 */ double getDeceleration(); diff --git a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegrator.java b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegrator.java index bc6eae181ff..8f0544cea2c 100644 --- a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegrator.java +++ b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegrator.java @@ -1,14 +1,21 @@ package fr.sncf.osrd.envelope_sim; +import static fr.sncf.osrd.envelope_sim.PhysicsRollingStock.getGradientAcceleration; +import static fr.sncf.osrd.envelope_sim.PhysicsRollingStock.getMaxEffort; + import com.google.common.collect.RangeMap; +import fr.sncf.osrd.envelope_sim.etcs.BrakingType; /** - * An utility class to help simulate the train, using numerical integration. It's used when + * A utility class to help simulate the train, using numerical integration. It's used when * simulating the train, and it is passed to speed controllers so they can take decisions about what * action to make. Once speed controllers took a decision, this same class is used to compute the * next position and speed of the train. */ public final class TrainPhysicsIntegrator { + // Gravity acceleration, in m/s² + public static final double GRAVITY_ACCELERATION = 9.81; + // A position delta lower than this value will be considered zero // Going back and forth with Distance and double (meters) often causes 1e-3 errors, // we need the tolerance to be higher than this @@ -25,18 +32,21 @@ public final class TrainPhysicsIntegrator { private final double directionSign; private final RangeMap tractiveEffortCurveMap; + private final BrakingType brakingType; private TrainPhysicsIntegrator( PhysicsRollingStock rollingStock, PhysicsPath path, Action action, double directionSign, - RangeMap tractiveEffortCurveMap) { + RangeMap tractiveEffortCurveMap, + BrakingType brakingType) { this.rollingStock = rollingStock; this.path = path; this.action = action; this.directionSign = directionSign; this.tractiveEffortCurveMap = tractiveEffortCurveMap; + this.brakingType = brakingType; } /** Simulates train movement */ @@ -46,8 +56,19 @@ public static IntegrationStep step( double initialSpeed, Action action, double directionSign) { + return step(context, initialLocation, initialSpeed, action, directionSign, BrakingType.CONSTANT); + } + + /** Simulates train movement */ + public static IntegrationStep step( + EnvelopeSimContext context, + double initialLocation, + double initialSpeed, + Action action, + double directionSign, + BrakingType brakingType) { var integrator = new TrainPhysicsIntegrator( - context.rollingStock, context.path, action, directionSign, context.tractiveEffortCurveMap); + context.rollingStock, context.path, action, directionSign, context.tractiveEffortCurveMap, brakingType); return integrator.step(context.timeStep, initialLocation, initialSpeed, directionSign); } @@ -65,16 +86,15 @@ private IntegrationStep step(double timeStep, double initialLocation, double ini } private IntegrationStep step(double timeStep, double position, double speed) { + if (action == Action.BRAKE) return newtonStep(timeStep, speed, getDeceleration(speed, position), directionSign); + double tractionForce = 0; var tractiveEffortCurve = tractiveEffortCurveMap.get(Math.min(Math.max(0, position), path.getLength())); assert tractiveEffortCurve != null; - double maxTractionForce = PhysicsRollingStock.getMaxEffort(speed, tractiveEffortCurve); + double maxTractionForce = getMaxEffort(speed, tractiveEffortCurve); double rollingResistance = rollingStock.getRollingResistance(speed); - double weightForce = getWeightForce(rollingStock, path, position); - - if (action == Action.ACCELERATE) tractionForce = maxTractionForce; - - boolean isBraking = (action == Action.BRAKE); + double averageGrade = getAverageGrade(rollingStock, path, position); + double weightForce = getWeightForce(rollingStock, averageGrade); if (action == Action.MAINTAIN) { tractionForce = rollingResistance - weightForce; @@ -82,20 +102,53 @@ private IntegrationStep step(double timeStep, double position, double speed) { else tractionForce = maxTractionForce; } - double acceleration = computeAcceleration( - rollingStock, rollingResistance, weightForce, speed, tractionForce, isBraking, directionSign); + if (action == Action.ACCELERATE) tractionForce = maxTractionForce; + double acceleration = + computeAcceleration(rollingStock, rollingResistance, weightForce, speed, tractionForce, directionSign); return newtonStep(timeStep, speed, acceleration, directionSign); } - /** Compute the weight force of a rolling stock at a given position on a given path */ - public static double getWeightForce(PhysicsRollingStock rollingStock, PhysicsPath path, double headPosition) { + private double getDeceleration(double speed, double position) { + assert (action == Action.BRAKE); + if (brakingType == BrakingType.CONSTANT) return rollingStock.getDeceleration(); + + var grade = getMinGrade(rollingStock, path, position); + var gradientAcceleration = getGradientAcceleration(grade); + return switch (brakingType) { + // See Subset referenced in RJSEtcsBrakeParams: §3.13.6.2.1.3. + case ETCS_EBD -> -rollingStock.getRJSEtcsBrakeParams().getSafeBrakingAcceleration(speed) + + gradientAcceleration; + // See Subset referenced in RJSEtcsBrakeParams: §3.13.6.3.1.3. + case ETCS_SBD -> -rollingStock.getRJSEtcsBrakeParams().getServiceBrakingAcceleration(speed) + + gradientAcceleration; + // See Subset referenced in RJSEtcsBrakeParams: §3.13.6.4.3. + case ETCS_GUI -> -rollingStock.getRJSEtcsBrakeParams().getNormalServiceBrakingAcceleration(speed) + + gradientAcceleration + + rollingStock.getRJSEtcsBrakeParams().getGradientAccelerationCorrection(grade, speed); + default -> throw new UnsupportedOperationException("Braking type not supported: " + brakingType); + }; + } + + /** Compute the average grade of a rolling stock at a given position on a given path in m/km */ + public static double getAverageGrade(PhysicsRollingStock rollingStock, PhysicsPath path, double headPosition) { var tailPosition = Math.min(Math.max(0, headPosition - rollingStock.getLength()), path.getLength()); headPosition = Math.min(Math.max(0, headPosition), path.getLength()); - var averageGrade = path.getAverageGrade(tailPosition, headPosition); - // get an angle from a meter per km elevation difference + return path.getAverageGrade(tailPosition, headPosition); + } + + /** Compute the weight force of a rolling stock at a given position on a given path */ + public static double getWeightForce(PhysicsRollingStock rollingStock, double grade) { + // get an angle from a m/km elevation difference // the curve's radius is taken into account in meanTrainGrade - var angle = Math.atan(averageGrade / 1000.0); // from m/km to m/m - return -rollingStock.getMass() * 9.81 * Math.sin(angle); + var angle = Math.atan(grade / 1000.0); // from m/km to m/m + return -rollingStock.getMass() * GRAVITY_ACCELERATION * Math.sin(angle); + } + + /** Compute the min grade of a rolling stock at a given position on a given path in m/km */ + public static double getMinGrade(PhysicsRollingStock rollingStock, PhysicsPath path, double headPosition) { + var tailPosition = Math.min(Math.max(0, headPosition - rollingStock.getLength()), path.getLength()); + headPosition = Math.min(Math.max(0, headPosition), path.getLength()); + return path.getMinGrade(tailPosition, headPosition); } /** @@ -107,15 +160,10 @@ public static double computeAcceleration( double weightForce, double currentSpeed, double tractionForce, - boolean isBraking, double directionSign) { assert tractionForce >= 0.; - if (isBraking) { - return rollingStock.getDeceleration(); - } - if (currentSpeed == 0 && directionSign > 0) { // If we are stopped and if the forces are not enough to compensate the opposite force, // the rolling resistance and braking force don't apply and the speed stays at 0 diff --git a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/overlays/EnvelopeDeceleration.java b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/overlays/EnvelopeDeceleration.java index ff5e8164a31..e9c56c3cb84 100644 --- a/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/overlays/EnvelopeDeceleration.java +++ b/core/envelope-sim/src/main/java/fr/sncf/osrd/envelope_sim/overlays/EnvelopeDeceleration.java @@ -4,6 +4,7 @@ import fr.sncf.osrd.envelope_sim.Action; import fr.sncf.osrd.envelope_sim.EnvelopeSimContext; import fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator; +import fr.sncf.osrd.envelope_sim.etcs.BrakingType; public class EnvelopeDeceleration { /** Generate a deceleration curve overlay */ @@ -12,15 +13,25 @@ public static void decelerate( double startPosition, double startSpeed, InteractiveEnvelopePartConsumer consumer, - double direction) { + double direction, + BrakingType brakingType) { if (!consumer.initEnvelopePart(startPosition, startSpeed, direction)) return; double position = startPosition; double speed = startSpeed; while (true) { - var step = TrainPhysicsIntegrator.step(context, position, speed, Action.BRAKE, direction); + var step = TrainPhysicsIntegrator.step(context, position, speed, Action.BRAKE, direction, brakingType); position += step.positionDelta; speed = step.endSpeed; if (!consumer.addStep(position, speed, step.timeDelta)) break; } } + + public static void decelerate( + EnvelopeSimContext context, + double startPosition, + double startSpeed, + InteractiveEnvelopePartConsumer consumer, + double direction) { + decelerate(context, startPosition, startSpeed, consumer, direction, BrakingType.CONSTANT); + } } diff --git a/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/Constants.kt b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/Constants.kt new file mode 100644 index 00000000000..c6e4fbfc97b --- /dev/null +++ b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/Constants.kt @@ -0,0 +1,6 @@ +package fr.sncf.osrd.envelope_sim.etcs + +/** See Subset referenced in ETCSBrakingSimulator: table in Appendix A.3.1. */ +const val tDriver = 4.0 // s +const val mRotatingMax = 15.0 // % +const val mRotatingMin = 2.0 // % diff --git a/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingCurves.kt b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingCurves.kt new file mode 100644 index 00000000000..60fabadc724 --- /dev/null +++ b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingCurves.kt @@ -0,0 +1,263 @@ +package fr.sncf.osrd.envelope_sim.etcs + +import fr.sncf.osrd.envelope.Envelope +import fr.sncf.osrd.envelope.OverlayEnvelopeBuilder +import fr.sncf.osrd.envelope.part.ConstrainedEnvelopePartBuilder +import fr.sncf.osrd.envelope.part.EnvelopePart +import fr.sncf.osrd.envelope.part.EnvelopePartBuilder +import fr.sncf.osrd.envelope.part.constraints.EnvelopeConstraint +import fr.sncf.osrd.envelope.part.constraints.EnvelopePartConstraintType +import fr.sncf.osrd.envelope.part.constraints.PositionConstraint +import fr.sncf.osrd.envelope.part.constraints.SpeedConstraint +import fr.sncf.osrd.envelope_sim.* +import fr.sncf.osrd.envelope_sim.overlays.EnvelopeDeceleration +import kotlin.math.max +import kotlin.math.min + +enum class BrakingCurveType { + EBD, // Emergency Brake Deceleration + EBI, // Emergency Brake Intervention + SBD, // Service Brake Deceleration + SBI, // Service Brake Intervention + GUI, // Guidance + PS, // Permitted Speed + IND // Indication +} + +enum class BrakingType { + CONSTANT, + ETCS_EBD, + ETCS_SBD, + ETCS_GUI +} + +/** Compute braking curves at every end of authority. */ +fun addBrakingCurvesAtEOAs( + envelope: Envelope, + context: EnvelopeSimContext, + endsOfAuthority: Collection +): Envelope { + val sortedEndsOfAuthority = endsOfAuthority.sortedBy { it.offsetEOA } + var beginPos = 0.0 + val builder = OverlayEnvelopeBuilder.forward(envelope) + for (endOfAuthority in sortedEndsOfAuthority) { + val targetPosition = endOfAuthority.offsetEOA.distance.meters + val targetSpeed = 0.0 + val overhead = + Envelope.make( + EnvelopePart.generateTimes( + listOf(EnvelopeProfile.CONSTANT_SPEED), + doubleArrayOf(beginPos, targetPosition), + doubleArrayOf(envelope.maxSpeed, envelope.maxSpeed) + ) + ) + val guiCurve = + computeBrakingCurve( + context, + overhead, + beginPos, + targetPosition, + targetSpeed, + BrakingType.ETCS_GUI + ) + val sbdCurve = + computeBrakingCurve( + context, + overhead, + beginPos, + targetPosition, + targetSpeed, + BrakingType.ETCS_SBD + ) + val fullIndicationCurve = + computeIndicationBrakingCurveFromRef( + context, + sbdCurve, + BrakingCurveType.SBD, + guiCurve, + beginPos + ) + val indicationCurve = keepBrakingCurveUnderOverlay(fullIndicationCurve, envelope, beginPos) + builder.addPart(indicationCurve) + beginPos = targetPosition + } + return builder.build() +} + +/** Compute braking curve: used to compute EBD, SBD or GUI. */ +private fun computeBrakingCurve( + context: EnvelopeSimContext, + envelope: Envelope, + beginPos: Double, + targetPosition: Double, + targetSpeed: Double, + brakingType: BrakingType +): EnvelopePart { + // If the stopPosition is below begin position, the input is invalid + // If the stopPosition is after the end of the path, the input is invalid except if it is an + // SVL, i.e. the target speed is 0 and the curve to compute is an EBD + if ( + targetPosition <= beginPos || + (targetPosition > context.path.length && + targetSpeed == 0.0 && + brakingType != BrakingType.ETCS_EBD) + ) + throw RuntimeException( + String.format( + "Trying to compute ETCS braking curve from out of bounds ERTMS end/limit of authority: %s", + targetPosition + ) + ) + val partBuilder = EnvelopePartBuilder() + partBuilder.setAttr(EnvelopeProfile.BRAKING) + val overlayBuilder = + ConstrainedEnvelopePartBuilder( + partBuilder, + PositionConstraint(beginPos, targetPosition), + SpeedConstraint(0.0, EnvelopePartConstraintType.FLOOR), + EnvelopeConstraint(envelope, EnvelopePartConstraintType.CEILING) + ) + EnvelopeDeceleration.decelerate( + context, + targetPosition, + targetSpeed, + overlayBuilder, + -1.0, + brakingType + ) + val brakingCurve = partBuilder.build() + return brakingCurve +} + +/** + * Compute Indication curve: EBI/SBD -> SBI -> PS -> IND. See Subset referenced in + * ETCSBrakingSimulator: figures 45 and 46. + */ +private fun computeIndicationBrakingCurveFromRef( + context: EnvelopeSimContext, + refBrakingCurve: EnvelopePart, + refBrakingCurveType: BrakingCurveType, + guiCurve: EnvelopePart, + beginPos: Double +): EnvelopePart { + assert(refBrakingCurve.endPos > beginPos) + val rollingStock = context.rollingStock + val tBs = + when (refBrakingCurveType) { + BrakingCurveType.EBI -> rollingStock.rjsEtcsBrakeParams.tBs2 + BrakingCurveType.SBD -> rollingStock.rjsEtcsBrakeParams.tBs1 + else -> + throw IllegalArgumentException( + "Expected EBI or SBD reference braking curve type, found: $refBrakingCurveType" + ) + } + + val reversedNewPos = ArrayList() + val reversedNewSpeeds = ArrayList() + for (i in refBrakingCurve.pointCount() - 1 downTo 0) { + val reversedNewPosIndex = refBrakingCurve.pointCount() - 1 - i + val speed = refBrakingCurve.getPointSpeed(i) + val sbiPosition = getSbiPosition(refBrakingCurve.getPointPos(i), speed, tBs) + val permittedSpeedPosition = getPermittedSpeedPosition(sbiPosition, speed) + val adjustedPermittedSpeedPosition = + getAdjustedPermittedSpeedPosition( + permittedSpeedPosition, + guiCurve.interpolatePosition(speed) + ) + val indicationPosition = getIndicationPosition(adjustedPermittedSpeedPosition, speed, tBs) + if (indicationPosition >= beginPos) { + reversedNewPos.add(indicationPosition) + reversedNewSpeeds.add(speed) + } else { + assert(reversedNewPosIndex > 0 && reversedNewPos[reversedNewPosIndex - 1] > beginPos) + // Interpolate to begin position if reaching a position before it + val prevPos = reversedNewPos[reversedNewPosIndex - 1] + val prevSpeed = reversedNewSpeeds[reversedNewPosIndex - 1] + val speedAtBeginPos = + prevSpeed + + (beginPos - prevPos) * (speed - prevSpeed) / (indicationPosition - prevPos) + reversedNewPos.add(beginPos) + reversedNewSpeeds.add(speedAtBeginPos) + break + } + } + + val nbPoints = reversedNewPos.size + val newPosArray = DoubleArray(nbPoints) + val newSpeedsArray = DoubleArray(nbPoints) + for (i in newPosArray.indices) { + newPosArray[i] = reversedNewPos[nbPoints - 1 - i] + newSpeedsArray[i] = reversedNewSpeeds[nbPoints - 1 - i] + } + val brakingCurve = + EnvelopePart.generateTimes(listOf(EnvelopeProfile.BRAKING), newPosArray, newSpeedsArray) + + return brakingCurve +} + +/** + * Keep the part of the full braking curve which is located underneath the overlay and intersects + * with it or with begin position. + */ +private fun keepBrakingCurveUnderOverlay( + fullBrakingCurve: EnvelopePart, + overlay: Envelope, + beginPos: Double +): EnvelopePart { + assert(fullBrakingCurve.maxSpeed >= overlay.minSpeed) + assert(fullBrakingCurve.endPos > beginPos) + val positions = fullBrakingCurve.clonePositions() + val speeds = fullBrakingCurve.cloneSpeeds() + val timeDeltas = fullBrakingCurve.cloneTimes() + val nbPoints = fullBrakingCurve.pointCount() + val partBuilder = EnvelopePartBuilder() + partBuilder.setAttr(EnvelopeProfile.BRAKING) + val overlayBuilder = + ConstrainedEnvelopePartBuilder( + partBuilder, + PositionConstraint(beginPos, overlay.endPos), + SpeedConstraint(0.0, EnvelopePartConstraintType.FLOOR), + EnvelopeConstraint(overlay, EnvelopePartConstraintType.CEILING) + ) + overlayBuilder.initEnvelopePart(positions[nbPoints - 1], speeds[nbPoints - 1], -1.0) + for (i in fullBrakingCurve.pointCount() - 2 downTo 0) { + if (!overlayBuilder.addStep(positions[i], speeds[i], timeDeltas[i])) break + } + return partBuilder.build() +} + +/** See Subset referenced in ETCSBrakingSimulator: §3.13.9.3.3.1 and §3.13.9.3.3.2. */ +private fun getSbiPosition(ebiOrSbdPosition: Double, speed: Double, tbs: Double): Double { + return getPreviousPosition(ebiOrSbdPosition, speed, tbs) +} + +/** See Subset referenced in ETCSBrakingSimulator: §3.13.9.3.5.1. */ +private fun getPermittedSpeedPosition(sbiPosition: Double, speed: Double): Double { + return getPreviousPosition(sbiPosition, speed, tDriver) +} + +/** See Subset referenced in ETCSBrakingSimulator: §3.13.9.3.5.4. */ +private fun getAdjustedPermittedSpeedPosition( + permittedSpeedPosition: Double, + guiPosition: Double = Double.POSITIVE_INFINITY +): Double { + return min(permittedSpeedPosition, guiPosition) +} + +/** See Subset referenced in ETCSBrakingSimulator: §3.13.9.3.6.1 and §3.13.9.3.6.2. */ +private fun getIndicationPosition( + permittedSpeedPosition: Double, + speed: Double, + tBs: Double +): Double { + val tIndication = max((0.8 * tBs), 5.0) + tDriver + return getPreviousPosition(permittedSpeedPosition, speed, tIndication) +} + +private fun getPreviousPosition(position: Double, speed: Double, elapsedTime: Double): Double { + return getPreviousPosition(position, speed * elapsedTime) +} + +private fun getPreviousPosition(position: Double, elapsedDistance: Double): Double { + return position - elapsedDistance +} diff --git a/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingSimulator.kt b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingSimulator.kt index 43800a7b64e..e73ed9d69a7 100644 --- a/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingSimulator.kt +++ b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/etcs/ETCSBrakingSimulator.kt @@ -1,21 +1,28 @@ package fr.sncf.osrd.envelope_sim.etcs import fr.sncf.osrd.envelope.Envelope -import fr.sncf.osrd.envelope_sim.PhysicsPath -import fr.sncf.osrd.envelope_sim.PhysicsRollingStock +import fr.sncf.osrd.envelope_sim.EnvelopeSimContext import fr.sncf.osrd.sim_infra.api.Path import fr.sncf.osrd.sim_infra.api.TravelledPath import fr.sncf.osrd.utils.units.Offset +/** + * In charge of computing and adding the ETCS braking curves. Formulas are found in `SUBSET-026-3 + * v400.pdf` from the file at + * https://www.era.europa.eu/system/files/2023-09/index004_-_SUBSET-026_v400.zip + */ interface ETCSBrakingSimulator { - val path: PhysicsPath - val rollingStock: PhysicsRollingStock - val timeStep: Double - - /** Compute the ETCS braking envelope from the MRSP, for each LOA and EOA. */ - fun addETCSBrakingParts( - mrsp: Envelope, - limitsOfAuthority: Collection, + val context: EnvelopeSimContext + + /** Compute the ETCS braking envelope for each LOA. */ + fun addSlowdownBrakingCurves( + envelope: Envelope, + limitsOfAuthority: Collection + ): Envelope + + /** Compute the ETCS braking envelope for each EOA. */ + fun addStopBrakingCurves( + envelope: Envelope, endsOfAuthority: Collection ): Envelope } @@ -38,18 +45,22 @@ data class EndOfAuthority( } } -class ETCSBrakingSimulatorImpl( - override val path: PhysicsPath, - override val rollingStock: PhysicsRollingStock, - override val timeStep: Double -) : ETCSBrakingSimulator { +class ETCSBrakingSimulatorImpl(override val context: EnvelopeSimContext) : ETCSBrakingSimulator { + override fun addSlowdownBrakingCurves( + envelope: Envelope, + limitsOfAuthority: Collection + ): Envelope { + if (limitsOfAuthority.isEmpty()) return envelope + // TODO: implement braking at LOAs CORRECTLY + return envelope + } - override fun addETCSBrakingParts( - mrsp: Envelope, - limitsOfAuthority: Collection, + override fun addStopBrakingCurves( + envelope: Envelope, endsOfAuthority: Collection ): Envelope { - if (limitsOfAuthority.isEmpty() && endsOfAuthority.isEmpty()) return mrsp - TODO("Not yet implemented") + if (endsOfAuthority.isEmpty()) return envelope + // TODO: compute curve at SvL and intersect + return addBrakingCurvesAtEOAs(envelope, context, endsOfAuthority) } } diff --git a/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/pipelines/MaxSpeedEnvelope.kt b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/pipelines/MaxSpeedEnvelope.kt index fefeed9d85e..6c76f2d7378 100644 --- a/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/pipelines/MaxSpeedEnvelope.kt +++ b/core/envelope-sim/src/main/kotlin/fr/sncf/osrd/envelope_sim/pipelines/MaxSpeedEnvelope.kt @@ -122,7 +122,7 @@ object MaxSpeedEnvelope { } cursor.nextPart() } - return etcsSimulator.addETCSBrakingParts(envelope, limitsOfAuthority, listOf()) + return etcsSimulator.addSlowdownBrakingCurves(envelope, limitsOfAuthority) } /** Generate braking curves overlay at every stop position */ @@ -179,7 +179,7 @@ object MaxSpeedEnvelope { stops .filter { it.isETCS } .map { EndOfAuthority(Offset(it.offset.meters), getDangerPoint(context, it)) } - return simulator.addETCSBrakingParts(envelope, listOf(), endsOfAuthority) + return simulator.addStopBrakingCurves(envelope, endsOfAuthority) } /** @@ -218,8 +218,7 @@ object MaxSpeedEnvelope { /** Generate a max speed envelope given a mrsp */ @JvmStatic fun from(context: EnvelopeSimContext, stopPositions: DoubleArray, mrsp: Envelope): Envelope { - val etcsSimulator = - ETCSBrakingSimulatorImpl(context.path, context.rollingStock, context.timeStep) + val etcsSimulator = ETCSBrakingSimulatorImpl(context) var maxSpeedEnvelope = addSlowdownBrakingCurves(etcsSimulator, context, mrsp) maxSpeedEnvelope = addStopBrakingCurves(etcsSimulator, context, stopPositions, maxSpeedEnvelope) diff --git a/core/envelope-sim/src/test/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegratorTest.java b/core/envelope-sim/src/test/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegratorTest.java index 1a8816cc521..68c2893fb83 100644 --- a/core/envelope-sim/src/test/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegratorTest.java +++ b/core/envelope-sim/src/test/java/fr/sncf/osrd/envelope_sim/TrainPhysicsIntegratorTest.java @@ -1,8 +1,7 @@ package fr.sncf.osrd.envelope_sim; import static fr.sncf.osrd.envelope_sim.SimpleContextBuilder.makeSimpleContext; -import static fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator.getWeightForce; -import static fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator.newtonStep; +import static fr.sncf.osrd.envelope_sim.TrainPhysicsIntegrator.*; import static fr.sncf.osrd.envelope_sim.allowances.mareco_impl.CoastingGenerator.coastFromBeginning; import static org.junit.jupiter.api.Assertions.*; @@ -90,9 +89,10 @@ public void testAccelerateAndCoast() { // make a huge traction effort double rollingResistance = testRollingStock.getRollingResistance(speed); - double weightForce = getWeightForce(testRollingStock, testPath, position); + double grade = getAverageGrade(testRollingStock, testPath, position); + double weightForce = getWeightForce(testRollingStock, grade); var acceleration = TrainPhysicsIntegrator.computeAcceleration( - testRollingStock, rollingResistance, weightForce, speed, 500000.0, false, +1); + testRollingStock, rollingResistance, weightForce, speed, 500000.0, +1); var step = newtonStep(TIME_STEP, speed, acceleration, +1); position += step.positionDelta; speed = step.endSpeed; diff --git a/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/FlatPath.java b/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/FlatPath.java index 4553e6f6eed..9710b694012 100644 --- a/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/FlatPath.java +++ b/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/FlatPath.java @@ -18,4 +18,9 @@ public double getLength() { public double getAverageGrade(double begin, double end) { return slope; } + + @Override + public double getMinGrade(double begin, double end) { + return slope; + } } diff --git a/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/SimpleRollingStock.java b/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/SimpleRollingStock.java index 0ba35134e63..ac546adfce7 100644 --- a/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/SimpleRollingStock.java +++ b/core/envelope-sim/src/testFixtures/java/fr/sncf/osrd/envelope_sim/SimpleRollingStock.java @@ -3,6 +3,7 @@ import com.google.common.collect.ImmutableRangeMap; import com.google.common.collect.Range; import edu.umd.cs.findbugs.annotations.SuppressFBWarnings; +import fr.sncf.osrd.railjson.schema.rollingstock.RJSEtcsBrakeParams; import java.util.ArrayList; @SuppressFBWarnings({"URF_UNREAD_PUBLIC_OR_PROTECTED_FIELD"}) @@ -90,6 +91,11 @@ public double getRollingResistanceDeriv(double speed) { return B + 2 * C * speed; } + @Override + public RJSEtcsBrakeParams getRJSEtcsBrakeParams() { + throw new UnsupportedOperationException("To be implemented"); + } + @Override public double getDeceleration() { return -constGamma; diff --git a/core/osrd-geom/build.gradle b/core/osrd-geom/build.gradle index 58cb8f08186..00bca43d060 100644 --- a/core/osrd-geom/build.gradle +++ b/core/osrd-geom/build.gradle @@ -1,4 +1,5 @@ plugins { + alias(libs.plugins.kotlin.jvm) id 'java' id 'jacoco' } diff --git a/core/osrd-railjson/build.gradle b/core/osrd-railjson/build.gradle index 50fcea14c79..819c03fc00b 100644 --- a/core/osrd-railjson/build.gradle +++ b/core/osrd-railjson/build.gradle @@ -1,4 +1,5 @@ plugins { + alias(libs.plugins.kotlin.jvm) id 'java' id 'jacoco' } @@ -20,8 +21,8 @@ dependencies { implementation libs.hppc // JSON parsing - implementation libs.moshi - implementation libs.moshi.adapters + api libs.moshi + api libs.moshi.adapters implementation libs.guava diff --git a/core/osrd-railjson/src/main/java/fr/sncf/osrd/railjson/schema/rollingstock/RJSEtcsBrakeParams.java b/core/osrd-railjson/src/main/java/fr/sncf/osrd/railjson/schema/rollingstock/RJSEtcsBrakeParams.java index c9be800c6dd..872e05dbb4a 100644 --- a/core/osrd-railjson/src/main/java/fr/sncf/osrd/railjson/schema/rollingstock/RJSEtcsBrakeParams.java +++ b/core/osrd-railjson/src/main/java/fr/sncf/osrd/railjson/schema/rollingstock/RJSEtcsBrakeParams.java @@ -9,39 +9,42 @@ */ public class RJSEtcsBrakeParams { + /** National Default Value: Available Adhesion. Found in Subset Appendix A.3.2 table. */ + private static final double mNvavadh = 0.0; + // A_brake_emergency: the emergency deceleration curve (values > 0 m/s²) @Json(name = "gamma_emergency") - public RJSSpeedIntervalValueCurve gammaEmergency; + private RJSSpeedIntervalValueCurve gammaEmergency; // A_brake_service: the full service deceleration curve (values > 0 m/s²) @Json(name = "gamma_service") - public RJSSpeedIntervalValueCurve gammaService; + private RJSSpeedIntervalValueCurve gammaService; // A_brake_normal_service: the normal service deceleration curve used to compute guidance curve (values > 0 m/s²) @Json(name = "gamma_normal_service") - public RJSSpeedIntervalValueCurve gammaNormalService; + private RJSSpeedIntervalValueCurve gammaNormalService; // Kdry_rst: the rolling stock deceleration correction factors for dry rails // Boundaries should be the same as gammaEmergency // Values (no unit) should be contained in [0, 1] @Json(name = "k_dry") - public RJSSpeedIntervalValueCurve kDry; + private RJSSpeedIntervalValueCurve kDry; // Kwet_rst: the rolling stock deceleration correction factors for wet rails // Boundaries should be the same as gammaEmergency // Values (no unit) should be contained in [0, 1] @Json(name = "k_wet") - public RJSSpeedIntervalValueCurve kWet; + private RJSSpeedIntervalValueCurve kWet; // Kn+: the correction acceleration factor on normal service deceleration in positive gradients // Values (in m/s²) should be contained in [0, 10] @Json(name = "k_n_pos") - public RJSSpeedIntervalValueCurve kNPos; + private RJSSpeedIntervalValueCurve kNPos; // Kn-: the correction acceleration factor on normal service deceleration in negative gradients // Values (in m/s²) should be contained in [0, 10] @Json(name = "k_n_neg") - public RJSSpeedIntervalValueCurve kNNeg; + private RJSSpeedIntervalValueCurve kNNeg; // T_traction_cut_off: time delay in s from the traction cut-off command to the moment the acceleration due to // traction is zero @@ -60,6 +63,49 @@ public class RJSEtcsBrakeParams { @Json(name = "t_be") public double tBe; + /** See Subset §3.13.6.2.1.4. */ + public double getSafeBrakingAcceleration(double speed) { + var aBrakeEmergency = getEmergencyBrakingDeceleration(speed); + var kDry = getRollingStockCorrectionFactorDry(speed); + var kWet = getRollingStockCorrectionFactorWet(speed); + return kDry * (kWet + mNvavadh * (1 - kWet)) * aBrakeEmergency; + } + + private double getEmergencyBrakingDeceleration(double speed) { + return gammaEmergency.getValue(speed); + } + + /** + * Corresponds to the correction factor of the emergency brake deceleration on dry tracks. + * The confidence level mNvebcl is the confidence level that the corresponding deceleration can be reached, + * but does not impact the calculation of kDry. See Subset §3.13.6.2.1.7. + */ + private double getRollingStockCorrectionFactorDry(double speed) { + return kDry.getValue(speed); + } + + /** Corresponds to the correction factor of the emergency brake deceleration on wet tracks. */ + private double getRollingStockCorrectionFactorWet(double speed) { + return kWet.getValue(speed); + } + + public double getServiceBrakingAcceleration(double speed) { + return gammaService.getValue(speed); + } + + public double getNormalServiceBrakingAcceleration(double speed) { + return gammaNormalService.getValue(speed); + } + + /** + * Gradient acceleration correction using on-board correction factors kN+ and kN-. + * See Subset, §3.13.6.4.2 and §3.13.6.4.3. + */ + public double getGradientAccelerationCorrection(double grade, double speed) { + var k = grade >= 0 ? kNPos.getValue(speed) : kNNeg.getValue(speed); + return -k * grade / 1000; + } + public static final class RJSSpeedIntervalValueCurve { // Speed in m/s (sorted ascending). External bounds are implicit to [0, rolling_stock.max_speed] public double[] boundaries; diff --git a/core/src/main/java/fr/sncf/osrd/standalone_sim/StandaloneSim.java b/core/src/main/java/fr/sncf/osrd/standalone_sim/StandaloneSim.java index bda637e8736..ff776e4f37a 100644 --- a/core/src/main/java/fr/sncf/osrd/standalone_sim/StandaloneSim.java +++ b/core/src/main/java/fr/sncf/osrd/standalone_sim/StandaloneSim.java @@ -92,7 +92,7 @@ public static StandaloneSimResult run( // extend the // neutral sections (with time to lower/raise pantograph...) context = context.updateCurves(rollingStock.addNeutralSystemTimes( - electrificationMap, trainSchedule.comfort, maxSpeedEnvelope, context.tractiveEffortCurveMap)); + electrificationMap, trainSchedule.comfort, context.tractiveEffortCurveMap)); var envelope = MaxEffortEnvelope.from(context, trainSchedule.initialSpeed, maxSpeedEnvelope); var simResultTrain = diff --git a/core/src/main/java/fr/sncf/osrd/train/RollingStock.java b/core/src/main/java/fr/sncf/osrd/train/RollingStock.java index e3f346138df..b40c4352d95 100644 --- a/core/src/main/java/fr/sncf/osrd/train/RollingStock.java +++ b/core/src/main/java/fr/sncf/osrd/train/RollingStock.java @@ -5,7 +5,6 @@ import com.google.common.collect.RangeMap; import com.google.common.collect.TreeRangeMap; import edu.umd.cs.findbugs.annotations.SuppressFBWarnings; -import fr.sncf.osrd.envelope.Envelope; import fr.sncf.osrd.envelope_sim.PhysicsRollingStock; import fr.sncf.osrd.envelope_sim.electrification.Electrification; import fr.sncf.osrd.envelope_sim.electrification.Electrified; @@ -123,6 +122,12 @@ public double getRollingResistanceDeriv(double speed) { return B + 2 * C * speed; } + @Override + public RJSEtcsBrakeParams getRJSEtcsBrakeParams() { + assert etcsBrakeParams != null; + return etcsBrakeParams; + } + public record ModeEffortCurves( boolean isElectric, TractiveEffortPoint[] defaultCurve, ConditionalEffortCurve[] curves) {} @@ -164,7 +169,7 @@ private boolean shouldCoast(Neutral n, Comfort comfort) { * Returns the tractive effort curve that matches best, along with the electrification * conditions that matched */ - protected CurveAndCondition findTractiveEffortCurve(Comfort comfort, Electrification electrification) { + private CurveAndCondition findTractiveEffortCurve(Comfort comfort, Electrification electrification) { if (electrification instanceof Neutral n) { if (shouldCoast(n, comfort)) { return new CurveAndCondition(COASTING_CURVE, new InfraConditions(null, null, null)); @@ -220,12 +225,10 @@ public CurvesAndConditions mapTractiveEffortCurves( * * @param electrificationMap The map of electrification conditions to use * @param comfort The comfort level to get the curves for - * @param maxSpeedEnvelope The maxSpeedEnvelope used to extend the neutral sections */ public RangeMap addNeutralSystemTimes( RangeMap electrificationMap, Comfort comfort, - Envelope maxSpeedEnvelope, RangeMap curvesUsed) { TreeRangeMap newCurves = TreeRangeMap.create(); @@ -260,7 +263,7 @@ public boolean isThermal() { } /** Return whether this rolling stock's has an electric mode */ - public final boolean isElectric() { + public boolean isElectric() { return modes.values().stream().anyMatch(ModeEffortCurves::isElectric); } diff --git a/core/src/main/kotlin/fr/sncf/osrd/standalone_sim/StandaloneSimulation.kt b/core/src/main/kotlin/fr/sncf/osrd/standalone_sim/StandaloneSimulation.kt index 9374d9c4cac..2ab10f84534 100644 --- a/core/src/main/kotlin/fr/sncf/osrd/standalone_sim/StandaloneSimulation.kt +++ b/core/src/main/kotlin/fr/sncf/osrd/standalone_sim/StandaloneSimulation.kt @@ -113,7 +113,6 @@ fun runStandaloneSimulation( rollingStock.addNeutralSystemTimes( electrificationMap, comfort, - maxSpeedEnvelope, context.tractiveEffortCurveMap ) )