I wrote this function to minimize the distances between orders assigned by the model by penalizing any deviations from the centroid of each route.

Optapy works like this:
1- Takes your problem facts and entities as classes
2- Updates the routes for the vehicle using the setters and getters methods you define:

In the vehicle class:

    @planning_list_variable(Customer, ['customer_range'])
    def get_customer_list(self):
        return self.customer_list

    def set_customer_list(self, customer_list):
        self.customer_list = customer_list

Referencing it in the routing solution:

@planning_solution
class VehicleRoutingSolution:
    """
    The VehicleRoutingSolution class represents both the problem and the solution
    in the vehicle routing domain. It stores references to all the problem facts
    (locations, depots, customers) and planning entities (vehicles) that define the problem.
    
    Attributes:
        name (str): The name of the solution.
        location_list (list of Location): A list of all locations involved in the routing.
        depot_list (list of Depot): A list of depots where vehicles start and end their routes.
        vehicle_list (list of Vehicle): A list of all vehicles used in the routing problem.
        customer_list (list of Customer): A list of all customers to be served by the vehicles.
        south_west_corner (Location): The southwestern corner of the bounding box for visualization.
        north_east_corner (Location): The northeastern corner of the bounding box for visualization.
        score (HardSoftScore, optional): The score of the solution, reflecting the quality of the solution.
    """

    def __init__(self, name, location_list, depot_list, vehicle_list, customer_list,
                 south_west_corner, north_east_corner, score=None):
        self.name = name
        self.location_list = location_list
        self.depot_list = depot_list
        self.vehicle_list = vehicle_list
        self.customer_list = customer_list
        self.south_west_corner = south_west_corner
        self.north_east_corner = north_east_corner
        self.score = score

    @planning_entity_collection_property(Vehicle)
    def get_vehicle_list(self):
        return self.vehicle_list

    @problem_fact_collection_property(Customer)
    @value_range_provider('customer_range', value_range_type=list)
    def get_customer_list(self):
        return self.customer_list

So, I tried using a constraint:

def penalize_deviation_from_centroid(vehicle):
    """
    Calculate a penalty based on the deviation of customer locations from the centroid.
    Ensure that it returns a numeric value (int/float) for penalizing deviations.
    We used ecludian distance here to reduce API calls.
    
    Args:
        vehicle (Vehicle): The vehicle with assigned customers.
    
    Returns:
        float: The total penalty based on deviations from the centroid.
    """
    if not vehicle.customer_list or len(vehicle.customer_list) == 0:
        return 0.0  # No customers, no penalty

    # Step 2: Calculate the centroid (average latitude and longitude)
    def calculate_centroid(customers):
        avg_lat = sum(float(customer.location.latitude) for customer in customers) / len(customers)
        avg_lon = sum(float(customer.location.longitude) for customer in customers) / len(customers)
        return Location(avg_lat, avg_lon)  # Return a Location object

    def euclidean_distance(loc1, loc2):
        """ Compute the Euclidean distance between two locations based on latitude and longitude. """
        return ((loc1.latitude - loc2.latitude) ** 2 + (loc1.longitude - loc2.longitude) ** 2) ** 0.5

    centroid = calculate_centroid(vehicle.customer_list)
    total_penalty = sum(euclidean_distance(customer.location, centroid) for customer in vehicle.customer_list)

    return total_penalty

def spherical_routes(constraint_factory):
    """
    Encourage routes to be circular by penalizing distance from the centroid of assigned locations.
    
    Args:
        constraint_factory (ConstraintFactory): The factory to create constraints.

    Returns:
        Constraint: The constraint penalizing large deviations from the route centroid.
    """
    return constraint_factory 
        .for_each(Vehicle) 
        .penalize("Penalize deviation from centroid", HardSoftScore.ofSoft(1_000),
                  lambda vehicle: int(penalize_deviation_from_centroid(vehicle) * WEIGHTS['spherical_routes']))

And then applying:
# Add to constraint provider
from optapy import constraint_provider

@constraint_provider
def vehicle_routing_constraints(constraint_factory):
    """
    Define all the constraints for the vehicle routing problem.

    This function combines all hard constraints related to vehicle capacity, CBM, weight, skills, and working time.

    Args:
        constraint_factory (ConstraintFactory): The factory to create constraints.

    Returns:
        list: A list of constraints that will be used to guide the solver.
    """
    return [
        # Hard constraints
        vehicle_cbm(constraint_factory), # Constraint Weight is 10 Million
        vehicle_weight(constraint_factory), # Constraint Weight is 1,000
        # vehicle_skill_constraint(constraint_factory),
        
        # # Soft constraints
        vehicle_time_constraint(constraint_factory), # Soft Weight of 1 Million
        total_vehicle_cost(constraint_factory), # Soft Weight is 1,000
        spherical_routes(constraint_factory), # # Soft Weight is 1,000
    ]

And then, the solving:

# Step 1: Setup the solver manager with the appropriate config
solver_config = optapy.config.solver.SolverConfig()
solver_config 
    .withSolutionClass(VehicleRoutingSolution) 
    .withEntityClasses(Vehicle) 
    .withConstraintProviderClass(vehicle_routing_constraints) 
    .withTerminationSpentLimit(Duration.ofSeconds(800))

# Step 2: Create the solver manager
solver_manager = solver_manager_create(solver_config)

# # Create the initial solution for the solver
solution = VehicleRoutingSolution(
    name="Vehicle Routing Problem with Actual Data",
    location_list=locations,
    depot_list=depots,
    vehicle_list=vehicles,
    customer_list=customers,
    south_west_corner=Location(29.990707246305476, 31.229210746581806),
    north_east_corner=Location(30.024396202211875, 31.262640488654238)
)

# Step 3: Solve the problem and get the solver job
SINGLETON_ID = 1  # A unique problem ID (can be any number)
solver_job = solver_manager.solve(SINGLETON_ID, lambda _: solution)

# Step 4: Get the best solution from the solver job
best_solution = solver_job.getFinalBestSolution()

I get this error:
org.optaplanner.jpyinterpreter.types.errors.org.optaplanner.jpyinterpreter.types.errors.AttributeError: org.optaplanner.jpyinterpreter.types.errors.AttributeError: object '<class list>' does not have attribute '__bool__'

Please Note: When I remove the spherical_routes constraint, the model runs just fine.

I tried so many times validating the output of the function and it’s always a number but it’s giving the same issue.

It’s trying for some reason to compare a list to number and lists don’t have __bool__ attribute to allow that operation.