This page is for an old version of Hexaly Optimizer. We recommend that you update your version and read the documentation for the latest stable release.

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.png

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 LocalSolver 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=]
/********** aircraft_landing.lsp **********/

use io;
use ls;

/* 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-1] {
        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][0..nbPlanes-1] = 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-1] <- int(earliestTime[p], targetTime[p]);

    // Landing time for each plane
    landingTime <- array(0..nbPlanes-1, (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-1] {
        constraint landingTime[p] >= max[previousPlane in 0..p-1](landingTime[previousPlane] + separationTime[landingOrder[previousPlane]][landingOrder[p]]);
    }

    // Cost for each plane
    for [p in 0..nbPlanes-1] {
        local planeIndex <- landingOrder[p];
        constraint landingTime[p] <= latestTime[planeIndex];
        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-1] (cost[p]);
    minimize totalCost;
}

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

/* Write the solution in a file following 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_10_5/bin/python
python aircraft_landing.py instances/airland1.txt
########## aircraft_landing.py ##########

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 = []
    target_time = []
    latest_time = []
    earliness_cost = []
    tardiness_cost = []
    separation_time = []

    for p in range(nb_planes):
        next(file_it) # Skip appearanceTime values
        earliest_time.append(int(next(file_it)))
        target_time.append(int(next(file_it)))
        latest_time.append(int(next(file_it)))
        earliness_cost.append(float(next(file_it)))
        tardiness_cost.append(float(next(file_it)))
        separation_time.append([None] * nb_planes)

        for pp in range(nb_planes):
            separation_time[p][pp] = int(next(file_it))
    return (nb_planes, earliest_time, target_time, latest_time, earliness_cost, tardiness_cost, separation_time)

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

def main(instance_file, output_file, time_limit):
    nb_planes, earliest_time, target_time, latest_time, earliness_cost, tardiness_cost, separation_time = 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_array = model.array(target_time)
        latest_time_array = model.array(latest_time)
        earliness_cost_array = model.array(earliness_cost)
        tardiness_cost_array = model.array(tardiness_cost)
        separation_time_array = model.array(separation_time)

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

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

        # 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_array, 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_array[plane_index])

            # Cost for each plane
            difference_to_target_time = abs(landing_time[p] - target_time_array[plane_index])
            unit_cost = earliness_cost_array[plane_index] if landing_time[p] < target_time_array[plane_index] else tardiness_cost_array[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
        #
        if len(sys.argv) >= 4:
            ls.param.time_limit = int(sys.argv[3])
        else:
            ls.param.time_limit = 100

        ls.solve()

        # Write the solution in a file following the following format:
        # - 1st line: value of the objective;
        # - 2nd line: for each position p, index of plane at position p.
        if len(sys.argv) >= 3:
            with open(sys.argv[2], '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 60
    main(instance_file, output_file, time_limit)
Compilation / Execution (Windows)
cl /EHsc aircraft_landing.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver105.lib
aircraft_landing instances\airland1.txt
Compilation / Execution (Linux)
g++ aircraft_landing.cpp -I/opt/localsolver_10_5/include -llocalsolver105 -lpthread -o aircraft_landing
./aircraft_landing instances/airland1.txt
/********** aircraft_landing.cpp **********/

#include <iostream>
#include <fstream>
#include <vector>
#include <string.h>
#include "localsolver.h"

using namespace localsolver;
using namespace std;

class AircraftLanding {
private:

    // Data from the problem
    int nbPlanes;
    int freezeTime;
    vector<lsint> appearanceTime;
    vector<lsint> earliestTime;
    vector<lsint> targetTime;
    vector<lsint> latestTime;
    vector<lsdouble> earlinessCost;
    vector<lsdouble> latenessCost;
    vector<vector<int>> separationTime;

    // LocalSolver
    LocalSolver localsolver;

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

    // 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 >> freezeTime;
        appearanceTime.resize(nbPlanes);
        earliestTime.resize(nbPlanes);
        targetTime.resize(nbPlanes);
        latestTime.resize(nbPlanes);
        earlinessCost.resize(nbPlanes);
        latenessCost.resize(nbPlanes);
        separationTime.resize(nbPlanes, vector<int>(nbPlanes));
        preferredTime.resize(nbPlanes);

        for (int i = 0; i < nbPlanes; ++i) {
            infile >> appearanceTime[i];
            infile >> earliestTime[i];
            infile >> targetTime[i];
            infile >> latestTime[i];
            infile >> earlinessCost[i];
            infile >> latenessCost[i];
            for (int j = 0; j < nbPlanes; ++j) {
                infile >> separationTime[i][j];
            }
        }
        infile.close();
    }

    LSExpression getMinLandingTime(LSExpression p, LSExpression prev, LSModel model, LSExpression separationTimeArray) {
        return model.iif(p > 0,
                         prev + model.at(separationTimeArray, 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 an "at" operator
        LSExpression targetTimeArray = model.array(targetTime.begin(), targetTime.end());
        LSExpression latestTimeArray = model.array(latestTime.begin(), latestTime.end());
        LSExpression earlinessCostArray = model.array(earlinessCost.begin(), earlinessCost.end());
        LSExpression latenessCostArray = model.array(latenessCost.begin(), latenessCost.end());
        LSExpression separationTimeArray = model.array();
        for (int i = 0; i < nbPlanes; i++) {
            LSExpression row = model.array(separationTime[i].begin(), separationTime[i].end());
            separationTimeArray.addOperand(row);
        }

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

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

        // 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(separationTimeArray, 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] <= latestTimeArray[planeIndex]);

            // Cost for each plane
            LSExpression unitCost = model.iif(landingTime[p] < targetTimeArray[planeIndex],
                                              earlinessCostArray[planeIndex],
                                              latenessCostArray[planeIndex]);
            LSExpression differenceToTargetTime = model.abs(landingTime[p] - targetTimeArray[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 */
    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
/********** AircraftLanding.cs **********/

using System;
using System.IO;
using System.Globalization;
using localsolver;

public class AircraftLanding : IDisposable
{
    // Data from the problem
    private int nbPlanes;
    private int[] earliestTime;
    private int[] targetTime;
    private int[] latestTime;
    private float[] earlinessCost;
    private float[] latenessCost;
    private int[,] separationTime;

    // LocalSolver
    private readonly LocalSolver localsolver;

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

    // 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]);

            earliestTime = new int[nbPlanes];
            targetTime = new int[nbPlanes];
            latestTime = new int[nbPlanes];
            earlinessCost = new float[nbPlanes];
            latenessCost = new float[nbPlanes];
            separationTime = new int[nbPlanes, nbPlanes];

            for (int p = 0; p < nbPlanes; ++p)
            {
                string[] secondLineSplitted = input.ReadLine().Split();
                earliestTime[p] = int.Parse(secondLineSplitted[2]);
                targetTime[p] = int.Parse(secondLineSplitted[3]);
                latestTime[p] = int.Parse(secondLineSplitted[4]);
                earlinessCost[p] = float.Parse(secondLineSplitted[5], CultureInfo.InvariantCulture);
                latenessCost[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)
                        {
                            separationTime[p, pp] = int.Parse(lineSplitted[i]);
                            pp++;
                        }
                    }
                }
            }
        }
    }

    private LSExpression GetMinLandingTime(LSExpression p, LSExpression prev, LSModel model, LSExpression separationTimeArray)
    {
        return model.If(p>0,
                        prev + model.At(separationTimeArray, 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 in order to be able to access them with "at" operators
        LSExpression targetTimeArray = model.Array(targetTime);
        LSExpression latestTimeArray = model.Array(latestTime);
        LSExpression earlinessCostArray = model.Array(earlinessCost);
        LSExpression latenessCostArray = model.Array(latenessCost);
        LSExpression separationTimeArray = model.Array(separationTime);

        // Int variables: preferred time for each plane
        preferredTime = new LSExpression[nbPlanes];
        for (int p = 0; p < nbPlanes; ++p)
        {
            preferredTime[p] = model.Int(earliestTime[p], targetTime[p]);
        }
        LSExpression preferredTimeArray = model.Array(preferredTime);

        // Landing time for each plane
        LSExpression landingTimeSelector = model.LambdaFunction((p, prev) =>
                        model.Max(preferredTimeArray[landingOrder[p]], GetMinLandingTime(p, prev, model, separationTimeArray)));

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

        // 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(separationTimeArray, 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] <= latestTimeArray[planeIndex]);

            // Cost for each plane
            LSExpression unitCost = model.If(landingTime[p] < targetTimeArray[planeIndex],
                                             earlinessCostArray[planeIndex],
                                             latenessCostArray[planeIndex]);
            LSExpression differenceToTargetTime = model.Abs(landingTime[p] - targetTimeArray[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 */
    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_10_5/bin/localsolver.jar
java -cp /opt/localsolver_10_5/bin/localsolver.jar:. AircraftLanding instances/airland1.txt
/********** AircraftLanding.java **********/

import java.util.*;
import java.io.*;
import localsolver.*;

public class AircraftLanding {

    // Data from the problem
    private int nbPlanes;
    private int[] earliestTime;
    private int[] targetTime;
    private int[] latestTime;
    private float[] earlinessCost;
    private float[] latenessCost;
    private int[][] separationTime;

    // LocalSolver
    private final LocalSolver localsolver;

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

    // 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

            earliestTime = new int[nbPlanes];
            targetTime = new int[nbPlanes];
            latestTime = new int[nbPlanes];
            earlinessCost = new float[nbPlanes];
            latenessCost = new float[nbPlanes];
            separationTime = new int[nbPlanes][nbPlanes];

            for (int p = 0; p < nbPlanes; p++) {
                input.nextInt(); // Skip appearanceTime values
                earliestTime[p] = input.nextInt();
                targetTime[p] = input.nextInt();
                latestTime[p] = input.nextInt();
                earlinessCost[p] = Float.parseFloat(input.next());
                latenessCost[p] = Float.parseFloat(input.next());
                for (int pp = 0; pp < nbPlanes; pp++) {
                    separationTime[p][pp] = input.nextInt();
                }
            }
        }
    }

    private LSExpression getMinLandingTime(LSExpression p, LSExpression prev, LSExpression separationTimeArray, 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(separationTimeArray, 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 in order to be able to access them with "at" operators
        LSExpression targetTimeArray = model.array(targetTime);
        LSExpression latestTimeArray = model.array(latestTime);
        LSExpression earlinessCostArray = model.array(earlinessCost);
        LSExpression latenessCostArray = model.array(latenessCost);
        LSExpression separationTimeArray = model.array(separationTime);

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

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

        // 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(separationTimeArray, 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(latestTimeArray, planeIndex)));

            // Cost for each plane
            LSExpression unitCost = model.iif(model.lt(model.at(landingTime, p), model.at(targetTimeArray, planeIndex)),
                                              model.at(earlinessCostArray, planeIndex),
                                              model.at(latenessCostArray, planeIndex));
            LSExpression differenceToTargetTime = model.abs(model.sub(model.at(landingTime, p), model.at(targetTimeArray, 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 */
    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);
        }
    }
}