My model does not work right, can someone help me?

Hello, currently i am trying to automize validating my solution of a Capacitated Vehicle Routing Problem of an Ant colony optimization algorithm. The code runs and gives me a solution, however, when printing the solution it does not seem to be the correct solution.

The problem consists of a CVRP with an unknown number of vehicles which can either be, a car and a tractor. With different capacities, car: 1000, tractor: 5000. They also both have different distance matrices. The model is infeasible, and I am not quite sure why. I hope someone can point me in the right direction

The matrices are: {'tractor': [[0, 0, 5385, 0], [0, 0, 5385, 0], [5906, 5906, 0, 5906], [0, 0, 5385, 0]], 'car': [[0, 0, 5015, 0], [0, 0, 5015, 0], [4657, 4657, 0, 4657], [0, 0, 5015, 0]]}

The capacities are: {'car': 1000, 'tractor': 5000}

the vehicle types: ['tractor', 'car'].

And lastly my code is:

import numpy as np

from gurobipy import Model, GRB, quicksum

import pandas as pd

# Global variables

# Gurobi license location (Not used in the code)

GURIBO_LICENSE = r"C:\Users\jkort\OneDrive\Bureaublad\Bachelor Thesis Jordi Kortekaas\Code\Guribo License\gurobi.lic"

# Distance matrix defining distances between clients

# Number of vehicles available for deliveries

# Function to define variables for the model

def define_variables(distance_matrix_dict, num_vehicles, capacity_dict, demand_list, vehicle_types):

    n = len(demand_list) #number of clients

    N = [i for i in range(1,n)] # set of clients

    V = [i for i in range(n)] # set of clients plus the depot (depot is 0)

    A = [(i,j) for i in V for j in V if i != j] # set of arcs connecting different locations

    # cost associated with each arc (i,j), i.e., the distance from i to j

    c = {(i,j,k): distance_matrix_dict[vehicle_types[k]][i][j] for i,j in A for k in range(len(vehicle_types))}

   

    # capacity of vehicles

    Q = {k: capacity_dict[vehicle_types[k]] for k in range(len(vehicle_types))}

    # demand at each client

    q = {i: demand_list[i] for i in N}

    # initialize the model

    mdl = Model("CVRP")

   

    # decision variables indicating whether arc (i,j) is used by vehicle k

   

    x = mdl.addVars(A, range(len(vehicle_types)), vtype=GRB.BINARY)

    print(x)

    # decision variables indicating the load of vehicle k after visiting client i

    u = mdl.addVars(N, range(len(vehicle_types)), vtype=GRB.CONTINUOUS)

   

    # set objective: minimize total distance traveled by all vehicles

    mdl.modelSense = GRB.MINIMIZE

    mdl.setObjective(quicksum(x[i,j,k]*c[i,j,k] for (i,j) in A for k in range(len(vehicle_types))))

    # Constraints:

    # Each client is visited exactly once by one vehicle

    mdl.addConstrs(quicksum(x[i, j, k] for j in V if j != i for k in range(len(vehicle_types))) == 1 for i in N)

    mdl.addConstrs(quicksum(x[i, j, k] for i in V if i != j for k in range(len(vehicle_types))) == 1 for j in N)

    # Load of vehicle k after visiting client j is greater than load after visiting client i, if i is visited before j by vehicle k

    mdl.addConstrs((x[i, j, k] == 1) >> (u[i, k]+q[j] == u[j, k]) for i, j in A if i != 0 and j != 0 for k in range(len(vehicle_types)))

    # Load of vehicle k after visiting each client i is greater or equal to the demand at i

    mdl.addConstrs(u[i, k] >= q[i] for i in N for k in range(len(vehicle_types)))

    # Load of vehicle k after visiting each client i is less or equal to the capacity of the vehicle

    mdl.addConstrs(u[i, k] <= Q[k] for i in N for k in range(len(vehicle_types)))

    # Each vehicle must leave the depot

    mdl.addConstrs(quicksum(x[0, i, k] for i in V if i != 0) == 1 for k in range(len(vehicle_types)))

    # Ensure that the same vehicle is used when arriving at client i and when leaving from it

    mdl.addConstrs(quicksum(x[i,j,k] for j in V if j != i) == quicksum(x[j,i,k] for j in V if j != i) for i in N for k in range(len(vehicle_types)))

   

    # Set the optimization tolerance and time limit

    mdl.Params.MIPGap = 0.1

    mdl.Params.TimeLimit = 30  # seconds

   

    # Write model to a file

    mdl.write("myLP.lp")

    # Solve the model

    mdl.optimize()

   

    # Print model details

    print(mdl)

    # do IIS if the model is infeasible

    # initialize the total distance

    total_distance = 0

    # initialize an empty DataFrame to store the route order

    df_route_order = pd.DataFrame()

    # Print solution

    for k in range(len(vehicle_types)):

        solution_list = [0]

        while True:

            added = False  # We will use this flag to detect when we can't find a next node

            for j in V:

                if j not in solution_list and x[solution_list[-1], j, k].X > 0.5:

                   

                    total_distance += c[solution_list[-1], j,k]  # add the distance to total_distance

                    solution_list.append(j)

                    added = True

                    break

            if not added or solution_list[-1] == 0:  # If we didn't add any new node, or we're back at the depot, the route for this vehicle is done

                break

        print(f"The route for the{vehicle_types[k]}: ", solution_list, f"The distance for the{vehicle_types[k]}: ",total_distance)

        df_route_order[f'vehicle_{k}'] = pd.Series(solution_list)  # add the solution_list to the DataFrame

    return total_distance, df_route_order

def main(distance_matrix_dict, num_vehicles, capacity_dict, demand_list, vehicle_types):

    total_distance, df_route_order = define_variables(distance_matrix_dict, num_vehicles, capacity_dict, demand_list, vehicle_types)

    return total_distance, df_route_order

if __name__ == "main":

    main(distance_matrix, num_vehicles, capacity, demand)

Riley Clement

Gurobi Staff

August 06, 2023 23:51

Hi Jordi,

Obviously we can't run your model without the data, but I note your code includes the following comment

# do IIS if the model is infeasible

Have you tried doing this? This will likely be the fastest way of discovering what is wrong.

- Riley

Jordi Kortekaas

August 07, 2023 11:01
Edited

Hi Riley,

I have, and it is not necessarily that it does not try to optimize the objective function. It does not find a solution, as the solution count is 0. I cannot find what is wrong. The variables of input are:

distance_matrix_dict = {'tractor': [[0, 0, 0, 460, 460, 157, 184, 0, 194, 722, 120, 290, 369, 165], [0, 0, 0, 460, 460, 157, 184, 0, 194, 722, 120, 290, 369, 165], [0, 0, 0, 460, 460, 157, 184, 0, 194, 722, 120, 290, 369, 165], [453, 453, 453, 0, 0, 546, 566, 453, 628, 1136, 480, 556, 164, 477], [453, 453, 453, 0, 0, 546, 566, 453, 628, 1136, 480, 556, 164, 477], [156, 156, 156, 542, 542, 0, 53, 156, 135, 808, 276, 217, 416, 322], [178, 178, 178, 567, 567, 48, 0, 178, 157, 760, 298, 189, 441, 344], [0, 0, 0, 460, 460, 157, 184, 0, 194, 722, 120, 290, 369, 165], [194, 194, 194, 626, 626, 135, 163, 194, 0, 757, 314, 317, 499, 360], [728, 728, 728, 1147, 1147, 808, 759, 728, 761, 0, 676, 936, 1069, 718], [120, 120, 120, 485, 485, 277, 305, 120, 314, 674, 0, 410, 408, 64], [292, 292, 292, 557, 557, 217, 184, 292, 316, 932, 412, 0, 431, 458], [374, 374, 374, 164, 164, 426, 446, 374, 507, 1057, 401, 435, 0, 398], [165, 165, 165, 482, 482, 322, 349, 165, 359, 714, 63, 455, 404, 0]],

'car': [[0, 0, 0, 359, 359, 109, 123, 0, 160, 434, 80, 236, 263, 119], [0, 0, 0, 359, 359, 109, 123, 0, 160, 434, 80, 236, 263, 119], [0, 0, 0, 359, 359, 109, 123, 0, 160, 434, 80, 236, 263, 119], [354, 354, 354, 0, 0, 420, 378, 354, 499, 761, 372, 418, 164, 378], [354, 354, 354, 0, 0, 420, 378, 354, 499, 761, 372, 418, 164, 378], [108, 108, 108, 417, 417, 0, 40, 108, 135, 419, 188, 187, 305, 227], [119, 119, 119, 378, 378, 36, 0, 119, 146, 383, 199, 144, 266, 238], [0, 0, 0, 359, 359, 109, 123, 0, 160, 434, 80, 236, 263, 119], [177, 177, 177, 498, 498, 135, 138, 177, 0, 439, 146, 286, 387, 187], [436, 436, 436, 767, 767, 419, 383, 436, 433, 0, 402, 519, 690, 428], [80, 80, 80, 375, 375, 189, 203, 80, 146, 402, 0, 316, 298, 50], [231, 231, 231, 418, 418, 187, 141, 231, 286, 516, 311, 0, 306, 350], [275, 275, 275, 164, 164, 311, 269, 275, 390, 683, 293, 309, 0, 300], [118, 118, 118, 381, 381, 227, 242, 118, 187, 427, 50, 355, 304, 0]]}

num_vehicles = 4

capacity_dict = {'car': 1100, 'tractor': 5000}

demand_list = [0, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750]

vehicle_types = ['tractor', 'car']

August 07, 2023 11:51

I don't think there is a feasible solution with this data. If I understand correctly, every client is visited by either the tractor or the car, and the vehicle's load increases by the demand at that client, but must remain below the capacity of the vehicle.

The total demand is 9750, while the total capacity of the vehicles is 6100 - so it seems that the problem cannot be satisfied but it is taking Gurobi a long time to prove infeasibility.

If you were to change the car capacity to 5000, so that the total capacity exceeds the demand, then you will find that solutions are found rather readily.

My model does not work right, can someone help me?

Comments

Didn't find what you were looking for?