diff --git a/docs/vrp_variants.rst b/docs/vrp_variants.rst
index 5118e91..92a23f9 100644
--- a/docs/vrp_variants.rst
+++ b/docs/vrp_variants.rst
@@ -192,6 +192,8 @@ where each item of the list is the maximum load per vehicle type. For example, i
 Note how the dimensions of ``load_capacity`` and ``cost`` are consistent: each list must have as many items as vehicle types, and the
 order of the items of the ``load_capacity`` list is consistent with the order of the ``cost`` list on every edge of the graph.
   
+Once the problem is solved, the type of vehicle per route can be queried with ``prob.best_routes_type``.  
+
 	
 VRP options
 ~~~~~~~~~~~
@@ -254,7 +256,20 @@ will add to the total travel cost each time a node is dropped. For example, if t
 	>>> prob.drop_penalty = 1000
 	
 This problem is sometimes referred to as the `capacitated profitable tour problem` or the `prize collecting tour problem.`
+
+Minimizing the global time span
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+It is possible to modify the objective function in order to solve a min-max problem. More specifically, the total time span can be minimized
+by setting the ``minimize_global_span`` to ``True``. Of course this assumes edges have a ``time`` argument:
+
+.. code-block:: python
+
+	>>> prob.minimize_global_span = True
 	
+.. note::
+
+   This may lead to poor computation times.
 	
 Other VRPs
 ~~~~~~~~~~
diff --git a/examples/cvrp_drop.py b/examples/cvrp_drop.py
index 037c011..7e732d2 100644
--- a/examples/cvrp_drop.py
+++ b/examples/cvrp_drop.py
@@ -17,14 +17,20 @@
 
 if __name__ == "__main__":
 
-    prob = VehicleRoutingProblem(G,
-                                 load_capacity=15,
-                                 drop_penalty=1000,
-                                 num_vehicles=4)
-    prob.solve()
+    prob = VehicleRoutingProblem(G, load_capacity=15, drop_penalty=1000, num_vehicles=4)
+    prob.solve(
+        preassignments=[  # locking these routes should yield prob.best_value == 7936
+            # [9, 14, 16],
+            # [12, 11, 4, 3, 1],
+            # [7, 13],
+            # [8, 10, 2, 5],
+        ],
+    )
     print(prob.best_value)
     print(prob.best_routes)
     print(prob.best_routes_cost)
     print(prob.best_routes_load)
     print(prob.node_load)
-    assert prob.best_value == 7548
+    assert prob.best_value == 8096
+
+    # why doesn't vrpy find 7936 ?
diff --git a/tests/test_issue110.py b/tests/test_issue110.py
new file mode 100644
index 0000000..e57723a
--- /dev/null
+++ b/tests/test_issue110.py
@@ -0,0 +1,19 @@
+from networkx import DiGraph
+from vrpy import VehicleRoutingProblem
+
+
+class TestIssue110:
+    def setup(self):
+        G = DiGraph()
+        G.add_edge("Source", 1, cost=[1, 2])
+        G.add_edge("Source", 2, cost=[2, 4])
+        G.add_edge(1, "Sink", cost=[0, 0])
+        G.add_edge(2, "Sink", cost=[2, 4])
+        G.add_edge(1, 2, cost=[1, 2])
+        G.nodes[1]["demand"] = 13
+        G.nodes[2]["demand"] = 13
+        self.prob = VehicleRoutingProblem(G, mixed_fleet=True, load_capacity=[10, 15])
+
+    def test_node_load(self):
+        self.prob.solve()
+        assert self.prob.best_routes_type == {1: 1, 2: 1}
diff --git a/tests/test_toy.py b/tests/test_toy.py
index ace24fb..8e23278 100644
--- a/tests/test_toy.py
+++ b/tests/test_toy.py
@@ -7,7 +7,6 @@
 
 
 class TestsToy:
-
     def setup(self):
         """
         Creates a toy graph.
@@ -58,10 +57,7 @@ def test_cspy_stops_capacity_duration(self):
         """Tests column generation procedure on toy graph
         with stop, capacity and duration constraints
         """
-        prob = VehicleRoutingProblem(self.G,
-                                     num_stops=3,
-                                     load_capacity=10,
-                                     duration=62)
+        prob = VehicleRoutingProblem(self.G, num_stops=3, load_capacity=10, duration=62)
         prob.solve(exact=False)
         assert prob.best_value == 85
         assert set(prob.best_routes_duration.values()) == {41, 62}
@@ -82,8 +78,7 @@ def test_cspy_stops_time_windows(self):
         assert prob.arrival_time[1]["Sink"] in [41, 62]
 
     def test_cspy_schedule(self):
-        """Tests if final schedule is time-window feasible
-        """
+        """Tests if final schedule is time-window feasible"""
         prob = VehicleRoutingProblem(
             self.G,
             num_stops=3,
@@ -93,13 +88,13 @@ def test_cspy_schedule(self):
         # Check departure times
         for k1, v1 in prob.departure_time.items():
             for k2, v2 in v1.items():
-                assert (self.G.nodes[k2]["lower"] <= v2)
-                assert (v2 <= self.G.nodes[k2]["upper"])
+                assert self.G.nodes[k2]["lower"] <= v2
+                assert v2 <= self.G.nodes[k2]["upper"]
         # Check arrival times
         for k1, v1 in prob.arrival_time.items():
             for k2, v2 in v1.items():
-                assert (self.G.nodes[k2]["lower"] <= v2)
-                assert (v2 <= self.G.nodes[k2]["upper"])
+                assert self.G.nodes[k2]["lower"] <= v2
+                assert v2 <= self.G.nodes[k2]["upper"]
 
     ###############
     # subsolve lp #
@@ -151,13 +146,13 @@ def test_LP_schedule(self):
         # Check departure times
         for k1, v1 in prob.departure_time.items():
             for k2, v2 in v1.items():
-                assert (self.G.nodes[k2]["lower"] <= v2)
-                assert (v2 <= self.G.nodes[k2]["upper"])
+                assert self.G.nodes[k2]["lower"] <= v2
+                assert v2 <= self.G.nodes[k2]["upper"]
         # Check arrival times
         for k1, v1 in prob.arrival_time.items():
             for k2, v2 in v1.items():
-                assert (self.G.nodes[k2]["lower"] <= v2)
-                assert (v2 <= self.G.nodes[k2]["upper"])
+                assert self.G.nodes[k2]["lower"] <= v2
+                assert v2 <= self.G.nodes[k2]["upper"]
 
     def test_LP_stops_elementarity(self):
         """Tests column generation procedure on toy graph"""
@@ -188,11 +183,9 @@ def test_clarke_wright(self):
     #########
 
     def test_all(self):
-        prob = VehicleRoutingProblem(self.G,
-                                     num_stops=3,
-                                     time_windows=True,
-                                     duration=63,
-                                     load_capacity=10)
+        prob = VehicleRoutingProblem(
+            self.G, num_stops=3, time_windows=True, duration=63, load_capacity=10
+        )
         prob.solve(cspy=False)
         lp_best = prob.best_value
         prob.solve(cspy=True)
@@ -217,9 +210,7 @@ def test_knapsack(self):
 
     def test_pricing_strategies(self):
         sol = []
-        for strategy in [
-                "Exact", "BestPaths", "BestEdges1", "BestEdges2", "Hyper"
-        ]:
+        for strategy in ["Exact", "BestPaths", "BestEdges1", "BestEdges2", "Hyper"]:
             prob = VehicleRoutingProblem(self.G, num_stops=4)
             prob.solve(pricing_strategy=strategy)
             sol.append(prob.best_value)
@@ -294,22 +285,25 @@ def test_fixed_cost(self):
         assert set(prob.best_routes_cost.values()) == {30 + 100, 40 + 100}
 
     def test_drop_nodes(self):
-        prob = VehicleRoutingProblem(self.G,
-                                     num_stops=3,
-                                     num_vehicles=1,
-                                     drop_penalty=100)
+        prob = VehicleRoutingProblem(
+            self.G, num_stops=3, num_vehicles=1, drop_penalty=100
+        )
         prob.solve()
         assert prob.best_value == 240
         assert prob.best_routes == {1: ["Source", 1, 2, 3, "Sink"]}
 
     def test_num_vehicles_use_all(self):
-        prob = VehicleRoutingProblem(self.G,
-                                     num_stops=3,
-                                     num_vehicles=2,
-                                     use_all_vehicles=True,
-                                     drop_penalty=100)
+        prob = VehicleRoutingProblem(
+            self.G, num_stops=3, num_vehicles=2, use_all_vehicles=True, drop_penalty=100
+        )
         prob.solve()
         assert len(prob.best_routes) == 2
+        prob.num_vehicles = 3
+        prob.solve()
+        assert len(prob.best_routes) == 3
+        prob.num_vehicles = 4
+        prob.solve()
+        assert len(prob.best_routes) == 4
 
     def test_periodic(self):
         self.G.nodes[2]["frequency"] = 2
@@ -322,10 +316,7 @@ def test_periodic(self):
                 frequency += 1
         assert frequency == 2
         assert prob.schedule[0] in [[1], [1, 2]]
-        prob = VehicleRoutingProblem(self.G,
-                                     num_stops=2,
-                                     periodic=2,
-                                     num_vehicles=1)
+        prob = VehicleRoutingProblem(self.G, num_stops=2, periodic=2, num_vehicles=1)
         prob.solve()
         assert prob.schedule == {}
 
diff --git a/vrpy/master_solve_pulp.py b/vrpy/master_solve_pulp.py
index f43f138..c9988c3 100644
--- a/vrpy/master_solve_pulp.py
+++ b/vrpy/master_solve_pulp.py
@@ -46,6 +46,7 @@ def solve(self, relax, time_limit):
         logger.debug("master problem relax %s" % relax)
         logger.debug("Status: %s" % pulp.LpStatus[self.prob.status])
         logger.debug("Objective: %s" % pulp.value(self.prob.objective))
+        # self.prob.writeLP("master.lp")
 
         if pulp.LpStatus[self.prob.status] != "Optimal":
             raise Exception("problem " + str(pulp.LpStatus[self.prob.status]))
@@ -194,8 +195,6 @@ def _solve(self, relax: bool, time_limit: Optional[int]):
                 if "artificial_bound_" in var.name:
                     var.upBound = 0
                     var.lowBound = 0
-            # self.prob.writeLP("master.lp")
-            # print(self.prob)
         # Solve with appropriate solver
         if self.solver == "cbc":
             self.prob.solve(
@@ -210,7 +209,8 @@ def _solve(self, relax: bool, time_limit: Optional[int]):
                 pulp.CPLEX_CMD(
                     msg=False,
                     timeLimit=time_limit,
-                    options=["set lpmethod 4", "set barrier crossover -1"],
+                    # options=["set lpmethod 4", "set barrier crossover -1"], # set barrier crossover -1 is deprecated
+                    options=["set lpmethod 4", "set solutiontype 2"],
                 )
             )
         elif self.solver == "gurobi":
@@ -360,7 +360,7 @@ def _add_drop_variables(self):
         drop[v] takes value 1 if and only if node v is dropped.
         """
         for node in self.G.nodes():
-            if self.G.nodes[node]["demand"] > 0 and node != "Source":
+            if node not in ["Source", "Sink"]:
                 self.drop[node] = pulp.LpVariable(
                     "drop_%s" % node,
                     lowBound=0,
diff --git a/vrpy/masterproblem.py b/vrpy/masterproblem.py
index 8b19186..1ddde88 100644
--- a/vrpy/masterproblem.py
+++ b/vrpy/masterproblem.py
@@ -7,7 +7,7 @@ class _MasterProblemBase:
         routes (list): Current routes/variables/columns.
         drop_penalty (int, optional): Value of penalty if node is dropped. Defaults to None.
         num_vehicles (int, optional): Maximum number of vehicles. Defaults to None.
-        use_all_vehicles (bool, optional): True if all vehicles specified by num_vehicles should be used. Defaults to False
+        use_all_vehicles (bool, optional): True if all vehicles specified by num_vehicles should be used. Defaults to False.
         periodic (bool, optional): True if vertices are to be visited periodically. Defaults to False.
         minimize_global_span (bool, optional): True if global span (maximum distance) is minimized. Defaults to False.
         solver (str): Name of solver to use.
diff --git a/vrpy/vrp.py b/vrpy/vrp.py
index bdd1b3a..2d27e39 100644
--- a/vrpy/vrp.py
+++ b/vrpy/vrp.py
@@ -281,23 +281,25 @@ def best_routes_cost(self):
         """Returns dict with route ids as keys and route costs as values."""
         cost = {}
         for route in self.best_routes:
-            edges = list(
-                zip(self.best_routes[route][:-1], self.best_routes[route][1:]))
+            edges = list(zip(self.best_routes[route][:-1], self.best_routes[route][1:]))
             k = self._best_routes_vehicle_type[route]
-            cost[route] = sum(
-                self._H.edges[i, j]["cost"][k] for (i, j) in edges)
+            cost[route] = sum(self._H.edges[i, j]["cost"][k] for (i, j) in edges)
         return cost
 
     @property
     def best_routes_load(self):
         """Returns dict with route ids as keys and route loads as values."""
         load = {}
-        if (not self.load_capacity or self.distribution_collection or
-                self.pickup_delivery):
+        if (
+            not self.load_capacity
+            or self.distribution_collection
+            or self.pickup_delivery
+        ):
             return load
         for route in self.best_routes:
             load[route] = sum(
-                self._H.nodes[v]["demand"] for v in self.best_routes[route])
+                self._H.nodes[v]["demand"] for v in self.best_routes[route]
+            )
         return load
 
     @property
@@ -310,8 +312,11 @@ def node_load(self):
         If simultaneous distribution and collection, load refers to truck load when leaving the node.
         """
         load = {}
-        if (not self.load_capacity and not self.pickup_delivery and
-                not self.distribution_collection):
+        if (
+            not self.load_capacity
+            and not self.pickup_delivery
+            and not self.distribution_collection
+        ):
             return load
         for i in self.best_routes:
             load[i] = {}
@@ -323,8 +328,9 @@ def node_load(self):
                 amount += self._H.nodes[v]["demand"]
                 if self.distribution_collection:
                     # in this case, load = Q + collect - demand
-                    amount += (self._H.nodes[v]["collect"] -
-                               2 * self._H.nodes[v]["demand"])
+                    amount += (
+                        self._H.nodes[v]["collect"] - 2 * self._H.nodes[v]["demand"]
+                    )
                 load[i][v] = amount
         return load
 
@@ -335,14 +341,13 @@ def best_routes_duration(self):
         if not self.duration and not self.time_windows:
             return duration
         for route in self.best_routes:
-            edges = list(
-                zip(self.best_routes[route][:-1], self.best_routes[route][1:]))
+            edges = list(zip(self.best_routes[route][:-1], self.best_routes[route][1:]))
             # Travel times
-            duration[route] = sum(
-                self._H.edges[i, j]["time"] for (i, j) in edges)
+            duration[route] = sum(self._H.edges[i, j]["time"] for (i, j) in edges)
             # Service times
-            duration[route] += sum(self._H.nodes[v]["service_time"]
-                                   for v in self.best_routes[route])
+            duration[route] += sum(
+                self._H.nodes[v]["service_time"] for v in self.best_routes[route]
+            )
 
         return duration
 
@@ -362,10 +367,15 @@ def arrival_time(self):
             for j in range(1, len(route)):
                 tail = route[j - 1]
                 head = route[j]
-                arrival[i][head] = max(
-                    arrival[i][tail] + self._H.nodes[tail]["service_time"] +
-                    self._H.edges[tail, head]["time"],
-                    self._H.nodes[head]["lower"])
+                arrival[i][head] = min(
+                    max(
+                        arrival[i][tail]
+                        + self._H.nodes[tail]["service_time"]
+                        + self._H.edges[tail, head]["time"],
+                        self._H.nodes[head]["lower"],
+                    ),
+                    self._H.nodes[head]["upper"],
+                )
         return arrival
 
     @property
@@ -387,7 +397,8 @@ def departure_time(self):
                 head = route[j]
                 departure[i][head] = min(
                     arrival[i][head] + self._H.nodes[head]["service_time"],
-                    self._H.nodes[head]["upper"])
+                    self._H.nodes[head]["upper"],
+                )
         return departure
 
     @property
@@ -431,11 +442,12 @@ def _pre_solve(self):
             G=self.G,
         )
         # Check feasibility
-        check_feasibility(load_capacity=self.load_capacity,
-                          G=self.G,
-                          duration=self.duration)
+        check_feasibility(
+            load_capacity=self.load_capacity, G=self.G, duration=self.duration
+        )
         self.num_vehicles, self._num_vehicles_schedule = check_periodic_num_vehicles(
-            periodic=self.periodic, num_vehicles=self.num_vehicles)
+            periodic=self.periodic, num_vehicles=self.num_vehicles
+        )
         # Lock preassigned routes
         if self._preassignments:
             check_preassignments(self._preassignments, self.G)
@@ -444,11 +456,15 @@ def _pre_solve(self):
         self._prune_graph()
         # Compute upper bound on number of stops as knapsack problem
         if self.load_capacity and not self.pickup_delivery:
-            num_stops = get_num_stops_upper_bound(self.G, self._max_capacity,
-                                                  self.num_stops,
-                                                  self.distribution_collection)
+            num_stops = get_num_stops_upper_bound(
+                self.G, self._max_capacity, self.num_stops, self.distribution_collection
+            )
             self.num_stops = num_stops
             logger.info("new upper bound : max num stops = %s" % num_stops)
+        # Set empty route cost if use_all_vehicles = True
+        if self.use_all_vehicles:
+            n_vehicle_types = len(self.G["Source"]["Sink"]["cost"])
+            self.G["Source"]["Sink"]["cost"] = n_vehicle_types * [1e10]
 
     def _initialize(self, solver):
         """Initialization with feasible solution."""
@@ -486,7 +502,8 @@ def _solve(self, dive, solver):
         elif len(self.G.nodes()) > 2:
             # Solve as MIP
             _, _ = self.masterproblem.solve(
-                relax=False, time_limit=self._get_time_remaining(mip=True))
+                relax=False, time_limit=self._get_time_remaining(mip=True)
+            )
             (
                 self._best_value,
                 self._best_routes_as_graphs,
@@ -498,13 +515,16 @@ def _column_generation(self):
             # Generate good columns
             self._find_columns()
             # Stop if time limit is passed
-            if (isinstance(self._get_time_remaining(), float) and
-                    self._get_time_remaining() == 0.0):
+            if (
+                isinstance(self._get_time_remaining(), float)
+                and self._get_time_remaining() == 0.0
+            ):
                 logger.info("time up !")
                 break
             # Stop if no improvement limit is passed or max iter exceeded
             if self._no_improvement > 1000 or (
-                    self._max_iter and self._iteration >= self._max_iter):
+                self._max_iter and self._iteration >= self._max_iter
+            ):
                 break
 
     def _find_columns(self):
@@ -512,38 +532,46 @@ def _find_columns(self):
         # Solve restricted relaxed master problem
         if self._dive:
             duals, relaxed_cost = self.masterproblem.solve_and_dive(
-                time_limit=self._get_time_remaining())
+                time_limit=self._get_time_remaining()
+            )
             if self.hyper_heuristic:
                 self.hyper_heuristic.init(relaxed_cost)
         else:
             duals, relaxed_cost = self.masterproblem.solve(
-                relax=True, time_limit=self._get_time_remaining())
+                relax=True, time_limit=self._get_time_remaining()
+            )
         logger.info("iteration %s, %.6s" % (self._iteration, relaxed_cost))
         pricing_strategy = self._get_next_pricing_strategy(relaxed_cost)
 
         # One subproblem per vehicle type
         for vehicle in range(self._vehicle_types):
             # Solve pricing problem with randomised greedy algorithm
-            if (self._greedy and not self.time_windows and
-                    not self.distribution_collection and
-                    not self.pickup_delivery and not self.minimize_global_span):
+            if (
+                self._greedy
+                and not self.time_windows
+                and not self.distribution_collection
+                and not self.pickup_delivery
+                and not self.minimize_global_span
+            ):
                 subproblem = self._def_subproblem(duals, vehicle, greedy=True)
                 self.routes, self._more_routes = subproblem.solve(n_runs=20)
                 # Add initial_routes
                 if self._more_routes:
-                    for r in (r for r in self.routes
-                              if r.graph["name"] not in self.masterproblem.y):
+                    for r in (
+                        r
+                        for r in self.routes
+                        if r.graph["name"] not in self.masterproblem.y
+                    ):
                         self.masterproblem.update(r)
 
             # Continue searching for columns
             self._more_routes = False
             if not self.minimize_global_span:
                 self._more_routes = self._solve_subproblem_with_heuristic(
-                    pricing_strategy=pricing_strategy,
-                    vehicle=vehicle,
-                    duals=duals)
+                    pricing_strategy=pricing_strategy, vehicle=vehicle, duals=duals
+                )
             else:
-                for route in self._routes[:len(self._routes)]:
+                for route in self._routes[: len(self._routes)]:
                     self._more_routes = self._solve_subproblem_with_heuristic(
                         pricing_strategy=pricing_strategy,
                         vehicle=vehicle,
@@ -580,39 +608,39 @@ def _solve_subproblem_with_heuristic(
         more_columns = False
         if self._pricing_strategy == "Hyper":
             if pricing_strategy == "BestPaths":
-                more_columns = self._attempt_solve_best_paths(vehicle=vehicle,
-                                                              duals=duals,
-                                                              route=route)
+                more_columns = self._attempt_solve_best_paths(
+                    vehicle=vehicle, duals=duals, route=route
+                )
             elif pricing_strategy == "BestEdges1":
-                more_columns = self._attempt_solve_best_edges1(vehicle=vehicle,
-                                                               duals=duals,
-                                                               route=route)
+                more_columns = self._attempt_solve_best_edges1(
+                    vehicle=vehicle, duals=duals, route=route
+                )
             elif pricing_strategy == "BestEdges2":
-                more_columns = self._attempt_solve_best_edges2(vehicle=vehicle,
-                                                               duals=duals,
-                                                               route=route)
+                more_columns = self._attempt_solve_best_edges2(
+                    vehicle=vehicle, duals=duals, route=route
+                )
             elif pricing_strategy == "Exact":
-                more_columns = self._attempt_solve_exact(vehicle=vehicle,
-                                                         duals=duals,
-                                                         route=route)
+                more_columns = self._attempt_solve_exact(
+                    vehicle=vehicle, duals=duals, route=route
+                )
         # old approach
         else:
             if pricing_strategy == "BestPaths":
-                more_columns = self._attempt_solve_best_paths(vehicle=vehicle,
-                                                              duals=duals,
-                                                              route=route)
+                more_columns = self._attempt_solve_best_paths(
+                    vehicle=vehicle, duals=duals, route=route
+                )
             elif pricing_strategy == "BestEdges1":
-                more_columns = self._attempt_solve_best_edges1(vehicle=vehicle,
-                                                               duals=duals,
-                                                               route=route)
+                more_columns = self._attempt_solve_best_edges1(
+                    vehicle=vehicle, duals=duals, route=route
+                )
             elif pricing_strategy == "BestEdges2":
-                more_columns = self._attempt_solve_best_edges2(vehicle=vehicle,
-                                                               duals=duals,
-                                                               route=route)
+                more_columns = self._attempt_solve_best_edges2(
+                    vehicle=vehicle, duals=duals, route=route
+                )
             if pricing_strategy == "Exact" or not more_columns:
-                more_columns = self._attempt_solve_exact(vehicle=vehicle,
-                                                         duals=duals,
-                                                         route=route)
+                more_columns = self._attempt_solve_exact(
+                    vehicle=vehicle, duals=duals, route=route
+                )
         return more_columns
 
     def _attempt_solve_best_paths(self, vehicle=None, duals=None, route=None):
@@ -626,7 +654,8 @@ def _attempt_solve_best_paths(self, vehicle=None, duals=None, route=None):
                 k_shortest_paths,
             )
             self.routes, self._more_routes = subproblem.solve(
-                self._get_time_remaining())
+                self._get_time_remaining()
+            )
             more_columns = self._more_routes
             if more_columns:
                 break
@@ -645,7 +674,8 @@ def _attempt_solve_best_edges1(self, vehicle=None, duals=None, route=None):
                 alpha,
             )
             self.routes, self._more_routes = subproblem.solve(
-                self._get_time_remaining(),)
+                self._get_time_remaining(),
+            )
             more_columns = self._more_routes
             if more_columns:
                 break
@@ -677,15 +707,16 @@ def _attempt_solve_best_edges2(self, vehicle=None, duals=None, route=None):
 
     def _attempt_solve_exact(self, vehicle=None, duals=None, route=None):
         subproblem = self._def_subproblem(duals, vehicle, route)
-        self.routes, self._more_routes = subproblem.solve(
-            self._get_time_remaining())
+        self.routes, self._more_routes = subproblem.solve(self._get_time_remaining())
         return self._more_routes
 
     def _get_next_pricing_strategy(self, relaxed_cost):
         "Return the appropriate pricing strategy based on input parameters"
         pricing_strategy = None
-        if (self._pricing_strategy == "Hyper" and
-                self._no_improvement != self._run_exact):
+        if (
+            self._pricing_strategy == "Hyper"
+            and self._no_improvement != self._run_exact
+        ):
             self._no_improvement_iteration = self._iteration
             if self._iteration == 0:
                 pricing_strategy = "BestPaths"
@@ -702,17 +733,16 @@ def _get_next_pricing_strategy(self, relaxed_cost):
         return pricing_strategy
 
     def _update_hyper_heuristic(self, relaxed_cost: float):
-        best_paths, best_paths_freq = self.masterproblem.get_heuristic_distribution(
-        )
+        best_paths, best_paths_freq = self.masterproblem.get_heuristic_distribution()
         self.hyper_heuristic.current_performance(
             new_objective_value=relaxed_cost,
             produced_column=self._more_routes,
             active_columns=best_paths_freq,
         )
         self.hyper_heuristic.move_acceptance()
-        self.hyper_heuristic.update_parameters(self._iteration,
-                                               self._no_improvement,
-                                               self._no_improvement_iteration)
+        self.hyper_heuristic.update_parameters(
+            self._iteration, self._no_improvement, self._no_improvement_iteration
+        )
 
     def _get_time_remaining(self, mip: bool = False):
         """
@@ -818,12 +848,14 @@ def _solve_with_clarke_wright(self):
         )
         alg.run()
         self._best_routes = dict(
-            zip([r for r in range(len(alg.best_routes))], alg.best_routes))
+            zip([r for r in range(len(alg.best_routes))], alg.best_routes)
+        )
         self._best_routes_vehicle_type = dict(
             zip(
                 [r for r in range(len(alg.best_routes))],
                 [0 for r in range(len(alg.best_routes))],
-            ))
+            )
+        )
 
     # Preprocessing #
 
@@ -835,9 +867,13 @@ def _get_initial_solution(self):
         """
         self._initial_routes = []
         # Run Clarke & Wright if possible
-        if (not self.time_windows and not self.pickup_delivery and
-                not self.distribution_collection and not self.mixed_fleet and
-                not self.periodic):
+        if (
+            not self.time_windows
+            and not self.pickup_delivery
+            and not self.distribution_collection
+            and not self.mixed_fleet
+            and not self.periodic
+        ):
             best_value = 1e10
             best_num_vehicles = 1e10
             for alpha in [x / 10 for x in range(1, 20)]:
@@ -858,15 +894,17 @@ def _get_initial_solution(self):
                     best_value = alg.best_value
                     best_num_vehicles = len(alg.best_routes)
             logger.info(
-                "Clarke & Wright solution found with value %s and %s vehicles" %
-                (best_value, best_num_vehicles))
+                "Clarke & Wright solution found with value %s and %s vehicles"
+                % (best_value, best_num_vehicles)
+            )
 
             # Run greedy algorithm if possible
-            alg = _Greedy(self.G, self.load_capacity, self.num_stops,
-                          self.duration)
+            alg = _Greedy(self.G, self.load_capacity, self.num_stops, self.duration)
             alg.run()
-            logger.info("Greedy solution found with value %s and %s vehicles" %
-                        (alg.best_value, len(alg.best_routes)))
+            logger.info(
+                "Greedy solution found with value %s and %s vehicles"
+                % (alg.best_value, len(alg.best_routes))
+            )
             self._initial_routes += alg.best_routes
 
             # TO DO: Run heuristics from VeRyPy for the CVRP
@@ -876,9 +914,7 @@ def _get_initial_solution(self):
         elif self.pickup_delivery:
             for v in self.G.nodes():
                 if "request" in self.G.nodes[v]:
-                    obvious_route = [
-                        "Source", v, self.G.nodes[v]["request"], "Sink"
-                    ]
+                    obvious_route = ["Source", v, self.G.nodes[v]["request"], "Sink"]
                     edges = [
                         ("Source", v),
                         (v, self.G.nodes[v]["request"]),
@@ -897,7 +933,7 @@ def _get_initial_solution(self):
     def _convert_initial_routes_to_digraphs(self):
         """
         Converts list of initial routes to list of Digraphs.
-        By default, initial routes are computed with vehicle type 0 (the first one in the list).
+        By default, initial routes are computed with the first feasible vehicle type.
         """
         self._routes = []
         self._routes_with_node = {}
@@ -910,7 +946,20 @@ def _convert_initial_routes_to_digraphs(self):
                 G.add_edge(i, j, cost=edge_cost)
                 total_cost += edge_cost
             G.graph["cost"] = total_cost
-            G.graph["vehicle_type"] = 0
+            if self.load_capacity:
+                for k in range(len(self.load_capacity)):
+                    if (
+                        sum(self.G.nodes[v]["demand"] for v in r)
+                        <= self.load_capacity[k]
+                    ):
+                        G.graph["vehicle_type"] = k
+                        break
+            else:
+                G.graph["vehicle_type"] = 0
+            if "vehicle_type" not in G.graph:
+                raise ValueError(
+                    "Could not find initial feasible solution. Check loading capacities."
+                )
             self._routes.append(G)
             for v in r[1:-1]:
                 if v in self._routes_with_node:
@@ -951,8 +1000,10 @@ def _add_fixed_costs(self):
     def _remove_infeasible_arcs_capacities(self):
         infeasible_arcs = []
         for (i, j) in self.G.edges():
-            if (self.G.nodes[i]["demand"] + self.G.nodes[j]["demand"] >
-                    self._max_capacity):
+            if (
+                self.G.nodes[i]["demand"] + self.G.nodes[j]["demand"]
+                > self._max_capacity
+            ):
                 infeasible_arcs.append((i, j))
         self.G.remove_edges_from(infeasible_arcs)
 
@@ -971,14 +1022,13 @@ def _remove_infeasible_arcs_time_windows(self):
                     # earliest time is coming straight from depot
                     self.G.nodes[v]["lower"] = max(
                         self.G.nodes[v]["lower"],
-                        self.G.nodes["Source"]["lower"] +
-                        self.G.edges["Source", v]["time"],
+                        self.G.nodes["Source"]["lower"]
+                        + self.G.edges["Source", v]["time"],
                     )
                     # Latest time is going straight to depot
                     self.G.nodes[v]["upper"] = min(
                         self.G.nodes[v]["upper"],
-                        self.G.nodes["Sink"]["upper"] -
-                        self.G.edges[v, "Sink"]["time"],
+                        self.G.nodes["Sink"]["upper"] - self.G.edges[v, "Sink"]["time"],
                     )
         self.G.remove_edges_from(infeasible_arcs)
 
@@ -1005,18 +1055,19 @@ def _set_zero_attributes(self):
 
         for v in self.G.nodes():
             for attribute in [
-                    "demand",
-                    "collect",
-                    "service_time",
-                    "lower",
-                    "upper",
+                "demand",
+                "collect",
+                "service_time",
+                "lower",
+                "upper",
             ]:
                 if attribute not in self.G.nodes[v]:
                     self.G.nodes[v][attribute] = 0
             # Ignore demand at Source/Sink
             if v in ["Source", "Sink"] and self.G.nodes[v]["demand"] > 0:
-                logger.warning("Demand %s at node %s is ignored." %
-                               (self.G.nodes[v]["demand"], v))
+                logger.warning(
+                    "Demand %s at node %s is ignored." % (self.G.nodes[v]["demand"], v)
+                )
                 self.G.nodes[v]["demand"] = 0
 
             # Set frequency = 1 if missing
@@ -1036,9 +1087,11 @@ def _readjust_sink_time_windows(self):
 
         if self.G.nodes["Sink"]["upper"] == 0:
             self.G.nodes["Sink"]["upper"] = max(
-                self.G.nodes[u]["upper"] + self.G.nodes[u]["service_time"] +
-                self.G.edges[u, "Sink"]["time"]
-                for u in self.G.predecessors("Sink"))
+                self.G.nodes[u]["upper"]
+                + self.G.nodes[u]["service_time"]
+                + self.G.edges[u, "Sink"]["time"]
+                for u in self.G.predecessors("Sink")
+            )
 
     def _update_dummy_attributes(self):
         """Adds dummy attributes on nodes and edges if missing."""
@@ -1067,8 +1120,7 @@ def _best_routes_as_node_lists(self):
         for route in self._best_routes_as_graphs:
             node_list = shortest_path(route, "Source", "Sink")
             self._best_routes[route_id] = node_list
-            self._best_routes_vehicle_type[route_id] = route.graph[
-                "vehicle_type"]
+            self._best_routes_vehicle_type[route_id] = route.graph["vehicle_type"]
             route_id += 1
         # Merge with preassigned complete routes
         for route in self._preassignments:
@@ -1078,8 +1130,7 @@ def _best_routes_as_node_lists(self):
                 best_cost = 1e10
                 for k in range(self._vehicle_types):
                     # If different vehicles, the cheapest feasible one is accounted for
-                    cost = sum(
-                        self._H.edges[i, j]["cost"][k] for (i, j) in edges)
+                    cost = sum(self._H.edges[i, j]["cost"][k] for (i, j) in edges)
                     load = sum(self._H.nodes[i]["demand"] for i in route)
                     if cost < best_cost:
                         if self.load_capacity: