# Copyright (c) 2026, Artelys (https://www.artelys.com)
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
# SPDX-License-Identifier: MPL-2.0
#
#
# @author Léo Rossignol <leo.rossignol at artelys.com>
# @author Alexis Godefroy <alexis.godefroy at artelys.com>
#
import os

import gzip
import numpy as np
import pandas as pd
import pathlib
import shutil

from emerge_res_integration_tool import config, constants

try:
    from knitro import KN_OBJGOAL_MINIMIZE
    from knitro.numpy import (
        KN_add_cons,
        KN_set_con_lobnds,
        KN_set_con_upbnds,
        KN_add_con_linear_struct,
        KN_new,
        KN_set_param,
        KN_PARAM_OUTDIR,
        KN_PARAM_OUTLEV,
        KN_PARAM_OUTMODE,
        KN_PARAM_OUT_CSVINFO,
        KN_PARAM_FINDIFF_RELSTEPSIZE,
        KN_PARAM_FINDIFF_ESTNOISE,
        KN_PARAM_NLP_ALGORITHM,
        KN_PARAM_GRADOPT,
        KN_PARAM_HESSOPT,
        KN_PARAM_OPTTOLABS,
        KN_PARAM_FEASTOL,
        KN_set_obj_goal,
        KN_add_vars,
        KN_set_var_lobnds,
        KN_set_var_upbnds,
        KN_set_var_primal_init_values,
        KN_RC_EVALFC,
        KN_get_obj_value,
        KN_RC_USER_TERMINATION,
        KN_add_eval_callback,
        KN_set_newpt_callback,
        KN_PARAM_NUMTHREADS,
        KN_solve,
        KN_get_solution,
        KN_free,
    )
except ImportError as exc:
    raise ImportError(
        "Knitro Python bindings are required to run the optimizer"
    ) from exc


class KnitroBatteryOptimizer:
    """
    Builds and solves the battery setpoint optimisation using Knitro.
    """

    def __init__(
        self,
        network_model,
        metrix_simulator,
        timeseries_manager,
        battery_fleet,
        output_dir: pathlib.Path,
    ):
        self.network_model = network_model
        self.metrix_simulator = metrix_simulator
        self.timeseries_manager = timeseries_manager
        self.battery_fleet = battery_fleet
        self.output_dir = output_dir
        self.horizon = timeseries_manager.horizon
        self.num_batteries = len(self.battery_fleet.batteries)
        self.num_p_vars = self.num_batteries * self.horizon
        self.num_e0_vars = self.num_batteries
        self.num_vars = self.num_p_vars + self.num_e0_vars
        self.final_solution = None
        self.power_scale = max(bat.pmax for bat in self.battery_fleet.batteries)  # Scaling factor for optim variables
        self._baseline_adequacy_costs = None
        self._baseline_redispatching_costs = None
        self._baseline_res_production = None
        self._final_res_production = None

    def set_baseline_run(self, run_metrix_normal: pd.DataFrame):
        """
        Store baseline adequacy/redispatching costs and RES production from a Metrix
        run.
        """
        if run_metrix_normal is None:
            self._baseline_adequacy_costs = None
            self._baseline_redispatching_costs = None
            self._baseline_res_production = None
            return

        # Build generator cost mapping consistent with the optimisation cost
        battery_ids = self.network_model.get_battery_ids()
        hydro_ids = self.network_model.get_hydro_ids()
        thermal_ids = self.network_model.get_thermal_ids()
        solar_ids = self.network_model.get_solar_ids()
        wind_ids = self.network_model.get_wind_ids()

        generator_cost_map = {}
        for gen_id in battery_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_BATTERY_UP, constants.ADEQUACY_BATTERY_DOWN)
        for gen_id in hydro_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_HYDRO_UP, constants.ADEQUACY_HYDRO_DOWN)
        for gen_id in thermal_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_THERMAL_UP, constants.ADEQUACY_THERMAL_DOWN)
        for gen_id in solar_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_SOLAR_UP, constants.ADEQUACY_SOLAR_DOWN)
        for gen_id in wind_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_WIND_UP, constants.ADEQUACY_WIND_DOWN)

        adequacy = self.metrix_simulator.get_adequacy_costs(generator_cost_map)
        redispatch = self.metrix_simulator.get_redispatching_costs()

        self._baseline_adequacy_costs = np.asarray(adequacy, dtype=float)[: self.horizon]
        self._baseline_redispatching_costs = np.asarray(redispatch, dtype=float)[: self.horizon]
        self._baseline_res_production = self._compute_res_production_from_results()

    def _p_var_index(self, bat_index: int, timestep: int) -> int:
        return bat_index * self.horizon + timestep

    def _e0_var_index(self, bat_index: int) -> int:
        return self.num_p_vars + bat_index

    def _compute_res_production_from_results(self) -> np.ndarray:
        """
        Compute total renewable generation (solar + wind + hydro) per timestep
        from the latest Metrix results file.
        """
        # Ensure results.csv exists (MetrixSimulation handles extraction if needed)
        try:
            self.metrix_simulator.extract_metrix_gz_results()
        except Exception:
            # If extraction fails, try reading directly; if that also fails, return zeros.
            config.logger.warning("Could not extract gzipped Metrix results, trying direct CSV read.")

        try:
            results_df = pd.read_csv(self.metrix_simulator.csv_results_file_path, sep=";")
        except Exception:
            config.logger.error("Unable to read Metrix results for RES production calculation.")
            return np.zeros(self.horizon, dtype=float)

        solar_ids = self.network_model.get_solar_ids()
        wind_ids = self.network_model.get_wind_ids()
        hydro_ids = self.network_model.get_hydro_ids()
        
        work_df = pd.DataFrame(index=range(self.horizon))
        timeseries_df = self.timeseries_manager.timeseries_df.reset_index(drop=True)
        results_df = results_df.reset_index(drop=True)

        # Calculating renewable production
        for gen_id in list(solar_ids) + list(wind_ids) + list(hydro_ids):
            col_base = gen_id
            col_adequacy = f"INIT_BAL_GEN_{gen_id}"
            col_redispatch = f"GEN_{gen_id}"

            if col_base in timeseries_df.columns:
                work_df[col_base] = timeseries_df[col_base].iloc[: self.horizon].values
            
            if col_adequacy in results_df.columns:
                work_df[col_adequacy] = results_df[col_adequacy].iloc[: self.horizon].values
            if col_redispatch in results_df.columns:
                work_df[col_redispatch] = results_df[col_redispatch].iloc[: self.horizon].values

        if work_df.empty or work_df.shape[1] == 0:
            return np.zeros(self.horizon, dtype=float)

        return work_df.sum(axis=1).to_numpy(dtype=float)

    def _initial_guess(self, use_custom_guess=False, custom_setpoints_mw=None, custom_e0_mwh=None):
        x0 = np.zeros(self.num_vars, dtype=float)
        
        if use_custom_guess:
            if custom_setpoints_mw is None or custom_e0_mwh is None:
                config.logger.warning("Custom guess requested but values not provided, using heuristic")
                use_custom_guess = False
            elif len(custom_setpoints_mw) != self.horizon:
                raise ValueError(
                    f"Custom setpoints has {len(custom_setpoints_mw)} values, "
                    f"expected {self.horizon}"
                )
            else:
                # Apply custom guess to all batteries
                for b_idx, bat in enumerate(self.battery_fleet.batteries):
                    profile_scaled = np.array(custom_setpoints_mw) / self.power_scale
                    start = self._p_var_index(b_idx, 0)
                    x0[start : start + self.horizon] = profile_scaled
                    x0[self._e0_var_index(b_idx)] = custom_e0_mwh / self.power_scale
                
                config.logger.info("Using custom initial guess")
        
        # Heuristic (par défaut ou si custom a échoué)
        if not use_custom_guess:
            rng = np.random.default_rng(seed=42)
            
            for b_idx, bat in enumerate(self.battery_fleet.batteries):
                profile = np.zeros(self.horizon)
                night_hours = range(0, self.horizon // 4)
                day_hours = range(self.horizon // 2, 3 * self.horizon // 4)

                charge_power = -0.2 * bat.pmax / self.power_scale
                discharge_power = 0.2 * bat.pmax / self.power_scale

                for t in night_hours:
                    profile[t] = charge_power
                for t in day_hours:
                    profile[t] = discharge_power

                noise = rng.uniform(
                    low=-0.1 * bat.pmax / self.power_scale,
                    high=0.1 * bat.pmax / self.power_scale,
                    size=self.horizon,
                )
                profile += noise
                start = self._p_var_index(b_idx, 0)
                x0[start : start + self.horizon] = profile
                soc0 = self.timeseries_manager.get_initial_soc(bat.id)
                x0[self._e0_var_index(b_idx)] = soc0 * bat.emax / self.power_scale
            
            config.logger.info("Using heuristic initial guess")
        
        return x0

    def _build_linear_constraints(self, kc):
        """
        Add linear SoC constraints:
        0 <= E0 - sum_{t1<=t} P(t1) <= Emax
        """
        # Initialization
        num_cons_per_bat = self.horizon
        total_cons = num_cons_per_bat * self.num_batteries
        # Compute constraints bounds per battery
        con_loBnds = []
        con_upBnds = []
        for b_idx, bat in enumerate(self.battery_fleet.batteries):
            for t in range(self.horizon):
                con_loBnds.append(0.0)
                con_upBnds.append(bat.emax / self.power_scale)
        # Add constraints to the Knitro optim instance
        KN_add_cons(kc, total_cons)
        KN_set_con_lobnds(kc, cLoBnds=con_loBnds)
        KN_set_con_upbnds(kc, cUpBnds=con_upBnds)
        # Set linear structure for each constraint
        c_idx = 0
        for b_idx, bat in enumerate(self.battery_fleet.batteries):
            e0_idx = self._e0_var_index(b_idx)

            for t in range(self.horizon):
                var_indices = [e0_idx]
                coefs = [1.0]
                for tau in range(t + 1):
                    var_indices.append(self._p_var_index(b_idx, tau))
                    coefs.append(-1.0)
                KN_add_con_linear_struct(kc, c_idx, var_indices, coefs)
                c_idx += 1

    def _evaluate_cost(self, x):
        """
        Apply setpoints, run Metrix once, and return the objective value.
        """
        # Get the dictionary of setpoints' lists for all batteries
        battery_setpoints_scaled = self.extract_setpoints(x)
        battery_setpoints_mw = {
            bat_id: profile * self.power_scale
            for bat_id, profile in battery_setpoints_scaled.items()
        }
        # Extract E0
        e0_values = {}
        for b_idx, bat in enumerate(self.battery_fleet.batteries):
            e0_idx = self._e0_var_index(b_idx)
            e0_values[bat.id] = x[e0_idx] * self.power_scale
        # Apply the new setpoints to the timeseries to prepare next iteration
        self.timeseries_manager.timeseries_df = (
            self.timeseries_manager.initial_df.copy()
        )

        self.timeseries_manager.apply_battery_setpoints(battery_setpoints_mw)
        #config.logger.debug(
        #    f"New timeseries df: {self.timeseries_manager.timeseries_df}"
        #)
        # Run PowSyBl-Metrix on new timeseries
        config.logger.info("Starting PowSyBl-Metrix run")
        return_code, std_err = self.metrix_simulator.run_metrix()
        config.logger.info(f"Metrix returned with exit code {return_code}")
        #config.logger.debug(std_err)
        # Compute adequacy costs: map generators IDs to their cost columns (up/down),
        # then process csv results file accordingly
        battery_ids = self.network_model.get_battery_ids()
        hydro_ids = self.network_model.get_hydro_ids()
        thermal_ids = self.network_model.get_thermal_ids()
        solar_ids = self.network_model.get_solar_ids()
        wind_ids = self.network_model.get_wind_ids()
        generator_cost_map = {}
        for gen_id in battery_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_BATTERY_UP, constants.ADEQUACY_BATTERY_DOWN)
        for gen_id in hydro_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_HYDRO_UP, constants.ADEQUACY_HYDRO_DOWN)
        for gen_id in thermal_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_THERMAL_UP, constants.ADEQUACY_THERMAL_DOWN)
        for gen_id in solar_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_SOLAR_UP, constants.ADEQUACY_SOLAR_DOWN)
        for gen_id in wind_ids:
            generator_cost_map[gen_id] = (constants.ADEQUACY_WIND_UP, constants.ADEQUACY_WIND_DOWN)
        adequacy_costs = self.metrix_simulator.get_adequacy_costs(generator_cost_map)
        config.logger.debug(f"Assessed adequacy costs: {adequacy_costs}")
        # Compute redispatching costs (directly extracted from csv results file)
        redispatching_costs = self.metrix_simulator.get_redispatching_costs()
        config.logger.debug(f"Assessed redispatching costs: {redispatching_costs}")
        
        # Check errors in Metrix run
        results_path = pathlib.Path(self.metrix_simulator.csv_results_file_path)
        gz_path = pathlib.Path(str(results_path) + ".gz")

        if not results_path.exists():
            if gz_path.exists():
                with gzip.open(gz_path, "rb") as f_in, open(results_path, "wb") as f_out:
                    shutil.copyfileobj(f_in, f_out)
            else:
                raise FileNotFoundError(
                    f"Metrix results file not found at {results_path} (or {gz_path})"
                )

        run_metrix_df = pd.read_csv(results_path, sep=";")

        # Timesteps with error
        error_indices = run_metrix_df.index[run_metrix_df["ERROR_CODE"] == 1].tolist()

        if error_indices:
            config.logger.warning(
                f"Metrix failed on timesteps {error_indices} -> applying penalty costs"
            )
            for idx in error_indices:
                adequacy_costs.iloc[idx] = 500_000
                redispatching_costs.iloc[idx] = 500_000
        
        total_cost = float((adequacy_costs + redispatching_costs).sum())
        # Update the list of best solutions found by Knitro
        self.current_evaluation = {
            "cost": total_cost,
            "setpoints": battery_setpoints_mw,
            "e0_values": e0_values,
            "adequacy_costs": adequacy_costs.tolist(),  # Par timestep
            "redispatching_costs": redispatching_costs.tolist(),  # Par timestep
        }
        # Return objective function value = total cost (AD + REDISP)
        config.logger.info(f"Assessed total cost: {total_cost}")
        return total_cost

    def extract_setpoints(self, x):
        """
        Convert a solution vector into a dictionary {bat_id: profile}.
        """
        x = np.asarray(x)
        # Initialize empty dictionary of battery setpoints
        battery_setpoints = {}
        for b_idx, bat in enumerate(self.battery_fleet.batteries):
            # Append an element to the dictionary with battery ID and
            # list of setpoints on the considered time horizon
            start = self._p_var_index(b_idx, 0)
            end = start + self.horizon
            battery_setpoints[bat.id] = x[start:end]
        return battery_setpoints

    def _save_best_solutions(self, filepath: str):
        """
        Persist the best solutions to a CSV for inspection.
        """
        if self.final_solution is None:
            config.logger.warning("No final solution to save")
            return
        
        os.makedirs(os.path.dirname(filepath), exist_ok=True)
        
        adequacy = np.asarray(self.final_solution["adequacy_costs"], dtype=float)
        redispatch = np.asarray(self.final_solution["redispatching_costs"], dtype=float)
        per_step_cost = adequacy + redispatch
        total_cumulated_cost = np.cumsum(per_step_cost)

        if (
            self._baseline_adequacy_costs is not None
            and self._baseline_redispatching_costs is not None
        ):
            baseline_total = (
                np.asarray(self._baseline_adequacy_costs, dtype=float)
                + np.asarray(self._baseline_redispatching_costs, dtype=float)
            )
        else:
            baseline_total = np.zeros(self.horizon, dtype=float)

        per_step_saved = baseline_total - per_step_cost
        total_cumulated_saved = np.cumsum(per_step_saved)

        if (
            self._baseline_res_production is not None
            and self._final_res_production is not None
        ):
            res_delta = (
                np.asarray(self._final_res_production, dtype=float)
                - np.asarray(self._baseline_res_production, dtype=float)
            )
            total_cumulated_avoided_curtailment = np.cumsum(res_delta)
        else:
            total_cumulated_avoided_curtailment = np.zeros(self.horizon, dtype=float)

        data = {"Timestep": list(range(self.horizon))}

        # Columns per battery
        for bat in self.battery_fleet.batteries:
            bat_id = bat.id
            profile_mw = np.asarray(
                self.final_solution["setpoints"][bat_id], dtype=float
            )[: self.horizon]
            e0 = float(self.final_solution["e0_values"][bat_id])
            soc_mwh = e0 - np.cumsum(profile_mw)
            soc = soc_mwh / float(bat.emax) if bat.emax > 0 else np.zeros(self.horizon)

            data[f"{bat_id}_E0 (MWh)"] = [e0] * self.horizon
            data[f"{bat_id}_profile (MW)"] = profile_mw
            data[f"{bat_id}_SOC"] = soc

        # Aggregated columns
        data[constants.OUTPUT_CUM_COST] = total_cumulated_cost
        data[constants.OUTPUT_ADEQUACY_COST] = adequacy
        data[constants.OUTPUT_REDISP_COST] = redispatch
        data[constants.OUTPUT_INITIAL_COST] = baseline_total
        data[constants.OUTPUT_CUM_SAVED_COST] = total_cumulated_saved
        data[constants.OUTPUT_CUM_CURTAILMENT] = total_cumulated_avoided_curtailment

        pd.DataFrame(data).to_csv(filepath, sep=";", index=False)

    def solve(self, enable_multistart=False, findiff_relstepsize=None, use_custom_guess=False, custom_setpoints_mw=None, custom_e0_mwh=None):
        # 1. Initialization steps -----------------------------------------------------
        # Folder path for outputs
        log_dir = self.output_dir / "knitro_logs"
        os.makedirs(log_dir, exist_ok=True)
        # Initialize dictionary for metrics about optim process
        iteration_data = {
            "best_obj": float("inf"),
            "iterations_without_improvement": 0,
            "max_stagnant_iterations": 500,  # stop if no significant improvement
            "improvement_threshold": 1e-1,  # minimal improvement (0.1%)
            "total_iterations": 0,
            "total_evals": 0,
        }

        # 2. Create and parameterize Knitro optim instance ----------------------------
        kc = KN_new()
        # # Logging, Knitro outputs parameters
        KN_set_param(kc, KN_PARAM_OUTDIR, log_dir)
        KN_set_param(kc, KN_PARAM_OUTLEV, 6)
        KN_set_param(kc, KN_PARAM_OUTMODE, 1)
        KN_set_param(kc, KN_PARAM_OUT_CSVINFO, 1)
        # # Mathematical parameters
        if findiff_relstepsize is not None:
            KN_set_param(kc, KN_PARAM_FINDIFF_RELSTEPSIZE, findiff_relstepsize)
            config.logger.info(f"Using FINDIFF_RELSTEPSIZE = {findiff_relstepsize}")
        else:
            KN_set_param(kc, KN_PARAM_FINDIFF_RELSTEPSIZE, 1e-2)
        # KN_set_param(kc, KN_PARAM_FINDIFF_ESTNOISE, 1) # parameter also helping stepsizes of derivatives, but I find less satisfactory solutions
        KN_set_param(kc, KN_PARAM_NLP_ALGORITHM, 4)  # better if black box use of Knitro
        #KN_set_param(kc, KN_PARAM_OPTTOLABS, 1e-4)  # stop parameter
        #KN_set_param(kc, KN_PARAM_FEASTOL, 1e-6)  # stop parameter
        KN_set_param(kc, KN_PARAM_GRADOPT, 2)
        KN_set_param(kc, KN_PARAM_HESSOPT, 3)
        # # Multistart, nb of threads
        if enable_multistart:
            KN_set_param(kc, "ms_enable", 1)
            KN_set_param(kc, "ms_maxsolves", 3)
        else:
            KN_set_param(kc, "ms_enable", 0)
        KN_set_param(kc, KN_PARAM_NUMTHREADS, 1)
        # # Set problem as minimization
        KN_set_obj_goal(kc, KN_OBJGOAL_MINIMIZE)

        # 3. Add variables (nb and bounds) to the instance ----------------------------
        # # - For each battery, n setpoints (where n is the number of timesteps) P_i
        # #   (1 <= i <= n) scaled with Pmax (unit: None).
        # # - In addition, for each battery, an initial charge e0 also scaled with
        # #   Pmax (unit: h).
        var_loBnds = []
        var_upBnds = []
        for b_idx, bat in enumerate(self.battery_fleet.batteries):
            for t in range(self.horizon):
                var_loBnds.append(-bat.pmax / self.power_scale)
                var_upBnds.append(bat.pmax / self.power_scale)
            var_loBnds.append(0.0)
            var_upBnds.append(bat.emax / self.power_scale)
        KN_add_vars(kc, self.num_vars)
        KN_set_var_lobnds(kc, xLoBnds=var_loBnds)
        KN_set_var_upbnds(kc, xUpBnds=var_upBnds)

        # 4. Set initial guess (starting point) ---------------------------------------
        x0 = self._initial_guess(
            use_custom_guess=use_custom_guess,
            custom_setpoints_mw=custom_setpoints_mw,
            custom_e0_mwh=custom_e0_mwh
        )
        KN_set_var_primal_init_values(kc, xInitVals=x0)

        # 5. Add the set of linear constraints to the instance ------------------------
        # # For each battery, there are n constraints (where n is the number of
        # # timesteps), all of them can be expressed with a linear relation between
        # # optim variables:
        # # - For each battery, for each t_i in [t_1, t_n]:
        # #   0 <= E0 - sum_{t_k<=t_i} P(t_k) <= Emax
        self._build_linear_constraints(kc)

        # 6. Define objective function ------------------------------------------------
        # # For this, we have to define the eval callback which contains Metrix calls
        def callback_eval(kc, cb, evalRequest, evalResult, userParams):
            if evalRequest.type != KN_RC_EVALFC:
                return 0
            x = evalRequest.x
            # Compute adequacy + redispatching costs for given set of setpoints (x)
            obj_value = self._evaluate_cost(x)
            evalResult.obj = obj_value
            # Update counter of function evaluations
            iteration_data["total_evals"] += 1
            return 0

        # 7. Define an additional user-defined callback -------------------------------
        # # This one is used to compute metrics about optimization process after each
        # # new solution estimate: tracking of the optimization iterations, stop
        # # process if user-defined criteria are met
        def callback_newpoint(kc, x, lambda_, userParams):
            iteration_data["total_iterations"] += 1
            try:
                current_obj = KN_get_obj_value(kc)
            except:
                return 0
            if iteration_data["total_iterations"] < 15:
                return 0
            # Compute and log improvement compared to the best value of objective func
            improvement = iteration_data["best_obj"] - current_obj
            if improvement > iteration_data["improvement_threshold"]:
                config.logger.info(
                    f"[Knitro Iter {iteration_data['total_iterations']}] "
                    f"Obj = {current_obj:.2f} | "
                    f"Improvement: {improvement:.2f}"
                )
                iteration_data["best_obj"] = current_obj
                iteration_data["iterations_without_improvement"] = 0
            else:
                iteration_data["iterations_without_improvement"] += 1
                config.logger.info(
                    f"[Knitro Iter {iteration_data['total_iterations']}] "
                    f"Obj = {current_obj:.2f} | "
                    f"No improvement ({iteration_data['iterations_without_improvement']}/{iteration_data['max_stagnant_iterations']})"
                )
            # Stop process if there is no improvement acording to user-defined criteria
            if (
                iteration_data["iterations_without_improvement"]
                >= iteration_data["max_stagnant_iterations"]
            ):
                config.logger.warning(
                    f"[Knitro] Stopping: {iteration_data['max_stagnant_iterations']} "
                    f"iterations without improvement (best: {iteration_data['best_obj']:.2f})"
                )
                return KN_RC_USER_TERMINATION
            return 0

        # 8. Add callbacks to the Knitro optim instance, solve ------------------------
        KN_add_eval_callback(kc, evalObj=True, funcCallback=callback_eval)
        KN_set_newpt_callback(kc, callback_newpoint)
        nStatus = KN_solve(kc)
        if nStatus == KN_RC_USER_TERMINATION:
            config.logger.info(
                f"Knitro process aborted after {iteration_data['total_iterations']} iterations"
                f" ({iteration_data['total_evals']} total Metrix calls)"
            )
        nStatus, objSol, x, lambda_ = KN_get_solution(kc)

        final_cost = self._evaluate_cost(x)
        final_adequacy = np.asarray(self.current_evaluation["adequacy_costs"], dtype=float)
        final_redispatch = np.asarray(self.current_evaluation["redispatching_costs"], dtype=float)
        final_total_cumulated_cost = float(np.cumsum(final_adequacy + final_redispatch)[-1])
        # RES production for the final (optimised) run
        self._final_res_production = self._compute_res_production_from_results()
        self.final_solution = {
            "status": nStatus,
            "objective": final_total_cumulated_cost,
            "x": x,
            "lambda": lambda_,
            "setpoints": self.current_evaluation["setpoints"],
            "e0_values": self.current_evaluation["e0_values"],
            "adequacy_costs": self.current_evaluation["adequacy_costs"],
            "redispatching_costs": self.current_evaluation["redispatching_costs"],
        }
        
        # # Free Knitro optim instance
        KN_free(kc)
        # # Save best solutions in a csv file
        self._save_best_solutions(self.output_dir / constants.BEST_RESULTS_PATH)

        return self.final_solution
