diff --git a/README.md b/README.md
index 5089d03..9bb0a4c 100644
--- a/README.md
+++ b/README.md
@@ -18,7 +18,7 @@ If you don't want to use Poetry, you can use install all packages listed in the
 You can use the `benchmark.py` script to solve instances. Here's an example:
 
 ``` shell
-poetry run python benchmark.py instances/*.json \
+poetry run python benchmark.py instances/n6-idx0-distribution0-travel0-serv1.json \
 --num_procs 4 \
 --weight_travel 1 \
 --weight_idle 2.5 \
@@ -28,6 +28,26 @@ poetry run python benchmark.py instances/*.json \
 --max_runtime 1
 ```
 
-This command solves all instances, four in parallel at a time. The objective weights are 1, 2.5 and 10 for the travel, idle and waiting time contributions, respectively. The selected algorithm is LNS with seed 1 and runs for one second.
+This command solves an instance with six clients with high service time variance. 
+The objective weights are 1, 2.5 and 10 for the travel, idle and waiting time, respectively. The selected algorithm is LNS with seed 1 and runs for one second.
 
 All results presented in our paper can be found in this repository. In particular, the raw results data is stored in `data/` and the notebooks that generate tables and figures can be found in `notebooks/`.
+
+
+## Citation
+
+Please consider citing our paper if this repository has been useful to you:
+
+``` bibtex
+@misc{Bekker2023,
+  title = {A Queueing-Based Approach for Integrated Routing and Appointment Scheduling},
+  author = {Bekker, Ren{\'e} and Bharti, Bharti and Lan, Leon and Mandjes, Michel},
+  year = {2023},
+  month = dec,
+  number = {arXiv:2312.02715},
+  eprint = {2312.02715},
+  primaryclass = {math},
+  publisher = {arXiv},
+  doi = {10.48550/arXiv.2312.02715},
+}
+```
diff --git a/benchmark.py b/benchmark.py
index 7722395..e2d131d 100644
--- a/benchmark.py
+++ b/benchmark.py
@@ -64,7 +64,7 @@ def parse_args():
     parser.add_argument("--weight_idle", type=float, default=2.5)
     parser.add_argument("--weight_wait", type=int, default=10)
     parser.add_argument("--fixed_weights", action="store_true")
-
+    parser.add_argument("--initialize_lns", action="store_true")
     parser.add_argument("--distance_scaling", type=float, default=1.0)
     parser.add_argument("--num_scenarios", type=int, default=100)
 
@@ -207,6 +207,7 @@ def solve(
     weight_wait: int,
     fixed_weights: bool,
     distance_scaling: float,
+    initialize_lns: bool,
     benchmark_runtime: Optional[bool],
     sol_dir: Optional[str],
     plot_dir: Optional[str],
@@ -242,15 +243,14 @@ def solve(
         kwargs["max_runtime"] = runtimes_by_size[data.dimension - 1]
 
     if algorithm == "lns":
-        visits = double_orientation_tsp(
-            data, data, cost_evaluator, **kwargs
-        ).solution.visits
+        extra = {}
+
+        if initialize_lns:
+            res = double_orientation_tsp(data, data, cost_evaluator, **kwargs)
+            extra["initial_visits"] = res.solution.visits
+
         result = large_neighborhood_search(
-            seed,
-            data,
-            cost_evaluator,
-            initial_visits=visits,
-            **kwargs,
+            seed, data, cost_evaluator, **extra, **kwargs
         )
     elif algorithm == "tsp":
         result = tsp(seed, data, cost_evaluator, **kwargs)