Wind Farm Layout Optimization Problem (WFLOP)

import hexaly.optimizer
import json
import math
import sys

def main(instance_file, output_file, time_limit):
    """
    Creates and solves the model.
    """
    # Read the input data
    degrees, nb_degrees, probabilities, inflow_speed, thrust_coeff,  dyn_coeff, \
            radius, nb_turbines, turbine_diameter, turbine_nominal_power, min_distance, \
            optimal_speed, min_speed, max_speed, disc_delta = read_data(instance_file)
    ## Additional required inputs
    points_x, points_y = build_discretization(radius, disc_delta)
    nb_locations = len(points_x)
    incompatibility_sets = compute_incompatibilities(points_x, points_y, min_distance)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        # Declare the optimization model
        model = optimizer.model

        # Decision Variable : List of booleans (one for each Wind Turbine)
        chosen_turbines = [model.bool() for location in range(nb_locations)]

        # Constraint : We have to choose nT locations on the nb_locations available positions
        model.constraint(sum(chosen_turbines[location] for location in range(nb_locations)) == nb_turbines)
        # Constraint : Two WT need to be spaced by at least d_min meters
        for location in range(nb_locations):
            for incompatible_location in incompatibility_sets[location]:
                model.constraint(chosen_turbines[location] + chosen_turbines[incompatible_location] <= 1)

        # Relative wind and power at each location for every possible wind angle
        turbine_wind = [[0 for angle in range(nb_degrees)] for location in range(nb_locations)]
        turbine_power = [[0 for angle in range(nb_degrees)] for location in range(nb_locations)]
        for location_1 in range(nb_locations):
            wt_1 = [points_x[location_1], points_y[location_1]]
            for angle in range(nb_degrees):
                # Wind at location i for a wind incidence angle j
                w_loss_sum = sum(chosen_turbines[location_2]
                        * (wake_loss(wt_1, [points_x[location_2], points_y[location_2]],
                                degrees[angle], dyn_coeff, turbine_diameter, thrust_coeff))**2
                                for location_2 in range(nb_locations))
                turbine_wind[location_1][angle] =  inflow_speed * (chosen_turbines[location_1]
                        - model.sqrt(w_loss_sum))

                # Power of the wind turbine of a WT at i with a wind incidence angle j
                polynomial_case = (min_speed <= turbine_wind[location_1][angle]) \
                        * (turbine_wind[location_1][angle] < optimal_speed)
                polynomial_value = ((turbine_wind[location_1][angle] - min_speed) \
                        / (optimal_speed - min_speed))**3
                constant_case = (optimal_speed <= turbine_wind[location_1][angle]) \
                        * (turbine_wind[location_1][angle] < max_speed)
                turbine_power[location_1][angle] = polynomial_case \
                        * turbine_nominal_power * polynomial_value \
                        +  constant_case * turbine_nominal_power


        # Objective : we maximize the Annual Energy Production (in Wh)
        # To calculate the AEP, we multiply the total power produced by the number of hours in one year
        objective = 8760 * sum(probabilities[angle] * turbine_power[location][angle]
                for location in range(nb_locations) for angle in range(nb_degrees))
        model.maximize(objective)

        # Finalize model and solve
        model.close()
        optimizer.param.time_limit = time_limit
        optimizer.solve()

        feasible_list = ["HxSolutionStatus.FEASIBLE", "HxSolutionStatus.OPTIMAL"]
        if (str(optimizer.solution.status) in feasible_list):
            # Store the solution in a map
            points = []
            for location in range(nb_locations):
                if (chosen_turbines[location].value == 1):
                    points.append((points_x[location], points_y[location]))

            # If asked, we print the result in the output file
            if (output_file is not None):
                with open(output_file, "w") as f:
                    print("Solution written in file ", output_file)
                    # Print every important info of the resolution
                    f.write("*************************************\n")
                    f.write("*** WIND FARM LAYOUT OPTIMIZATION ***\n")
                    f.write("*************************************\n\n")

                    f.write("Annual Energy Production : " + str(objective.value)
                            + " Wh (" + str(round(objective.value / 1e9 * 100) / 100)
                            + " GWh).\n\n")
                    f.write("Locations of the Wind Turbines : \n")

                    for turbine in range(nb_turbines):
                        point = points[turbine]
                        f.write("   - (" + str(point[0]) + ", " + str(point[1]) + ")\n")

        else:
            print("No feasible solution have been found within the allowed time.")

def read_data(instance_file):
    """
    Reads the input files of the problem.
    """
    # Opening the .json file
    with open(instance_file, 'r') as file:
        data = json.load(file)

    # Wind information
    degrees = data["wind"]["degrees"]     # Discretization of the degrees rose
    nb_degrees = len(degrees)             # Number of discretizations
    probabilities = data["wind"]["probs"] # Probabilities of wind in every direction
    inflow_speed = data["wind"]["speed"]  # Inflow wind speed
    thrust_coeff = data["wind"]["cT"]     # Thrust coefficient 
    dyn_coeff = data["wind"]["kY"]        # Dynamic coefficient

    # Information about the problem
    radius = data["radius"]             # Radius of the perimeter to fullfil
    nb_turbines = data["nT"]            # Number of Wind Turbines
    turbine_diameter = data["D"]                   # Diameter of each Turbine
    turbine_nominal_power = data["p_rated"] * 1e6  # Power of the Turbine
    min_distance = 2 * turbine_diameter            # Minimal distance between two turbines
    optimal_speed = data["u_rated"]     # Optimal wind speed (maximum efficiency)
    min_speed = data["u_cut_in"]        # Minimal wind speed producing energy
    max_speed = data["u_cut_out"]       # Maximal wind speed producing energy

    disc_delta = 1.7 * turbine_diameter # Distance between two points in the discretization

    return degrees, nb_degrees, probabilities, inflow_speed, thrust_coeff, \
            dyn_coeff, radius, nb_turbines, turbine_diameter, turbine_nominal_power, \
            min_distance, optimal_speed, min_speed, max_speed, disc_delta

def build_discretization(radius, disc_delta):
    """
    Builds the discretization in a circular field of radius {radius} and with a
    distance between points of {disc_delta}.

    The discretization follows the one described in 'eawe' paper, i.e. building
    a regular uniform mesh inside the circle, and a point every degree on the
    frontier of the circle.
    """
    points_x = []
    points_y = []
    # Maximum number of points in any direction from the center
    max_onedir_points = int(radius / disc_delta + 1)
    # Interior points
    for c_x in range(max_onedir_points):
        for side_x in [-1, 1]:
            if (not(c_x == 0 and side_x == 1)):
                new_x = side_x * c_x * disc_delta
                for c_y in range(max_onedir_points):
                    for side_y in [-1, 1]:
                        new_y = side_y * c_y * disc_delta
                        if (math.sqrt(new_x**2 + new_y**2) < radius \
                                and not(c_y == 0 and side_y == 1)):
                            points_x.append(new_x)
                            points_y.append(new_y)

    # Points on the border of the circle
    for deg in range(360):
        points_x.append(radius * math.sin(deg * math.pi / 180))
        points_y.append(radius * math.cos(deg * math.pi / 180))

    return points_x, points_y

def compute_incompatibilities(points_x, points_y, min_distance):
    """
    Compute incompatibility sets for every location on the field.

    N_i[i] := {l in 0...nb_locations : l != i && ||l - i|| < min_distance}
    """
    # Number of available points
    nb_locations = len(points_x)
    incompatibilities = [[] for location in range(nb_locations)]
    # Fill the set for every location
    for location_1 in range(nb_locations):
        for location_2 in range(location_1+1, nb_locations):
            if (math.sqrt((points_x[location_1] - points_x[location_2])**2
                    + (points_y[location_1] - points_y[location_2])**2) < min_distance):
                incompatibilities[location_1].append(location_2)
                incompatibilities[location_2].append(location_1)

    return incompatibilities

def distances(location_1, location_2, theta):
    """
    This function calculates the distance between two points, considering a
    frame of reference with angle theta.
    """
    # Rotation of the reference
    theta_deg = 270 - theta
    theta_rad = theta_deg * math.pi / 180
    cos_wind = math.cos(-theta_rad)
    sin_wind = math.sin(-theta_rad)

    # Change the reference of the coordinates
    location_2_x = (location_2[0] * cos_wind) - (location_2[1] * sin_wind)
    location_2_y = (location_2[0] * sin_wind) + (location_2[1] * cos_wind)

    location_1_x = (location_1[0] * cos_wind) - (location_1[1] * sin_wind)
    location_1_y = (location_1[0] * sin_wind) + (location_1[1] * cos_wind)

    return {"parallel":location_1_x - location_2_x, "perpendicular":location_1_y - location_2_y}

def wake_loss(turbine_1, turbine_2, theta, dyn_coeff, turbine_diameter, thrust_coeff):
    """
    Wake loss evaluated at i, caused by a Wind Turbine at l, considering a wind
    of angle {theta}. 
    This loss is a simplified version of Bastankhah's Gaussian model.
    """
    # We calculate parallel and perpendicular distances
    d = distances(turbine_1, turbine_2, theta)
    # If x_i - x_l < 0, the wake loss is zero
    if (d["parallel"] > 0):
        # Standard deviation of the wake deficit
        s_y = dyn_coeff * d["parallel"] + turbine_diameter / math.sqrt(8)
        # Separation of the product in two components
        coeff = thrust_coeff / (8 * (s_y / turbine_diameter)**2)
        pdt_1 = 1 - math.sqrt(1 - coeff)
        pdt_2 = math.exp(-1/2 * (d["perpendicular"] / s_y)**2)

        return pdt_1 * pdt_2

    return 0

# Main
if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python wind_farm_layout_optimization.py inputFile "
              + "[solFile] [timeLimit]")
        sys.exit(1)
    # Get arguments
    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
    # Create the model and solve it
    main(instance_file, output_file, time_limit)

#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <fstream>
#include <iostream>
#include <numeric>
#include <vector>
#include "libs/lib-cpp/json.hpp"
#include <map>
#include <sstream>
#include <vector>
#include <string>
#include <any>

using namespace hexaly;
using namespace std;
using json = nlohmann::json;

class WindFarmLayoutOptimization {
public:

    // Wind information 
    vector<double> degrees;  // Discretization of the degrees rose
    int nbDegrees;                // Number of discretizations
    vector<double> probabilities; // Probabilities of wind in every direction
    double inflowSpeed;           // Inflow wind speed
    double thrustCoeff;      // Thrust coefficient 
    double dynCoeff;         // Dynamic coefficient

    // Information about the problem
    double radius;           // Radius of the perimeter to fullfil
    int nbTurbines;          // Number of Wind Turbines
    double turbineDiameter;     // Diameter of each Turbine
    double turbineNominalPower; // Power of the Turbine
    double minDistance;         // Minimal distance between two turbines
    double optimalSpeed;     // Optimal wind speed (maximum efficiency)
    double minSpeed;         // Minimal wind speed producing energy
    double maxSpeed;         // Maximal wind speed producing energy 

    // Discretization data
    double discDelta;                 // Distance between two points 
    vector<double> pointsX;           // Available locations (X-coord)
    vector<double> pointsY;           // Available locations (Y-coord)
    int nbLocations;                  // Number of available locations
    map<int, vector<int>> incompatibilitySets; // Incompatibility sets

    // Hexaly Optimizer
    HexalyOptimizer optimizer;
    // Decision variable
    std::vector<HxExpression> chosenTurbines;

    // Wind at each location
    std::vector<std::vector<HxExpression>> turbineWind;
    // Power of each wind turbine
    std::vector<std::vector<HxExpression>> turbinePower;

    // Objective function
    HxExpression objective;

    double PI = 3.141592653589793;

    /* Create and solve the model */
    void solve(int timelimit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Decision Variable : List of booleans (one for each Wind Turbine)
        chosenTurbines.resize(nbLocations);
        for (int location = 0; location < nbLocations; ++location) {
            chosenTurbines[location] = model.boolVar();
        }

        // Constraint : We have to choose nbTurbines locations on the nbLocations available positions
        HxExpression chosenTurbinesSum = model.sum();
        for (int location = 0; location < nbLocations; location++) {
            chosenTurbinesSum.addOperand(chosenTurbines[location]);
        }
        model.constraint(chosenTurbinesSum == nbTurbines);

        // Constraint : Two WT need to be spaced by at least dMin meters
        for (int location = 0; location < nbLocations; location++) {
            for (int incompatibleLocation : incompatibilitySets[location]) {
                model.constraint(chosenTurbines[location] + chosenTurbines[incompatibleLocation] <= 1);
            }
        }

        // Wind and power at each location for every possible wind angle
        turbineWind.resize(nbLocations);
        turbinePower.resize(nbLocations);
        for (int location = 0; location < nbLocations; ++location) {
            turbineWind[location].resize(nbDegrees);
            turbinePower[location].resize(nbDegrees);
        }
        for (int location1 = 0; location1 < nbLocations; location1++) {
            vector<double> wt1;
            wt1.push_back(pointsX[location1]);
            wt1.push_back(pointsY[location1]);

            for (int angle = 0; angle < nbDegrees; angle++) {
                HxExpression sumWakeLossesSqd = model.sum();

                for (int location2 = 0; location2 < nbLocations; location2++) {
                    vector<double> wt2;
                    wt2.push_back(pointsX[location2]);
                    wt2.push_back(pointsY[location2]);

                    double wakeL = wakeLoss(wt1, wt2, degrees[angle]);

                    sumWakeLossesSqd.addOperand(
                            model.pow(model.prod(wakeL, chosenTurbines[location2]), 2));
                }

                // Wind at location i for a wind incidence angle j
                turbineWind[location1][angle] = model.prod(inflowSpeed, model.sub(
                        chosenTurbines[location1], model.sqrt(sumWakeLossesSqd)));
                // Power of the wind turbine of a WT at i with a wind incidence angle j
                HxExpression polynomialCase = model.prod(
                        model.leq(minSpeed, turbineWind[location1][angle]),
                        model.lt(turbineWind[location1][angle], optimalSpeed));
                HxExpression polynomialValue = model.pow(
                        model.div(model.sub(turbineWind[location1][angle], minSpeed),
                                optimalSpeed - minSpeed), 3);
                HxExpression constantCase = model.prod(
                        model.leq(optimalSpeed, turbineWind[location1][angle]),
                        model.lt(turbineWind[location1][angle], maxSpeed));
                turbinePower[location1][angle] = model.sum(
                        model.prod(turbineNominalPower, polynomialCase, polynomialValue),
                        model.prod(turbineNominalPower, constantCase));
            }
        }

        // Objective : we maximize the Annual Energy Production (in Wh)
        objective = model.sum();
        for (int location = 0; location < nbLocations; location++) {
            for (int angle = 0; angle < nbDegrees; angle++) {
                // To calculate the AEP, we multiply the total power produced 
                // by the number of hours in one year
                objective.addOperand(8760 * probabilities[angle] * turbinePower[location][angle]);
            }
        }
        model.maximize(objective);

        // Finalize model and solve
        model.close();
        optimizer.getParam().setTimeLimit(timelimit);
        optimizer.solve();
    }

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

        json data;
        infile >> data;
        infile.close();

        // Wind information
        nbDegrees = data["wind"]["J"];
        inflowSpeed = data["wind"]["speed"];
        thrustCoeff = data["wind"]["cT"];
        dynCoeff = data["wind"]["kY"];

        vector<double> djson = data["wind"]["degrees"];
        vector<double> pjson = data["wind"]["probs"];

        for (int i = 0; i < djson.size(); i++) {
            degrees.push_back(djson[i]);
            probabilities.push_back(pjson[i]);
        }

        // Information about the problem
        radius = data["radius"];
        nbTurbines = data["nT"];
        turbineDiameter = data["D"];
        turbineNominalPower= (double)data["p_rated"] * 1e6;
        minDistance = 2.0 * turbineDiameter;
        optimalSpeed = data["u_rated"];
        minSpeed = data["u_cut_in"];
        maxSpeed = data["u_cut_out"];

        discDelta = turbineDiameter * 1.7;

        // Additional required inputs
        buildDiscretization();
        computeIncompatibilities();
    }

    /*
     *  Build the discretization in a circular field of radius {radius} and
     *  with a distance between points of {discDelta}.
     *
     *  The discretization follows the one described in 'eawe' paper, i.e.
     *  building a regular uniform mesh inside the circle, and a point every
     *  degree on the frontier of the circle.
     */
    void buildDiscretization() {
        // Maximum number of points in any direction from the center
        int maxOnedirPoints = round(radius / discDelta) + 1;
        vector<int> sides = {-1, 1};
        // Interior points
        for (int cX = 0; cX < maxOnedirPoints; cX++) {
            for (int sideX : sides) {
                if (!(cX == 0 && sideX == 1)) {
                    double newX = sideX * cX * discDelta;
                    for (int cY = 0; cY < maxOnedirPoints; cY++) {
                        for (int sideY : sides) {
                            double newY = sideY * cY * discDelta;
                            if (sqrt(pow(newX, 2) + pow(newY, 2)) < radius
                                    && !(cY == 0 && sideY == 1)) {
                                pointsX.push_back(newX);
                                pointsY.push_back(newY);
                            }
                        }
                    } 
                }
            }
        }
        // Points on the border of the circle
        for (int deg = 0; deg < 360; deg++) {
            double radians = deg * PI / 180;
            double x = radius * sin(radians);
            double y = radius * cos(radians);
            pointsX.push_back(x);
            pointsY.push_back(y);
        }
    }

    /*
     * Compute incompatibility sets for every location on the field.
     *
     * N_i[i] := {l in 0...nbTurbines : l != i && ||l - i|| < minDistance}
     */
    void computeIncompatibilities() {
        nbLocations = pointsX.size();

        for (int location1 = 0; location1 < nbLocations; location1++) {
            for (int location2 = location1+1; location2 < nbLocations; location2++) {
                if (sqrt(pow(pointsX[location1] - pointsX[location2], 2) 
                        + pow(pointsY[location1] - pointsY[location2], 2)) < minDistance) {
                    incompatibilitySets[location1].push_back(location2);
                    incompatibilitySets[location2].push_back(location1);
                }
            }
        }
    }

    /*
     * Calculate the distance between two points, considering a frame of
     * reference with angle theta.
     */
    map<string, double> distances(vector<double> location1,
            vector<double> location2,
            double theta) {
        map<string, double> res;
        // Rotation of the reference
        double thetaDeg = 270 - theta;
        double thetaRad = thetaDeg * PI / 180;

        double cosWind = cos(-thetaRad);
        double sinWind = sin(-thetaRad);

        // Change the reference of the coordinates
        double location2X = (location2[0] * cosWind) - (location2[1] * sinWind);
        double location2Y = (location2[0] * sinWind) + (location2[1] * cosWind);

        double location1X = (location1[0] * cosWind) - (location1[1] * sinWind);
        double location1Y = (location1[0] * sinWind) + (location1[1] * cosWind);

        res["parallel"] = location1X - location2X;
        res["perpendicular"] =  location1Y - location2Y;

        return res;

    }

    /*
     * Wake loss evaluated at i, caused by a Wind Turbine at l, considering
     * a wind of angle {theta}.
     * This loss is a simplified version of Bastankhah's Gaussian model.
     */
    double wakeLoss(vector<double> turbine1,
            vector<double> turbine2,
            double theta) {
        // We calculate parallel and perpendicular distances
        map<string, double> d = distances(turbine1, turbine2, theta);
        // If x_i - x_l < 0, the wake loss is zero
        if (d["parallel"] > 0) {
            // Standard deviation of the wake deficit
            double sY = dynCoeff * d["parallel"] + turbineDiameter / sqrt(8);
            // Separation of the product in two components
            double coeff =  thrustCoeff / (8.0 * pow(sY / turbineDiameter, 2));
            double pdt1 = 1.0 - sqrt(1.0 - coeff);
            double pdt2 = exp(-0.5 * pow(d["perpendicular"] / sY, 2));

            return pdt1 * pdt2;
        }
        return 0;
    }
 
    /* 
     * Write the solution in the given file, giving :
     *    - The instance
     *    - The timelimit
     *    - The Annual Energy Production (objective function, in Wh and GWh)
     *    - Location (x,y) of every turbine
     */
    void writeSolution(const string& fileName) {
        if (optimizer.getSolution().getStatus() == 2
                || optimizer.getSolution().getStatus() == 3) {

            vector<vector<double>> points;

            for (int i = 0; i < nbLocations; i++) {
                if (chosenTurbines[i].getValue() == 1) {
                    vector<double> point = {pointsX[i], pointsY[i]};
                    points.push_back(point);
                }
            }

            // If asked, we print the result in the output file
            ofstream outfile;
            outfile.exceptions(ofstream::failbit | ofstream::badbit);
            outfile.open(fileName.c_str());

            printf("Solution written in file %s\n\n", fileName.c_str());

            outfile << "*************************************\n";
            outfile << "*** WIND FARM LAYOUT OPTIMIZATION ***\n";
            outfile << "*************************************\n\n";

            outfile << "Annual Energy Production : " << 
                    objective.getDoubleValue() << " Wh"; 
            outfile << " (" << round(objective.getDoubleValue() / 1e9 * 100) / 100 <<
                    " GWh).\n\n";
            outfile << "Locations of the Wind Turbines : \n";

            for (int turbine = 0; turbine < nbTurbines; turbine++) {
                outfile << "   - (" << points[turbine][0] << ", " << points[turbine][1] << ")\n";
            }
        } else {
            printf("No feasible solution have been found within the allowed time.\n");
        }
        
    }
};

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

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

    try {
        WindFarmLayoutOptimization model;

        model.readData(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;
    }
}

using System;
using System.IO;
using System.Linq;
using System.Collections.Generic;
using Hexaly.Optimizer;

using Newtonsoft.Json.Linq;
using Newtonsoft.Json;

public class WindFarmLayoutOptimization : IDisposable
{
    // Wind information
    private double[] degrees;       // Discretization of the degrees rose
    private int nbDegrees;          // Number of discretizations
    private double[] probabilities; // Probabilities of wind in every direction
    private double inflowSpeed;     // Inflow wind speed
    private double thrustCoeff;     // Thrust coefficient
    private double dynCoeff;        // Dynamic coefficient

    // Information about the problem
    private double radius;          // Radius of the perimeter to fullfil
    private int nbTurbines;         // Number of Wind Turbines
    private double turbineDiameter;     // Diameter of each Turbine
    private double turbineNominalPower; // Power of the Turbine
    private double minDistance;         // Minimal distance between two turbines
    private double optimalSpeed;    // Optimal wind speed (maximum efficiency)
    private double minSpeed;        // Minimal wind speed producing energy
    private double maxSpeed;        // Maximal wind speed producing energy

    // Discretization data
    private double discDelta;       // Distance between two points
    List<double> pointsX;           // Available locations (X-coord)
    List<double> pointsY;           // Available locations (Y-coord)
    int nbLocations;                // Number of available locations
    Dictionary<int, List<int>> incompatibilitySets; // Incompatibility sets

    // Hexaly Optimizer
    HexalyOptimizer optimizer;
    // Decision variable
    HxExpression[] chosenTurbines;

    // Wind at each location
    HxExpression[,] turbineWind;
    // Power of each wind turbine
    HxExpression[,] turbinePower;

    // Objective function
    HxExpression objective;

    public WindFarmLayoutOptimization()
    {
        optimizer = new HexalyOptimizer();
    }

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

    /* Create and solve the model */
    private void Solve(int timelimit)
    {
        // Declare the optimization model
        HxModel model = optimizer.GetModel();

        // Decision Variable : List of booleans (one for each Wind Turbine)
        chosenTurbines = new HxExpression[nbLocations];
        for (int location = 0; location < nbLocations; ++location)
        {
            chosenTurbines[location] = model.Bool();
        }

        // Constraint : We have to choose nbTurbines locations on the nbLocations available positions
        model.Constraint(model.Eq(model.Sum(chosenTurbines), nbTurbines));

        // Constraint : Two WT need to be spaced by at least minDistance meters
        for (int location = 0; location < nbLocations; location++)
        {
            foreach (int incompatibleLocation in incompatibilitySets[location])
            {
                model.Constraint(model.Leq(
                        model.Sum(chosenTurbines[location], chosenTurbines[incompatibleLocation]), 1));
            }
        }

        // Wind and power at each location for every possible wind angle
        turbineWind = new HxExpression[nbLocations, nbDegrees];
        turbinePower = new HxExpression[nbLocations, nbDegrees];
        for (int location1 = 0; location1 < nbLocations; location1++)
        {
            List<double> wt1 = new List<double>();
            wt1.Add(pointsX[location1]);
            wt1.Add(pointsY[location1]);

            for (int angle = 0; angle < nbDegrees; angle++)
            {
                HxExpression sumWakeLossesSqd = model.Sum();

                for (int location2 = 0; location2 < nbLocations; location2++)
                {
                    List<double> wt2 = new List<double>();
                    wt2.Add(pointsX[location2]);
                    wt2.Add(pointsY[location2]);

                    double wakeLoss = WakeLoss(wt1, wt2, degrees[angle]);
                    sumWakeLossesSqd.AddOperand(
                            model.Pow(model.Prod(wakeLoss, chosenTurbines[location2]), 2));
                }
                // Wind at location i for a wind incidence angle j
                turbineWind[location1, angle] = model.Prod(
                        inflowSpeed, model.Sub(chosenTurbines[location1], model.Sqrt(sumWakeLossesSqd)));
                // Power of the wind turbine of a WT at i with a wind incidence angle j
                HxExpression polynomialValue = model.Pow(model.Div(
                        model.Sub(turbineWind[location1, angle], minSpeed), optimalSpeed - minSpeed), 3);
                HxExpression polynomialCase = model.Prod(
                        model.Leq(minSpeed, turbineWind[location1, angle]),
                        model.Lt(turbineWind[location1, angle], optimalSpeed));
                HxExpression constantCase = model.Prod(
                        model.Leq(optimalSpeed, turbineWind[location1, angle]),
                        model.Lt(turbineWind[location1, angle], maxSpeed));

                turbinePower[location1, angle] = model.Sum(
                        model.Prod(turbineNominalPower, model.Prod(polynomialCase, polynomialValue)),
                        model.Prod(turbineNominalPower, constantCase));
            }
        }

        // Objective : we maximize the Annual Energy Production (in Wh)
        objective = model.Sum();
        for (int location = 0; location < nbLocations; location++)
        {
            for (int angle = 0; angle < nbDegrees; angle++)
            {
                // To calculate the AEP, we multiply the total power produced 
                // by the number of hours in one year
                objective.AddOperand(8760 * model.Prod(probabilities[angle], turbinePower[location, angle]));
            }
        }
        model.Maximize(objective);

        // Finalize model and solve
        model.Close();
        optimizer.GetParam().SetTimeLimit(timelimit);
        optimizer.Solve();
    }

    /* Read instance data */
    private void ReadData(string fileName)
    {
        JObject data = JObject.Parse(File.ReadAllText(fileName));
        JObject wind = (JObject)data["wind"];

        // Wind Information
        nbDegrees = (int)wind["J"];
        inflowSpeed = (double)wind["speed"];
        thrustCoeff = (double)wind["cT"];
        dynCoeff = (double)wind["kY"];
        degrees = wind["degrees"].ToObject<double[]>();
        probabilities = wind["probs"].ToObject<double[]>();

        // Information on the problem
        radius = (double)data["radius"];
        nbTurbines = (int)data["nT"];
        turbineDiameter = (double)data["D"];
        turbineNominalPower = (double)data["p_rated"] * 1e6;
        minDistance = 2.0 * turbineDiameter;
        optimalSpeed = (double)data["u_rated"];
        minSpeed = (double)data["u_cut_in"];
        maxSpeed = (double)data["u_cut_out"];

        discDelta = turbineDiameter * 1.7;

        // Additional required inputs
        BuildDiscretization();
        ComputeIncompatibilities();
    }

    /*
     *  Build the discretization in a circular field of radius {radius} and
     *  with a distance between points of {discDelta}.
     *
     *  The discretization follows the one described in 'eawe' paper, i.e.
     *  building a regular uniform mesh inside the circle, and a point every
     *  degree on the frontier of the circle.
     */
    private void BuildDiscretization()
    {
        // Maximum number of points in any direction from the center
        int maxOnedirPoints = (int)Math.Round(radius / discDelta) + 1;
        List<int> sides = new List<int>() {-1, 1};

        pointsX = new List<double>();
        pointsY = new List<double>();

        // Interior points
        for (int cX = 0; cX < maxOnedirPoints; cX++)
        {
            foreach (int sideX in sides)
            {
                if (!(cX == 0 && sideX == 1))
                {
                    double newX = sideX * cX * discDelta;
                    for (int cY = 0; cY < maxOnedirPoints; cY++)
                    {
                        foreach (int sideY in sides)
                        {
                            double newY = sideY * cY * discDelta;
                            if (Math.Sqrt(Math.Pow(newX, 2)
                                    + Math.Pow(newY, 2)) < radius
                                    && !(cY == 0 && sideY == 1))
                            {
                                pointsX.Add(newX);
                                pointsY.Add(newY);
                            }
                        }
                    }
                }
            }
        }
        // Points on the border of the circle
        for (int deg = 0; deg < 360; deg++)
        {
            double radians = deg * Math.PI / 180;
            double x = radius * Math.Sin(radians);
            double y = radius * Math.Cos(radians);
            pointsX.Add(x);
            pointsY.Add(y);
        }
    }

    /*
     * Compute incompatibility sets for every location on the field.
     *
     * N_i[i] := {l in 0...nbTurbines : l != i && ||l - i|| < minDistance}
     */
    private void ComputeIncompatibilities()
    {
        nbLocations = pointsX.Count;
        incompatibilitySets = new Dictionary<int, List<int>>();

        for (int location1 = 0; location1 < nbLocations; location1++)
        {
            for (int location2 = location1 + 1; location2 < nbLocations; location2++)
            {
                if (Math.Sqrt(Math.Pow(pointsX[location1] - pointsX[location2], 2)
                        + Math.Pow(pointsY[location1] - pointsY[location2], 2)) < minDistance)
                {
                    if (!incompatibilitySets.ContainsKey(location1))
                        incompatibilitySets[location1] = new List<int>();

                    if (!incompatibilitySets.ContainsKey(location2))
                        incompatibilitySets[location2] = new List<int>();

                    incompatibilitySets[location1].Add(location2);
                    incompatibilitySets[location2].Add(location1);
                }
            }
        }
    }

    /*
     * Calculate the distance between two points, considering a frame of
     * reference with angle theta.
     */
    private Dictionary<string, double> Distances(List<double> location1,
           List<double> location2,
            double theta)
    {
        Dictionary<string, double> res = new Dictionary<string, double>();
        // Rotation of the reference
        double thetaDeg = 270 - theta;
        double thetaRad = thetaDeg * Math.PI / 180;

        double cosWind = Math.Cos(-thetaRad);
        double sinWind = Math.Sin(-thetaRad);

        // Change the reference of the coordinates
        double location2X = (location2[0] * cosWind) - (location2[1] * sinWind);
        double location2Y = (location2[0] * sinWind) + (location2[1] * cosWind);

        double location1X = (location1[0] * cosWind) - (location1[1] * sinWind);
        double location1Y = (location1[0] * sinWind) + (location1[1] * cosWind);

        // Store the values in a "dictionary"
        res["parallel"] = location1X - location2X;
        res["perpendicular"] = location1Y - location2Y;

        return res;

    }

    /*
     * Wake loss evaluated at i, caused by a Wind Turbine at l, considering
     * a wind of angle {theta}.
     * This loss is a simplified version of Bastankhah's Gaussian model.
     */
    private double WakeLoss(List<double> turbine1,
            List<double> turbine2,
            double theta)
    {
        // We calculate parallel and perpendicular distances
        Dictionary<string, double> d = Distances(turbine1, turbine2, theta);
        // If x_i - x_l < 0, the wake loss is zero
        if (d["parallel"] > 0)
        {
            // Standard deviation of the wake deficit
            double sY = dynCoeff * d["parallel"] + turbineDiameter / Math.Sqrt(8.0);
            // Separation of the product in two components
            double coeff = thrustCoeff / (8.0 * Math.Pow(sY / turbineDiameter, 2));
            double pdt1 = 1.0 - Math.Sqrt(1.0 - coeff);
            double pdt2 = Math.Exp(-0.5 * Math.Pow(d["perpendicular"] / sY, 2));

            return pdt1 * pdt2;
        }
        return 0;
    }

    /* 
     * Write the solution in the given file, giving :
     *    - The Annual Energy Production (objective function, in Wh and GWh)
     *    - Location (x,y) of every turbine
     */
    private void WriteSolution(string fileName)
    {
        if (optimizer.GetSolution().GetStatus() == HxSolutionStatus.Feasible
                || optimizer.GetSolution().GetStatus() == HxSolutionStatus.Optimal)
        {
            List<List<double>> points = new List<List<double>>();
            for (int location = 0; location < nbLocations; location++)
            {
                if (chosenTurbines[location].GetValue() == 1)
                {
                    List<double> point = new List<double>() {pointsX[location], pointsY[location]};
                    points.Add(point);
                }
            }

            // If asked, we print the result in the output file
            using StreamWriter output = new StreamWriter(fileName);
            Console.WriteLine("Solution written in file " + fileName);

            output.WriteLine("*************************************");
            output.WriteLine("*** WIND FARM LAYOUT OPTIMIZATION ***");
            output.WriteLine("*************************************");
            output.WriteLine();

            output.WriteLine("Annual Energy Production : "
                    + objective.GetDoubleValue() + " Wh ("
                    + Math.Round((double)objective.GetDoubleValue() / 1e9, 2) + " GWh).");
            output.WriteLine();
            output.WriteLine("Locations of the Wind Turbines : ");

            for (int turbine = 0; turbine < nbTurbines; turbine++)
            {
                output.WriteLine("   - (" + points[turbine][0] + ", " + points[turbine][1] + ")");
            }
        }
        else
        {
            Console.WriteLine("No feasible solution have been found within the allowed time.");
        }

    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: WindFarmLayoutOptimization 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] : "60";

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

import java.util.*;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import java.io.*;
import com.hexaly.optimizer.*;

import com.google.gson.JsonObject;
import com.google.gson.JsonArray;
import com.google.gson.JsonParser;
import com.google.gson.stream.JsonReader;


import java.lang.Math;
import java.nio.file.Files;
import java.nio.file.Paths;

public class WindFarmLayoutOptimization {
    // Wind information 
    private double[] degrees;        // Discretization of the degrees rose
    private int nbDegrees;           // Number of discretizations
    private double[] probabilities;  // Probabilities of wind in every direction
    private double inflowSpeed;      // Inflow wind speed
    private double thrustCoeff;      // Thrust coefficient 
    private double dynCoeff;         // Dynamic coefficient 

    // Information about the problem
    private double radius;           // Radius of the perimeter to fullfil
    private int nbTurbines;          // Number of Wind Turbines
    private double turbineDiameter;     // Diameter of each Turbine
    private double turbineNominalPower; // Power of the Turbine
    private double minDistance;         // Minimal distance between two turbines
    private double optimalSpeed;     // Optimal wind speed (maximum efficiency)
    private double minSpeed;         // Minimal wind speed producing energy
    private double maxSpeed;         // Maximal wind speed producing energy

    // Discretization data
    private double discDelta;               // Distance between two points
    List<Double> pointsX;                   // Available locations (X-coord)
    List<Double> pointsY;                   // Available locations (Y-coord)
    private int nbLocations;                // Number of available locations
    private Map<Integer, List<Integer>> incompatibilitySets; // Incompatibility sets

    // Hexaly Optimizer
    final HexalyOptimizer optimizer;
    // Decision variable
    private HxExpression[] chosenTurbines;

    // Wind at each location
    private HxExpression[][] turbineWind;
    // Power of each wind turbine
    private HxExpression[][] turbinePower;

    // Objective function
    private HxExpression objective;

    public WindFarmLayoutOptimization(HexalyOptimizer optimizer) throws IOException {
        this.optimizer = optimizer;
    }

    /* Create and solve the model */
    public void solve(int timelimit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Decision Variable : List of booleans (one for each Wind Turbine)
        chosenTurbines = new HxExpression[nbLocations];
        for (int location = 0; location < nbLocations; location++) {
            chosenTurbines[location] = model.boolVar();
        }
        // Constraint : We have to choose nbTurbines locations on the nbLocations available positions
        model.constraint(model.eq(model.sum(chosenTurbines), nbTurbines));

        // Constraint : Two WT need to be spaced by at least d_min meters
        for (int location = 0; location < nbLocations; location++) {
            for (int incompatibleLocation : incompatibilitySets.get(location)) {
                model.constraint(model.leq(model.sum(chosenTurbines[location], chosenTurbines[incompatibleLocation]), 1));
            }
        }
        // Wind and power at each location for every possible wind angle
        turbineWind = new HxExpression[nbLocations][nbDegrees];
        turbinePower = new HxExpression[nbLocations][nbDegrees];

        for (int location1 = 0; location1 < nbLocations; location1++) {
            List<Double> wt1 = new ArrayList<Double>();
            wt1.add(pointsX.get(location1));
            wt1.add(pointsY.get(location1));

            for (int angle = 0; angle < nbDegrees; angle++) {
                HxExpression sumWakeLossesSqd = model.sum();

                for (int location2 = 0; location2 < nbLocations; location2++) {
                    List<Double> wt2 = new ArrayList<Double>();
                    wt2.add(pointsX.get(location2));
                    wt2.add(pointsY.get(location2));

                    double wakeLoss = wakeLoss(wt1, wt2, degrees[angle]);
                    sumWakeLossesSqd.addOperand(
                            model.prod(chosenTurbines[location2], Math.pow(wakeLoss, 2)));
                }
                // Wind at location i for a wind incidence angle j
                turbineWind[location1][angle] = model.prod(inflowSpeed,
                        model.sub(chosenTurbines[location1], model.sqrt(sumWakeLossesSqd)));

                // Power of the wind turbine of a WT at i with a wind incidence angle j
                HxExpression polynomialCase = model.prod(
                        model.leq(minSpeed, turbineWind[location1][angle]),
                        model.lt(turbineWind[location1][angle], optimalSpeed));
                HxExpression polynomialValue = model.pow(
                        model.div(model.sub(turbineWind[location1][angle], minSpeed),
                        optimalSpeed - minSpeed), 3);
                HxExpression constantCase = model.prod(
                        model.leq(optimalSpeed, turbineWind[location1][angle]),
                        model.lt(turbineWind[location1][angle], maxSpeed));
                turbinePower[location1][angle] = model.sum(
                        model.prod(turbineNominalPower,
                                model.prod(polynomialCase, polynomialValue)),
                        model.prod(turbineNominalPower, constantCase));
            }
        }
        // Objective : we maximize the Annual Energy Production (in Wh)
        // To calculate the AEP, we multiply the total power produced by the number of hours in one year
        objective = model.sum();
        for (int location = 0; location < nbLocations; location++) {
            for (int angle = 0; angle < nbDegrees; angle++) {
                objective.addOperand(model.prod(8760, 
                        model.prod(probabilities[angle], turbinePower[location][angle])));
            }
        }
        model.maximize(objective);

        // Finalize model and solve
        model.close();
        optimizer.getParam().setTimeLimit(timelimit);
        optimizer.solve();
    }

    /* Read instance data */
    public void readData(String fileName) throws IOException {
        JsonReader reader = new JsonReader(
                new InputStreamReader(new FileInputStream(fileName)));
        JsonObject data = JsonParser.parseReader(reader).getAsJsonObject();

        // Wind information
        JsonObject wind = data.getAsJsonObject("wind");
        nbDegrees = wind.get("J").getAsInt();
        inflowSpeed = wind.get("speed").getAsDouble();
        thrustCoeff = wind.get("cT").getAsDouble();
        dynCoeff = wind.get("kY").getAsDouble();

        JsonArray degreesArray = (JsonArray) wind.get("degrees");
        JsonArray probsArray = (JsonArray) wind.get("probs");

        degrees = new double[nbDegrees];
        probabilities = new double[nbDegrees];
        for (int angle = 0; angle < nbDegrees; angle++) {
            degrees[angle] = degreesArray.get(angle).getAsDouble();
            probabilities[angle] = probsArray.get(angle).getAsDouble();
        }

        // Information about the problem
        radius = data.get("radius").getAsDouble();
        nbTurbines = data.get("nT").getAsInt();
        turbineDiameter = data.get("D").getAsDouble();
        turbineNominalPower = data.get("p_rated").getAsDouble() * 1e6;
        minDistance = 2.0 * turbineDiameter;
        optimalSpeed = data.get("u_rated").getAsDouble();
        minSpeed = data.get("u_cut_in").getAsDouble();
        maxSpeed = data.get("u_cut_out").getAsDouble();

        discDelta = turbineDiameter * 1.7;
        // Additional required inputs
        buildDiscretization();
        computeIncompatibilities();
    }

    /*
     *  Build the discretization in a circular field of radius {radius} and
     *  with a distance between points of {discDelta}.
     *
     *  The discretization follows the one described in 'eawe' paper, i.e.
     *  building a regular uniform mesh inside the circle, and a point every
     *  degree on the frontier of the circle.
     */
    public void buildDiscretization() {
        // Maximum number of points in any direction from the center
        int maxOnedirPoints = (int)Math.round(radius / discDelta) + 1;
        int[] sides = {-1, 1};

        pointsX = new ArrayList<Double>();
        pointsY = new ArrayList<Double>();

        // Interior points
        for (int cX = 0; cX < maxOnedirPoints; cX++) {
            for (int sideX : sides) {
                if (!(cX == 0 && sideX == 1)) {
                    double newX = sideX * cX * discDelta;
                    for (int cY = 0; cY < maxOnedirPoints; cY++) {
                        for (int sideY : sides) {
                            double newY = sideY * cY * discDelta;
                            if (Math.sqrt(Math.pow(newX, 2) 
                                    + Math.pow(newY, 2)) < radius
                                    && !(cY == 0 && sideY == 1)) {
                                pointsX.add(newX);
                                pointsY.add(newY);
                            }
                        }
                    } 
                }
            }
        }
        // Points on the border of the circle
        for (int deg = 0; deg < 360; deg++) {
            double radians = Math.toRadians(deg);
            double x = radius * Math.sin(radians);
            double y = radius * Math.cos(radians);
            pointsX.add(x);
            pointsY.add(y);
        }
    }

    /*
     * Compute incompatibility sets for every location on the field.
     *
     * N_i[i] := {l in 0...nbLocations : l != i && ||l - i|| < minDistance}
     */
    public void computeIncompatibilities() {
        nbLocations = pointsX.size();

        incompatibilitySets = new HashMap<>();
        for (int location = 0; location < nbLocations; location++) {
            incompatibilitySets.put(location, new ArrayList<>());
        }

        for (int location1 = 0; location1 < nbLocations; location1++) {
            for (int location2 = location1 + 1; location2 < nbLocations; location2++) {
                if (Math.sqrt(Math.pow(pointsX.get(location1) - pointsX.get(location2), 2)
                        + Math.pow(pointsY.get(location1) - pointsY.get(location2), 2)) < minDistance) {
                    incompatibilitySets.get(location1).add(location2);
                    incompatibilitySets.get(location2).add(location1);
                }
            }
        }
    }

    /*
     * Calculate the distance between two points, considering a frame of
     * reference with angle theta.
     */
    public Map<String, Double> distances(List<Double> location1,
            List<Double> location2, 
            double theta) {
        Map<String, Double> res = new HashMap<>();
        // Rotation of the reference
        double thetaDeg = 270 - theta;
        double thetaRad = Math.toRadians(thetaDeg);

        double cosWind = Math.cos(-thetaRad);
        double sinWind = Math.sin(-thetaRad);

        // Change the reference of the coordinates
        double location2X = (location2.get(0) * cosWind) - (location2.get(1) * sinWind);
        double location2Y = (location2.get(0) * sinWind) + (location2.get(1) * cosWind);

        double location1X = (location1.get(0) * cosWind) - (location1.get(1) * sinWind);
        double location1Y = (location1.get(0) * sinWind) + (location1.get(1) * cosWind);

        res.put("parallel", location1X - location2X);
        res.put("perpendicular", location1Y - location2Y);

        return res;

    }

    /*
     * Wake loss evaluated at i, caused by a Wind Turbine at l, considering
     * a wind of angle {theta}.
     * This loss is a simplified version of Bastankhah's Gaussian model.
     */
    public double wakeLoss(List<Double> turbine1,
            List<Double> turbine2,
            double theta) {
        // We calculate parallel and perpendicular distances
        Map<String, Double> d = distances(turbine1, turbine2, theta);
        // If x_i - x_l < 0, the wake loss is zero
        if (d.get("parallel") > 0) {
            // Standard deviation of the wake deficit
            double sY = dynCoeff * d.get("parallel") + turbineDiameter / Math.sqrt(8.0);
            // Separation of the product in two components
            double coeff = thrustCoeff / (8.0 * Math.pow(sY / turbineDiameter, 2));
            double pdt1 = 1.0 - Math.sqrt(1.0 - coeff);
            double pdt2 = Math.exp(-0.5 * Math.pow(d.get("perpendicular") / sY, 2));

            return pdt1 * pdt2;
        }
        return 0;
    }

    /* 
     * Write the solution in the given file, giving :
     *    - The instance 
     *    - The timeLimit
     *    - The Annual Energy Production (objective function, in Wh and GWh)
     *    - Location (x,y) of every turbine
     */
    public void writeSolution(String fileName) throws IOException {
        if (optimizer.getSolution().getStatus() == HxSolutionStatus.Feasible
                || optimizer.getSolution().getStatus() == HxSolutionStatus.Optimal) {
            List<List<Double>> points = new ArrayList<>();
            for (int location = 0; location < nbLocations; location++) {
                if (chosenTurbines[location].getValue() == 1) {
                    points.add(Arrays.asList(pointsX.get(location), pointsY.get(location)));
                }
            }
            // If asked, we print the result in the output file
            try (PrintWriter output = new PrintWriter(fileName)) {
                System.out.println("Solution written in file " + fileName);

                output.print("*************************************\n");
                output.print("*** WIND FARM LAYOUT OPTIMIZATION ***\n");
                output.print("*************************************\n\n");

                double r_AEP = Math.round(objective.getDoubleValue() / 1e9 * 100.0) / 100.0;
                output.print("Annual Energy Production : " + objective.getDoubleValue()
                        + " Wh (" + r_AEP + " GWh).\n\n");
                output.print("Locations of the Wind Turbines : \n");

                for (int turbine = 0; turbine < nbTurbines; turbine++) {
                    output.print("   - (" + points.get(turbine).get(0) + ", "
                            + points.get(turbine).get(1) + ")\n");
                }
            } 
        } else {
            System.err.println("No feasible solution have been found within the allowed time.");
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java WindFarmLayoutOptimization inputFile"
                    + " [solFile] [timeLimit]");
            System.exit(1);
        }
        // Read the arguments
        String instanceFile = args[0];
        String outputFile = args.length > 1 ? args[1] : null;
        String strTimeLimit = args.length > 2 ? args[2] : "60";

        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            WindFarmLayoutOptimization model = new WindFarmLayoutOptimization(optimizer);

            model.readData(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));

            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            // If an exception occurs, we catch it, print it and give its StackTrace
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}