2022-04-28 17:29:24 +02:00
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __author__ = "Max Kannenberg"
# __copyright__ = "2020-2022"
# __license__ = "ISC"
2022-05-04 16:34:17 +02:00
# Calculate the running time of a train run on a path with special settings with information from the corresponding YAML files with the file paths `trainDirectory`, `pathDirectory`, `settingsDirectory`.
2022-04-28 17:29:24 +02:00
# calculate a train run focussing on using the minimum possible running time
2022-08-17 22:30:32 +02:00
function calculateMinimumRunningTime ( CSs :: Vector { Dict } , settings :: Settings , train :: Train )
2022-06-22 11:11:44 +02:00
startingPoint = SupportPoint ( )
startingPoint [ :s ] = CSs [ 1 ] [ :s_entry ]
2022-04-28 17:29:24 +02:00
calculateForces! ( startingPoint , CSs , 1 , " default " , train , settings . massModel ) # traction effort and resisting forces (in N)
2022-06-22 11:11:44 +02:00
drivingCourse :: Vector { Dict } = [ startingPoint ] # List of support points
2022-04-28 17:29:24 +02:00
for csId in 1 : length ( CSs )
CS = CSs [ csId ]
2022-06-22 11:11:44 +02:00
# determine the different flags for switching between the states for creating moving phases
2022-06-05 15:41:28 +02:00
s_braking = brakingDistance ( drivingCourse [ end ] [ :v ] , CS [ :v_exit ] , train . a_braking , settings . approxLevel )
2022-08-17 22:26:46 +02:00
calculateForces! ( drivingCourse [ end ] , CSs , csId , " default " , train , settings . massModel ) # tractive effort and resisting forces (in N)
2022-04-28 17:29:24 +02:00
previousSpeedLimitReached = false
stateFlags = Dict ( :endOfCSReached => drivingCourse [ end ] [ :s ] > CS [ :s_exit ] ,
:brakingStartReached => drivingCourse [ end ] [ :s ] + s_braking == CS [ :s_exit ] ,
:tractionDeficit => drivingCourse [ end ] [ :F_T ] < drivingCourse [ end ] [ :F_R ] , # or add another flag for equal forces?
:resistingForceNegative => drivingCourse [ end ] [ :F_R ] < 0.0 ,
2022-08-17 22:26:46 +02:00
:previousSpeedLimitReached => false ,
2022-04-28 17:29:24 +02:00
:speedLimitReached => drivingCourse [ end ] [ :v ] > CS [ :v_limit ] ,
:error => false )
2022-08-17 22:26:46 +02:00
# determine the behavior sections for this characteristic section. It has to be at least one of those BS: "breakFree", "clearing", "accelerating", "cruising", "diminishing", "coasting", "braking" or "halt")
while ! stateFlags [ :endOfCSReached ] # s < s_exit
if stateFlags [ :error ]
error ( " ERROR in calc in CS " , csId , " : BS= " , drivingCourse [ end ] [ :behavior ] , " s= " , drivingCourse [ end ] [ :s ] , " s_braking= " , s_braking , " v_limit= " , CS [ :v_limit ] , " v= " , drivingCourse [ end ] [ :v ] , " v_exit= " , CS [ :v_exit ] , " with the flags: endOfCS: " , stateFlags [ :endOfCSReached ] , " brakingStart: " , stateFlags [ :brakingStartReached ] , " F_T<F_R: " , stateFlags [ :tractionDeficit ] , " F_R<0: " , stateFlags [ :resistingForceNegative ] , " v_previousLimit: " , stateFlags [ :previousSpeedLimitReached ] , " v_limit: " , stateFlags [ :speedLimitReached ] , " error: " , stateFlags [ :error ] )
end
2022-04-28 17:29:24 +02:00
2022-08-17 22:26:46 +02:00
if ! stateFlags [ :brakingStartReached ] # s+s_braking < s_exit
if ! stateFlags [ :tractionDeficit ]
if drivingCourse [ end ] [ :F_T ] > drivingCourse [ end ] [ :F_R ] && drivingCourse [ end ] [ :v ] == 0.0
( drivingCourse , stateFlags ) = addBreakFreeSection! ( drivingCourse , stateFlags , CSs , csId , settings , train )
elseif stateFlags [ :previousSpeedLimitReached ]
( drivingCourse , stateFlags ) = addClearingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train )
elseif drivingCourse [ end ] [ :F_T ] > drivingCourse [ end ] [ :F_R ] && ! stateFlags [ :speedLimitReached ]
( drivingCourse , stateFlags ) = addAcceleratingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train )
elseif drivingCourse [ end ] [ :F_T ] == drivingCourse [ end ] [ :F_R ] && ! stateFlags [ :speedLimitReached ]
# cruise only one step
if settings . stepVariable == :distance
s_cruising = settings . stepSize
elseif settings . stepVariable == time
s_cruising = Δs_with_Δt ( settings . stepSize , drivingCourse [ end ] [ :a ] , drivingCourse [ end ] [ :v ] )
elseif settings . stepVariable == velocity
s_cruising = train . length / ( 10.0 ) # TODO which step size should be used?
end
( drivingCourse , stateFlags ) = addCruisingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train , " cruising " , s_cruising )
elseif drivingCourse [ end ] [ :F_R ] < 0 && stateFlags [ :speedLimitReached ]
s_braking = brakingDistance ( drivingCourse [ end ] [ :v ] , CS [ :v_exit ] , train . a_braking , settings . approxLevel )
s_cruising = CS [ :s_exit ] - drivingCourse [ end ] [ :s ] - s_braking
if s_cruising > 0.0
( drivingCourse , stateFlags ) = addCruisingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train , " downhillBraking " , s_cruising )
else
stateFlags [ :brakingStartReached ] = true
end
elseif drivingCourse [ end ] [ :F_T ] == drivingCourse [ end ] [ :F_R ] || stateFlags [ :speedLimitReached ]
s_braking = brakingDistance ( drivingCourse [ end ] [ :v ] , CS [ :v_exit ] , train . a_braking , settings . approxLevel )
s_cruising = CS [ :s_exit ] - drivingCourse [ end ] [ :s ] - s_braking
if s_cruising > 1 / 10 ^ ( settings . approxLevel ) # TODO: define another minimum cruising length?
( drivingCourse , stateFlags ) = addCruisingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train , " cruising " , s_cruising )
else
stateFlags [ :brakingStartReached ] = true
end
2022-04-28 17:29:24 +02:00
else
2022-08-17 22:26:46 +02:00
error ( )
2022-04-28 17:29:24 +02:00
end
2022-08-17 22:26:46 +02:00
elseif stateFlags [ :tractionDeficit ]
( drivingCourse , stateFlags ) = addDiminishingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train )
2022-04-28 17:29:24 +02:00
else
2022-08-17 22:26:46 +02:00
error ( )
2022-04-28 17:29:24 +02:00
end
2022-08-17 22:26:46 +02:00
else #if !stateFlags[:endOfCSReached] # s < s_exit
( drivingCourse , stateFlags ) = addBrakingSection! ( drivingCourse , stateFlags , CSs , csId , settings , train )
#else
# error()
2022-04-28 17:29:24 +02:00
end
2022-08-12 13:24:27 +02:00
2022-08-17 22:26:46 +02:00
if CS [ :s_exit ] - drivingCourse [ end ] [ :s ] < 1 / 10 ^ ( settings . approxLevel )
drivingCourse [ end ] [ :s ] = CS [ :s_exit ] # round s up to CS[:s_exit]
2022-08-12 13:24:27 +02:00
2022-08-17 22:26:46 +02:00
# set state flag
stateFlags [ :endOfCSReached ] = true
end
2022-08-12 13:24:27 +02:00
end
2022-08-17 22:26:46 +02:00
#if s == s_exit
# halt
#end
2022-04-28 17:29:24 +02:00
2022-07-12 16:59:45 +02:00
# for testing: # TODO
2022-04-28 17:29:24 +02:00
if drivingCourse [ end ] [ :s ] != CS [ :s_exit ]
println ( " ERROR: In CS " , csId , " the train run ends at s= " , drivingCourse [ end ] [ :s ] , " and not s_exit= " , CS [ :s_exit ] )
end
if drivingCourse [ end ] [ :v ] > CS [ :v_exit ]
println ( " ERROR: In CS " , csId , " the train run ends with v= " , drivingCourse [ end ] [ :v ] , " and not with v_exit= " , CS [ :v_exit ] )
end
end #for
2022-08-17 22:26:46 +02:00
drivingCourse = addHalt! ( drivingCourse , CSs , length ( CSs ) , settings , train )
2022-04-28 17:29:24 +02:00
2022-08-17 22:30:32 +02:00
return drivingCourse
2022-04-28 17:29:24 +02:00
end #function calculateMinimumRunningTime
2022-06-03 17:24:16 +02:00
2022-06-05 15:25:12 +02:00
2022-06-03 17:24:16 +02:00
"""
calculateTractiveEffort ( v , tractiveEffortVelocityPairs )
Calculate the trains tractive effort with the ` tractiveEffortVelocityPairs ` dependend on the velocity ` v ` .
...
# Arguments
- ` v::AbstractFloat ` : the current velocity in m / s .
- ` tractiveEffortVelocityPairs::Array{} ` : the trains pairs for velocity in m / s and tractive effort in N as one array containing an array for each pair .
...
# Examples
``` julia-repl
julia > calculateTractiveEffort ( 20.0 , [ ( 0.0 , 180000 ) , ( 20.0 , 100000 ) , ( 40.0 , 60000 ) , ( 60.0 , 40000 ) , ( 80.0 , 30000 ) ] )
100000
julia > calculateTractiveEffort ( 30.0 , [ ( 0.0 , 180000 ) , ( 20.0 , 100000 ) , ( 40.0 , 60000 ) , ( 60.0 , 40000 ) , ( 80.0 , 30000 ) ] )
80000
```
"""
function calculateTractiveEffort ( v :: AbstractFloat , tractiveEffortVelocityPairs :: Array { } )
2022-06-05 15:25:12 +02:00
if v < 0.0
#println("v=",v)
return 0.0
end
2022-06-03 17:24:16 +02:00
for row in 1 : length ( tractiveEffortVelocityPairs )
nextPair = tractiveEffortVelocityPairs [ row ]
if nextPair [ 1 ] == v
return nextPair [ 2 ]
elseif nextPair [ 1 ] > v
# interpolate for a straight line between the two surrounding points with the formula: F=(v-v_(row-1))*(F_row-F_(row-1))/(v_row-v_(row-1))+F_(row-1)
previousPair = tractiveEffortVelocityPairs [ row - 1 ]
F_T_interpolation = ( v - previousPair [ 1 ] ) * ( nextPair [ 2 ] - previousPair [ 2 ] ) / ( nextPair [ 1 ] - previousPair [ 1 ] ) + previousPair [ 2 ]
return F_T_interpolation
end #if
end #for
# if v gets higher than the velocities in tractiveEffortVelocityPairs the last tractive effort will be used
# TODO: also an extrapolation could be used
return tractiveEffortVelocityPairs [ end ] [ 2 ]
end #function calculateTractiveEffort
2022-06-05 15:25:12 +02:00
2022-06-03 17:24:16 +02:00
"""
calculate and return the path resistance dependend on the trains position and mass model
"""
function calculatePathResistance ( CSs :: Vector { Dict } , csId :: Integer , s :: Real , massModel , train :: Train )
if massModel == :mass_point
2022-06-05 15:41:28 +02:00
pathResistance = forceFromCoefficient ( CSs [ csId ] [ :r_path ] , train . m_train_full )
2022-06-03 17:24:16 +02:00
elseif massModel == :homogeneous_strip
pathResistance = 0.0
s_rear = s - train . length # position of the rear of the train
while csId > 0 && s_rear < CSs [ csId ] [ :s_exit ]
2022-06-05 15:41:28 +02:00
pathResistance = pathResistance + ( min ( s , CSs [ csId ] [ :s_exit ] ) - max ( s_rear , CSs [ csId ] [ :s_entry ] ) ) / train . length * forceFromCoefficient ( CSs [ csId ] [ :r_path ] , train . m_train_full )
2022-06-03 17:24:16 +02:00
csId = csId - 1
if csId == 0
2022-07-12 16:59:45 +02:00
# TODO: currently for values < s_trainrun_start the values of s_trainrun_start will be used
2022-06-05 15:41:28 +02:00
return pathResistance + ( CSs [ 1 ] [ :s_entry ] - s_rear ) / train . length * forceFromCoefficient ( CSs [ 1 ] [ :r_path ] , train . m_train_full )
2022-06-03 17:24:16 +02:00
end #if
end #while
end #if
return pathResistance
end #function calculatePathResistance
2022-06-05 15:25:12 +02:00
2022-06-03 17:24:16 +02:00
"""
2022-06-22 11:11:44 +02:00
calculate and return tractive and resisting forces for a support point
2022-06-03 17:24:16 +02:00
"""
2022-06-22 11:11:44 +02:00
function calculateForces! ( supportPoint :: Dict , CSs :: Vector { Dict } , csId :: Integer , bsType :: String , train :: Train , massModel )
2022-06-03 17:24:16 +02:00
# calculate resisting forces
2022-06-22 11:11:44 +02:00
supportPoint [ :R_traction ] = tractionUnitResistance ( supportPoint [ :v ] , train )
2022-06-05 15:25:12 +02:00
if train . transportType == :freight
2022-06-22 11:11:44 +02:00
supportPoint [ :R_wagons ] = freightWagonsResistance ( supportPoint [ :v ] , train )
2022-06-05 15:25:12 +02:00
elseif train . transportType == :passenger
2022-06-22 11:11:44 +02:00
supportPoint [ :R_wagons ] = passengerWagonsResistance ( supportPoint [ :v ] , train )
2022-06-05 15:25:12 +02:00
end
2022-06-22 11:11:44 +02:00
supportPoint [ :R_train ] = supportPoint [ :R_traction ] + supportPoint [ :R_wagons ]
supportPoint [ :R_path ] = calculatePathResistance ( CSs , csId , supportPoint [ :s ] , massModel , train )
supportPoint [ :F_R ] = supportPoint [ :R_train ] + supportPoint [ :R_path ]
2022-06-03 17:24:16 +02:00
# calculate tractive effort
2022-06-05 15:25:12 +02:00
if bsType == " braking " || bsType == " coasting " || bsType == " halt "
2022-06-22 11:11:44 +02:00
supportPoint [ :F_T ] = 0.0
2022-06-03 17:24:16 +02:00
elseif bsType == " cruising "
2022-06-22 11:11:44 +02:00
supportPoint [ :F_T ] = min ( max ( 0.0 , supportPoint [ :F_R ] ) , calculateTractiveEffort ( supportPoint [ :v ] , train . tractiveEffort ) )
2022-06-03 17:24:16 +02:00
else # bsType == "accelerating" || bsType == "diminishing" || 'default'
2022-06-22 11:11:44 +02:00
supportPoint [ :F_T ] = calculateTractiveEffort ( supportPoint [ :v ] , train . tractiveEffort )
2022-06-03 17:24:16 +02:00
end
2022-06-22 11:11:44 +02:00
return supportPoint
2022-06-03 17:24:16 +02:00
end #function calculateForces!
"""
TODO
"""
function moveAStep ( previousPoint :: Dict , stepVariable :: Symbol , stepSize :: Real , csId :: Integer )
# stepSize is the currentStepSize depending on the accessing function
# TODO: csId is only for error messages. Should it be removed?
#= 0 8 / 3 1 T O D O : H o w t o c h e c k i f t h e t r a i n s t o p p s d u r i n g t h i s s t e p ? I s h o u l d t h r o w a n e r r o r m y s e l f t h a t I c a t c h i n h i g h e r h i e r a r c h i e s . =#
2022-06-22 11:11:44 +02:00
# create the next support point
newPoint = SupportPoint ( )
2022-06-03 17:24:16 +02:00
# calculate s, t, v, E
if stepVariable == :distance # distance step method
2022-06-30 18:42:35 +02:00
Δs = stepSize # step size (in m)
2022-06-03 17:24:16 +02:00
if previousPoint [ :a ] == 0.0
if previousPoint [ :v ] == 0.0
error ( " ERROR: The train tries to cruise at v=0.0 m/s at s= " , previousPoint [ :s ] , " in CS " , csId , " . " )
end
2022-06-30 18:42:35 +02:00
Δt = Δt_with_constant_v ( Δs , previousPoint [ :v ] ) # step size (in s)
Δv = 0.0 # step size (in m/s)
2022-06-03 17:24:16 +02:00
else
2022-06-05 15:41:28 +02:00
# check if the parts of the following square roots will be <0.0 in the functions Δt_with_Δs and Δv_with_Δs
2022-06-30 18:42:35 +02:00
squareRootPartIsNegative = ( previousPoint [ :v ] / previousPoint [ :a ] ) ^ 2 + 2 * Δs / previousPoint [ :a ] < 0.0 || previousPoint [ :v ] ^ 2 + 2 * Δs * previousPoint [ :a ] < 0.0
2022-06-03 17:24:16 +02:00
if previousPoint [ :a ] < 0.0 && squareRootPartIsNegative
error ( " ERROR: The train stops during the accelerating section in CS " , csId , " because the tractive effort is lower than the resistant forces. " ,
" Before the stop the last point has the values s= " , previousPoint [ :s ] , " m, v= " , previousPoint [ :v ] , " m/s, a= " , previousPoint [ :a ] , " m/s^2, " ,
" F_T= " , previousPoint [ :F_T ] , " N, R_traction= " , previousPoint [ :R_traction ] , " N, R_wagons= " , previousPoint [ :R_wagons ] , " N, R_path= " , previousPoint [ :R_path ] , " N. " )
end
2022-06-30 18:42:35 +02:00
Δt = Δt_with_Δs ( Δs , previousPoint [ :a ] , previousPoint [ :v ] ) # step size (in s)
Δv = Δv_with_Δs ( Δs , previousPoint [ :a ] , previousPoint [ :v ] ) # step size (in m/s)
2022-06-03 17:24:16 +02:00
end
elseif stepVariable == :time # time step method
2022-06-30 18:42:35 +02:00
Δt = stepSize # step size (in s)
Δs = Δs_with_Δt ( Δt , previousPoint [ :a ] , previousPoint [ :v ] ) # step size (in m)
Δv = Δv_with_Δt ( Δt , previousPoint [ :a ] ) # step size (in m/s)
2022-06-03 17:24:16 +02:00
elseif stepVariable == :velocity # velocity step method
if previousPoint [ :a ] == 0.0
if previousPoint [ :v ] == 0.0
error ( " ERROR: The train tries to cruise at v=0.0 m/s at s= " , previousPoint [ :s ] , " in CS " , csId , " . " )
end
2022-06-30 18:42:35 +02:00
Δs = stepSize # step size (in m)
2022-06-03 17:24:16 +02:00
# TODO what is the best default step size for constant v? define Δs or Δt?
2022-06-30 18:42:35 +02:00
Δt = Δt_with_constant_v ( Δs , previousPoint [ :v ] ) # step size (in s)
Δv = 0.0 # step size (in m/s)
2022-06-03 17:24:16 +02:00
else
2022-06-30 18:42:35 +02:00
Δv = stepSize * sign ( previousPoint [ :a ] ) # step size (in m/s)
Δs = Δs_with_Δv ( Δv , previousPoint [ :a ] , previousPoint [ :v ] ) # step size (in m)
Δt = Δt_with_Δv ( Δv , previousPoint [ :a ] ) # step size (in s)
2022-06-03 17:24:16 +02:00
end
end #if
2022-06-30 18:42:35 +02:00
newPoint [ :s ] = previousPoint [ :s ] + Δs # position (in m)
newPoint [ :t ] = previousPoint [ :t ] + Δt # point in time (in s)
newPoint [ :v ] = previousPoint [ :v ] + Δv # velocity (in m/s)
2022-06-03 17:24:16 +02:00
return newPoint
end #function moveAStep
2022-06-05 15:25:12 +02:00
2022-06-03 17:24:16 +02:00
"""
# if the rear of the train is still located in a former characteristic section it has to be checked if its speed limit can be kept
"""
2022-08-17 12:18:27 +02:00
function getLowestSpeedLimit ( CSs :: Vector { Dict } , csWithTrainHeadId :: Integer , s :: Real , trainLength :: Real )
2022-06-03 17:24:16 +02:00
v_limit = CSs [ csWithTrainHeadId ] [ :v_limit ]
s_exit = CSs [ csWithTrainHeadId ] [ :s_exit ]
if csWithTrainHeadId > 1 && s - trainLength < CSs [ csWithTrainHeadId ] [ :s_entry ]
formerCsId = csWithTrainHeadId - 1
while formerCsId > 0 && s - trainLength < CSs [ formerCsId ] [ :s_exit ]
2022-07-12 16:59:45 +02:00
if CSs [ formerCsId ] [ :v_limit ] < v_limit # TODO: is the position of the train's rear < s_trainrun_start, v_limit of the first CS is used
2022-06-03 17:24:16 +02:00
v_limit = CSs [ formerCsId ] [ :v_limit ]
s_exit = CSs [ formerCsId ] [ :s_exit ]
end
formerCsId = formerCsId - 1
end
end
2022-08-17 12:18:27 +02:00
lowestSpeedLimit = Dict ( :v => v_limit , :s_end => s_exit + trainLength )
return lowestSpeedLimit
end #function getLowestSpeedLimit
2022-06-03 17:24:16 +02:00
2022-06-05 15:25:12 +02:00
2022-06-03 17:24:16 +02:00
"""
2022-06-05 15:25:12 +02:00
TODO
2022-06-03 17:24:16 +02:00
"""
2022-08-19 19:51:02 +02:00
function getNextPointOfInterest ( pointsOfInterest :: Vector { NamedTuple } , s :: Real )
2022-06-05 15:25:12 +02:00
for POI in pointsOfInterest
if POI [ 1 ] > s
return POI
2022-06-03 17:24:16 +02:00
end
end
error ( " ERROR in getNextPointOfInterest: There is no POI higher than s= " , s , " m. " )
end #function getNextPointOfInterest
2022-07-12 16:59:45 +02:00
## create vectors with the moving section's points of interest and with the characteristic sections with secured braking and accelerating behavior
2022-06-03 17:24:16 +02:00
function determineCharacteristics ( path :: Path , train :: Train , settings :: Settings )
2022-07-12 16:59:45 +02:00
# determine the positions of the points of interest depending on the interesting part of the train (front/rear) and the train's length
2022-08-19 19:51:02 +02:00
pointsOfInterest = NamedTuple [ ]
2022-07-12 16:59:45 +02:00
if ! isempty ( path . poi )
for POI in path . poi
s_poi = POI [ :station ]
if POI [ :measure ] == " rear "
s_poi += train . length
end
2022-08-19 19:51:02 +02:00
push! ( pointsOfInterest , ( s = s_poi , label = POI [ :label ] ) )
2022-07-12 16:59:45 +02:00
end
2022-08-19 19:51:02 +02:00
sort! ( pointsOfInterest , by = x -> x [ :s ] )
2022-07-12 16:59:45 +02:00
end
characteristicSections = CharacteristicSections ( path , train . v_limit , train . length , pointsOfInterest )
characteristicSections = secureBrakingBehavior! ( characteristicSections , train . a_braking , settings . approxLevel )
2022-06-03 17:24:16 +02:00
2022-07-12 16:59:45 +02:00
return ( characteristicSections , pointsOfInterest )
2022-06-03 17:24:16 +02:00
end #function determineCharacteristics