Capacitated Arc Routing Problem (CARP)¶

#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;
    }
}