Simulating traffic around a lock#
In this notebook, we simulate a lock on a network which randomly generated vessels have to pass. We add a pre-coded complex lock object on the graph. Vessels are locked together if they can fit inside the lock, and arrive within the clustering time window.
0. Import libraries#
# package(s) used for creating and geo-locating the graph
import networkx as nx
import pyproj
from shapely.geometry import Point, LineString
from shapely.ops import transform
# package(s) related to the simulation (creating the vessel, running the simulation)
import datetime
import simpy
import opentnsim
from opentnsim.core.logutils import logbook2eventtable
from opentnsim.core.plotutils import generate_vessel_gantt_chart
from scipy.stats import norm, uniform, expon
# import of modules important for locking
from opentnsim.lock import lock as lock_module
from opentnsim.vessel_traffic_service import vessel_traffic_service as vessel_traffic_service_module
from opentnsim.graph import mixins as graph_module
# package(s) needed for inspecting the output
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
print("This notebook is executed with OpenTNSim version {}".format(opentnsim.__version__))
This notebook is executed with OpenTNSim version 1.3.7
1. Define object classes#
# make your preferred Vessel class out of available mix-ins.
Vessel = type(
"Vessel",
(
lock_module.PassesLockComplex, # allows to interact with a lock
opentnsim.core.Identifiable, # allows to give the object a name and a random ID,
opentnsim.core.Movable, # allows the object to move, with a fixed speed, while logging this activity
opentnsim.core.VesselProperties, # allows vessel to have dimensions, namely a length (L), width (B), and draught (T)
opentnsim.core.ExtraMetadata, # allow additional information, such as an arrival time (required for passing a lock)
opentnsim.mixins.HasMultiDiGraph, # allow to operate on a graph that can include parallel edges from and to the same nodes
opentnsim.output.HasOutput, # allow additional output to be stored
),
{}
)
2. Create graph#
# define reference systems
wgs84eqd = pyproj.CRS('4087')
wgs84rad = pyproj.CRS('4326')
# define transformer functions
wgs84eqd_to_wgs84rad = pyproj.transformer.Transformer.from_crs(wgs84eqd,wgs84rad,always_xy=True).transform #equidistant wgs84 to radial wgs84
wgs84rad_to_wgs84eqd = pyproj.transformer.Transformer.from_crs(wgs84rad,wgs84eqd,always_xy=True).transform #radial wgs84 to equidistant wgs84
# create a directed graph
graph = nx.DiGraph()
# add nodes
graph.add_node('-1',geometry=transform(wgs84eqd_to_wgs84rad,Point(-350600,0)))
graph.add_node('0',geometry=transform(wgs84eqd_to_wgs84rad,Point(-5000,0)))
graph.add_node('1',geometry=transform(wgs84eqd_to_wgs84rad,Point(5000,0)))
graph.add_node('+1',geometry=transform(wgs84eqd_to_wgs84rad,Point(350600,0)))
# add edges
graph.add_edge('-1','0', weight=1)
graph.add_edge('0','-1', weight=1)
graph.add_edge('0','1', weight=1)
graph.add_edge('1','0', weight=1)
graph.add_edge('1','+1', weight=1)
graph.add_edge('+1','1', weight=1)
graph_module.plot_graph(graph)
3. Run simulation#
def generate_vessel(
env,
name,
start_node,
end_node,
arrival_time,
vessel_speed=4,
vessel_length=100,
vessel_beam=20,
vessel_draft=10,
vessel_type="tanker"):
"""
Creates and returns a Vessel object with a computed route through the environment graph.
Parameters:
----------
env : Environment
The simulation environment containing the graph and other context.
name : str
Human readabile identifier for the vessel.
start_node : str or int
The starting node in the graph (converted to string).
end_node : str or int
The destination node in the graph (converted to string).
arrival_time : pd.Timestamp
The scheduled arrival time of the vessel at the start node.
vessel_speed : float, optional
Speed of the vessel in knots or simulation units (default is 4).
vessel_length : float, optional
Length of the vessel in meters (default is 100).
vessel_beam : float, optional
Beam (width) of the vessel in meters (default is 20).
vessel_draft : float, optional
Draught (depth below waterline) of the vessel in meters (default is 10).
vessel_type : str, optional
Type of vessel (e.g., "tanker", "cargo", "container") (default is "tanker").
Returns:
-------
Vessel or None
A Vessel object initialized with the given parameters and route.
Returns None if no valid path exists between start_node and end_node.
"""
# Ensure nodes are strings
start_node = str(start_node)
end_node = str(end_node)
try:
route = nx.dijkstra_path(env.graph, start_node, end_node)
except nx.NetworkXNoPath:
print(f"⚠️ No path from {start_node} to {end_node}. Vessel {name} not created.")
return None
geometry = env.graph.nodes[start_node]['geometry']
data_vessel = {
"env": env,
"name": name,
"geometry": geometry,
"route": route,
"v": vessel_speed,
"L": vessel_length,
"B": vessel_beam,
"T": vessel_draft,
"type": vessel_type,
"arrival_time": arrival_time,
}
vessel = Vessel(**data_vessel)
return vessel
def generate_vessels_with_distributions(
env,
num_vessels,
start_time,
arrival_dist_up=None,
arrival_dist_down=None,
seed_up=None,
seed_down=None):
"""
Generates a list of vessels with interarrival times drawn from specified distributions
for upward and downward directions. Supports independent seeding for reproducibility.
Parameters
----------
env : Environment
The simulation environment containing the graph and vessel context.
num_vessels : int
Total number of vessels to generate. Vessels alternate between up and down directions.
start_time : pd.Timestamp
The initial timestamp from which vessel arrivals begin.
arrival_dist_up : callable, optional
A function returning interarrival times (in minutes) for upward-moving vessels.
If None, defaults to an exponential distribution with mean 20 minutes.
arrival_dist_down : callable, optional
A function returning interarrival times (in minutes) for downward-moving vessels.
If None, defaults to an exponential distribution with mean 20 minutes.
seed_up : int or None, optional
Seed for the random number generator used in upward direction.
seed_down : int or None, optional
Seed for the random number generator used in downward direction.
Returns
-------
list of Vessel
A list of Vessel objects with assigned routes and arrival times.
Vessels for which no valid path exists are skipped.
"""
vessels = []
# Create independent random generators
rng_up = np.random.default_rng(seed_up)
rng_down = np.random.default_rng(seed_down)
# Default to exponential distribution with mean 20 minutes
if arrival_dist_up is None:
arrival_dist_up = lambda: rng_up.exponential(scale=20)
if arrival_dist_down is None:
arrival_dist_down = lambda: rng_down.exponential(scale=20)
up_time = start_time
down_time = start_time
for i in range(num_vessels):
if i % 2 == 0:
# Upward direction: -1 → +1
start_node, end_node = "-1", "+1"
delta_minutes = arrival_dist_up()
arrival_time = up_time + pd.Timedelta(minutes=delta_minutes)
up_time = arrival_time
else:
# Downward direction: +1 → -1
start_node, end_node = "+1", "-1"
delta_minutes = arrival_dist_down()
arrival_time = down_time + pd.Timedelta(minutes=delta_minutes)
down_time = arrival_time
vessel = generate_vessel(
env=env,
name=f"Vessel {i + 1}",
start_node=start_node,
end_node=end_node,
arrival_time=arrival_time
)
if vessel:
vessels.append(vessel)
return vessels
def mission(env, vessel):
"""
Method that defines the mission of the vessel.
In this case:
keep moving along the path until its end point is reached
"""
while True:
yield from vessel.move()
if vessel.geometry == nx.get_node_attributes(env.graph, "geometry")[vessel.route[-1]]:
break
# start simpy environment
simulation_start = datetime.datetime(2025, 1, 1, 0, 0, 0)
env = simpy.Environment(initial_time=simulation_start.timestamp())
env.epoch = simulation_start
# add graph to environment
env.graph = graph
# add components important for locking to the environment
env.vessel_traffic_service = vessel_traffic_service_module.VesselTrafficService(graph=graph)
lock = lock_module.IsLockComplex(
env=env,
name='Lock',
node_open='0',
node_A = '0',
node_B = '1',
distance_lock_doors_A_to_waiting_area_A = 4800,
distance_lock_doors_B_to_waiting_area_B = 4800,
distance_from_start_node_to_lock_doors_A = 4800,
distance_from_end_node_to_lock_doors_B = 4800,
lock_length = 400,
lock_width = 50,
lock_depth = 15,
levelling_time = 300,
sailing_distance_to_crossing_point = 1800,
doors_opening_time= 300,
doors_closing_time= 300,
speed_reduction_factor_lock_chamber=0.5,
sailing_in_time_gap_through_doors = 300,
sailing_in_speed_sea = 1.5,
sailing_in_speed_canal = 1.5,
sailing_out_time_gap_through_doors = 120,
sailing_time_before_opening_lock_doors = 600,
sailing_time_before_closing_lock_doors = 120,
registration_nodes = ['-1','+1'],
)
start_time = pd.Timestamp('2025-01-01 00:00:00')
# Example using numpy
# Create independent random generators
rng_up = np.random.default_rng(123)
rng_down = np.random.default_rng(456)
# Exponential distributions with different means (scale = mean)
arrival_dist_up = lambda: rng_up.exponential(scale=5) # mean 30 minutes
arrival_dist_down = lambda: rng_down.exponential(scale=5) # mean 15 minutes
# Example using scipy.stats and other distributions
# # Create independent random generators
# rng_up = np.random.default_rng(123)
# rng_down = np.random.default_rng(456)
# # Exponential distributions using SciPy
# arrival_dist_up = lambda: expon(scale=30).rvs(random_state=rng_up)
# arrival_dist_down = lambda: expon(scale=15).rvs(random_state=rng_down)
vessels = generate_vessels_with_distributions(
env=env,
num_vessels=11, #40
start_time=start_time,
arrival_dist_up=arrival_dist_up,
arrival_dist_down=arrival_dist_down,
seed_up=123,
seed_down=456
)
for vessel in vessels:
env.process(mission(env, vessel))
env.run()
4. Inspect output#
# load the logbook data into a dataframe
lock_df = pd.DataFrame.from_dict(lock.lock_chamber.logbook)
print("'{}' logbook data:".format(lock.name))
print('')
display(lock_df)
'Lock' logbook data:
| Message | Timestamp | Value | Geometry | |
|---|---|---|---|---|
| 0 | Lock doors closing start | 2025-01-01 00:02:27.392402 | {} | 0 |
| 1 | Lock doors closing stop | 2025-01-01 00:07:27.392402 | {} | 0 |
| 2 | Lock chamber converting start | 2025-01-01 00:07:27.392402 | {} | 0 |
| 3 | Lock chamber converting stop | 2025-01-01 00:12:27.392402 | {} | 1 |
| 4 | Lock doors opening start | 2025-01-01 00:12:27.392402 | {} | 1 |
| 5 | Lock doors opening stop | 2025-01-01 00:17:27.392402 | {} | 1 |
| 6 | Lock doors closing start | 2025-01-02 00:38:56.161300 | {} | 1 |
| 7 | Lock doors closing stop | 2025-01-02 00:43:56.161300 | {} | 1 |
| 8 | Lock chamber converting start | 2025-01-02 00:43:56.161300 | {} | 1 |
| 9 | Lock chamber converting stop | 2025-01-02 00:48:56.161300 | {} | 0 |
| 10 | Lock doors opening start | 2025-01-02 00:48:56.161300 | {} | 0 |
| 11 | Lock doors opening stop | 2025-01-02 00:53:56.161300 | {} | 0 |
| 12 | Lock doors closing start | 2025-01-02 01:35:28.000000 | {} | 0 |
| 13 | Lock doors closing stop | 2025-01-02 01:40:28.000000 | {} | 0 |
| 14 | Lock chamber converting start | 2025-01-02 01:40:28.000000 | {} | 0 |
| 15 | Lock chamber converting stop | 2025-01-02 01:45:28.000000 | {} | 1 |
| 16 | Lock doors opening start | 2025-01-02 01:45:28.000000 | {} | 1 |
| 17 | Lock doors opening stop | 2025-01-02 01:50:28.000000 | {} | 1 |
| 18 | Lock doors closing start | 2025-01-02 02:21:11.064110 | {} | 1 |
| 19 | Lock doors closing stop | 2025-01-02 02:26:11.064110 | {} | 1 |
| 20 | Lock chamber converting start | 2025-01-02 02:26:11.064110 | {} | 1 |
| 21 | Lock chamber converting stop | 2025-01-02 02:31:11.064110 | {} | 0 |
| 22 | Lock doors opening start | 2025-01-02 02:31:11.064110 | {} | 0 |
| 23 | Lock doors opening stop | 2025-01-02 02:36:11.064110 | {} | 0 |
| 24 | Lock doors closing start | 2025-01-02 03:01:43.000000 | {} | 0 |
| 25 | Lock doors closing stop | 2025-01-02 03:06:43.000000 | {} | 0 |
| 26 | Lock chamber converting start | 2025-01-02 03:06:43.000000 | {} | 0 |
| 27 | Lock chamber converting stop | 2025-01-02 03:11:43.000000 | {} | 1 |
| 28 | Lock doors opening start | 2025-01-02 03:11:43.000000 | {} | 1 |
| 29 | Lock doors opening stop | 2025-01-02 03:16:43.000000 | {} | 1 |
Gantt chart of event table#
df_eventtable = opentnsim.core.logutils.logbook2eventtable([*vessels, lock.lock_chamber])
fig = generate_vessel_gantt_chart(df_eventtable)
Time-distance diagram of vessels passing the lock and planning info#
def cm_to_pixels(cm):
return cm * 37.8 # Set figure height to 10 cmfig.update_layout(height=cm_to_pixels(10))
# We can plot the time-distance diagram
fig = lock.create_time_distance_plot(vessels = vessels,
xlimmin = -6050,
xlimmax = 6050,
ylimmin = pd.Timestamp('2025-01-01 22:00:00'),
ylimmax = pd.Timestamp('2025-01-02 09:00:00'),
method='Plotly')
fig.update_layout(height=cm_to_pixels(20))
Vessel delays: individual delays and overall average#
delays = []
for vessel in vessels:
vessel_df = pd.DataFrame(vessel.logbook)
waiting_stop = vessel_df[vessel_df.Message == "Waiting stop"]
if not waiting_stop.empty:
delay = (waiting_stop.Timestamp-vessel.metadata["arrival_time"]).iloc[0]
else:
delay = pd.Timedelta(seconds=0)
delays.append(delay)
print(f"The average vessel delay is {np.round(np.average(delays).total_seconds()/60,1)} minutes")
The average vessel delay is 0.0 minutes
Simulated intensity (vessels per hour)#
‘get_vessels_during_leveling’:
Identify locking cycles (by looking at the lock logbook)
From the vessel list identify which vessels were in the lock during that locking cycle
‘calculate_cycle_looptimes’:
Using the info derived from ‘get_vessels_during_leveling’ calculate looptimes
‘calculate_detailed_cycle_time’:
using the info from ‘get_vessels_during_leveling’ and ‘calculate_cycle_looptimes’ calculate detailed cycle times
leveling_cycles = opentnsim.lock.logutils.get_vessels_during_leveling(lock.lock_chamber, vessels)
looptimes_df = opentnsim.lock.logutils.calculate_cycle_looptimes(leveling_cycles, vessels)
Tc_df = opentnsim.lock.logutils.calculate_detailed_cycle_time(lock.lock_chamber, vessels, leveling_cycles)
display(pd.DataFrame(leveling_cycles))
display(looptimes_df)
display(Tc_df)
| leveling_start | leveling_stop | vessels_present | |
|---|---|---|---|
| 0 | 2025-01-01 00:07:27.392402 | 2025-01-01 00:12:27.392402 | [] |
| 1 | 2025-01-02 00:43:56.161300 | 2025-01-02 00:48:56.161300 | [Vessel 2, Vessel 4, Vessel 6, Vessel 8] |
| 2 | 2025-01-02 01:40:28.000000 | 2025-01-02 01:45:28.000000 | [Vessel 1, Vessel 3, Vessel 5, Vessel 7] |
| 3 | 2025-01-02 02:26:11.064110 | 2025-01-02 02:31:11.064110 | [Vessel 10] |
| 4 | 2025-01-02 03:06:43.000000 | 2025-01-02 03:11:43.000000 | [Vessel 9, Vessel 11] |
| cycle | looptime_seconds | |
|---|---|---|
| 0 | 1 | 0.000000 |
| 1 | 2 | NaN |
| 2 | 3 | 899.999999 |
| 3 | 4 | 899.910762 |
| 4 | 5 | 899.999999 |
| t_l_up | sum_t_i_up | T_close_up | T_waterlevel_up | T_open_up | sum_t_u_up | t_l_down | sum_t_i_down | T_close_down | T_waterlevel_down | T_open_down | sum_t_u_down | Tc_seconds | up_vessels | down_vessels | I_s | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.000000 | 0.000000 | 300.0 | 300.0 | 300.0 | 0.000000 | 0.000000 | 988.768898 | 300.0 | 300.0 | 300.0 | 602.980563 | 3391.749461 | [] | [Vessel 2, Vessel 4, Vessel 6, Vessel 8] | 4.245597 |
| 1 | 899.999999 | 988.768899 | 300.0 | 300.0 | 300.0 | 602.980561 | 899.910762 | 340.172787 | 300.0 | 300.0 | 300.0 | 48.596112 | 5580.429120 | [Vessel 1, Vessel 3, Vessel 5, Vessel 7] | [Vessel 10] | 3.225558 |
for index, row in Tc_df.iterrows():
print('Locking cycle {} has an intensity of {:.2f} vessels per hour'.format(index+1, row['I_s']))
Locking cycle 1 has an intensity of 4.25 vessels per hour
Locking cycle 2 has an intensity of 3.23 vessels per hour
Estimated capacity (vessels per hour)#
n_max = 4
vessel_speed_outside_of_lock = 4
# Part III, Ch3, Eq. 3.2 (NB: de helft van de looptime wordt hier effectief geimplementeerd door de sailing to lock te berekenen)
t_sailing_to_lock = lock.sailing_distance_to_crossing_point/vessel_speed_outside_of_lock
T_entering = t_sailing_to_lock + (n_max-1)*lock.sailing_in_time_gap_through_doors + 50/2
# Part III, Ch3, Eq. 3.3
T_operation = lock.doors_closing_time + lock.lock_chamber.levelling_time + lock.doors_opening_time
# Part III, Ch3, Eq. 3.4 (NB: de helft van de looptime wordt hier effectief geimplementeerd door de sailing out of lock te berekenen)
t_sailing_out_of_lock = lock.sailing_distance_to_crossing_point/vessel_speed_outside_of_lock
T_exiting = t_sailing_out_of_lock + (n_max-1)*lock.sailing_out_time_gap_through_doors + 350/2
# Part III, Ch3, Eq. 3.1
T_locking = T_entering + T_operation + T_exiting
T_c = 2 * T_locking
C_s = 2*n_max / (T_c/3600)
print(f"The capacity of the lock is {np.round(C_s,1)} vessels per hour")
The capacity of the lock is 4.4 vessels per hour
lock.vessel_planning.iloc[-1]
id fa18838f-3930-4406-8e7e-e3a2a828def0
node_from 1
node_to 0
lock_chamber Lock
L 100
B 20
T 10
operation_index 3
time_of_registration 2025-01-01 00:21:11.223921
time_of_acceptance 2025-01-01 00:21:11.223921
time_potential_lock_door_opening_stop 2025-01-02 02:05:30.891323172
time_arrival_at_waiting_area 2025-01-02 00:21:11.223921
time_arrival_at_lineup_area NaN
time_lock_passing_start 2025-01-02 02:08:00.891323172
time_lock_entry_start 2025-01-02 02:15:30.891323172
time_lock_entry_stop 2025-01-02 02:21:11.064109642
time_lock_departure_start 2025-01-02 02:36:11.064109642
time_lock_departure_stop 2025-01-02 02:36:59.660221994
time_lock_passing_stop 2025-01-02 02:44:29.660221994
time_potential_lock_door_closure_start 2025-01-02 02:21:11.064109642
direction 1.0
delay 0 days 01:34:19.667402172
Name: 10, dtype: object
lock.operation_planning.iloc[2]
node_from 0
node_to 1
direction 0
lock_chamber Lock
vessels [<__main__.Vessel object at 0x7f3d4d3d9a90>, <...
capacity_L 0
capacity_B -30
time_potential_lock_door_opening_stop 2025-01-02 01:08:59.141862586
time_operation_start 2025-01-02 01:11:29.141862586
time_entry_start 2025-01-02 01:18:59.141862586
time_entry_stop 2025-01-02 01:35:27.910761408
time_door_closing_start 2025-01-02 01:35:27.910761408
time_door_closing_stop 2025-01-02 01:40:27.910761408
time_levelling_start 2025-01-02 01:40:27.910761408
time_levelling_stop 2025-01-02 01:45:27.910761408
time_door_opening_start 2025-01-02 01:45:27.910761408
time_door_opening_stop 2025-01-02 01:50:27.910761408
time_departure_start 2025-01-02 01:50:27.910761408
time_departure_stop 2025-01-02 02:00:30.891323172
time_operation_stop 2025-01-02 02:08:00.891323172
time_potential_lock_door_closure_start 2025-01-02 02:02:30.891323172
wlev_A NaN
wlev_B NaN
maximum_individual_delay 0 days 01:06:48.646226938
total_delay 0 days 01:58:00.050114584
status unavailable
Name: 2, dtype: object
lock.operation_planning.iloc[3]
node_from 1
node_to 0
direction 1
lock_chamber Lock
vessels [<__main__.Vessel object at 0x7f3d4d40de50>]
capacity_L 300
capacity_B 30
time_potential_lock_door_opening_stop 2025-01-02 01:55:30.891323172
time_operation_start 2025-01-02 01:55:30.891323172
time_entry_start 2025-01-02 02:15:30.891323172
time_entry_stop 2025-01-02 02:21:11.064109642
time_door_closing_start 2025-01-02 02:21:11.064109642
time_door_closing_stop 2025-01-02 02:26:11.064109642
time_levelling_start 2025-01-02 02:26:11.064109642
time_levelling_stop 2025-01-02 02:31:11.064109642
time_door_opening_start 2025-01-02 02:31:11.064109642
time_door_opening_stop 2025-01-02 02:36:11.064109642
time_departure_start 2025-01-02 02:36:11.064109642
time_departure_stop 2025-01-02 02:36:59.660221994
time_operation_stop 2025-01-02 02:44:29.660221994
time_potential_lock_door_closure_start 2025-01-02 02:38:59.660221994
wlev_A NaN
wlev_B NaN
maximum_individual_delay 0 days 01:34:19.667402172
total_delay 0 days 01:34:19.667402172
status unavailable
Name: 3, dtype: object