diff --git a/vrp-core/src/construction/heuristics/metrics.rs b/vrp-core/src/construction/heuristics/metrics.rs
index 981df54a5..e733e5fb3 100644
--- a/vrp-core/src/construction/heuristics/metrics.rs
+++ b/vrp-core/src/construction/heuristics/metrics.rs
@@ -5,9 +5,11 @@ mod metrics_test;
 use crate::construction::enablers::{TotalDistanceTourState, TotalDurationTourState, WaitingTimeActivityState};
 use crate::construction::features::MaxVehicleLoadTourState;
 use crate::construction::heuristics::{InsertionContext, RouteContext, RouteState};
+use crate::models::common::Distance;
 use crate::models::problem::{TransportCost, TravelTime};
 use rosomaxa::algorithms::math::*;
 use rosomaxa::prelude::*;
+use rosomaxa::utils::parallel_collect;
 use std::cmp::Ordering;
 
 /// Gets max load variance in tours.
@@ -180,11 +182,9 @@ pub fn get_distance_gravity_mean(insertion_ctx: &InsertionContext) -> Float {
     let profile = insertion_ctx.solution.routes.first().map(|route_ctx| &route_ctx.route().actor.vehicle.profile);
 
     if let Some(profile) = profile {
-        let medoids = insertion_ctx
-            .solution
-            .routes
-            .iter()
-            .filter_map(|route_ctx| get_medoid(route_ctx, transport))
+        let medoids = parallel_collect(&insertion_ctx.solution.routes, |route_ctx| get_medoid(route_ctx, transport))
+            .into_iter()
+            .flatten()
             .collect::<Vec<_>>();
 
         let mut distances = Vec::with_capacity(medoids.len() * 2);
@@ -212,46 +212,53 @@ pub fn group_routes_by_proximity(insertion_ctx: &InsertionContext) -> RouteProxi
     let profile = &solution.routes.first().map(|route_ctx| &route_ctx.route().actor.vehicle.profile)?;
     let transport = insertion_ctx.problem.transport.as_ref();
 
-    let indexed_medoids = solution
-        .routes
-        .iter()
-        .enumerate()
-        .map(|(idx, route_ctx)| (idx, get_medoid(route_ctx, transport)))
-        .collect::<Vec<_>>();
+    let indexed_route_clusters =
+        parallel_collect(&solution.routes, |route_ctx| get_approx_clusters(route_ctx, transport))
+            .into_iter()
+            .enumerate()
+            .collect::<Vec<_>>();
 
-    Some(
-        indexed_medoids
+    Some(parallel_collect(&indexed_route_clusters, |(outer_idx, outer_clusters)| {
+        let mut route_distances = indexed_route_clusters
             .iter()
-            .map(|(outer_idx, outer_medoid)| {
-                let mut route_distances = indexed_medoids
-                    .iter()
-                    .filter(move |(inner_idx, _)| *outer_idx != *inner_idx)
-                    .map(move |(inner_idx, inner_medoid)| {
-                        let distance = match (outer_medoid, inner_medoid) {
-                            (Some(outer_medoid), Some(inner_medoid)) => {
-                                let distance = transport.distance_approx(profile, *outer_medoid, *inner_medoid);
-                                if distance < 0. {
-                                    None
-                                } else {
-                                    Some(distance)
-                                }
-                            }
-                            _ => None,
-                        };
-                        (*inner_idx, distance)
-                    })
-                    .collect::<Vec<_>>();
-
-                route_distances.sort_unstable_by(|(_, a_distance), (_, b_distance)| match (a_distance, b_distance) {
-                    (Some(a_distance), Some(b_distance)) => compare_floats(*a_distance, *b_distance),
-                    (Some(_), None) => Ordering::Less,
-                    _ => Ordering::Greater,
-                });
-
-                route_distances
+            .filter(move |(inner_idx, _)| *outer_idx != *inner_idx)
+            .map(move |(inner_idx, inner_clusters)| {
+                let distance = match (outer_clusters, inner_clusters) {
+                    (Some(outer_clusters), Some(inner_clusters)) => {
+                        // get a sum of distances between all pairs of clusters
+                        let pair_distance = outer_clusters
+                            .iter()
+                            .flat_map(|outer| inner_clusters.iter().map(move |inner| (inner, outer)))
+                            .map(|(&o, &i)| {
+                                transport.distance_approx(profile, o, i).max(0.)
+                                    + transport.distance_approx(profile, i, o).max(0.)
+                            })
+                            .sum::<Distance>()
+                            / 2.;
+
+                        let total_pairs = outer_clusters.len() * inner_clusters.len();
+                        if total_pairs == 0 {
+                            None
+                        } else {
+                            // get average distance between clusters
+                            Some(pair_distance / total_pairs as Float)
+                        }
+                    }
+                    _ => None,
+                };
+
+                (*inner_idx, distance)
             })
-            .collect::<Vec<_>>(),
-    )
+            .collect::<Vec<_>>();
+
+        route_distances.sort_unstable_by(|(_, a_distance), (_, b_distance)| match (a_distance, b_distance) {
+            (Some(a_distance), Some(b_distance)) => compare_floats(*a_distance, *b_distance),
+            (Some(_), None) => Ordering::Less,
+            _ => Ordering::Greater,
+        });
+
+        route_distances
+    }))
 }
 
 fn get_values_from_route_state<'a>(
@@ -281,3 +288,32 @@ fn get_medoid(route_ctx: &RouteContext, transport: &(dyn TransportCost)) -> Opti
         .min_by(|(sum_a, _), (sum_b, _)| compare_floats(*sum_a, *sum_b))
         .map(|(_, location)| location)
 }
+
+fn get_approx_clusters(route_ctx: &RouteContext, transport: &(dyn TransportCost)) -> Option<Vec<usize>> {
+    const CLUSTER_SIZE_RATIO: Float = 0.2;
+
+    let distance_threshold = route_ctx.state().get_total_distance().copied()? * CLUSTER_SIZE_RATIO;
+    let profile = &route_ctx.route().actor.vehicle.profile;
+
+    Some(
+        route_ctx
+            .route()
+            .tour
+            .all_activities()
+            .map(|activity| activity.place.location)
+            .fold((Distance::default(), None, Vec::default()), |(mut acc_distance, prev, mut clusters), current| {
+                if let Some(prev) = prev {
+                    acc_distance += transport.distance_approx(profile, prev, current).max(0.);
+                    if acc_distance > distance_threshold {
+                        clusters.push(current);
+                        acc_distance = Distance::default();
+                    }
+                } else {
+                    clusters.push(current)
+                }
+
+                (acc_distance, Some(current), clusters)
+            })
+            .2,
+    )
+}