Aircraft Landing

Principles learned

  • Use a lambda expression to define a recursive array

  • Use a list variable to model a sequencing problem

  • Access a multi-dimensional array with an “at” operator

  • Use ternary conditions

Problem

../_images/aircraft_landing.svg

In the aircraft landing problem, landing times of a set of planes have to be scheduled. Each plane can land in a predetermined time window around a target time. The objective is to minimize the “cost penalty” due to landings before or after target times. A separation time has to be respected between two successive planes.

Download the example


Data

The instances provided come from the OR-Library, J.E. Beasly.

The format of the data files is as follows:

  • First line: number of planes and freeze time.

  • From the second line, for each plane i:

    • Appearance time

    • Earliest time of landing

    • Target time of landing

    • Latest time of landing

    • Penalty cost per unit of time for landing before target

    • Penalty cost per unit of time for landing after target

    • For each plane j: separation time required after i lands before j can land

Only the static problem is presented here, so the appearance time variables and the freeze time variable are not used.

Program

This Hexaly model defines a landing order as a list variable. The i-th element of the list corresponds to the index of the i-th plane to land. To ensure that all planes are scheduled, the size of the list variable is constrained to be equal to the number of planes.

Another decision variable is introduced: the preferred landing time for each plane. A plane should land after the target time only if it has to wait after the landing of the preceding plane. Thus, the preferred landing time of a plane is between the earliest landing time and the target time.

The landing time of each plane is defined by a recursive function as the maximum between the preferred landing time and the first time the plane can land, respecting the separation time after the preceding plane.

A last constraint is added to ensure that the landing times respect the separation time with all the previous planes, and not only with the one immediately preceding.

The objective is to minimize the penalty cost due to the earliness or the lateness of each plane. To compute this cost for each plane, a ternary condition is used to select either the earliness cost or the lateness cost depending on the difference between the landing time and the target time.

Execution:
localsolver aircraft_landing.lsp inFileName=instances/airland1.txt [lsTimeLimit=] [solFileName=]
use io;

/* Read instance data */
function input() {
    local usage = "Usage: localsolver aircraft_landing.lsp "
            + "inFileName=inputFile [lsTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;
    local inFile = io.openRead(inFileName);

    nbPlanes = inFile.readInt();
    inFile.readInt(); // Skip freezeTime value

    for [p in 0...nbPlanes] {
        inFile.readInt(); // Skip appearanceTime values
        earliestTime[p] = inFile.readInt();
        targetTime[p] = inFile.readInt();
        latestTime[p] = inFile.readInt();
        earlinessCost[p] = inFile.readDouble();
        tardinessCost[p] = inFile.readDouble();
        separationTime[p][j in 0...nbPlanes] = inFile.readInt();
    }
    inFile.close();
}

/* Declare the optimization model */
function model() {
    // A list variable: landingOrder[i] is the index of the ith plane to land
    landingOrder <- list(nbPlanes);

    // All planes must be scheduled
    constraint count(landingOrder) == nbPlanes;

    // Int variable: preferred landing time for each plane
    preferredTime[p in 0...nbPlanes] <- int(earliestTime[p], targetTime[p]);

    // Landing time for each plane
    landingTime <- array(0...nbPlanes, (p, prev) => max(preferredTime[landingOrder[p]],
            p > 0 ? prev + separationTime[landingOrder[p - 1]][landingOrder[p]] : 0));

    // Landing times must respect the separation time with every previous plane
    for [p in 1...nbPlanes] {
        constraint landingTime[p] >= max[previousPlane in 0...p](landingTime[previousPlane]
                + separationTime[landingOrder[previousPlane]][landingOrder[p]]);
    }

    for [p in 0...nbPlanes] {
        local planeIndex <- landingOrder[p];
        
        // Constraint on latest landing time
        constraint landingTime[p] <= latestTime[planeIndex];

        // Cost for each plane
        cost[p] <- (landingTime[p] < targetTime[planeIndex] ?
                earlinessCost[planeIndex] :
                tardinessCost[planeIndex]) * abs(landingTime[p] - targetTime[planeIndex]);
    }

    // Minimize the total cost
    totalCost <- sum[p in 0...nbPlanes] (cost[p]);
    minimize totalCost;
}

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

/* Write the solution in a file with the following format:
 * - 1st line: value of the objective;
 * - 2nd line: for each position p, index of plane at position p. */
function output() {
    if (solFileName == nil) return;
    local solFile = io.openWrite(solFileName);
    solFile.println(totalCost.isUndefined() ? "(invalid)" : totalCost.value);
    for [p in landingOrder.value]
        solFile.print(p, " ");
    solFile.println();
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python aircraft_landing.py instances\airland1.txt
Execution (Linux)
export PYTHONPATH=/opt/localsolver_13_0/bin/python
python aircraft_landing.py instances/airland1.txt
import localsolver
import sys


def read_elem(filename):
    with open(filename) as f:
        return [str(elem) for elem in f.read().split()]


#
# Read instance data
#
def read_instance(instance_file):
    file_it = iter(read_elem(instance_file))
    nb_planes = int(next(file_it))
    next(file_it)  # Skip freezeTime value
    earliest_time_data = []
    target_time_data = []
    latest_time_data = []
    earliness_cost_data = []
    tardiness_cost_data = []
    separation_time_data = []

    for p in range(nb_planes):
        next(file_it)  # Skip appearanceTime values
        earliest_time_data.append(int(next(file_it)))
        target_time_data.append(int(next(file_it)))
        latest_time_data.append(int(next(file_it)))
        earliness_cost_data.append(float(next(file_it)))
        tardiness_cost_data.append(float(next(file_it)))
        separation_time_data.append([None] * nb_planes)

        for pp in range(nb_planes):
            separation_time_data[p][pp] = int(next(file_it))

    return nb_planes, earliest_time_data, target_time_data, latest_time_data, \
        earliness_cost_data, tardiness_cost_data, separation_time_data


def get_min_landing_time(p, prev, model, separation_time, landing_order):
    return model.iif(
        p > 0,
        prev + model.at(separation_time, landing_order[p - 1], landing_order[p]),
        0)


def main(instance_file, output_file, time_limit):
    nb_planes, earliest_time_data, target_time_data, latest_time_data, \
        earliness_cost_data, tardiness_cost_data, separation_time_data = \
        read_instance(instance_file)

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

        # A list variable: landingOrder[i] is the index of the ith plane to land
        landing_order = model.list(nb_planes)

        # All planes must be scheduled
        model.constraint(model.count(landing_order) == nb_planes)

        # Create LocalSolver arrays to be able to access them with an "at" operator
        target_time = model.array(target_time_data)
        latest_time = model.array(latest_time_data)
        earliness_cost = model.array(earliness_cost_data)
        tardiness_cost = model.array(tardiness_cost_data)
        separation_time = model.array(separation_time_data)

        # Int variable: preferred landing time for each plane
        preferred_time_vars = [model.int(earliest_time_data[p], target_time_data[p])
                               for p in range(nb_planes)]
        preferred_time = model.array(preferred_time_vars)

        # Landing time for each plane
        landing_time_lambda = model.lambda_function(
            lambda p, prev:
                model.max(
                    preferred_time[landing_order[p]],
                    get_min_landing_time(p, prev, model, separation_time, landing_order)))
        landing_time = model.array(model.range(0, nb_planes), landing_time_lambda)

        # Landing times must respect the separation time with every previous plane
        for p in range(1, nb_planes):
            last_separation_end = [
                landing_time[previous_plane]
                + model.at(
                    separation_time,
                    landing_order[previous_plane],
                    landing_order[p])
                for previous_plane in range(p)]
            model.constraint(landing_time[p] >= model.max(last_separation_end))

        total_cost = model.sum()
        for p in range(nb_planes):
            plane_index = landing_order[p]

            # Constraint on latest landing time
            model.constraint(landing_time[p] <= latest_time[plane_index])

            # Cost for each plane
            difference_to_target_time = abs(landing_time[p] - target_time[plane_index])
            unit_cost = model.iif(
                landing_time[p] < target_time[plane_index],
                earliness_cost[plane_index],
                tardiness_cost[plane_index])
            total_cost.add_operand(unit_cost * difference_to_target_time)

        # Minimize the total cost
        model.minimize(total_cost)

        model.close()

        # Parameterize the solver
        ls.param.time_limit = time_limit

        ls.solve()

        #
        # Write the solution in a file with the following format:
        # - 1st line: value of the objective;
        # - 2nd line: for each position p, index of plane at position p.
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("%d\n" % total_cost.value)
                for p in landing_order.value:
                    f.write("%d " % p)
                f.write("\n")


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

    instance_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) >= 3 else None
    time_limit = int(sys.argv[3]) if len(sys.argv) >= 4 else 20
    main(instance_file, output_file, time_limit)
Compilation / Execution (Windows)
cl /EHsc aircraft_landing.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver130.lib
aircraft_landing instances\airland1.txt
Compilation / Execution (Linux)
g++ aircraft_landing.cpp -I/opt/localsolver_13_0/include -llocalsolver130 -lpthread -o aircraft_landing
./aircraft_landing instances/airland1.txt
#include "localsolver.h"
#include <fstream>
#include <iostream>
#include <string.h>
#include <vector>

using namespace localsolver;
using namespace std;

class AircraftLanding {
private:
    // Data from the problem
    int nbPlanes;
    int tmp;
    vector<int> earliestTimeData;
    vector<int> targetTimeData;
    vector<int> latestTimeData;
    vector<double> earlinessCostData;
    vector<double> latenessCostData;
    vector<vector<int>> separationTimeData;

    // LocalSolver
    LocalSolver localsolver;

    // Decision variables
    LSExpression landingOrder;
    vector<LSExpression> preferredTimeVars;

    // Landing time for each plane
    LSExpression landingTime;

    // Objective
    LSExpression totalCost;

public:
    /* Read instance data */
    void readInstance(const string& fileName) {
        ifstream infile;
        infile.exceptions(ifstream::failbit | ifstream::badbit);
        infile.open(fileName.c_str());

        infile >> nbPlanes;
        infile >> tmp; // Skip freezeTime value
        earliestTimeData.resize(nbPlanes);
        targetTimeData.resize(nbPlanes);
        latestTimeData.resize(nbPlanes);
        earlinessCostData.resize(nbPlanes);
        latenessCostData.resize(nbPlanes);
        separationTimeData.resize(nbPlanes, vector<int>(nbPlanes));
        preferredTimeVars.resize(nbPlanes);

        for (int i = 0; i < nbPlanes; ++i) {
            infile >> tmp; // Skip appearanceTime values
            infile >> earliestTimeData[i];
            infile >> targetTimeData[i];
            infile >> latestTimeData[i];
            infile >> earlinessCostData[i];
            infile >> latenessCostData[i];
            for (int j = 0; j < nbPlanes; ++j) {
                infile >> separationTimeData[i][j];
            }
        }
        infile.close();
    }

    LSExpression getMinLandingTime(LSExpression p, LSExpression prev, LSModel model, LSExpression separationTime) {
        return model.iif(p > 0, prev + model.at(separationTime, landingOrder[p - 1], landingOrder[p]), 0);
    }

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

        // A list variable: landingOrder[i] is the index of the ith plane to land
        landingOrder = model.listVar(nbPlanes);

        // All planes must be scheduled
        model.constraint(model.count(landingOrder) == nbPlanes);

        // Create LocalSolver arrays to be able to access them with "at" operators
        LSExpression targetTime = model.array(targetTimeData.begin(), targetTimeData.end());
        LSExpression latestTime = model.array(latestTimeData.begin(), latestTimeData.end());
        LSExpression earlinessCost = model.array(earlinessCostData.begin(), earlinessCostData.end());
        LSExpression latenessCost = model.array(latenessCostData.begin(), latenessCostData.end());
        LSExpression separationTime = model.array();
        for (int i = 0; i < nbPlanes; ++i) {
            LSExpression row = model.array(separationTimeData[i].begin(), separationTimeData[i].end());
            separationTime.addOperand(row);
        }

        // Int variables: preferred time for each plane
        for (int p = 0; p < nbPlanes; ++p) {
            preferredTimeVars[p] = model.intVar(earliestTimeData[p], targetTimeData[p]);
        }
        LSExpression preferredTime = model.array(preferredTimeVars.begin(), preferredTimeVars.end());

        // Landing time for each plane
        LSExpression landingTimeLambda = model.createLambdaFunction([&](LSExpression p, LSExpression prev) {
            return model.max(preferredTime[landingOrder[p]], getMinLandingTime(p, prev, model, separationTime));
        });
        landingTime = model.array(model.range(0, nbPlanes), landingTimeLambda);

        // Landing times must respect the separation time with every previous plane
        for (int p = 1; p < nbPlanes; ++p) {
            LSExpression lastSeparationEnd = model.max();
            for (int previousPlane = 0; previousPlane < p; ++previousPlane) {
                lastSeparationEnd.addOperand(landingTime[previousPlane] +
                                             model.at(separationTime, landingOrder[previousPlane], landingOrder[p]));
            }
            model.constraint(landingTime[p] >= lastSeparationEnd);
        }

        totalCost = model.sum();
        for (int p = 0; p < nbPlanes; ++p) {
            // Constraint on latest landing time
            LSExpression planeIndex = landingOrder[p];
            model.constraint(landingTime[p] <= latestTime[planeIndex]);

            // Cost for each plane
            LSExpression unitCost =
                model.iif(landingTime[p] < targetTime[planeIndex], earlinessCost[planeIndex], latenessCost[planeIndex]);
            LSExpression differenceToTargetTime = model.abs(landingTime[p] - targetTime[planeIndex]);
            totalCost.addOperand(unitCost * differenceToTargetTime);
        }

        // Minimize the total cost
        model.minimize(totalCost);

        model.close();

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

        localsolver.solve();
    }

    /* Write the solution in a file with the following format:
     * - 1st line: value of the objective;
     * - 2nd line: for each position p, index of plane at position p. */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());
        outfile << totalCost.getDoubleValue() << endl;
        LSCollection landingOrderCollection = landingOrder.getCollectionValue();
        for (int i = 0; i < nbPlanes; ++i) {
            outfile << landingOrderCollection[i] << " ";
        }
        outfile << endl;
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: aircraft_landing inputFile [outputFile] [timeLimit]" << endl;
        return 1;
    }
    const char* instanceFile = argv[1];
    const char* solFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "20";
    try {
        AircraftLanding model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        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 AircraftLanding.cs /reference:localsolvernet.dll
AircraftLanding instances\airland1.txt
using System;
using System.IO;
using System.Globalization;
using localsolver;

public class AircraftLanding : IDisposable
{
    // Data from the problem
    private int nbPlanes;
    private int[] earliestTimeData;
    private int[] targetTimeData;
    private int[] latestTimeData;
    private float[] earlinessCostData;
    private float[] latenessCostData;
    private int[,] separationTimeData;

    // LocalSolver
    private readonly LocalSolver localsolver;

    // Decision variables
    private LSExpression landingOrder;
    private LSExpression[] preferredTimeVars;

    // Landing time for each plane
    private LSExpression landingTime;

    // Objective
    private LSExpression totalCost;

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

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

    /* Read instance data */
    private void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] firstLineSplitted = input.ReadLine().Split();
            nbPlanes = int.Parse(firstLineSplitted[1]);

            earliestTimeData = new int[nbPlanes];
            targetTimeData = new int[nbPlanes];
            latestTimeData = new int[nbPlanes];
            earlinessCostData = new float[nbPlanes];
            latenessCostData = new float[nbPlanes];
            separationTimeData = new int[nbPlanes, nbPlanes];

            for (int p = 0; p < nbPlanes; ++p)
            {
                string[] secondLineSplitted = input.ReadLine().Split();
                earliestTimeData[p] = int.Parse(secondLineSplitted[2]);
                targetTimeData[p] = int.Parse(secondLineSplitted[3]);
                latestTimeData[p] = int.Parse(secondLineSplitted[4]);
                earlinessCostData[p] = float.Parse(
                    secondLineSplitted[5],
                    CultureInfo.InvariantCulture
                );
                latenessCostData[p] = float.Parse(
                    secondLineSplitted[6],
                    CultureInfo.InvariantCulture
                );
                int pp = 0;

                while (pp < nbPlanes)
                {
                    string[] lineSplitted = input.ReadLine().Split(' ');
                    for (int i = 0; i < lineSplitted.Length; ++i)
                    {
                        if (lineSplitted[i].Length > 0)
                        {
                            separationTimeData[p, pp] = int.Parse(lineSplitted[i]);
                            pp++;
                        }
                    }
                }
            }
        }
    }

    private LSExpression GetMinLandingTime(
        LSExpression p,
        LSExpression prev,
        LSModel model,
        LSExpression separationTime
    )
    {
        return model.If(
            p > 0,
            prev + model.At(separationTime, landingOrder[p - 1], landingOrder[p]),
            0
        );
    }

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

        // A list variable: landingOrder[i] is the index of the ith plane to land
        landingOrder = model.List(nbPlanes);

        // All planes must be scheduled
        model.Constraint(model.Count(landingOrder) == nbPlanes);

        // Create LocalSolver arrays to be able to access them with "at" operators
        LSExpression targetTime = model.Array(targetTimeData);
        LSExpression latestTime = model.Array(latestTimeData);
        LSExpression earlinessCost = model.Array(earlinessCostData);
        LSExpression latenessCost = model.Array(latenessCostData);
        LSExpression separationTime = model.Array(separationTimeData);

        // Int variables: preferred time for each plane
        preferredTimeVars = new LSExpression[nbPlanes];
        for (int p = 0; p < nbPlanes; ++p)
            preferredTimeVars[p] = model.Int(earliestTimeData[p], targetTimeData[p]);
        LSExpression preferredTime = model.Array(preferredTimeVars);

        // Landing time for each plane
        LSExpression landingTimeLambda = model.LambdaFunction(
            (p, prev) =>
                model.Max(
                    preferredTime[landingOrder[p]],
                    GetMinLandingTime(p, prev, model, separationTime)
                )
        );

        landingTime = model.Array(model.Range(0, nbPlanes), landingTimeLambda);

        // Landing times must respect the separation time with every previous plane
        for (int p = 1; p < nbPlanes; ++p)
        {
            LSExpression lastSeparationEnd = model.Max();
            for (int previousPlane = 0; previousPlane < p; ++previousPlane)
            {
                lastSeparationEnd.AddOperand(
                    landingTime[previousPlane]
                        + model.At(separationTime, landingOrder[previousPlane], landingOrder[p])
                );
            }
            model.Constraint(landingTime[p] >= lastSeparationEnd);
        }

        totalCost = model.Sum();
        for (int p = 0; p < nbPlanes; ++p)
        {
            LSExpression planeIndex = landingOrder[p];

            // Constraint on latest landing time
            model.Constraint(landingTime[p] <= latestTime[planeIndex]);

            // Cost for each plane
            LSExpression unitCost = model.If(
                landingTime[p] < targetTime[planeIndex],
                earlinessCost[planeIndex],
                latenessCost[planeIndex]
            );
            LSExpression differenceToTargetTime = model.Abs(
                landingTime[p] - targetTime[planeIndex]
            );
            totalCost.AddOperand(unitCost * differenceToTargetTime);
        }

        // Minimize the total cost
        model.Minimize(totalCost);

        model.Close();

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

        localsolver.Solve();
    }

    /* Write the solution in a file with the following format:
    * - 1st line: value of the objective;
    * - 2nd line: for each position p, index of plane at position p. */
    private void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(totalCost.GetDoubleValue());
            LSCollection landingOrderCollection = landingOrder.GetCollectionValue();
            for (int i = 0; i < nbPlanes; ++i)
                output.Write(landingOrderCollection.Get(i) + " ");
            output.WriteLine();
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: AircraftLanding 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 (AircraftLanding model = new AircraftLanding())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac AircraftLanding.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. AircraftLanding instances\airland1.txt
Compilation / Execution (Linux)
javac AircraftLanding.java -cp /opt/localsolver_13_0/bin/localsolver.jar
java -cp /opt/localsolver_13_0/bin/localsolver.jar:. AircraftLanding instances/airland1.txt
import java.util.*;
import java.io.*;
import localsolver.*;

public class AircraftLanding {

    // Data from the problem
    private int nbPlanes;
    private int[] earliestTimeData;
    private int[] targetTimeData;
    private int[] latestTimeData;
    private float[] earlinessCostData;
    private float[] latenessCostData;
    private int[][] separationTimeData;

    // LocalSolver
    private final LocalSolver localsolver;

    // Decision variables
    private LSExpression landingOrder;
    private LSExpression[] preferredTimeVars;

    // Landing time for each plane
    private LSExpression landingTime;

    // Objective
    private LSExpression totalCost;

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

    /* Read instance data */
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            nbPlanes = input.nextInt();
            input.nextInt(); // Skip freezeTime value

            earliestTimeData = new int[nbPlanes];
            targetTimeData = new int[nbPlanes];
            latestTimeData = new int[nbPlanes];
            earlinessCostData = new float[nbPlanes];
            latenessCostData = new float[nbPlanes];
            separationTimeData = new int[nbPlanes][nbPlanes];

            for (int p = 0; p < nbPlanes; ++p) {
                input.nextInt(); // Skip appearanceTime values
                earliestTimeData[p] = input.nextInt();
                targetTimeData[p] = input.nextInt();
                latestTimeData[p] = input.nextInt();
                earlinessCostData[p] = Float.parseFloat(input.next());
                latenessCostData[p] = Float.parseFloat(input.next());
                for (int pp = 0; pp < nbPlanes; ++pp) {
                    separationTimeData[p][pp] = input.nextInt();
                }
            }
        }
    }

    private LSExpression getMinLandingTime(LSExpression p, LSExpression prev, LSExpression separationTime,
        LSModel model) {
        LSExpression planeIndex = model.at(landingOrder, p);
        LSExpression previousPlaneIndex = model.at(landingOrder, model.sub(p, 1));
        return model.iif(model.gt(p, 0), model.sum(prev, model.at(separationTime, previousPlaneIndex, planeIndex)), 0);
    }

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

        // A list variable: landingOrder[i] is the index of the ith plane to land
        landingOrder = model.listVar(nbPlanes);

        // All planes must be scheduled
        model.constraint(model.eq(model.count(landingOrder), nbPlanes));

        // Create LocalSolver arrays to be able to access them with "at" operators
        LSExpression targetTime = model.array(targetTimeData);
        LSExpression latestTime = model.array(latestTimeData);
        LSExpression earlinessCost = model.array(earlinessCostData);
        LSExpression latenessCost = model.array(latenessCostData);
        LSExpression separationTime = model.array(separationTimeData);

        // Int variables: preferred time for each plane
        preferredTimeVars = new LSExpression[nbPlanes];
        for (int p = 0; p < nbPlanes; ++p) {
            preferredTimeVars[p] = model.intVar(earliestTimeData[p], targetTimeData[p]);
        }
        LSExpression preferredTime = model.array(preferredTimeVars);

        // Landing time for each plane
        LSExpression landingTimeLambda = model
            .lambdaFunction((p, prev) -> model.max(model.at(preferredTime, model.at(landingOrder, p)),
                getMinLandingTime(p, prev, separationTime, model)));
        landingTime = model.array(model.range(0, nbPlanes), landingTimeLambda);

        // Landing times must respect the separation time with every previous plane
        for (int p = 1; p < nbPlanes; ++p) {
            LSExpression lastSeparationEnd = model.max();
            for (int previousPlane = 0; previousPlane < p; ++previousPlane) {
                lastSeparationEnd.addOperand(model.sum(model.at(landingTime, previousPlane),
                    model.at(separationTime, model.at(landingOrder, previousPlane), model.at(landingOrder, p))));
            }
            model.constraint(model.geq(model.at(landingTime, p), lastSeparationEnd));
        }

        totalCost = model.sum();
        for (int p = 0; p < nbPlanes; ++p) {
            LSExpression planeIndex = model.at(landingOrder, p);

            // Constraint on latest landing time
            model.addConstraint(model.leq(model.at(landingTime, p), model.at(latestTime, planeIndex)));

            // Cost for each plane
            LSExpression unitCost = model.iif(model.lt(model.at(landingTime, p), model.at(targetTime, planeIndex)),
                model.at(earlinessCost, planeIndex), model.at(latenessCost, planeIndex));
            LSExpression differenceToTargetTime = model
                .abs(model.sub(model.at(landingTime, p), model.at(targetTime, planeIndex)));
            totalCost.addOperand(model.prod(unitCost, differenceToTargetTime));
        }

        // Minimize the total cost
        model.minimize(totalCost);

        model.close();

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

        localsolver.solve();
    }

    /*
     * Write the solution in a file with the following format:
     * - 1st line: value of the objective;
     * - 2nd line: for each position p, index of plane at position p.
     */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(new FileWriter(fileName))) {
            output.println(totalCost.getDoubleValue());
            LSCollection landingOrderCollection = landingOrder.getCollectionValue();
            for (int i = 0; i < nbPlanes; ++i) {
                output.print(landingOrderCollection.get(i) + " ");
            }
            output.println();
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: AircraftLanding inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }

        String instanceFile = args[0];
        String outputFile = args.length > 1 ? args[1] : null;
        String strTimeLimit = args.length > 2 ? args[2] : "20";
        try (LocalSolver localsolver = new LocalSolver()) {
            AircraftLanding model = new AircraftLanding(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);
        }
    }
}