Capacitated Arc Routing Problem (CARP)

Principles learned

  • Add multiple list decision variables

  • Use a lambda expression to compute a sum with a variable number of terms

  • Add a disjoint constraint

  • Define a sequence of expressions

  • Use of “contains” and “at” operators

Problem

../_images/capacitatedarcrouting.svg

In the capacitated arc routing problem (CARP), a fleet of delivery vehicles with uniform capacity must service edges with known demand. The vehicles start and end their routes at a common depot. Each edge can only be served by one vehicle. The objective is to assign a sequence of edges to serve to each truck of the fleet minimizing the total distance traveled such that all edges are served and the total demand served by each truck does not exceed its capacity.

Download the example


Data

The instances provided are those that can be downloaded from the DIMACS website.

The format of the data files is as follows:

  • The number of nodes follows the keyword VERTICES.

  • The number of required edges (edges with positive demand) follows the keyword ARISTAS_REQ.

  • The number of non-required edges (edges with zero demand) follows the keyword ARISTAS_NOREQ.

  • The number of vehicles available follows the keyword VEHICULOS.

  • The capacity of a truck follows the keyword CAPACIDAD.

  • After the keyword LISTA_ARISTAS_REQ, for each required edge is listed its demand and its cost.

  • After the keyword LISTA_ARISTAS_NOREQ, for each non-required edge is listed its cost.

  • The index of the depot node follows the keyword DEPOSITO.

Program

This Hexaly model defines a route for each truck k as a list variable (edgesSequences[k]). It corresponds to the sequence of edges visited and serviced. To ensure that all edges must be served by at most one truck, all the list variables are constrained to form a disjunction, thanks to the “disjoint” operator.

The quantity delivered on each visit is the demand on the edge of this visit. This expression is just defined with an “at” operator to access the array requiredEdgesDemand. The total quantity delivered by each truck is computed with a function to apply the sum operator over all edges visited (and serviced) by a truck. Note that the number of terms in this sum varies automatically with the size of the list. This quantity is constrained to be lower than the truck capacity.

For each truck, the distance travelled from the visit n-1 to the visit n is accessed with an “at” operator to the multi-dimensional array distanceBetweenEdges, with the first index. Here again we use a function to sum distances along each route.

Execution:
localsolver capacitated_arc_routing.lsp inFileName=instances/egl-e1-A.dat [lsTimeLimit=] [solFileName=]
use io;

/* Read instance data. The input files follow the format of the DIMACS challenge */
function input() {
    usage = "Usage: localsolver capacitated_arc_routing.lsp "
            + "inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;
    readData();
}

/* Declare the optimization model */
function model() {
    // Sequence of edges visited and "serviced" by each truck
    // A sequence can contain an edge in one direction or its reverse
    edgesSequences[k in 0...nbTrucks] <- list(2 * nbRequiredEdges);

    // An edge must be serviced by at most one truck
    constraint disjoint(edgesSequences);

    // An edge can be travelled in both directions but its demand must be
    // satisfied only once
    for [i in 0...nbRequiredEdges] {
        constraint contains(edgesSequences, 2 * i) + contains(edgesSequences, 2 * i + 1) == 1;
    }

    for [k in 0...nbTrucks] {
        local sequence <- edgesSequences[k];
        local c <- count(sequence);

        // Quantity in each truck
        routeQuantity <- sum(sequence, j => requiredEdgesDemands[j]);
        // Capacity constraint : a truck must not exceed its capacity
        constraint routeQuantity <= truckCapacity;

        // Distance travelled by each truck
        routeDistance[k] <- sum(1...c, i => requiredEdgesCosts[sequence[i]] +
                distanceBetweenEdges[sequence[i - 1]][sequence[i]]) +
                (c > 0 ? requiredEdgesCosts[sequence[0]] + distanceFromDepot[sequence[0]] +
                distanceToDepot[sequence[c - 1]] : 0);
    }

    // Total distance travelled
    totalDistance <- sum[k in 0...nbTrucks](routeDistance[k]);

    // Objective: minimize the distance travelled
    minimize totalDistance;
}

/* Parametrize the solver */
function param() {
    if (lsTimeLimit == nil) lsTimeLimit = 20;
}

/* Write the solution in a file with the following format:
 * - total distance
 * - number of routes
 * - for each truck, the edges visited */
function output() {
    if (solFileName == nil) return;
    local outfile = io.openWrite(solFileName);
    outfile.println("Objective function value : ", totalDistance.value);
    outfile.println("Number of routes : ", nbTrucks);
    for [k in 0...nbTrucks] {
        outfile.print("Sequence of truck ", k + 1, ": ");
        local sequence = edgesSequences[k].value;
        for [i in 0...sequence.count()] {
            outfile.print("(", requiredEdgesOrigin[sequence[i]], 
                    ", ", requiredEdgesDestination[sequence[i]], ")  ");
        }
        outfile.println();
    }
}

function readData() {
    local inFile = io.openRead(inFileName);
    for [_ in 0...2] inFile.readln();
    nbNodes = inFile.readln().trim().split(":")[1].toInt();
    nbRequiredEdges = inFile.readln().trim().split(":")[1].toInt();
    nbNotRequiredEdges = inFile.readln().trim().split(":")[1].toInt();
    nbTrucks = inFile.readln().trim().split(":")[1].toInt();
    truckCapacity = inFile.readln().trim().split(":")[1].toInt();
    for [_ in 0...2] inFile.readln();

    inFile.readln();
    requiredNodes[0] = -1;
    for [i in 0...nbNodes][j in 0...nbNodes] nodeNeighbors[i][j] = 0;
    for [i in 0...nbRequiredEdges] {
        local elements = inFile.readln().trim().split("   ");
        local firstNode = elements[0].trim().split(",")[0].trim().split()[1].toInt();
        local secondNode = elements[0].trim().split(",")[1].trim().replace(")", "").toInt();
        local cost = elements[1].trim().split()[1].toInt();
        local demand = elements[2].trim().split()[1].toInt();
        // Even indices store direct edges, and odd indices store reverse edges
        requiredEdgesCosts[2 * i] = cost;
        requiredEdgesDemands[2 * i] = demand;
        requiredEdgesCosts[2 * i + 1] = cost;
        requiredEdgesDemands[2 * i + 1] = demand;
        requiredEdgesOrigin[2 * i] = firstNode;
        requiredEdgesDestination[2 * i] = secondNode;
        requiredEdgesOrigin[2 * i + 1] = secondNode;
        requiredEdgesDestination[2 * i + 1] = firstNode;
        addRequiredNodes(firstNode);
        addRequiredNodes(secondNode);
        nodeNeighbors[firstNode - 1][secondNode - 1] = cost;
        nodeNeighbors[secondNode - 1][firstNode - 1] = cost;
    }

    if (nbNotRequiredEdges > 0) {
        inFile.readln();
        for [i in 0...nbNotRequiredEdges] {
            local elements = inFile.readln().trim().split("   ");
            local firstNode = elements[0].trim().split(",")[0].trim().split()[1].toInt();
            local secondNode = elements[0].trim().split(",")[1].trim().replace(")", "").toInt();
            local cost = elements[1].trim().split()[1].toInt();
            nodeNeighbors[firstNode - 1][secondNode - 1] = cost;
            nodeNeighbors[secondNode - 1][firstNode - 1] = cost;
        }
    }
    depotNode = inFile.readln().trim().split(":")[1].toInt();
    nbRequiredNodes = requiredNodes.count();
    for [i in 0...nbNodes] {
        shortestPath[i] = 0;
    }
    findRequiredPaths(shortestPath);
    findDistanceBetweenEdges();
    findDistanceToDepot();
    findDistanceFromDepot(shortestPath);
}

/* Finds the shortest path from one node "origin" to all the other nodes of the graph
   thanks to the Dijkstra's algorithm */
function minDistance(shortestPath, sptSet) {
    minValue = +inf;
    minIndex = -1;
    for [i in 0...nbNodes] {
        if (sptSet[i] == false && shortestPath[i] <= minValue) {
            minValue = shortestPath[i];
            minIndex = i;
        }
    }
    return minIndex;
}

function shortestPathFinder(origin, shortestPath) {
    local positiveInfinity =  +inf;
    for [i in 0...nbNodes] {
        shortestPath[i] = positiveInfinity;
        sptSet[i] = false;
    }
    shortestPath[origin - 1] = 0;
    for [_ in 0...nbNodes] {
        currentNode = minDistance(shortestPath, sptSet);
        sptSet[currentNode] = true;
        currentNeighbors = nodeNeighbors[currentNode];
        for [neighbor in 0...nbNodes] {
            if (currentNeighbors[neighbor] != 0) {
                distance = currentNeighbors[neighbor];
                if (!sptSet[neighbor] &&
                        shortestPath[currentNode] + distance < shortestPath[neighbor]) {
                    shortestPath[neighbor] = shortestPath[currentNode] + distance;
                }
            }
        }
    }
}

function findRequiredPaths(shortestPath) {
    for [i in 0...nbRequiredNodes] {
        node = requiredNodes[i];
        shortestPathFinder(node, shortestPath);
        for [j in 0...nbNodes] {
            requiredDistance[i][j] = shortestPath[j];
        }
    }
}

function findDistanceBetweenEdges() {
    for [i in 0...2*nbRequiredEdges] {
        for [j in 0...2*nbRequiredEdges] {
            if (requiredEdgesDestination[i] == requiredEdgesOrigin[j]) {
                distanceBetweenEdges[i][j] = 0;
            } else {
                for [k in 0...nbRequiredNodes] {
                    if (requiredNodes[k] == requiredEdgesDestination[i]) {
                        distanceBetweenEdges[i][j] = requiredDistance[k][requiredEdgesOrigin[j] - 1];
                    }
                }
            }
        }
    }
}

function findDistanceToDepot() {
    for [i in 0...2*nbRequiredEdges] {
        if (requiredEdgesDestination[i] == depotNode) {
            distanceToDepot[i] = 0;
        } else {
            for [k in 0...nbRequiredNodes] {
                if (requiredNodes[k] == requiredEdgesDestination[i]) {
                    distanceToDepot[i] = requiredDistance[k][depotNode - 1];
                }
            }
        }
    }
}

function findDistanceFromDepot(shortestPath) {
    for [i in 0...2*nbRequiredEdges] {
        if (depotNode == requiredEdgesOrigin[i]) {
            distanceFromDepot[i] = 0;
        } else {
            depotIsRequiredNode = false;
            for [k in 0...nbRequiredNodes] {
                if (requiredNodes[k] == depotNode) {
                    depotIsRequiredNode = true;
                    distanceFromDepot[i] = requiredDistance[k][requiredEdgesOrigin[i] - 1];
                }
            }
            if (!depotIsRequiredNode) {
                shortestPathFinder(depotNode, shortestPath);
                distanceFromDepot[i] = shortestPath[requiredEdgesOrigin[i] - 1];
            }
        }
    }
}

function addRequiredNodes(node) {
    local n = requiredNodes.count();
    local found = false;
    for [i in 0...n] {
        if (requiredNodes[i] == node) {
            found = true;
        }
    }
    if (!found) {
        if (requiredNodes[0] == -1) {
            requiredNodes[0] = node;
        } else {
            requiredNodes[n] = node;
        }
    }
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python capacitated_arc_routing.py instances\egl-e1-A.dat
Execution (Linux)
export PYTHONPATH=/opt/localsolver_13_0/bin/python
python capacitated_arc_routing.py instances/egl-e1-A.dat
import sys
import localsolver


def main(instance_file, output_file, str_time_limit):
    #
    # Read instance data
    #
    instance = CarpInstance(instance_file)

    with localsolver.LocalSolver() as ls:
        #
        # Declare the optimization model
        #
        model = ls.model

        # Sequence of edges visited and "serviced" by each truck
        edges_sequences_vars = [model.list(2 * instance.nb_required_edges)
                                for _ in range(instance.nb_trucks)]
        edges_sequences = model.array(edges_sequences_vars)

        # Create distance and cost arrays to be able to access it with an "at" operator
        demands = model.array(instance.demands_data)
        costs = model.array(instance.costs_data)
        dist_from_depot = model.array(instance.dist_from_depot_data)
        dist_to_depot = model.array(instance.dist_to_depot_data)
        edges_dist = model.array()
        for n in range(2 * instance.nb_required_edges):
            edges_dist.add_operand(model.array(instance.edges_dist_data[n]))

        # An edge must be serviced by at most one truck
        model.constraint(model.disjoint(edges_sequences))

        # An edge can be travelled in both directions but its demand must be 
        # satisfied only once
        for i in range(instance.nb_required_edges):
            model.constraint(
                model.contains(edges_sequences, 2 * i)
                + model.contains(edges_sequences, 2 * i + 1) 
                == 1)

        route_distances = [None] * instance.nb_trucks
        for k in range(instance.nb_trucks):
            sequence = edges_sequences_vars[k]
            c = model.count(sequence)

            # Quantity in each truck
            demand_lambda = model.lambda_function(lambda j: demands[j])
            route_quantity = model.sum(sequence, demand_lambda)
            # Capacity constraint : a truck must not exceed its capacity
            model.constraint(route_quantity <= instance.truck_capacity)

            # Distance travelled by each truck
            dist_lambda = model.lambda_function(
                lambda i:
                    costs[sequence[i]]
                    + model.at(
                        edges_dist,
                        sequence[i - 1],
                        sequence[i]))
            route_distances[k] = model.sum(model.range(1, c), dist_lambda) \
                + model.iif(
                    c > 0,
                    costs[sequence[0]] + dist_from_depot[sequence[0]] \
                        + dist_to_depot[sequence[c - 1]],
                    0)

        # Total distance travelled
        total_distance = model.sum(route_distances)

        # Objective: minimize the distance travelled
        model.minimize(total_distance)

        model.close()

        # Parameterize the solver
        ls.param.time_limit = int(str_time_limit)

        ls.solve()

        #
        # Write the solution in a file with the following format:
        #  - total distance
        #  - number of routes
        #  - for each truck, the edges visited
        #
        if output_file:
            with open(output_file, 'w') as f:
                f.write("Objective function value : {}\nNumber of routes : {}\n".format(
                    total_distance.value, instance.nb_trucks))
                for k in range(instance.nb_trucks):
                    f.write("Sequence of truck {}: ".format(k + 1))
                    sequence = edges_sequences_vars[k].value
                    c = len(sequence)
                    for i in range(c):
                        f.write("({}, {})  ".format(instance.origins_data[sequence[i]],
                                instance.destinations_data[sequence[i]]))
                    f.write("\n")


class CarpInstance:

    def read_elem(self, filename):
        with open(filename) as f:
            return [str(elem) for elem in f.read().strip().split("\n")]

    # The input files follow the format of the DIMACS challenge
    def __init__(self, filename):
        file_it = iter(self.read_elem(filename))
        for _ in range(2):
            next(file_it)
        nb_nodes = int(next(file_it).strip().split(":")[1])
        self.nb_required_edges = int(next(file_it).strip().split(":")[1])
        nb_not_required_edges = int(next(file_it).strip().split(":")[1])
        self.nb_trucks = int(next(file_it).strip().split(":")[1])
        self.truck_capacity = int(next(file_it).strip().split(":")[1])
        for _ in range(3):
            next(file_it)
        self.demands_data = list()
        self.costs_data = list()
        self.origins_data = list()
        self.destinations_data = list()
        required_nodes = list()
        node_neighbors = [([0] * nb_nodes) for _ in range(nb_nodes)]
        for _ in range(self.nb_required_edges):
            elements = next(file_it)
            edge = tuple(map(int, elements.strip().split("   ")[0][2:-1].strip().split(",")))
            cost = int(elements.strip().split("   ")[1].strip().split()[1])
            demand = int(elements.strip().split("   ")[2].strip().split()[1])
            for _ in range(2):
                self.costs_data.append(cost)
                self.demands_data.append(demand)
            # even indices store direct edges, and odd indices store reverse edges
            self.origins_data.append(edge[0])
            self.destinations_data.append(edge[1])
            self.origins_data.append(edge[1])
            self.destinations_data.append(edge[0])
            if edge[0] not in required_nodes:
                required_nodes.append(edge[0])
            if edge[1] not in required_nodes:
                required_nodes.append(edge[1])
            node_neighbors[edge[0] - 1][edge[1] - 1] = cost
            node_neighbors[edge[1] - 1][edge[0] - 1] = cost
        if nb_not_required_edges > 0:
            next(file_it)
            for _ in range(nb_not_required_edges):
                elements = next(file_it)
                edge = tuple(map(int, elements.strip().split("   ")[0][2:-1].strip().split(",")))
                cost = int(elements.strip().split("   ")[1].strip().split()[1])
                node_neighbors[edge[0] - 1][edge[1] - 1] = cost
                node_neighbors[edge[1] - 1][edge[0] - 1] = cost
        depot_node = int(next(file_it).strip().split(":")[1])
        # Finds the shortest path from one "required node" to another
        nb_required_nodes = len(required_nodes)
        required_distances = list()
        for node in required_nodes:
            paths = self.shortest_path_finder(node, nb_nodes, node_neighbors)
            required_distances.append(paths)
        # Since we can explore the edges in both directions, we will represent all possible
        # edges with an index
        self.edges_dist_data = None
        self.find_distance_between_edges(nb_required_nodes, required_nodes, required_distances)
        self.dist_to_depot_data = None
        self.find_distance_to_depot(nb_required_nodes, depot_node,
                                    required_nodes, required_distances)
        self.dist_from_depot_data = None
        self.find_distance_from_depot(nb_required_nodes, nb_nodes, depot_node,
                                      required_nodes, required_distances, node_neighbors)

    # Finds the shortest path from one node "origin" to all the other nodes of the graph
    # thanks to the Dijkstra's algorithm
    def min_distance(self, nb_nodes, shortest_path, sptSet):
        min = sys.maxsize
        for i in range(nb_nodes):
            if shortest_path[i] < min and sptSet[i] == False:
                min = shortest_path[i]
                min_index = i
        return min_index

    def shortest_path_finder(self, origin, nb_nodes, node_neighbors):
        shortest_path = [sys.maxsize] * nb_nodes
        shortest_path[origin - 1] = 0
        sptSet = [False] * nb_nodes
        for _ in range(nb_nodes):
            current_node = self.min_distance(nb_nodes, shortest_path, sptSet)
            sptSet[current_node] = True
            current_neighbors = node_neighbors[current_node]
            for neighbor in range(nb_nodes):
                if current_neighbors[neighbor] != 0:
                    distance = current_neighbors[neighbor]
                    if ((sptSet[neighbor] == False) and
                            (shortest_path[current_node] + distance < shortest_path[neighbor])):
                        shortest_path[neighbor] = distance + shortest_path[current_node]
        return shortest_path

    def find_distance_between_edges(self, nb_required_nodes, required_nodes, required_distances):
        self.edges_dist_data = [[None] * (2 * self.nb_required_edges)
                                for _ in range(2 * self.nb_required_edges)]
        for i in range(2 * self.nb_required_edges):
            for j in range(2 * self.nb_required_edges):
                if self.destinations_data[i] == self.origins_data[j]:
                    self.edges_dist_data[i][j] = 0
                else:
                    for k in range(nb_required_nodes):
                        if required_nodes[k] == self.destinations_data[i]:
                            self.edges_dist_data[i][j] = required_distances[k][
                                self.origins_data[j] - 1]

    def find_distance_to_depot(
            self, nb_required_nodes, depot_node, required_nodes, required_distances):
        self.dist_to_depot_data = [None] * (2 * self.nb_required_edges)
        for i in range(2 * self.nb_required_edges):
            if self.destinations_data[i] == depot_node:
                self.dist_to_depot_data[i] = 0
            else:
                for k in range(nb_required_nodes):
                    if required_nodes[k] == self.destinations_data[i]:
                        self.dist_to_depot_data[i] = required_distances[k][depot_node-1]

    def find_distance_from_depot(self, nb_required_nodes, nb_nodes, depot_node,
                                 required_nodes, required_distances, node_neighbors):
        self.dist_from_depot_data = [None] * (2 * self.nb_required_edges)
        for i in range(2 * self.nb_required_edges):
            if depot_node == self.origins_data[i]:
                self.dist_from_depot_data[i] = 0
            else:
                depot_is_required_node = False
                for k in range(nb_required_nodes):
                    if required_nodes[k] == depot_node:
                        depot_is_required_node = True
                        self.dist_from_depot_data[i] = required_distances[k][
                            self.origins_data[i] - 1]
                if not depot_is_required_node:
                    shortest_paths_from_depot = self.shortest_path_finder(
                        depot_node, nb_nodes, node_neighbors)
                    self.dist_from_depot_data[i] = shortest_paths_from_depot[self.origins_data[i] - 1]


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python capacitated_arc_routing.py input_file [output_file] [time_limit]")
        sys.exit(1)

    instance_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) > 2 else None
    str_time_limit = sys.argv[3] if len(sys.argv) > 3 else "20"

    main(instance_file, output_file, str_time_limit)
Compilation / Execution (Windows)
cl /EHsc capacitated_arc_routing.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver130.lib
capacitated_arc_routing instances\egl-e1-A.dat
Compilation / Execution (Linux)
g++ capacitated_arc_routing.cpp -I/opt/localsolver_13_0/include -llocalsolver130 -lpthread -o capacitated_arc_routing
./capacitated_arc_routing instances/egl-e1-A.dat
#include "localsolver.h"
#include <algorithm>
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace localsolver;
using namespace std;

class CapacitatedArcRouting {

public:
    // Data from the problem
    int nbTrucks;
    int truckCapacity;
    int nbNodes;
    int nbRequiredEdges;
    int nbRequiredNodes;
    int nbNotRequiredNodes;
    int depotNode;
    vector<int> demandsData;
    vector<int> costsData;
    vector<int> originsData;
    vector<int> destinationsData;
    vector<int> requiredNodes;
    vector<vector<int>> nodesDistData;
    vector<vector<int>> edgesDistData;
    vector<int> distToDepotData;
    vector<int> distFromDepotData;
    vector<vector<int>> nodeNeighborsData;

    // LocalSolver
    LocalSolver localsolver;

    // Decision variables
    vector<LSExpression> edgesSequencesVars;
    LSExpression totalDistance;

    /* Read instance data */
    void readInstance(const string& fileName) { readCarp(fileName); }

    void solve(int limit) {
        // Declare the optimization model
        LSModel model = localsolver.getModel();

        // Sequence of edges visited and "serviced" by each truck
        edgesSequencesVars.resize(nbTrucks);
        for (int i = 0; i < nbTrucks; ++i) {
            // A sequence can contain an edge in one direction or its reverse
            edgesSequencesVars[i] = model.listVar(2 * nbRequiredEdges);
        }
        LSExpression edgesSequences = model.array(edgesSequencesVars.begin(), edgesSequencesVars.end());

        // Create distance and cost arrays to be able to access it with an "at" operator
        LSExpression demands = model.array(demandsData.begin(), demandsData.end());
        LSExpression costs = model.array(costsData.begin(), costsData.end());
        LSExpression distFromDepot = model.array(distFromDepotData.begin(), distFromDepotData.end());
        LSExpression distToDepot = model.array(distToDepotData.begin(), distToDepotData.end());
        LSExpression edgesDist = model.array();
        for (int n = 0; n < 2 * nbRequiredEdges; ++n) {
            edgesDist.addOperand(model.array(edgesDistData[n].begin(), edgesDistData[n].end()));
        }

        // An edge must be serviced by at most one truck
        model.constraint(model.disjoint(edgesSequences));

        // An edge can be travelled in both directions but its demand must be satisfied only once
        for (int i = 0; i < nbRequiredEdges; ++i) {
            model.constraint(model.contains(edgesSequences, 2 * i) + model.contains(edgesSequences, 2 * i + 1) == 1);
        }

        vector<LSExpression> routeDistances(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k) {
            LSExpression sequence = edgesSequencesVars[k];
            LSExpression c = model.count(sequence);

            // Quantity in each truck
            LSExpression demandLambda =
                model.createLambdaFunction([&](LSExpression j) { return demands[j]; });
            LSExpression routeQuantity = model.sum(sequence, demandLambda);
            // Capacity constraint : a truck must not exceed its capacity
            model.constraint(routeQuantity <= truckCapacity);

            // Distance travelled by each truck
            LSExpression distLambda = model.createLambdaFunction(
                [&](LSExpression i) { return costs[sequence[i]] + model.at(edgesDist, sequence[i - 1], sequence[i]); });
            routeDistances[k] =
                model.sum(model.range(1, c), distLambda) +
                model.iif(c > 0, costs[sequence[0]] + distFromDepot[sequence[0]] + distToDepot[sequence[c - 1]], 0);
        }

        // Total distance travelled
        totalDistance = model.sum(routeDistances.begin(), routeDistances.end());

        // Objective : minimize the distance travelled
        model.minimize(totalDistance);

        model.close();

        // Parametrize the solver
        localsolver.getParam().setTimeLimit(limit);

        localsolver.solve();
    }

    /* Write the solution in a file with the following format:
     * - total distance
     * - number of routes
     * - for each truck, the edges visited */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << "Objective function value : " << totalDistance.getValue() << endl;
        outfile << "Number of routes : " << nbTrucks << endl;
        for (int k = 0; k < nbTrucks; ++k) {
            outfile << "Sequence of truck " << k + 1 << ": ";
            LSCollection sequence = edgesSequencesVars[k].getCollectionValue();
            int c = sequence.count();
            for (int i = 0; i < c; ++i) {
                outfile << "(" << originsData[sequence[i]] << ", " << destinationsData[sequence[i]] << ")  ";
            }
            outfile << endl;
        }
    }

private:
    // The input files follow the format of the DIMACS challenge
    void readCarp(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }

        string str;
        char* pch;
        char* line;
        int first = -1;
        int second = -1;
        for (int i = 0; i < 2; ++i)
            getline(infile, str);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbNodes = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbRequiredEdges = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbNotRequiredNodes = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbTrucks = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        truckCapacity = atoi(pch);
        for (int i = 0; i < 3; ++i)
            getline(infile, str);
        nodeNeighborsData.resize(nbNodes);
        for (int i = 0; i < nbNodes; ++i) {
            nodeNeighborsData[i].resize(nbNodes);
        }
        for (int i = 0; i < nbNodes; ++i) {
            for (int j = 0; j < nbNodes; ++j) {
                nodeNeighborsData[i][j] = 0;
            }
        }
        getline(infile, str);
        for (int i = 0; i < nbRequiredEdges; ++i) {
            extractEdge(str, pch, line, first, second, true);
            extractCost(str, pch, line, first, second, true);
            extractDemand(str, pch, line);
            getline(infile, str);
        }
        if (nbNotRequiredNodes > 0) {
            getline(infile, str);
            for (int i = 0; i < nbNotRequiredNodes; ++i) {
                extractEdge(str, pch, line, first, second, false);
                extractCost(str, pch, line, first, second, false);
                getline(infile, str);
            }
        }
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        depotNode = atoi(pch);
        nbRequiredNodes = requiredNodes.size();
        vector<int> shortestPath;
        shortestPath.resize(nbNodes);
        for (int i = 0; i < nbRequiredNodes; ++i) {
            shortestPathFinder(requiredNodes[i], shortestPath);
            for (int j = 0; j < nbNodes; ++j) {
                if (shortestPath[j] == std::numeric_limits<int>::max()) {
                }
            }
            nodesDistData.push_back(shortestPath);
        }
        findDistanceBetweenEdges();
        findDistanceToDepot();
        findDistanceFromDepot(shortestPath);
    }

    // Finds the shortest path from one node "origin" to all the other nodes of the graph
    // thanks to the Dijkstra's algorithm.
    int minDistance(vector<int> shortestPath, vector<bool> sptSet) {
        // Initializes min value
        int min = std::numeric_limits<int>::max(), minIndex = -1;
        for (int i = 0; i < nbNodes; ++i) {
            if (sptSet[i] == false && shortestPath[i] <= min) {
                min = shortestPath[i];
                minIndex = i;
            }
        }
        return minIndex;
    }

    void shortestPathFinder(int origin, vector<int>& shortestPath) {
        vector<int> currentNeighbors;
        int currentNode;
        int distance;
        vector<bool> sptSet;
        sptSet.resize(nbNodes);
        for (int i = 0; i < nbNodes; ++i) {
            shortestPath[i] = std::numeric_limits<int>::max();
            sptSet[i] = false;
        }
        shortestPath[origin - 1] = 0;
        for (int count = 0; count < nbNodes; ++count) {
            currentNode = minDistance(shortestPath, sptSet);
            sptSet[currentNode] = true;
            currentNeighbors = nodeNeighborsData[currentNode];
            for (int neighbor = 0; neighbor < nbNodes; ++neighbor) {
                if (currentNeighbors[neighbor] != 0) {
                    distance = currentNeighbors[neighbor];
                    if ((sptSet[neighbor] == false) &&
                        (shortestPath[currentNode] + distance < shortestPath[neighbor])) {
                        shortestPath[neighbor] = shortestPath[currentNode] + distance;
                    }
                }
            }
        }
    }

    void findDistanceBetweenEdges() {
        edgesDistData.resize(2 * nbRequiredEdges);
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            edgesDistData[i].resize(2 * nbRequiredEdges);
        }
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            for (int j = 0; j < 2 * nbRequiredEdges; ++j) {
                if (destinationsData[i] == originsData[j]) {
                    edgesDistData[i][j] = 0;
                } else {
                    for (int k = 0; k < nbRequiredNodes; ++k) {
                        if (requiredNodes[k] == destinationsData[i]) {
                            edgesDistData[i][j] = nodesDistData[k][originsData[j] - 1];
                        }
                    }
                }
            }
        }
    }

    void findDistanceToDepot() {
        distToDepotData.resize(2 * nbRequiredEdges);
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (destinationsData[i] == depotNode) {
                distToDepotData[i] = 0;
            } else {
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes[k] == destinationsData[i]) {
                        distToDepotData[i] = nodesDistData[k][depotNode - 1];
                    }
                }
            }
        }
    }

    void findDistanceFromDepot(vector<int>& shortestPath) {
        distFromDepotData.resize(2 * nbRequiredEdges);
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (depotNode == originsData[i]) {
                distFromDepotData[i] = 0;
            } else {
                bool depotIsRequiredNode = false;
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes[k] == depotNode) {
                        depotIsRequiredNode = true;
                        distFromDepotData[i] = nodesDistData[k][originsData[i] - 1];
                    }
                }
                if (!depotIsRequiredNode) {
                    shortestPathFinder(depotNode, shortestPath);
                    distFromDepotData[i] = shortestPath[originsData[i] - 1];
                }
            }
        }
    }

    void extractEdge(const string& str, char* pch, char* line, int& first, int& second, bool req) {
        line = strdup(str.c_str());
        pch = strtok(line, " ");
        if (strcmp(pch, "(") == 0) {
            pch = strtok(NULL, ",");
            first = atoi(pch);
            pch = strtok(NULL, ")");
            second = atoi(pch);
        }

        if (req == true) {
            // Even indices store direct edges, and odd indices store reverse edges
            originsData.push_back(first);
            destinationsData.push_back(second);
            originsData.push_back(second);
            destinationsData.push_back(first);
            if (find(requiredNodes.begin(), requiredNodes.end(), first) == requiredNodes.end()) {
                requiredNodes.push_back(first);
            }
            if (find(requiredNodes.begin(), requiredNodes.end(), second) == requiredNodes.end()) {
                requiredNodes.push_back(second);
            }
        }
    }

    void extractCost(const string& str, char* pch, char* line, int& first, int& second, bool req) {
        line = strdup(str.c_str());
        pch = strtok(line, "coste");
        int cost;
        pch = strtok(NULL, "coste");
        cost = atoi(pch);
        if (req == true) {
            for (int i = 0; i < 2; ++i)
                costsData.push_back(cost);
        }
        nodeNeighborsData[first - 1][second - 1] = cost;
        nodeNeighborsData[second - 1][first - 1] = cost;
    }

    void extractDemand(const string& str, char* pch, char* line) {
        line = strdup(str.c_str());
        pch = strtok(line, "demanda");
        int demand;
        pch = strtok(NULL, "demanda");
        demand = atoi(pch);
        for (int i = 0; i < 2; ++i)
            demandsData.push_back(demand);
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: capacitated_arc_routing inputFile [outputFile] [timeLimit]" << endl;
        return 1;
    }

    const char* instanceFile = argv[1];
    const char* outputFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "20";

    try {
        CapacitatedArcRouting model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (outputFile != NULL)
            model.writeSolution(outputFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %LS_HOME%\bin\localsolvernet.dll .
csc CapacitatedArcRouting.cs /reference:localsolvernet.dll
CapacitatedArcRouting instances\egl-e1-A.dat
using System;
using System.IO;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using localsolver;

public class CapacitatedArcRouting : IDisposable
{
    // Data from the problem
    int nbTrucks;
    int truckCapacity;
    int nbNodes;
    int nbRequiredEdges;
    int nbRequiredNodes;
    int nbNotRequiredEdges;
    int depotNode;
    int[] demandsData;
    int[] costsData;
    List<int> requiredNodes;
    int[] originsData;
    int[] destinationsData;
    int[] distToDepotData;
    int[] distFromDepotData;
    int[][] nodesDistData;
    int[][] edgesDistData;
    int[][] nodeNeighborsData;

    // LocalSolver
    LocalSolver localsolver;

    // Decision variables
    LSExpression[] edgesSequencesVars;
    LSExpression totalCost;

    public CapacitatedArcRouting()
    {
        localsolver = new LocalSolver();
    }

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputCarp(fileName);
    }

    public void Dispose()
    {
        if (localsolver != null)
            localsolver.Dispose();
    }

    void Solve(int limit)
    {
        // Declare the optimization model
        LSModel model = localsolver.GetModel();

        // Sequence of edges visited and "serviced" by each truck
        edgesSequencesVars = new LSExpression[nbTrucks];
        for (int i = 0; i < nbTrucks; ++i)
        {
            // A sequence can contain an edge in one direction or its reverse
            edgesSequencesVars[i] = model.List(2 * nbRequiredEdges);
        }
        LSExpression edgesSequences = model.Array(edgesSequencesVars);

        // Creates distance and cost arrays to be able to access it with an "at" operator
        LSExpression demands = model.Array(demandsData);
        LSExpression costs = model.Array(costsData);
        LSExpression distFromDepot = model.Array(distFromDepotData);
        LSExpression distToDepot = model.Array(distToDepotData);
        LSExpression edgesDist = model.Array(edgesDistData);

        // An edge must be serviced by at most one truck
        model.Constraint(model.Disjoint(edgesSequences));

        // An edge can be travelled in both directions but its demand must be satisfied only once
        for (int i = 0; i < nbRequiredEdges; ++i)
        {
            model.Constraint(
                model.Contains(edgesSequences, 2 * i) + model.Contains(edgesSequences, 2 * i + 1)
                    == 1
            );
        }

        LSExpression[] routeDistances = new LSExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; ++k)
        {
            LSExpression sequence = edgesSequencesVars[k];
            LSExpression c = model.Count(sequence);

            // Quantity in each truck
            LSExpression demandLambda = model.LambdaFunction(j => demands[j]);
            LSExpression routeQuantity = model.Sum(sequence, demandLambda);
            // Capacity constraint : a truck must not exceed its capacity
            model.Constraint(routeQuantity <= truckCapacity);

            // Distance travelled by each truck
            LSExpression distLambda = model.LambdaFunction(
                i => costs[sequence[i]] + edgesDist[sequence[i - 1], sequence[i]]
            );
            routeDistances[k] =
                model.Sum(model.Range(1, c), distLambda)
                + model.If(
                    c > 0,
                    costs[sequence[0]] + distFromDepot[sequence[0]] + distToDepot[sequence[c - 1]],
                    0
                );
        }

        // Total Distance travelled
        totalCost = model.Sum(routeDistances);

        // Objective : minimize the distance travelled
        model.Minimize(totalCost);

        model.Close();

        // Parametrize the solver
        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.Solve();
    }

    /* Write the solution in a file with the following format:
    - total distance
    - number of routes
    - for each truck, the edges visited */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("Objective function value : " + totalCost.GetValue());
            output.WriteLine("Number of routes : " + nbTrucks);
            for (int k = 0; k < nbTrucks; ++k)
            {
                output.Write("Sequence of truck " + (k + 1) + ": ");
                LSCollection sequence = edgesSequencesVars[k].GetCollectionValue();
                int c = sequence.Count();
                for (int i = 0; i < c; ++i)
                {
                    output.Write(
                        "("
                            + originsData[sequence[i]]
                            + ", "
                            + destinationsData[sequence[i]]
                            + ")  "
                    );
                }
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: CapacitatedArcRouting inputFile [solFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "20";

        using (CapacitatedArcRouting model = new CapacitatedArcRouting())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

    // The input files follow the format of the DIMACS challenge
    private void ReadInputCarp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            for (int i = 0; i < 2; ++i)
                input.ReadLine();
            splitted = input.ReadLine().Split(':');
            nbNodes = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            nbRequiredEdges = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            nbNotRequiredEdges = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            nbTrucks = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            truckCapacity = int.Parse(splitted[1]);
            for (int i = 0; i < 2; ++i)
                input.ReadLine();
            // Even indices store direct edges, and odd indices store reverse edges
            originsData = new int[2 * nbRequiredEdges];
            destinationsData = new int[2 * nbRequiredEdges];
            demandsData = new int[2 * nbRequiredEdges];
            costsData = new int[2 * nbRequiredEdges];
            nodeNeighborsData = new int[nbNodes][];
            for (int i = 0; i < nbNodes; ++i)
                nodeNeighborsData[i] = new int[nbNodes];
            for (int i = 0; i < nbNodes; ++i)
            {
                for (int j = 0; j < nbNodes; ++j)
                    nodeNeighborsData[i][j] = 0;
            }
            requiredNodes = new List<int>();
            input.ReadLine();
            for (int i = 0; i < nbRequiredEdges; ++i)
            {
                splitted = input.ReadLine().Split(' ');
                int firstNode = int.Parse(splitted[2].Remove(splitted[2].Length - 1));
                int secondNode = int.Parse(splitted[3].Remove(splitted[3].Length - 1));
                int cost = int.Parse(splitted[7]);
                int demand = int.Parse(splitted[11]);
                costsData[2 * i] = cost;
                demandsData[2 * i] = demand;
                costsData[2 * i + 1] = cost;
                demandsData[2 * i + 1] = demand;
                originsData[2 * i] = firstNode;
                destinationsData[2 * i] = secondNode;
                originsData[2 * i + 1] = secondNode;
                destinationsData[2 * i + 1] = firstNode;
                if (!requiredNodes.Contains(firstNode))
                    requiredNodes.Add(firstNode);
                if (!requiredNodes.Contains(secondNode))
                    requiredNodes.Add(secondNode);
                nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
            }
            if (nbNotRequiredEdges > 0)
            {
                input.ReadLine();
                for (int i = 0; i < nbNotRequiredEdges; ++i)
                {
                    splitted = input.ReadLine().Split(' ');
                    int firstNode = int.Parse(splitted[2].Remove(splitted[2].Length - 1));
                    int secondNode = int.Parse(splitted[3].Remove(splitted[3].Length - 1));
                    int cost = int.Parse(splitted[7]);
                    nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                    nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
                }
            }
            splitted = input.ReadLine().Split(':');
            depotNode = int.Parse(splitted[1]);
        }
        nbRequiredNodes = requiredNodes.Count;
        int[] shortestPath = new int[nbNodes];
        FindRequiredPaths(ref shortestPath);
        FindDistanceBetweenEdges();
        FindDistanceToDepot();
        FindDistanceFromDepot(ref shortestPath);
    }

    // Finds the shortest path from one node "origin" to all the other nodes of the graph
    // thanks to the Dijkstra's algorithm
    public int MinDistance(int[] shortestPath, bool[] sptSet)
    {
        // Initializes min value
        int min = int.MaxValue,
            minIndex = -1;
        for (int i = 0; i < nbNodes; ++i)
            if (sptSet[i] == false && shortestPath[i] <= min)
            {
                min = shortestPath[i];
                minIndex = i;
            }
        return minIndex;
    }

    public void ShortestPathFinder(int origin, ref int[] shortestPath)
    {
        int[] currentNeighbors;
        int currentNode,
            distance;
        bool[] sptSet = new bool[nbNodes];
        for (int i = 0; i < nbNodes; ++i)
        {
            shortestPath[i] = int.MaxValue;
            sptSet[i] = false;
        }
        shortestPath[origin - 1] = 0;
        for (int count = 0; count < nbNodes; ++count)
        {
            currentNode = MinDistance(shortestPath, sptSet);
            sptSet[currentNode] = true;
            currentNeighbors = nodeNeighborsData[currentNode];
            for (int neighbor = 0; neighbor < nbNodes; ++neighbor)
            {
                if (currentNeighbors[neighbor] != 0)
                {
                    distance = currentNeighbors[neighbor];
                    if (
                        !sptSet[neighbor]
                        && shortestPath[currentNode] + distance < shortestPath[neighbor]
                    )
                    {
                        shortestPath[neighbor] = shortestPath[currentNode] + distance;
                    }
                }
            }
        }
    }

    public void FindRequiredPaths(ref int[] shortestPath)
    {
        nodesDistData = new int[nbRequiredNodes][];
        for (int i = 0; i < nbRequiredNodes; ++i)
            nodesDistData[i] = new int[nbNodes];
        for (int i = 0; i < nbRequiredNodes; ++i)
        {
            int node = requiredNodes[i];
            ShortestPathFinder(node, ref shortestPath);
            for (int j = 0; j < nbNodes; ++j)
                nodesDistData[i][j] = shortestPath[j];
        }
    }

    public void FindDistanceBetweenEdges()
    {
        edgesDistData = new int[2 * nbRequiredEdges][];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
            edgesDistData[i] = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
        {
            for (int j = 0; j < 2 * nbRequiredEdges; ++j)
            {
                if (destinationsData[i] == originsData[j])
                    edgesDistData[i][j] = 0;
                else
                {
                    for (int k = 0; k < nbRequiredNodes; ++k)
                    {
                        if (requiredNodes[k] == destinationsData[i])
                            edgesDistData[i][j] = nodesDistData[k][originsData[j] - 1];
                    }
                }
            }
        }
    }

    public void FindDistanceToDepot()
    {
        distToDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
        {
            if (destinationsData[i] == depotNode)
                distToDepotData[i] = 0;
            else
            {
                for (int k = 0; k < nbRequiredNodes; ++k)
                {
                    if (requiredNodes[k] == destinationsData[i])
                        distToDepotData[i] = nodesDistData[k][depotNode - 1];
                }
            }
        }
    }

    public void FindDistanceFromDepot(ref int[] shortestPath)
    {
        distFromDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
        {
            if (depotNode == originsData[i])
                distFromDepotData[i] = 0;
            else
            {
                bool depotIsRequiredNode = false;
                for (int k = 0; k < nbRequiredNodes; ++k)
                {
                    if (requiredNodes[k] == depotNode)
                    {
                        depotIsRequiredNode = true;
                        distFromDepotData[i] = nodesDistData[k][originsData[i] - 1];
                    }
                }
                if (!depotIsRequiredNode)
                {
                    ShortestPathFinder(depotNode, ref shortestPath);
                    distFromDepotData[i] = shortestPath[originsData[i] - 1];
                }
            }
        }
    }
}
Compilation / Execution (Windows)
javac CapacitatedArcRouting.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. CapacitatedArcRouting instances\egl-e1-A.dat
Compilation / Execution (Linux)
javac CapacitatedArcRouting.java -cp /opt/localsolver_13_0/bin/localsolver.jar
java -cp /opt/localsolver_13_0/bin/localsolver.jar:. CapacitatedArcRouting instances/egl-e1-A.dat
import java.util.*;
import java.io.*;
import localsolver.*;

public class CapacitatedArcRouting {

    // Data from the problem
    int nbTrucks;
    int truckCapacity;
    int nbNodes;
    int nbRequiredEdges;
    int nbRequiredNodes;
    int nbNotRequiredEdges;
    int depotNode;
    int[] demandsData;
    int[] costsData;
    int[] originsData;
    int[] destinationsData;
    List<Integer> requiredNodes;
    int[] distToDepotData;
    int[] distFromDepotData;
    int[][] nodesDistData;
    int[][] edgesDistData;
    int[][] nodeNeighborsData;

    // LocalSolver
    private final LocalSolver localsolver;

    // Decision variables
    LSExpression[] edgesSequencesVars;
    LSExpression totalDistance;

    private CapacitatedArcRouting(LocalSolver localsolver) {
        this.localsolver = localsolver;
    }

    private void solve(int limit) {
        // Declare the optimization model
        LSModel model = localsolver.getModel();

        // Sequence of edges visited and "serviced" by each truck
        edgesSequencesVars = new LSExpression[nbTrucks];
        for (int i = 0; i < nbTrucks; ++i) {
            edgesSequencesVars[i] = model.listVar(2 * nbRequiredEdges);
        }
        LSExpression edgesSequences = model.array(edgesSequencesVars);

        // Create distance and cost arrays to be able to access it with an "at" operator
        LSExpression demands = model.array(demandsData);
        LSExpression costs = model.array(costsData);
        LSExpression distFromDepot = model.array(distFromDepotData);
        LSExpression distToDepot = model.array(distToDepotData);
        LSExpression edgesDist = model.array(edgesDistData);

        // An edge must be serviced by at most one truck
        model.constraint(model.disjoint(edgesSequences));

        // An edge can be travelled in both directions but its demand must be satisfied
        // only once
        for (int i = 0; i < nbRequiredEdges; ++i) {
            model.constraint(model
                .eq(model.sum(model.contains(edgesSequences, 2 * i), model.contains(edgesSequences, (2 * i) + 1)), 1));
        }
        LSExpression[] routeDistances = new LSExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; ++k) {
            LSExpression sequence = edgesSequencesVars[k];
            LSExpression c = model.count(sequence);

            // Quantity in each truck
            LSExpression demandLambda = model.lambdaFunction(j -> model.at(demands, j));
            LSExpression routeQuantity = model.sum(sequence, demandLambda);
            // Capacity constraint : a truck must not exceed its capacity
            model.constraint(model.leq(routeQuantity, truckCapacity));

            // Distance travelled by each truck
            LSExpression distLambda = model.lambdaFunction(i -> model.sum(model.at(costs, model.at(sequence, i)),
                model.at(edgesDist, model.at(sequence, model.sub(i, 1)), model.at(sequence, i))));
            routeDistances[k] = model.sum(model.sum(model.range(1, c), distLambda),
                model.iif(model.gt(c, 0),
                    model.sum(model.at(costs, model.at(sequence, 0)), model.at(distFromDepot, model.at(sequence, 0)),
                        model.at(distToDepot, model.at(sequence, model.sub(c, 1)))),
                    0));
        }

        // Total Distance travelled
        totalDistance = model.sum(routeDistances);

        // Objective : minimize the distance travelled
        model.minimize(totalDistance);

        model.close();

        // Parametrize the solver
        localsolver.getParam().setTimeLimit(limit);

        localsolver.solve();
    }

    // Write the solution in a file with the following format:
    // - total distance
    // - number of routes
    // - for each truck, the edges visited
    private void WriteSolution(String filename) throws IOException {
        try (PrintWriter output = new PrintWriter(filename)) {
            output.println("Objective function value : " + totalDistance.getValue());
            output.println("Number of routes : " + nbTrucks);
            for (int k = 0; k < nbTrucks; ++k) {
                output.print("Sequence of truck " + (k + 1) + ": ");
                LSCollection sequence = edgesSequencesVars[k].getCollectionValue();
                int c = sequence.count();
                for (int i = 0; i < c; ++i) {
                    output.print("(" + originsData[(int) sequence.get(i)] + ", "
                        + destinationsData[(int) sequence.get(i)] + ")  ");
                }
                output.println();
            }
        }
    }

    // The input files follow the format of the DIMACS challenge
    private void readInstance(String filename) throws IOException {
        try (Scanner input = new Scanner(new File(filename))) {
            String[] splitted;
            for (int i = 0; i < 2; ++i)
                input.nextLine();
            splitted = input.nextLine().split(":");
            nbNodes = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            nbRequiredEdges = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            nbNotRequiredEdges = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            nbTrucks = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            truckCapacity = Integer.parseInt(splitted[1].trim());
            for (int i = 0; i < 2; ++i)
                input.nextLine();
            requiredNodes = new ArrayList<Integer>();
            // Even indices store direct edges, and odd indices store reverse edges
            demandsData = new int[2 * nbRequiredEdges];
            costsData = new int[2 * nbRequiredEdges];
            originsData = new int[2 * nbRequiredEdges];
            destinationsData = new int[2 * nbRequiredEdges];
            nodeNeighborsData = new int[nbNodes][];
            for (int j = 0; j < nbNodes; ++j) {
                nodeNeighborsData[j] = new int[nbNodes];
            }
            input.nextLine();
            for (int i = 0; i < nbRequiredEdges; ++i) {
                splitted = input.nextLine().split(" ");
                int firstNode = Integer.parseInt(splitted[2].substring(0, splitted[2].length() - 1));
                int secondNode = Integer.parseInt(splitted[3].substring(0, splitted[3].length() - 1));
                int cost = Integer.parseInt(splitted[7]);
                int demand = Integer.parseInt(splitted[11]);
                costsData[2 * i] = cost;
                demandsData[2 * i] = demand;
                costsData[2 * i + 1] = cost;
                demandsData[2 * i + 1] = demand;
                originsData[2 * i] = firstNode;
                destinationsData[2 * i] = secondNode;
                originsData[2 * i + 1] = secondNode;
                destinationsData[2 * i + 1] = firstNode;
                if (!requiredNodes.contains(firstNode)) {
                    requiredNodes.add(firstNode);
                }
                if (!requiredNodes.contains(secondNode)) {
                    requiredNodes.add(secondNode);
                }
                nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
            }
            if (nbNotRequiredEdges > 0) {
                input.nextLine();
                for (int i = 0; i < nbNotRequiredEdges; ++i) {
                    splitted = input.nextLine().split(" ");
                    int firstNode = Integer.parseInt(splitted[2].substring(0, splitted[2].length() - 1));
                    int secondNode = Integer.parseInt(splitted[3].substring(0, splitted[3].length() - 1));
                    int cost = Integer.parseInt(splitted[7]);
                    nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                    nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
                }
            }
            splitted = input.nextLine().split(":");
            depotNode = Integer.parseInt(splitted[1].trim());
            nbRequiredNodes = requiredNodes.size();
            int[] shortestPath = new int[nbNodes];
            findRequiredPaths(shortestPath);
            findDistanceBetweenEdges();
            findDistanceToDepot();
            findDistanceFromDepot(shortestPath);
        }
    }

    // Finds the shortest path from one node "origin" to all the other nodes of the
    // graph
    // thanks to the Dijkstra's algorithm
    public int minDistance(int[] shortestPath, boolean[] sptSet) {
        int min = Integer.MAX_VALUE, minIndex = -1;
        for (int i = 0; i < nbNodes; ++i)
            if (sptSet[i] == false && shortestPath[i] <= min) {
                min = shortestPath[i];
                minIndex = i;
            }
        return minIndex;
    }

    public void shortestPathFinder(int origin, int[] shortestPath) {
        int[] currentNeighbors;
        int currentNode, distance;
        boolean[] sptSet = new boolean[nbNodes];
        for (int i = 0; i < nbNodes; ++i) {
            shortestPath[i] = Integer.MAX_VALUE;
            sptSet[i] = false;
        }
        shortestPath[origin - 1] = 0;
        for (int count = 0; count < nbNodes; ++count) {
            currentNode = minDistance(shortestPath, sptSet);
            sptSet[currentNode] = true;
            currentNeighbors = nodeNeighborsData[currentNode];
            for (int neighbor = 0; neighbor < nbNodes; ++neighbor) {
                if (currentNeighbors[neighbor] != 0) {
                    distance = currentNeighbors[neighbor];
                    if (!sptSet[neighbor] && shortestPath[currentNode] + distance < shortestPath[neighbor]) {
                        shortestPath[neighbor] = shortestPath[currentNode] + distance;
                    }
                }
            }
        }
    }

    public void findRequiredPaths(int[] shortestPath) {
        nodesDistData = new int[nbRequiredNodes][];
        for (int i = 0; i < nbRequiredNodes; ++i) {
            nodesDistData[i] = new int[nbNodes];
        }
        for (int i = 0; i < nbRequiredNodes; ++i) {
            int node = requiredNodes.get(i);
            shortestPathFinder(node, shortestPath);
            for (int j = 0; j < nbNodes; ++j) {
                nodesDistData[i][j] = shortestPath[j];
            }
        }
    }

    public void findDistanceBetweenEdges() {
        edgesDistData = new int[2 * nbRequiredEdges][];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            edgesDistData[i] = new int[2 * nbRequiredEdges];
        }
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            for (int j = 0; j < 2 * nbRequiredEdges; ++j) {
                if (destinationsData[i] == originsData[j]) {
                    edgesDistData[i][j] = 0;
                } else {
                    for (int k = 0; k < nbRequiredNodes; ++k) {
                        if (requiredNodes.get(k) == destinationsData[i]) {
                            edgesDistData[i][j] = nodesDistData[k][originsData[j] - 1];
                        }
                    }
                }
            }
        }
    }

    public void findDistanceToDepot() {
        distToDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (destinationsData[i] == depotNode) {
                distToDepotData[i] = 0;
            } else {
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes.get(k) == destinationsData[i]) {
                        distToDepotData[i] = nodesDistData[k][depotNode - 1];
                    }
                }
            }
        }
    }

    public void findDistanceFromDepot(int[] shortestPath) {
        distFromDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (depotNode == originsData[i]) {
                distFromDepotData[i] = 0;
            } else {
                boolean depotIsRequiredNode = false;
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes.get(k) == depotNode) {
                        depotIsRequiredNode = true;
                        distFromDepotData[i] = nodesDistData[k][originsData[i] - 1];
                    }
                }
                if (!depotIsRequiredNode) {
                    shortestPathFinder(depotNode, shortestPath);
                    distFromDepotData[i] = shortestPath[originsData[i] - 1];
                }
            }
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java CapacitatedArcRouting inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }
        try (LocalSolver localsolver = new LocalSolver()) {
            String instanceFile = args[0];
            String outputFile = args.length > 1 ? args[1] : null;
            String strTimeLimit = args.length > 2 ? args[2] : "20";

            CapacitatedArcRouting model = new CapacitatedArcRouting(localsolver);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.WriteSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}