test_closedloop_tube_main.py

#!/usr/bin/env python
# coding: utf-8
'''
### Biomethane production Optimal Control with the (online-tube) 'NMPC' feedback online controller (Biogoals.Twin tool)
**Authors:** Davide Carecci  
**Created:** May 2025  
**Last Modified:** April 23, 2026  
**Version:** 1.1  
**Institution:** Politecnico di Milano  
#### Revision History

| Date       | Version | Author(s)              | Description                                                                                                    |
|------------|---------|------------------------|----------------------------------------------------------------------------------------------------------------|
| 2025-05-01 | 1.0     | D. Carecci             | Initial implementation for experimental validation over the reactors in the BioTA lab (USM, Valparaiso, Chile) |
| 2026-04-23 | 1.1     | D. Carecci             | Code and folder cleaning and documentation for handle over a copy to A2A S.p.A                                 |

_Note_: This code is meant to be run on a Raspberry Pi connected to the real system in the BioTA lab (USM, Valparaiso, Chile) for the experimental validation of the tube-based NMPC controller. 
It reads the necessary input files, runs the model and EKF, and saves the results for further analysis and plotting. 
The code is structured to allow for modifications and extensions as needed for different test scenarios or model configurations. 

_Note_: 1. modify manually the path to the files to be loaded and for saving the results (e.g. 'read_all')
        2. Modify 'modelname', 'save_out' and 'final_challange' variables
        3. Output path to 'prepare_inputs' 
        4. the inorganic carbon concentration inside the feedstocks ('sic') (set this equal to the values computed by the Modelica model).

_Note_: to run NMPC library on RPi: 
        1. Comment the ```sys.path.insert(0, 'c:\modelonimpact-1.8.1\oct-dist\myenv\lib\site-packages')``` in 'AM2HNgb_xi.py' and 'NMPC_controller.py' files. 
        2. Uncomment the "In[15]" cell
        3. Modify the name of y_df_data_on and y_df_data_off files in the "In[6]" cell, to read the real measurements from the RPi instead of the simulated ones.
        4. Combine this code with the codes present in the '*Biogoals.Architecture*' repository.

_Note_: needed at start-up: 
        1. initial states and covariance matrix of openloop models and EKF (they must contain a row at Timestamp = start_timestamp model and EKF)
        2. z_star_kminus1_filename (it must contain a row at Timestamp = start_timestamp NMPC)
        3. y_previous_predictions_df (it must contain a row at Timestamp = start_timestamp NMPC)
        4. y_star_nom_all_df (only for plot)
        5. z_star_all_df (only for plot)
        6. v_star_all_df (only for plot)

_Note_: to run OFFLINE-TUBE NMPC formulation, just enforce "ancillary.z_star[:,0] - ancillary.z_star0" as constraint in the ancillary problem (see below, commented lines)
_Note_: to run CLASSICAL NMPC formulation:
        1. comment everything related to the "ancillary" instantiation. 
        2. in the "nominal_online" problem, comment everyting related to the "zstar0", especially "nominal_online.set_init_state(z_star0.values[0])"
        3. rename "nominal_online" as "classical"
        4. add classical.set_init_state(y_df_ekf[(y_df_ekf['Timestamp']==classical.integrator_parameters['start_timestamp'])].iloc[0,2:11].to_numpy(dtype=float))

_Note_: shall I use df.loc[:,col_name] = value to override col_name in a df becouse it returns a warning?
'''
# In[1]: IMPORT STANDARD LIBRARIES
import numpy as np
import logging
import os
import pandas as pd
from datetime import datetime, timedelta
import ast

#------------------------------------------------------------------------------------------------------------#
# In[2]: IMPORT CUSTOM LIBRARIES
from general_utils.read_file import*
from general_utils.save_df import*
from general_utils.create_dataframe import create_dataframe
from AM2HNgb_xi import*
from model_read_file import read_all
from system_AM2HNgb_xi import*
from output_AM2HNgb_xi import*
from forward_euler_AM2HNgb_xi import*
from discrete_system_AM2HNgb_xi import*
from general_utils.sample_and_create_df import sample_df
from general_utils.convert_u_to_pwm import*
# Functions from the 'NMPC' library neeeded for the main only
from general_utils.process_parameters import*
from EKF import*
from ekf_AM2HNgb_xi import*
from Parameter_estimation import*
from general_utils.sample_and_create_df import*
from transform_matrix_to_df import*
from NMPC_controller import*
from general_utils.FlexiblePlotter import*

#------------------------------------------------------------------------------------------------------------#
# In[3]: DECLARE THE MAIN FUNCTION (ALLOW TO CALL IT FROM OTHER SCRIPTS) AND SETUP LOGGER
def main(modelname):
    '''MAIN FUNCTION TO RUN THE TUBE-BASED NMPC ON RASPBERRY PI
    INPUTS: modelname (str): The name of the model to use (e.g. '1' or '2').
    OUTPUTS: Saves the control action and logs to specified files.
    '''
    # DECLARE SIMULATION META-OPTIONS
    save_out = True # If True, saves intermediate results externally. Needed to be true for proper closed-loop operations. However, if the code is modified accordingly, it can give flexibility to what is saved and what is not of interest.
    # modelname = '1' #or '2' # If the desire is to run the present code stand alone, delete main function and unindent the code below. Then uncomment this line and the lines at the bottom of the code.
    directory = os.getcwd() # Modify if this code from another directory with respect to where the input and output files are stored and willing to be saved
    testname = 'Test_closedloop' # If there is a specific test subfolder name of the <directory>, specify it here (e.g., '/Test_closedloop')
    #------------------------------------------------------------------------------------------------------------#
    # DECLARE LOGGER
    # Configure logging
    log_directory = os.path.join(directory, testname, "logs")
    os.makedirs(log_directory, exist_ok=True)
    log_file = os.path.join(log_directory, 'main_nmpc_logging.log')
    logger = logging.getLogger()
    # Clean up any existing handlers for this logger
    for handler in logger.handlers[:]:
        logger.removeHandler(handler)
        handler.close()
    # File handler for writing logs to a file
    file_handler = logging.FileHandler(log_file, mode='a')
    file_handler.setLevel(logging.INFO)
    file_formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
    file_handler.setFormatter(file_formatter)
    # Add handlers to the logger
    logger.addHandler(file_handler)
    #if not logger.handlers:  # Avoid duplicate handlers
    # # Console handler for displaying logs in the terminal
    # console_handler = logging.StreamHandler()
    # console_handler.setLevel(logging.INFO)
    # console_formatter = logging.Formatter('%(levelname)s - %(message)s')
    # console_handler.setFormatter(console_formatter)
    # # Add handlers to the logger
    # logger.addHandler(console_handler)
    logger.setLevel(logging.INFO)
    logger.info('############################################### START TUBE MAIN ###############################################')

    #------------------------------------------------------------------------------------------------------------#
    # READ RUN INDEX
    index_control_run_df = pd.read_csv(os.path.join(directory, testname, "index_control_run.txt"))
    index_control_run = index_control_run_df['IndexRun'].values[0]
    logger.info('###############################################\n')
    logger.info(f"Index control run: {index_control_run}")
    #------------------------------------------------------------------------------------------------------------#
    # READ ALL INPUTS
    maize, cowslurry, tomatosauce, d_flow, u_flow, init, theta = read_all(os.path.join(directory, testname), 
                                [], # Add path to PATH if needed
                                os.path.join(directory, testname, 'Input\\'), # Input path (where the most of the filenames below are located)
                                ['parameter_update.json', 'integrator.json', 'parameter.json', 'parameter_bounds.json'], # Integrator setting and process parameters for the model
                                [f'Manual_flowrates.csv', f'{modelname}_u_actual_ANCILLARY.csv', # Needed to set current loads in the model
                                'MaizeSilage_comp_raw.csv','CowSlurry_comp_raw.csv', 'TomatoSauce_comp_raw.csv', # Compositions of the co-feedstocks
                                f'{modelname}_Openloop_model.csv', f'{modelname}_EKF.csv', f'{modelname}_EKF_COV.csv', # Needed to initialize model and EKF
                                f'{modelname}_Openloop_model_ref.csv'],)
    integrator_parameters = theta['integrator_parameters']
    integrator_parameters['model']['start_timestamp'] = datetime.strptime(integrator_parameters['model']['start_timestamp'], "%Y-%m-%d %H:%M:%S")
    integrator_parameters['model']['end_timestamp'] = datetime.strptime(integrator_parameters['model']['end_timestamp'], "%Y-%m-%d %H:%M:%S")
    integrator_parameters['ekf']['start_timestamp'] = datetime.strptime(integrator_parameters['ekf']['start_timestamp'], "%Y-%m-%d %H:%M:%S")
    integrator_parameters['ekf']['end_timestamp'] = datetime.strptime(integrator_parameters['ekf']['end_timestamp'], "%Y-%m-%d %H:%M:%S")
    integrator_parameters['nmpc']['start_timestamp'] = datetime.strptime(integrator_parameters['nmpc']['start_timestamp'], "%Y-%m-%d %H:%M:%S")
    integrator_parameters['nmpc']['end_timestamp'] = datetime.strptime(integrator_parameters['nmpc']['end_timestamp'], "%Y-%m-%d %H:%M:%S")

    #------------------------------------------------------------------------------------------------------------#
    # In[4]: INSTANTIATE MODEL REFERENCE (REFERENCE PARAMETERS ARE GIVEN AS ATTRIBUTE)
    model = AM2HNgb_xi(   
        modelname=f'{modelname}_model', 
        system_func = system_AM2HNgb_xi,
        output_func = output_AM2HNgb_xi,
        discrete_system_func=discrete_system_AM2HNgb_xi,
        input_func = u,
        param_input_func = param_u,
        process_parameters = theta['params_ref'].copy(),
        integrator = forward_euler_AM2HNgb_xi,
        other_parameters = theta['parameters_other_ref'].copy(),
        init_state=[],
        logger = logger,
        integrator_parameters=integrator_parameters['model'].copy()
    )
    #------------------------------------------------------------------------------------------------------------#
    # PREPARE INPUTS FOR THE MODEL (default for the other blocks if not provided in instantiation)
    model_inputs = model.prepare_inputs(maize, cowslurry, tomatosauce, d_flow, u_flow, init['x_model_ref'], 
                                        mult_factors = [86400/300, 86400/300, 1], fill_option=['first_row','first_row','first_row'],
                                        save_to_csv=True, log=True, sic=[0,125,0], output_path = os.path.join(directory, testname, "Input"))
    #------------------------------------------------------------------------------------------------------------#
    # RUN MODEL WITH REFERENCE PARAMETERS (ESTIMATED FROM OFFLINE IDENTIFICATION)
    results = model.integrate(
                uncertain_param_names = [])
    outputs = model.evaluate_output(
                results,
                [])
    all = model.result_grouping(results,outputs)
    # SAMPLE AND SAVE OUT STATES AND OUTPUTS
    y_df_ref = sample_and_create_df(all, all.keys(), model.time_points, 1, model.integrator_parameters['start_timestamp'], False)
    if save_out == True:
        # In theory the duplicates are overridden, but evenutally save y_df_ref.iloc[1:,:11]!!
        # To save all and not only the states:
        save_df_with_check(y_df_ref.drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname, "Input", f'{modelname}_Openloop_model_ref.csv'), log=True)
    
    #------------------------------------------------------------------------------------------------------------#
    # In[5]: PREPARE INPUTS AND RUN OPENLOOP MODEL (WITH THE LAST ESTIMATED PARAMETERS FROM ONLINE IDENTIFICATION)
    # Modify parameters and initial states
    model.process_parameters = theta['params'].copy()
    model.other_parameters = theta['parameters_other'].copy()
    model_inputs = model.prepare_inputs(maize, cowslurry, tomatosauce, d_flow, u_flow, init['x_model'], 
                                        mult_factors = [86400/300, 86400/300, 1], fill_option=['first_row','first_row','first_row'],
                                        save_to_csv=True, log=True, sic=[0,125,0], output_path = os.path.join(directory, testname, "Input"))
    # RUN MODEL WITH LAST-AVAILABLE PARAMETERS
    results = model.integrate(
                uncertain_param_names = [])
    outputs = model.evaluate_output(
                results,
                [])
    all = model.result_grouping(results,outputs)
    # SAMPLE AND SAVE OUT STATES AND OUTPUTS
    y_df = sample_and_create_df(all, all.keys(), model.time_points, 1, model.integrator_parameters['start_timestamp'], False)
    if save_out == True:
        # In theory the duplicates are overridden, but evenutally save y_df.iloc[1:,:11]!!
        # To save all and not only the states:
        save_df_with_check(y_df.drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname, "Input", f'{modelname}_Openloop_model.csv'), log=True)

    #------------------------------------------------------------------------------------------------------------#
    # In[6]: READ PAST MEASUREMENTS DATA
    y_df_on_path = os.path.join(directory, testname, f'{modelname}_y_df_on.txt')
    y_df_off_path = os.path.join(directory, testname, f'{modelname}_y_df_off.txt')
    y_df_data_on = pd.read_csv(y_df_on_path, sep=',', header=0, parse_dates=['Timestamp'])
    y_df_data_off = pd.read_csv(y_df_off_path, sep=',', header=0, parse_dates=['Timestamp'])
    y_df_data_off_f = sample_df(y_df_data_off, 24*7) # Note: every 7 days 1 point of S2 measurement, to be realistic in parameter online identification!!
    # Compute some quantities 
    y_df_data_on['ch4_rate'] = y_df_data_on['gas_rate']*y_df_data_on['xch4_interp']
    y_df_data_on['co2_rate'] = y_df_data_on['gas_rate']*y_df_data_on['xco2_interp']
    y_df_data_on['co2ch4'] = y_df_data_on['xco2_interp']/y_df_data_on['xch4_interp']
    #------------------------------------------------------------------------------------------------------------#
    # READ THE LIST OF SIMULATED-MEASURED VARIABLES' COUPLES
    var_couples_list = read_and_split_txt(os.path.join(directory, testname, "var_couples_list.txt"), log=True)
    var_couples_list = [ast.literal_eval(f"[{item}]") for item in var_couples_list[0]]
    
    #------------------------------------------------------------------------------------------------------------#
    # In[7]: INSTANTIATE EKF
    kR = np.array([67.17/7, 0.00061*3, 0.00048*3]) #integrator_parameters['ekf']['kR'] # Note: nominally, use the ones specified in integrator_parameters['ekf']['kR'].
    kQ = np.array(integrator_parameters['ekf']['kQ'])
    R = np.eye(len(var_couples_list[:3]))*kR#*list(y_trajectory[0]) # Uncomment if you want to use the first row of the measurements as scaling factor for R, otherwise it is just kR (tuning parameter)
    Q = np.eye(model.init_state.shape[0])*kQ*(np.abs(np.array([2.03466,1.61127,0.3769,1.58906/10,1.32249/2.5,0.07159/10,11.8705*2.5,8.53828,0.54181])))
    ekf = EKF(model=model,
            integrator=ekf_AM2HNgb_xi,
            modelname=f'{modelname}_ekf',
            var_couples_list=var_couples_list[:3],
            Q = Q,
            R = R,
            integrator_parameters=integrator_parameters['ekf'].copy(),
            logger=logger)
    #------------------------------------------------------------------------------------------------------------#
    # PREPARE EKF INPUTS
    ekf_inputs = ekf.prepare_inputs(maize, cowslurry, tomatosauce, d_flow, u_flow, init['x_ekf'], init['P_ekf'],
                    mult_factors = [86400/300, 86400/300, 1], fill_option=['first_row','first_row','first_row'],
                        save_to_csv=True, log=True, sic=[0,125,0], output_path = os.path.join(directory, testname, "Input"))
    #------------------------------------------------------------------------------------------------------------#
    # RUN EKF
    y_traj = ekf.prepare_online_data(y_df_data_on,
                            log = True)
    results_ekf = ekf.integrate_ekf(y_traj,
                log = True)
    outputs_ekf = ekf.model.evaluate_output(
                results_ekf[0],
                [])
    all_ekf = ekf.model.result_grouping(results_ekf[0],outputs_ekf)
    # SAMPLE AND SAVE OUT STATES AND OUTPUTS
    y_df_ekf = sample_and_create_df(all_ekf, all_ekf.keys(), ekf.time_points, 1, ekf.integrator_parameters['start_timestamp'], False)
    if save_out == True:
        # In theory the duplicates are overridden, but evenutally save y_df_ekf.iloc[1:,:11]!!
        # To save all and not only the states:
        save_df_with_check(y_df_ekf.drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname,"Input", f'{modelname}_EKF.csv'), log=True)
    #------------------------------------------------------------------------------------------------------------#
    # EXTRACT DIAGONAL OF COVARIANCE MATRIX P 
    P_ekf = dict(zip(['P' + s for s in model.get_state_dict().keys()], np.array([results_ekf[1][:,i,i] for i in range(len(model.init_state))])))
    P_df_ekf = sample_and_create_df(P_ekf, P_ekf.keys(), ekf.time_points, 1, ekf.integrator_parameters['start_timestamp'], False)
    # SAVE OUT DIAGONAL COVARIANCE
    if save_out == True:
        save_df_with_check(P_df_ekf, os.path.join(directory, testname,"Input",f'{modelname}_EKF_diagCOV.csv'), log=True)
    if y_df_ekf['Timestamp'].all() != P_df_ekf['Timestamp'].all():
        raise ValueError('Timestamps do not match')
    combined_df = pd.concat([y_df_ekf, P_df_ekf.copy().drop(columns=['Timestamp'])], axis=1)
    # Compute states confidence intervals (95%CI)
    all_df_ekf = ekf.compute_confidence_intervals(combined_df, 1, 0.95)
    #------------------------------------------------------------------------------------------------------------#
    # SAMPLE AND SAVE OUT WHOLE COVARIANCE MATRIX P
    all_ekf_P = transform_matrix_to_df(results_ekf[1], list(model.get_state_dict().keys()))
    #all_ekf_P.insert(0, 'Timestamp', y_df_ekf['Timestamp'])
    all_ekf_P = convert_data_structure(all_ekf_P)
    y_df_Pekf = sample_and_create_df(all_ekf_P, all_ekf_P.keys(), ekf.time_points, 1, ekf.integrator_parameters['start_timestamp'], False)
    if save_out == True:
        # In theory the duplicates are overridden, but evenutally save y_df_Pekf.iloc[1:]!!
        save_df_with_check(y_df_Pekf, os.path.join(directory, testname, "Input", f'{modelname}_EKF_COV.csv'), log=True)
    
    #------------------------------------------------------------------------------------------------------------#
    # In[8]: COMPUTE OUTPUT ONE-STEP-AHEAD PREDICTION ERROR
    # One-step actully means one-control-interval ahead (multiple integration steps ahead, depending on the control interval and integration time step).
    # Compares real measurements with the predictions of the open-loop model initialized in the EKF estimates.
    y_previous_predictions_df = read_csv_file(os.path.join(directory, testname, "Output", f'{modelname}_Openloop_model_osaprediction.csv'), log=True)
    y_previous_predictions_df['qC_gb'] = y_previous_predictions_df['xC_gb']*y_previous_predictions_df['qTot']
    y_df_previouspred = y_previous_predictions_df.merge(y_df_data_on, on='Timestamp')
    y_df_previouspred = y_df_previouspred[(y_df_previouspred['Timestamp']==integrator_parameters['model']['end_timestamp'])]
    y_df_previouspred['qM_gb'] = y_df_previouspred['xM_gb']*y_df_previouspred['qTot']
    y_df_previouspred['co2ch4_gb'] = y_df_previouspred['xC_gb']/y_df_previouspred['xM_gb']
    for value in var_couples_list[:3]:
        logger.info(f'One-control interval ahead prediction errors are: {value[0]} : {y_df_previouspred[value[0]] - y_df_previouspred[value[1]]}')
    # Compute percentage error and store in a dictionary
    percentage_errors = {}
    for value in var_couples_list[:3]:
        error = (y_df_previouspred[value[0]] - y_df_previouspred[value[1]]) / y_df_previouspred[value[1]] * 100
        percentage_errors[value[0]] = error.values[0]
    # Create a DataFrame from the dictionary
    percentage_errors_df = pd.DataFrame([percentage_errors])
    percentage_errors_df.insert(0, 'Timestamp', pd.Timestamp(integrator_parameters['model']['end_timestamp']))
    # Save the DataFrame to a CSV file
    save_df_with_check(percentage_errors_df, os.path.join(directory, testname, "Output", f'{modelname}_OSAPredictionErrors_wrtdata.csv'), log=True)

    # Compares the pure EKF estimates with the predictions of the open-loop model initialized in the EKF estimates.
    # Gives a quantification about how much the EKF is 'deviating' from the open-loop model predictions
    # Compute absolute percentage error between all common columns in y_df_previouspred and y_df_ekf at the last timestamp (end of the control interval)
    common_columns = set(y_df_previouspred.columns).intersection(set(y_df_ekf.columns))
    timestamp_filter = y_df_previouspred['Timestamp'] == integrator_parameters['model']['end_timestamp']
    if timestamp_filter.any():
        row_previouspred = y_df_previouspred[timestamp_filter].iloc[0]
        row_ekf = y_df_ekf[y_df_ekf['Timestamp'] == integrator_parameters['model']['end_timestamp']].iloc[0]
        absolute_percentage_errors = {}
        for column in common_columns:
            if column != 'Timestamp':
                value_previouspred = row_previouspred[column]
                value_ekf = row_ekf[column]
                if value_ekf != 0:
                    error = abs((value_previouspred - value_ekf) / value_ekf) * 100
                    absolute_percentage_errors[column] = error
                else:
                    error = abs((value_previouspred - (value_ekf+1e-6)) / (value_ekf+1e-6)) * 100
                    absolute_percentage_errors[column] = error
        # Convert the errors to a DataFrame for better visualization
        absolute_percentage_errors_df = pd.DataFrame([absolute_percentage_errors])
        absolute_percentage_errors_df.insert(0, 'Timestamp', pd.Timestamp(integrator_parameters['model']['end_timestamp']))
        # Save the DataFrame to a CSV file
        save_df_with_check(absolute_percentage_errors_df, os.path.join(directory, testname, "Output", f'{modelname}_OSAPredictionErrors_wrtEKF.csv'), log=True)
    else:
        logger.warning(f"No matching row found for Timestamp == {integrator_parameters['model']['end_timestamp']}")
    
    #------------------------------------------------------------------------------------------------------------# 
    # In[9]: CHECK EKF ERROR ON VFA (S2) LAST POINT
    # Declare the feasibility bounds for process parameters
    bound_dict = theta['params_bounds'].copy()
    LB = [bound[0] * scale[1] for bound,scale in zip(bound_dict.values(), theta['params_ref'].values())]
    UB = [bound[1] * scale[1] for bound,scale in zip(bound_dict.values(), theta['params_ref'].values())]
    bnds = list([(lb,ub) for lb,ub in zip(LB,UB)])
    bound_scaled_dict = dict(zip(bound_dict.keys(),bnds))
    #------------------------------------------------------------------------------------------------------------#
    # Compute error on the last S2
    y_df_data_off_cut = y_df_data_off_f[(y_df_data_off_f['Timestamp'])<=model.integrator_parameters['end_timestamp']]       # Note: cut data to the end of the EKF
    x_df_ekf_cut = pd.concat([init['x_ekf'][(init['x_ekf']['Timestamp']<=ekf.integrator_parameters['start_timestamp'])], 
                            y_df_ekf.iloc[1:,:11].drop(columns='Xh')], axis=0, ignore_index=True)                           # Note: concatenate past states with EKF results at this run (avoid duplicate of initial state)
    last_ekf_s2_meas_estimate = x_df_ekf_cut.merge(y_df_data_off_cut, on='Timestamp').iloc[-1]['S2']                        
    error_ekf = np.abs(last_ekf_s2_meas_estimate - y_df_data_off_cut['TVFA'].iloc[-1])
    logger.info(f'EKF error on S2 on last point ({y_df_data_off_cut["Timestamp"].iloc[-1]}) is {error_ekf} mmol/L/day')
    #------------------------------------------------------------------------------------------------------------#
    # # CHECK IF PARAMETER UPDATE IS NEEDED (IF ERROR IS TOO HIGH AND THERE IS DATA IN THE WINDOW) 
    # # If first run, I want to update parameters also if there is no data in the window. In addition, delay in S2 measurement must be considered (TO BE IMPLEMENTED)
    # parameter_update_S2_threshold = 20 # mmol/L (to be tuned, maybe better to put a percentage error threshold with respect to the S2 value)
    # if error_ekf > parameter_update_S2_threshold and y_df_data_off_cut['Timestamp'].between(model.integrator_parameters['end_timestamp'] - timedelta(hours=integrator_parameters['nmpc']['control_interval']), model.integrator_parameters['end_timestamp']).any(): #METTERE PERCENTUALE
    #     #------------------------------------------------------------------------------------------------------------#
    #     # Update integrator parameters of the model (over the window of past data for parameter update)
    #     model.integrator_parameters['start_timestamp'] = model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght'])
    #     model_inputs = model.prepare_inputs(maize, cowslurry, tomatosauce, d_flow, u_flow, init['x_model'], sic=[0, 125,0],
    #                                         mult_factors = [86400/300, 86400/300, 1], fill_option=['first_row','first_row','first_row'], log=True, output_path = os.path.join(directory, testname, "Input"))
    #     logger.info(f'Parameter update now {model.integrator_parameters["end_timestamp"]} back to {model.integrator_parameters["start_timestamp"]}')
    #     #------------------------------------------------------------------------------------------------------------#
    #     # DO PARAMETER UPDATE
    #     estimator = Parameter_estimation(model, 
    #                                     f'{modelname}_estimator',
    #                                     var_couples_list,
    #                                     uncertain_param_names = ['k5','k6','K2I','Ks2','mumax2'],
    #                                     nominal_process_parameters=theta['params'].copy(),
    #                                     nominal_other_parameters=theta['parameters_other'].copy(),
    #                                     var_couples_list_offline = [var_couples_list[3]],
    #                                     y_df_data_on = y_df_data_on,
    #                                     y_df_data_off = y_df_data_off_f,
    #                                     bound_scaled_dict = bound_scaled_dict,
    #                                     logger=logger,
    #                                     )                                                           # Note: start and end timestamps if not given are the ones of 'model'
    #     results_update = parameter_update(estimator,
    #                     f'{modelname}_update',
    #                     model.u_values,
    #                     model.param_u_values,
    #                     model.delta_t,
    #                     model.time_points, 
    #                     'nelder-mead',
    #                     800,
    #                     False,
    #                     bool_run_sensitivity=True,
    #                     deltap=[0.3, -0.3],
    #                     sens_method='rr',
    #                     initcond=False,
    #                     sens_type='local',
    #                     rho_rate=30,
    #                     rho_start=10,
    #                     freezone=1,
    #                     active_refpenalty=1e-6,
    #                     logs = [True, True, False, True])
    #     #------------------------------------------------------------------------------------------------------------#
    #     # Check if I want to save this update or discard it. If save, save it to csv
    #     # Consideration 31.01.2025: I shall 'accept' the updated parameters only if EKF run with the updated parameters is better than the one with the nominal parameters, on the same window of data
    #     # Update json file with new estimated process parameters
    #     update_parameters('parameters_update.json', results_update[0], 'process_parameters', log=True)
    #     update_parameters('parameters_update.json', results_update[1], 'other_parameters', log=True)
    #     # Update the current model process parameters with new estimated ones
    #     model.process_parameters = results_update[0]
    #     model.other_parameters = results_update[1]
    #     # Save estimator status to csv
    #     if save_out == True:
    #         for name in estimator.uncertain_param_names:
    #             estim_status = pd.DataFrame([estimator.status[name]])
    #             save_df_with_check(estim_status, os.path.join(directory, testname, "Input", f'{modelname}_Openloop_model_{name}.csv'), log=True)
    #     # Save prior and posterior states to CSV (+ override 'Openloop_model.csv' with the posterior)
    #     y_df_prior = results_update[2]
    #     y_df_opt = results_update[3]
    #     if save_out == True:
    #         # In theory the duplicates are overridden, but evenutally save df.iloc[1:,:11]!!
    #         save_df_with_check(y_df_prior.drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname, "Input", f'{modelname}_Openloop_model_prior_run{index_control_run}.csv'), log=True)
    #         save_df_with_check(y_df_opt.drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname, "Input", f'{modelname}_Openloop_model_posterior_run{index_control_run}.csv'), log=True)
    #         save_df_with_check(y_df_opt.drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname, "Input",f'{modelname}_Openloop_model.csv'), log=True)
    #     #------------------------------------------------------------------------------------------------------------# 
    #     # Restart EKF in posterior (IN REALITY THERE IS DELAY FOR COMPUTATION! See notes 23.10.2024)
    #     last_post_s2_meas_estimate = y_df_opt.merge(y_df_data_off_cut, on='Timestamp').iloc[-1]['S2']
    #     error_post = np.abs(last_post_s2_meas_estimate - y_df_data_off_cut['TVFA'].iloc[-1])
    #     if error_ekf > error_post: # I reinit in posterior ONLY if the EKF estimate at the last available offline meas. point is worst than the one provided by the openloop posterior model
    #         # Consideratoion 31.01.2025: I shall 'accept' the updated parameters only if EKF run with the updated parameters is better than the one with the nominal parameters, on the same window of data
    #         if y_df_opt.loc[y_df_opt.index[-1], 'Timestamp'] == y_df_ekf.loc[y_df_ekf.index[-1], 'Timestamp']:  # Note: if 'y_df_opt' in locals()
    #             logger.info('Future initial state of EKF ovverride with posterior. Before it was: \n')
    #             logger.info(f'{y_df_ekf.loc[y_df_ekf.index[-1]]}')
    #             logger.info('Now it is: \n')
    #             logger.info(f'{y_df_opt.loc[y_df_opt.index[-1]]}')
    #             y_df_ekf.loc[y_df_ekf.index[-1]] = y_df_opt.loc[y_df_opt.index[-1]]                             # Note: override y_df_ekf with y_df_post, to initialize NMPC accordingly
    #         else:
    #             raise ValueError('EKF shall be re-initialized, but Timestamps do not match with posterior dataframe')
    #         # To save all and not only the states
    #         save_df_with_check(y_df_opt.iloc[-1:,:].drop(columns=['Xh','S2_comp_adm1']), os.path.join(directory, testname, "Input", f'{modelname}_EKF.csv'), log=True)
    #         logger.info('Future initial state of EKF ovverride with posterior, saved to csv')
    #         # Open point: if no EKF reinit in posterior...what to do? Model and EKF are not in sync anymore? Run in 'parallel'?
    #     #------------------------------------------------------------------------------------------------------------# 
    #     # Move back to the start integrator_parameters[model]...in this case shall be equal to the ones of EKF
    #     model.integrator_parameters['start_timestamp'] = ekf.integrator_parameters['start_timestamp']
    #     model.integrator_parameters['end_timestamp'] = ekf.integrator_parameters['end_timestamp']
    
    #------------------------------------------------------------------------------------------------------------# 
    # In[10]: LOAD SETPOINTS (FROM THE OFFLINE OPTIMIZATION MADE WITH AGRI-ACODM)
    y_ref = read_and_split_txt(os.path.join(directory, testname, "Input", 'Output_setpoints.txt'), log=True)
    y_ref = create_dataframe(['Qch4_ref', 'co2ch4_ref'], [pd.Series(y_ref[1]), pd.Series(y_ref[2])], datetime(2025,4,30,0,0))
    x_ref = read_and_split_txt(os.path.join(directory, testname, "Input", 'State_setpoints.txt'), log=True)
    x_ref = create_dataframe(['Xh_ref', 'X1_ref', 'X2_ref', 'S1_ref', 'S2_ref', 'csi_ref','xMgb_ref'], [pd.Series(x_ref[1]), 
                                                                                                        pd.Series(x_ref[2]),
                                                                                                        pd.Series(x_ref[3]),
                                                                                                        pd.Series(x_ref[4]),
                                                                                                        pd.Series(x_ref[5]),
                                                                                                        pd.Series(x_ref[6]),
                                                                                                        pd.Series(x_ref[7]),
                                                                                                    ], datetime(2025,4,30,10,0))
    # Convert units from L/h to mmol/L/day
    y_ref['Qch4_ref'] = udm_gas_conversion(y_ref['Qch4_ref'], model.other_parameters['T'], model.other_parameters['Pt'], model.other_parameters['V'], 'Lh')
    y_ref['S2_ref'] = x_ref['S2_ref']
    
    # #------------------------------------------------------------------------------------------------------------# 
    # In[11]: INSTANTIATE NMPC(s)
    integrator_parameters['nmpc']['slackness'] = integrator_parameters['nmpc']['slackness'] == 'true'
    end_nmpc = integrator_parameters['nmpc']['start_timestamp'] + timedelta(hours=integrator_parameters['nmpc']['Hp']*integrator_parameters['nmpc']['control_interval']) 
    if integrator_parameters['nmpc']['end_timestamp'] != end_nmpc:
        logging.warning('There is an inconsistency between the end timestamp of the NMPC and the computed one from the integrator parameters!')
    ancillary = NMPC(model=model,
                modelname=f'{modelname}',
                control_interval= integrator_parameters['nmpc']['control_interval'],                                                            # Note: number of hours between one MPC run and the next
                N = integrator_parameters['nmpc']['Hp']*integrator_parameters['nmpc']['control_interval']/24,                                   # Note: days of prediction horizon...must be coherent with start and end timestamps of integrator_parameters!!
                Nu = integrator_parameters['nmpc']['Hc']*integrator_parameters['nmpc']['control_interval']/24,                                  # Note: days of control horizon
                ref_df=y_ref,
                var_couples_list=var_couples_list[:3],
                integrator_parameters=integrator_parameters['nmpc'].copy(),
                state_scales=np.array([2.03466286,1.61127425,0.37690481,1.58906061,1.32249497,0.07159448,11.87046368,2.53827702,0.54180507]),   # Note: model.init_state
                output_scales=[26.1232906165, 0.8456824872417199, 11.87046368],
                input_scale=0.01,
                problemtype=integrator_parameters['nmpc']['problemtype'],                                                                       # Note: should be 'ANCILLARY' in online-tube NMPC formulation
                slackness= integrator_parameters['nmpc']['slackness'],                                                                          # Note: should be 'False' in online-tube NMPC formulation
                formulation= integrator_parameters['nmpc']['formulation'],                                                                      # Note: should be 'ancillaryonline' in online-tube NMPC formulation
                logger=logger)
    #------------------------------------------------------------------------------------------------------------# 
    if integrator_parameters['nmpc']['formulation'] == 'ancillaryonline':
        nominal_online = NMPC(model=model,
                modelname=f'{modelname}',
                control_interval= integrator_parameters['nmpc']['control_interval'],                                                            # Note: number of hours between one MPC run and the next
                N = integrator_parameters['nmpc']['Hp']*integrator_parameters['nmpc']['control_interval']/24,                                   # Note: days of prediction horizon...must be coherent with start and end timestamps of integrator_parameters!!
                Nu = integrator_parameters['nmpc']['Hc']*integrator_parameters['nmpc']['control_interval']/24,                                  # Note: days of control horizon
                ref_df=y_ref,
                var_couples_list=var_couples_list[:3],
                integrator_parameters=integrator_parameters['nmpc'].copy(),
                state_scales=np.array([2.03466286,1.61127425,0.37690481,1.58906061,1.32249497,0.07159448,11.87046368,2.53827702,0.54180507]),   # Note: model.init_state
                output_scales= [26.1232906165, 0.8456824872417199, 11.87046368],
                input_scale=0.01,
                problemtype='NOMINAL',                                                                                                          # Note: should be 'NOMINAL' in online-tube NMPC formulation
                slackness=True,                                                                                                                 # Note: should be 'True' in online-tube NMPC formulation
                formulation='nominal',                                                                                                          # Note: should be 'nominal' in online-tube NMPC formulation
                logger=logger)
    
    #------------------------------------------------------------------------------------------------------------# 
    # In[12]: DECLARE INDEX OF FILES TO BE READ FOR INITIAL GUESSES OF NMPC 
    index_run_read = 0 if index_control_run == 0 else index_control_run - 1 
    #------------------------------------------------------------------------------------------------------------# 
    # PREPARE NMPC INPUTS
    nmpc_inputs = ancillary.prepare_inputs(maize, cowslurry, tomatosauce, d_flow, u_flow,
                        mult_factors = [86400/300, 86400/300, 1], fill_option=['first_row','first_row','first_row'], 
                        log=True, sic=[0,125,0], output_path = os.path.join(directory, testname, "Input"))
    if integrator_parameters['nmpc']['formulation'] == 'ancillaryonline':
        nmpc_inputs_nominal = nominal_online.prepare_inputs(maize, cowslurry, tomatosauce, d_flow, u_flow,
                        mult_factors = [86400/300, 86400/300, 1], fill_option=['first_row','first_row','first_row'], 
                        log=True, sic=[0,125,0], output_path = os.path.join(directory, testname, "Input"))
    #------------------------------------------------------------------------------------------------------------# 
    # RUN NMPC (ONLINE-NOMINAL)
    if integrator_parameters['nmpc']['formulation'] == 'ancillaryonline':
        logging.info('Solve nominal..................................................................................................')
        nominal_online.model.other_parameters['observer'] = 'alternative_gb'
        nominal_online.model.set_casadi_functions()
        cost, prediction_timestamps = nominal_online.set_ext_cost_function(log = True)
        # Generate a new list with 1-hour intervals (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
        filled_timestamps = pd.date_range(start=min(prediction_timestamps), end=max(prediction_timestamps), freq="1H")
        filled_timestamps_list = filled_timestamps.tolist()
        nominal_online.add_constraints()
        #------------------------------------------------------------------------------------------------------------# 
        # Initialize nominal state trajectories with the last available value of z*
        z_star_kminus1_filename = os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_run{index_run_read}_z_star_ANCILLARY.csv')
        z_star_kminus1_df = read_csv_file(z_star_kminus1_filename, log=True)
        z_star0_df = z_star_kminus1_df[(z_star_kminus1_df['Timestamp']>=nominal_online.integrator_parameters['start_timestamp']) & (z_star_kminus1_df['Timestamp']<=nominal_online.integrator_parameters['start_timestamp'])]
        z_star0 = z_star0_df.drop(columns='Timestamp')
        nominal_online.set_init_state(z_star0.values[0])
        #------------------------------------------------------------------------------------------------------------# 
        # Warm-starting of the optimization with the previous run's optimal trajectories (if they exist)
        Uk_guess_filename = os.path.join(directory, testname, "Input", f'{nominal_online.modelname}_run{index_run_read}_u_star_{nominal_online.problemtype}.csv')
        Xk_guess_filename = os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_run{index_run_read}_x_star_{nominal_online.problemtype}.csv')
        if os.path.exists(Uk_guess_filename) and os.path.exists(Xk_guess_filename):
            Uk_guess_df = read_csv_file(Uk_guess_filename, log=True)
            # Ensure the 'Timestamp' column is in datetime format (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
            Uk_guess_df['Timestamp'] = pd.to_datetime(Uk_guess_df['Timestamp'])
            Uk_guess_df.set_index('Timestamp', inplace=True)
            # Create a new index with 1-hour intervals
            Uk_guess_df = Uk_guess_df.reindex(pd.date_range(start=Uk_guess_df.index.min(), end=Uk_guess_df.index.max(), freq='1H'))
            # Interpolate the missing values
            Uk_guess_df = Uk_guess_df.fillna(method='pad')
            Uk_guess_df.reset_index(inplace=True)
            Uk_guess_df.rename(columns={'index': 'Timestamp'}, inplace=True)
            Uk_guess_df.reset_index(inplace=True)
            # Remove first row
            Uk_guess_df = Uk_guess_df.drop(index=0) #                           # Note: remove first row (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
            Uk_guess_df = Uk_guess_df[(Uk_guess_df['Timestamp'] >= nominal_online.integrator_parameters['start_timestamp']) &
                                    (Uk_guess_df['Timestamp'] <= nominal_online.integrator_parameters['end_timestamp']-timedelta(hours=integrator_parameters['nmpc']['control_interval']))] #filter out start/end_timestamp
            # Sample Uk_guess_df every control_interval hours within the start and end range
            Uk_guess_df = sample_df(Uk_guess_df, integrator_parameters['nmpc']['control_interval'])
            Uk_guess_df.drop(columns=['index'], inplace=True, errors='ignore')
            #----------------------------------------------------------#
            Xk_guess_df = read_csv_file(Xk_guess_filename, log=True)
            # Remove first row
            Xk_guess_df.reset_index(inplace=True)
            Xk_guess_df = Xk_guess_df.drop(index=0)                             # Note: remove first row (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
            # Filter out start/end_timestamp and then sample every control_interval hours
            Xk_guess_df = Xk_guess_df[(Xk_guess_df['Timestamp'] >= nominal_online.integrator_parameters['start_timestamp']) &
                                        (Xk_guess_df['Timestamp'] <= nominal_online.integrator_parameters['end_timestamp']-timedelta(hours=integrator_parameters['nmpc']['control_interval']))] #filter out start/end_timestamp
            Xk_guess_df = sample_df(Xk_guess_df, integrator_parameters['nmpc']['control_interval'])
            Xk_guess_df.drop(columns=['index'], inplace=True, errors='ignore')
            #----------------------------------------------------------#
            Xk_guess = Xk_guess_df.drop(columns='Timestamp')                    # Note: check last timestamp of the guesses = start_timestamp of nmpc
            Uk_guess = Uk_guess_df.drop(columns='Timestamp')                    # Note: check last timestamp of the guesses = start_timestamp of nmpc
            Xk_guess = Xk_guess.values[:,:].T
            Uk_guess = Uk_guess.values[:,:].T
        else:
            Xk_guess = np.tile(nominal_online.init_state, (1, nominal_online.nintervals))
            Uk_guess = np.tile(100, (1, nominal_online.nintervals))
        # Set initial guesses
        initial_guesses = {
                        'eps1': np.tile(0, (9, 1)).T,                           # Note: save and read from csv TO BE IMPLEMENTED for slack variables
                        'eps2': np.tile(0, (9, 1)).T}                           # Note: save and read from csv TO BE IMPLEMENTED for slack variables
        initial_guesses['Xk'] = Xk_guess
        initial_guesses['Uk'] = Uk_guess
        nominal_online.set_init_decvar(initial_guesses)
        #------------------------------------------------------------------------------------------------------------# 
        # SET OPTI PARAMETERS: weights and optimization variables' bounds
        opti_param_values ={ 'weights_y': np.array(integrator_parameters['nmpc']['opti_parameters']['weights_y']),
                            'weights_du': integrator_parameters['nmpc']['opti_parameters']['weights_du'],
                            'bnds_u': np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_u']),
                            'bnds_du':  np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_du']),
                            'bnds_x': np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_x'])}
        if nominal_online.slackness:
            opti_param_values['weights_eps1'] = 100*integrator_parameters['nmpc']['opti_parameters']['weights_eps1']
            opti_param_values['weights_eps2'] = 100*integrator_parameters['nmpc']['opti_parameters']['weights_eps2']
        nominal_online.set_parameters(opti_param_values)
        #------------------------------------------------------------------------------------------------------------# 
        nominal_online.setup_solver(maxiter=300, print_level=5)
        #------------------------------------------------------------------------------------------------------------# 
        decvar_star, y_star, costs_star, constraints_star = nominal_online.run_optimization(log=True)
        #------------------------------------------------------------------------------------------------------------# 
        # SAVE OUT DECISION VARIABLES VALUES (TO BE USED AS WARM START FOR NEXT RUN)
        # STATES
        x_star = [decvar_star['x_star'][:,0]] #First element is the initial state
        for i in range(len(nominal_online.x_list)):
            if nominal_online.solution is not None:
                x_star.append(nominal_online.solution.value(nominal_online.x_list[i]))
            else:
                x_star.append(nominal_online.opti.debug.value(nominal_online.x_list[i]))
        x_star_df_nom = pd.DataFrame(np.array(x_star), columns=model.get_state_dict().keys())
        x_star_df_nom.insert(0, 'Timestamp', filled_timestamps_list)
        save_df_with_check(x_star_df_nom, os.path.join(directory, testname, "Output",f'{nominal_online.modelname}_z_star_NOMINAL.csv'), log=True)
        save_df_with_check(x_star_df_nom, os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_run{index_control_run}_x_star_{nominal_online.problemtype}.csv'),log=True)
        # CONTROL ACTION
        u_star_df_nom = pd.DataFrame(np.round(np.concatenate((decvar_star['u_star'],decvar_star['u_star'][-1:])),1).T, columns=['Uk'])                      # Note: duplicate last element
        u_star_df_nom.insert(0, 'Timestamp', prediction_timestamps)
        save_df_with_check(u_star_df_nom.iloc[:1], os.path.join(directory, testname, "Input", f'{nominal_online.modelname}_v_star_NOMINAL.csv'), log=True)
        save_df_with_check(u_star_df_nom, os.path.join(directory, testname, "Input", f'{nominal_online.modelname}_run{index_control_run}_u_star_{nominal_online.problemtype}.csv'), log=True)
        # OUTPUTS
        y_star = []
        for i in range(len(nominal_online.y_list)):
            if nominal_online.solution is not None:
                y_star.append(nominal_online.solution.value(nominal_online.y_list[i]))
            else:
                y_star.append(nominal_online.opti.debug.value(nominal_online.y_list[i]))
        y_star_df_nom = pd.DataFrame(np.array(y_star), columns=[value[0] for value in var_couples_list[:3]])
        y_star_df_nom.insert(0, 'Timestamp', filled_timestamps_list)
        y_star_df_nom['qM_gb'] = y_star_df_nom['xM_gb']*y_star_df_nom['qTot']
        y_star_df_nom['co2ch4_gb'] = y_star_df_nom['xC_gb']/y_star_df_nom['xM_gb']
        save_df_with_check(y_star_df_nom, os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_y_star_NOMINAL.csv'), log=True)
        save_df_with_check(y_star_df_nom, os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_run{index_control_run}_y_star_{nominal_online.problemtype}.csv'), log=True)
        # Save general controller status
        nmpc_status = nominal_online.status.copy()
        nmpc_status = pd.DataFrame([nmpc_status])
        save_df_with_check(nmpc_status, os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_nmpc_status_{nominal_online.problemtype}.csv'), log=True)
        # SLACK VARIABLES
        if nominal_online.slackness:
            eps1_star_df_nom = pd.DataFrame([decvar_star['eps1_star']], columns=model.get_state_dict().keys())
            eps2_star_df_nom = pd.DataFrame([decvar_star['eps2_star']], columns=model.get_state_dict().keys())
            eps1_star_df_nom.insert(0, 'Timestamp', filled_timestamps_list[0])
            eps2_star_df_nom.insert(0, 'Timestamp', filled_timestamps_list[0])
            save_df_with_check(eps1_star_df_nom, os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_eps1_star_NOMINAL.csv'), log=True)
            save_df_with_check(eps2_star_df_nom, os.path.join(directory, testname, "Output", f'{nominal_online.modelname}_eps2_star_NOMINAL.csv'), log=True)
    
    #------------------------------------------------------------------------------------------------------------#
    # In[13]: RUN NMPC (ANCILLARY)
    logging.info('Solve ancillary..................................................................................................')
    ancillary.model.other_parameters['observer'] = 'alternative_gb'
    ancillary.model.set_casadi_functions()
    #------------------------------------------------------------------------------------------------------------# 
    # IF NOT NOMINAL, READ z^* and \nu^* FROM FILE (cut between start and end timestamps)
    v_star_nom_filename = os.path.join(directory, testname, "Input", f'{ancillary.modelname}_run{index_control_run}_u_star_NOMINAL.csv')
    z_star_nom_filename = os.path.join(directory, testname, "Output",f'{ancillary.modelname}_run{index_control_run}_x_star_NOMINAL.csv')
    y_star_nom_filename = os.path.join(directory, testname, "Output", f'{ancillary.modelname}_run{index_control_run}_y_star_NOMINAL.csv')
    if os.path.exists(v_star_nom_filename) and os.path.exists(z_star_nom_filename):
        v_star_nom_df = read_csv_file(v_star_nom_filename, log=True)
        z_star_nom_df = read_csv_file(z_star_nom_filename, log=True)
        y_star_nom_df = read_csv_file(y_star_nom_filename, log=True)
        # Filter between start and end-step timestamps (NB end-step becouse it shall be given as it outputs the nominal problem and the shape of decision variables i.e. without values at end)
        v_star_nom_df = v_star_nom_df[(v_star_nom_df['Timestamp']>=ancillary.integrator_parameters['start_timestamp']) & (v_star_nom_df['Timestamp']<=ancillary.integrator_parameters['end_timestamp']-timedelta(hours=ancillary.integrator_parameters['control_interval']))]
        z_star_nom_df = z_star_nom_df[(z_star_nom_df['Timestamp']>=ancillary.integrator_parameters['start_timestamp']) & (z_star_nom_df['Timestamp']<=ancillary.integrator_parameters['end_timestamp']-timedelta(hours=ancillary.integrator_parameters['control_interval']))]
        y_star_nom_df = y_star_nom_df[(y_star_nom_df['Timestamp']>=ancillary.integrator_parameters['start_timestamp']) & (y_star_nom_df['Timestamp']<=ancillary.integrator_parameters['end_timestamp']-timedelta(hours=ancillary.integrator_parameters['control_interval']))]
        # Sample every control_interval hours within the start and end range
        v_star_nom_df = sample_df(v_star_nom_df, ancillary.integrator_parameters['control_interval'])
        v_star_nom_df.drop(columns=['index'], inplace=True, errors='ignore')
        y_star_nom_df = sample_df(y_star_nom_df, ancillary.integrator_parameters['control_interval'])
        y_star_nom_df.drop(columns=['index'], inplace=True, errors='ignore')
        z_star_nom_df = sample_df(z_star_nom_df, ancillary.integrator_parameters['control_interval'])
        z_star_nom_df.drop(columns=['index'], inplace=True, errors='ignore')
        # Check last timestamp of the guesses = start_timestamp of nmpc
        v_star_nom = v_star_nom_df.drop(columns='Timestamp')
        z_star_nom = z_star_nom_df.drop(columns='Timestamp')
        y_star_nom = y_star_nom_df.drop(columns=['Timestamp', 'qM_gb','co2ch4_gb'])
        z_star_nom = (z_star_nom.values).T
        y_star_nom = (y_star_nom.values).T
        v_star_nom = v_star_nom.values
        v_star_nom = np.reshape(v_star_nom,(1,max(v_star_nom.shape)))
    else:
        z_star_nom = np.tile(model.init_state, (1, ancillary.nintervals))
        y_star_nom = np.tile(np.array([49.8,0.54181,0.45819]), (1, ancillary.nintervals))
        v_star_nom = np.tile(100, (1, ancillary.nintervals))
    # Create dataframes for z_star_nom, y_star_nom, and v_star_nom (JUST TO PLOT VERY INITIAL GUESS)
    z_star_nom_df = pd.DataFrame(z_star_nom.T, columns=model.get_state_dict().keys())
    z_star_nom_df.insert(0, 'Timestamp', pd.date_range(start=ancillary.integrator_parameters['start_timestamp'], end=ancillary.integrator_parameters['end_timestamp'], freq=f'{ancillary.integrator_parameters["control_interval"]}H')[:-1])  # Note: exclude the last timestamp to match the shape
    y_star_nom_df = pd.DataFrame(y_star_nom.T, columns=[value[0] for value in var_couples_list[:3]])
    y_star_nom_df.insert(0, 'Timestamp', pd.date_range(start=ancillary.integrator_parameters['start_timestamp'], end=ancillary.integrator_parameters['end_timestamp'], freq=f'{ancillary.integrator_parameters["control_interval"]}H')[:-1])  # Note: exclude the last timestamp to match the shape
    y_star_nom_df['qM_gb'] = y_star_nom_df['xM_gb']*y_star_nom_df['qTot']
    y_star_nom_df['co2ch4_gb'] = y_star_nom_df['xC_gb']/y_star_nom_df['xM_gb']
    v_star_nom_df = pd.DataFrame(v_star_nom.T, columns=['Uk'])
    v_star_nom_df.insert(0, 'Timestamp', pd.date_range(start=ancillary.integrator_parameters['start_timestamp'], end=ancillary.integrator_parameters['end_timestamp'], freq=f'{ancillary.integrator_parameters["control_interval"]}H')[:-1])  # Note: exclude the last timestamp to match the shape
    ancillary.z_star = z_star_nom
    ancillary.v_star = v_star_nom
    #------------------------------------------------------------------------------------------------------------# 
    cost, prediction_timestamps = ancillary.set_ext_cost_function(log = True)
    # Generate a new list with 1-hour intervals (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
    filled_timestamps = pd.date_range(start=min(prediction_timestamps), end=max(prediction_timestamps), freq="1H")
    filled_timestamps_list = filled_timestamps.tolist()
    ancillary.add_constraints()
    #------------------------------------------------------------------------------------------------------------# 
    # # !!!! UNCOMMENT THE FOLLOWING 3 LINES TO ALLOW FOR OFFLINE-TUBE NMPC FORMULATION (KEEPING THE SAME CODE STRUCTURE AS THE ONLINE-TUBE) !!!!
    # init_zstar_anc = ancillary.z_star[:,0] - ancillary.z_star0
    # offline_tube = init_zstar_anc == 0
    # ancillary.additional_ext_contraint(constraint_name="offline_mode", constraint=offline_tube)
    #------------------------------------------------------------------------------------------------------------# 
    # SET INIT_MPC TO THE LAST EKF STATE ESTIMATE (OR TO THE LAST POSTERIOR IF EKF REINIT IN POSTERIOR). ACTUAL FEEDBACK.
    ancillary.set_init_state(y_df_ekf[(y_df_ekf['Timestamp']==ancillary.integrator_parameters['start_timestamp'])].iloc[0,2:11].to_numpy(dtype=float))
    #------------------------------------------------------------------------------------------------------------# 
    # READ OLD OPTIMAL RESULTS FOR WARM START-UP I.E. SET AS INITIAL GUESSES
    Uk_guess_filename = os.path.join(directory, testname, "Input", f'{ancillary.modelname}_run{index_run_read}_u_star_{ancillary.problemtype}.csv')
    Xk_guess_filename = os.path.join(directory, testname, "Output",f'{ancillary.modelname}_run{index_run_read}_x_star_{ancillary.problemtype}.csv')
    if os.path.exists(Uk_guess_filename) and os.path.exists(Xk_guess_filename):
        Uk_guess_df = read_csv_file(Uk_guess_filename, log=True)
        # Ensure the 'Timestamp' column is in datetime format (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
        Uk_guess_df['Timestamp'] = pd.to_datetime(Uk_guess_df['Timestamp'])
        Uk_guess_df.set_index('Timestamp', inplace=True)
        # Create a new index with 1-hour intervals
        Uk_guess_df = Uk_guess_df.reindex(pd.date_range(start=Uk_guess_df.index.min(), end=Uk_guess_df.index.max(), freq='1H'))
        # Interpolate the missing values
        Uk_guess_df = Uk_guess_df.fillna(method='pad')
        Uk_guess_df.reset_index(inplace=True)
        Uk_guess_df.rename(columns={'index': 'Timestamp'}, inplace=True)
        Uk_guess_df.reset_index(inplace=True)
        # Remove first row
        Uk_guess_df = Uk_guess_df.drop(index=0) #remove first row (added on 25.05.2025 to allow sporadic/asynchronous NMPC runs)
        Uk_guess_df = Uk_guess_df[(Uk_guess_df['Timestamp'] >= nominal_online.integrator_parameters['start_timestamp']) &
                                (Uk_guess_df['Timestamp'] <= nominal_online.integrator_parameters['end_timestamp']-timedelta(hours=integrator_parameters['nmpc']['control_interval']))]
        # Sample Uk_guess_df every control_interval hours within the start and end range
        Uk_guess_df = sample_df(Uk_guess_df, integrator_parameters['nmpc']['control_interval'])
        Uk_guess_df.drop(columns=['index'], inplace=True, errors='ignore')
        #----------------------------------------------------------#
        Xk_guess_df = read_csv_file(Xk_guess_filename, log=True)
        # Remove first row
        Xk_guess_df.reset_index(inplace=True)
        Xk_guess_df = Xk_guess_df.drop(index=0)
        # Filter out start/end_timestamp and then sample every control_interval hours
        Xk_guess_df = Xk_guess_df[(Xk_guess_df['Timestamp'] >= nominal_online.integrator_parameters['start_timestamp']) &
                                    (Xk_guess_df['Timestamp'] <= nominal_online.integrator_parameters['end_timestamp']-timedelta(hours=integrator_parameters['nmpc']['control_interval']))]
        Xk_guess_df = sample_df(Xk_guess_df, integrator_parameters['nmpc']['control_interval'])
        Xk_guess_df.drop(columns=['index'], inplace=True, errors='ignore')
        #----------------------------------------------------------#
        # Check last timestamp of the guesses = start_timestamp of nmpc
        Xk_guess = Xk_guess_df.drop(columns='Timestamp')
        Uk_guess = Uk_guess_df.drop(columns='Timestamp')
        Xk_guess = Xk_guess.values[:,:].T
        Uk_guess = Uk_guess.values[:,:].T
    else:
        Xk_guess = np.tile(ancillary.init_state, (1, ancillary.nintervals))
        Uk_guess = np.tile(100, (1, ancillary.nintervals))
        # Only for plotting purposes at index=0
        Uk_guess_df = pd.DataFrame(Uk_guess.T, columns=['Uk'])
        Uk_guess_df.insert(0, 'Timestamp', pd.date_range(start=ancillary.integrator_parameters['start_timestamp'], end=ancillary.integrator_parameters['end_timestamp'], freq=f'{ancillary.integrator_parameters["control_interval"]}H')[:-1])  # Note: exclude the last timestamp to match the shape
    # Set initial guesses
    initial_guesses = {
                    'eps1': np.tile(0, (9, 1)).T,                           # Note: save and read from csv TO BE IMPLEMENTED for slack variables
                    'eps2': np.tile(0, (9, 1)).T}                           # Note: save and read from csv TO BE IMPLEMENTED for slack variables
    initial_guesses['Xk'] = Xk_guess
    initial_guesses['Uk'] = Uk_guess
    # If online-tube, set also z_star0 as initial guess for z_star trajectory (to be read from csv)
    if ancillary.formulation == 'ancillaryonline':
        initial_guesses.update({'z_star0': ancillary.z_star[:,0]})
    ancillary.set_init_decvar(initial_guesses)
    #------------------------------------------------------------------------------------------------------------# 
    # SET OPTI PARAMETERS: weights and optimization variables' bounds
    opti_param_values ={ 'weights_y': np.array(integrator_parameters['nmpc']['opti_parameters']['weights_y']),
                        'weights_du': integrator_parameters['nmpc']['opti_parameters']['weights_du'],
                        'bnds_u': np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_u']),
                        'bnds_du':  np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_du']),
                        'bnds_x': np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_x_anc'])}
    opti_param_values['weights_yN'] = integrator_parameters['nmpc']['opti_parameters']['weights_yN']
    opti_param_values['weights_xstar'] = integrator_parameters['nmpc']['opti_parameters']['weights_xstar']
    opti_param_values['weights_ustar'] = integrator_parameters['nmpc']['opti_parameters']['weights_ustar']
    if ancillary.formulation == 'ancillaryonline':
        opti_param_values['bnds_z'] = np.array(integrator_parameters['nmpc']['opti_parameters']['bnds_z'])
    ancillary.set_parameters(opti_param_values)
    #------------------------------------------------------------------------------------------------------------# 
    ancillary.setup_solver(maxiter=300, print_level = 5)
    #------------------------------------------------------------------------------------------------------------# 
    decvar_star, y_star, costs_star, constraints_star = ancillary.run_optimization(log=True)
    #------------------------------------------------------------------------------------------------------------# 
    # SAVE OUT DECISION VARIABLES VALUES (TO BE USED AS WARM START FOR NEXT RUN)
    # STATES
    x_star = [decvar_star['x_star'][:,0]]                                                                                       # Note: first element is the initial state
    for i in range(len(ancillary.x_list)):
        if ancillary.solution is not None:
            x_star.append(ancillary.solution.value(ancillary.x_list[i]))
        else:
            x_star.append(ancillary.opti.debug.value(ancillary.x_list[i]))
    x_star_df = pd.DataFrame(np.array(x_star), columns=model.get_state_dict().keys())
    x_star_df.insert(0, 'Timestamp', filled_timestamps_list)
    save_df_with_check(x_star_df, os.path.join(directory, testname, "Output", f'{ancillary.modelname}_run{index_control_run}_x_star_{ancillary.problemtype}.csv'), log=True)
    # CONTROL ACTION
    u_star_df = pd.DataFrame(np.round(np.concatenate((decvar_star['u_star'],decvar_star['u_star'][-1:])),1).T, columns=['Uk'])  # Note: duplicate last element
    u_star_df.insert(0, 'Timestamp', prediction_timestamps)
    save_df_with_check(u_star_df.iloc[:1], os.path.join(directory, testname, "Input", f'{ancillary.modelname}_u_actual_{ancillary.problemtype}.csv'), log=True)
    save_df_with_check(u_star_df, os.path.join(directory, testname, "Input", f'{ancillary.modelname}_run{index_control_run}_u_star_{ancillary.problemtype}.csv'), log=True)
    # OUTPUTS
    y_star = []
    for i in range(len(ancillary.y_list)):
        if ancillary.solution is not None:
            y_star.append(ancillary.solution.value(ancillary.y_list[i]))
        else:
            y_star.append(ancillary.opti.debug.value(ancillary.y_list[i]))
    y_star_df = pd.DataFrame(np.array(y_star), columns=[value[0] for value in var_couples_list[:3]])
    y_star_df.insert(0, 'Timestamp', filled_timestamps_list)
    y_star_df['qM_gb'] = y_star_df['xM_gb']*y_star_df['qTot']
    y_star_df['co2ch4_gb'] = y_star_df['xC_gb']/y_star_df['xM_gb']
    save_df_with_check(y_star_df, os.path.join(directory, testname, "Output", f'{ancillary.modelname}_run{index_control_run}_y_star_{ancillary.problemtype}.csv'), log=True)
    # Save general controller status
    nmpc_status = ancillary.status.copy()
    nmpc_status = pd.DataFrame([nmpc_status])
    save_df_with_check(nmpc_status, os.path.join(directory, testname, "Output", f'{ancillary.modelname}_nmpc_status_{ancillary.problemtype}.csv'), log=True)
    # OSA-PREDICTION
    y_df_future = pd.merge(x_star_df, y_star_df, on='Timestamp', how='outer', suffixes=('', '_drop'))                           # Note: add suffixes to avoid duplication
    y_df_future = y_df_future.loc[:, ~y_df_future.columns.str.endswith('_drop')]                                                # Note: drop duplicate columns
    save_df_with_check(y_df_future.iloc[1:], os.path.join(directory, testname, "Output", f'{modelname}_Openloop_model_osaprediction.csv'), log=True)
    y_df_future['qC_gb'] = y_df_future['xC_gb']*y_df_future['qTot']                                                             # Note: only for plotting
    # IF ANCILLARY ONLINE, SAVE ALSO THE UPDATED z_0^* AND z^* VALUES
    if ancillary.formulation == 'ancillaryonline':
        # RE-OPTIMIZED STATES
        z_star_0 = np.array(decvar_star['z_star0'])
        z_star_0_df = pd.DataFrame(np.reshape(z_star_0, (1,z_star_0.shape[0])), columns=model.get_state_dict().keys())
        z_star_0_df.insert(0, 'Timestamp', ancillary.integrator_parameters['start_timestamp'])
        save_df_with_check(z_star_0_df, os.path.join(directory, testname, "Output", f'{ancillary.modelname}_run{index_control_run}_z_star0.csv'), log=True)
        z_star = decvar_star['z_star']
        z_star_df = pd.DataFrame(np.array(z_star), columns=model.get_state_dict().keys())
        z_star_df.insert(0, 'Timestamp', filled_timestamps_list)
        save_df_with_check(z_star_df, os.path.join(directory, testname, "Output", f'{ancillary.modelname}_run{index_control_run}_z_star_{ancillary.problemtype}.csv'), log=True)
        # RE-OPTIMIZED CONTROL ACTION: you fool! The control action is the same as the one in the nominal problem
        # RE-OPTIMIZED OUTPUTS
        y_star_anc = []
        for i in range(len(ancillary.y_star_list)):
            if ancillary.solution is not None:
                y_star_anc.append(ancillary.solution.value(ancillary.y_star_list[i]))
            else:
                y_star_anc.append(ancillary.opti.debug.value(ancillary.y_star_list[i]))
        y_star_anc_df = pd.DataFrame(np.array(y_star_anc), columns=[value[0] for value in var_couples_list[:3]])
        y_star_anc_df.insert(0, 'Timestamp', filled_timestamps_list)
        y_star_anc_df['qM_gb'] = y_star_anc_df['xM_gb']*y_star_anc_df['qTot']
        y_star_anc_df['co2ch4_gb'] = y_star_anc_df['xC_gb']/y_star_anc_df['xM_gb']
        save_df_with_check(y_star_anc_df, os.path.join(directory, testname, "Output", f'{ancillary.modelname}_run{index_control_run}_y_star_anc_{ancillary.problemtype}.csv'), log=True)
    
    #------------------------------------------------------------------------------------------------------------# 
    # In[14]: PLOT VS DATA
    # Create an instance of FlexiblePlotter
    plotter = FlexiblePlotter(default_figsize = (12, 12), logger=logger)
    # Calculate the x_limits
    end_time = pd.to_datetime(y_df_data_on['Timestamp'].iloc[-1]) + timedelta(hours=integrator_parameters['nmpc']['control_interval']*integrator_parameters['nmpc']['Hp'])
    start_time = model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght'])
    date_string = f'{start_time.strftime("%d%m%H")}-{end_time.strftime("%d%m%H")}'
    # Read the data from the CSV files
    df_ekf_past = init['x_ekf'][init['x_ekf']['Timestamp']<=ekf.integrator_parameters['start_timestamp']]
    df_model_past = init['x_model'][init['x_model']['Timestamp']<=model.integrator_parameters['start_timestamp']]
    df_Pekf_past = init['P_ekf'][init['P_ekf']['Timestamp']<=ekf.integrator_parameters['start_timestamp']]
    df_model_ref_past = init['x_model_ref'][init['x_model_ref']['Timestamp']<=model.integrator_parameters['start_timestamp']]
    df_data_past_off = y_df_data_off[y_df_data_off['Timestamp']<=model.integrator_parameters['end_timestamp']]
    df_data_past = y_df_data_on[y_df_data_on['Timestamp']<=model.integrator_parameters['end_timestamp']]
    z_star_all_df = read_csv_file(os.path.join(directory, testname, 'Output', f'{nominal_online.modelname}_z_star_NOMINAL.csv'))
    z_star_past_df = z_star_all_df[z_star_all_df['Timestamp']<=model.integrator_parameters['end_timestamp']]
    y_star_nom_all_df = read_csv_file(os.path.join(directory, testname, 'Output', f'{nominal_online.modelname}_y_star_NOMINAL.csv'))
    y_star_nom_past_df = y_star_nom_all_df[y_star_nom_all_df['Timestamp']<=model.integrator_parameters['end_timestamp']]
    df_u_past = u_flow[u_flow['Timestamp']<=model.integrator_parameters['start_timestamp']]
    v_star_all_df = read_csv_file(os.path.join(directory, testname, "Input", f'{nominal_online.modelname}_v_star_NOMINAL.csv'))
    v_star_past_df = v_star_all_df[v_star_all_df['Timestamp']<=model.integrator_parameters['start_timestamp']]
    # Extract a subset of columns from df_Pekf_past
    columns_to_extract = ['Timestamp', 'Xh1-Xh1', 'Xh2-Xh2', 'Xh3-Xh3', 'X1-X1', 'X2-X2', 'S1-S1', 'S2-S2', 'csi-csi', 'xM_gb-xM_gb']
    df_Pekf_past_subset = df_Pekf_past[columns_to_extract]
    # Rename the columns with names specified in a list
    new_column_names = ['Timestamp', 'PXh1', 'PXh2', 'PXh3', 'PX1', 'PX2', 'PS1', 'PS2', 'Pcsi', 'PxM_gb']
    df_Pekf_past_subset.columns = new_column_names
    #------------------------------------------------------------------------------------------------------------# 
    # INCREMENTAL MULTIPLE PLOT OF MAIN QUANTITIES
    data = {"data_on": y_df_data_on,
            "data_off": y_df_data_off,
            "AM2HN": y_df,
            "EKF": y_df_ekf,
            "oldpred":y_previous_predictions_df,
            "newpred":y_df_future,
            "df_ekf_past": df_ekf_past,
            "df_model_past": df_model_past,
            }
    # Define the variables to be plotted
    variable_groups = [
        [("data_on", "gas_rate"), ("AM2HN", "qTot"), ("EKF", "qTot"), ("oldpred", "qTot"), ("newpred", "qTot"), ("df_ekf_past", "qTot"), ("df_model_past", "qTot")],  # First subplot
        [("data_on", "ch4_rate"), ("AM2HN", "qM_gb"), ("EKF", "qM_gb"), ("oldpred", "qM_gb"), ("newpred", "qM_gb"), ("df_ekf_past", "qM_gb"), ("df_model_past", "qM_gb")],  # Second subplot
        [("data_on", "co2_rate"), ("AM2HN", "qC_gb"), ("EKF", "qC_gb"), ("oldpred", "qC_gb"), ("newpred", "qC_gb"), ("df_ekf_past", "qC_gb"), ("df_model_past", "qC_gb")],  # Third subplot
        [("data_on", "co2ch4"), ("AM2HN", "co2ch4_gb"), ("EKF", "co2ch4_gb"), ("oldpred", "co2ch4_gb"), ("newpred", "co2ch4_gb"), ("df_ekf_past", "co2ch4_gb"), ("df_model_past", "co2ch4_gb")],  # Second subplot
        [("data_off", "TVFA"), ("AM2HN", "S2"), ("EKF", "S2"), ("oldpred", "S2"), ("newpred", "S2"), ("df_ekf_past", "S2"), ("df_model_past", "S2")],  # Second subplot
    ]
    # Calculate the x_limits
    x_limits = [(start_time, end_time) for _ in range(len(variable_groups))]
    # Use the FlexiblePlotter to create the plot
    plotter.plot_grouped_variables_across_dfs(
        data=data,
        variable_groups=variable_groups,
        x_axis_format="datetime",
        x_limits=x_limits,
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,                                                            # Note: hours of tick interval x
        grid=True,
        show_plot=False,
        y_labels=[None for _ in range(len(variable_groups))],                   # Note: no y-axis labels for subplots
        plot_types=[
            ['line', 'line', 'line', 'scatter', 'scatter','line','line'],       # First subplot
            ['line', 'line', 'line', 'scatter', 'scatter','line','line'],       # Second subplot 
            ['line', 'line', 'line', 'scatter', 'scatter','line','line'],       # Third subplot 
            ['line', 'line', 'line', 'scatter', 'scatter','line','line'],       # Fourth subplot 
            ['scatter', 'line', 'line', 'scatter', 'scatter','line','line']     # Fifth subplot (data_off as scatter to distinguish from data_on)
        ],
        legend_param={"loc": "upper left"},
        colors={"df_ekf_past":["green" for _ in range(len(variable_groups[0]))], 
                "df_model_past":["orange" for _ in range(len(variable_groups[0]))]},
        alphas={"df_ekf_past": [0.5, 0.5, 0.5, 0.5], "df_model_past": [0.5, 0.5, 0.5, 0.5], "oldpred":[0.5, 0.5, 0.5, 0.5], "newpred":[0.5, 0.5, 0.5, 0.5]},
        y_limits=[(20, 140), (10, 70), (10, 70), (0.5, 1.5), 
                (0, 40)
                ],
    )
    output_path = os.path.join(directory, testname, 'plots')
    # Create the output directory if it does not exist
    os.makedirs(os.path.dirname(output_path), exist_ok=True)
    # Save the figure
    plt.savefig(os.path.join(output_path, f"{modelname}_{date_string}_model&ekf_ALT_withpred.png"), bbox_inches='tight', dpi=300)
    # Log the output path
    logger.info(f"Plot saved to {output_path}")
    #------------------------------------------------------------------------------------------------------------#
    # PLOT ALL STATES
    # Ensure numeric columns are cast to float
    numeric_columns = ['Xh1', 'Xh2', 'Xh3', 'X1', 'X2', 'S1', 'S2', 'csi','xM_gb']
    for col in numeric_columns:
        if col in y_df.columns:
            y_df[col] = pd.to_numeric(y_df[col], errors='coerce')
        if col in y_df_ekf.columns:
            y_df_ekf[col] = pd.to_numeric(y_df_ekf[col], errors='coerce')
    # Create an instance of FlexiblePlotter
    plotter = FlexiblePlotter(default_figsize = (12, int(2*len(numeric_columns))), logger=logger)
    # Define the variables to be plotted
    data = {"model":y_df,
            "ekf":y_df_ekf,
            "oldpred":y_previous_predictions_df,
            "newpred":y_df_future,
            "df_ekf_past": df_ekf_past,
            "df_model_past": df_model_past,
            }
    # Dynamically create variable_groups based on keys in data and numeric_columns
    variable_groups = []
    for col in numeric_columns:
        group = []
        for key in data.keys():
            if col in data[key].columns:
                group.append((key, col))
        if group:  # Only add non-empty groups
            variable_groups.append(group)
    # Calculate the x_limits: use the ones computed above
    x_limits = [(start_time, end_time) for _ in range(len(variable_groups))]
    # Use the FlexiblePlotter to create the plot
    plotter.plot_grouped_variables_across_dfs(
        data=data,
        variable_groups=variable_groups,
        x_axis_format="datetime",
        x_limits=x_limits,
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,
        grid=True,
        show_plot=False,
        y_labels=[None for _ in range(len(variable_groups))], 
        plot_types=[
        ['line', 'line', 'scatter', 'scatter','line','line'] for _ in range(len(variable_groups))], 
        legend_param={"loc": "upper left"},
        colors={"df_ekf_past":["green" for _ in range(len(variable_groups[0]))],
                "ekf":["green" for _ in range(len(variable_groups[0]))], 
                "model":["orange" for _ in range(len(variable_groups[0]))],
                "df_model_past":["orange" for _ in range(len(variable_groups[0]))]},
        alphas={"df_ekf_past": [0.5 for _ in range(len(variable_groups[0]))], "df_model_past": [0.5 for _ in range(len(variable_groups[0]))], 
        "oldpred":[0.5 for _ in range(len(variable_groups[0]))], "newpred":[0.5 for _ in range(len(variable_groups[0]))]},
    )
    # Save the figure
    plt.savefig(os.path.join(output_path, f"{modelname}_{date_string}_{'&'.join(data.keys())}_states_ALT_withpred.png"), bbox_inches='tight', dpi=300)
    # Log the output path
    logger.info(f"Plot saved to {output_path}")
    #------------------------------------------------------------------------------------------------------------# 
    # PLOT THE DIAGONAL OF THE COVARIANCE MATRIX P (scaled)
    # Concatenate P_df_ekf and df_Pekf_past_subset by rows
    P_df_combined = pd.concat([df_Pekf_past_subset, P_df_ekf], axis=0, ignore_index=True)
    P_ekf_scales = np.array([2.03466286,1.61127425,0.37690481,1.58906061,1.32249497,0.07159448,11.87046368,2.53827702,0.54180507])
    if 'Timestamp' in P_df_combined.columns:
        P_df_combined.drop(columns=['Timestamp'], inplace=True)
    P_ekf_scaled = P_df_combined.copy()
    for i, col in enumerate(P_ekf_scaled.columns):
        P_ekf_scaled[col] = np.sqrt(P_ekf_scaled[col]) / np.abs(P_ekf_scales[i])
    if 'Timestamp' not in P_ekf_scaled.columns:
        P_ekf_scaled.insert(0, 'Timestamp', pd.concat([df_Pekf_past_subset['Timestamp'], P_df_ekf['Timestamp']], ignore_index=True))
    numeric_columns = ['PXh1', 'PXh2', 'PXh3', 'PX1', 'PX2', 'PS1', 'PS2', 'Pcsi','PxM_gb']
    # Create an instance of FlexiblePlotter
    plotter = FlexiblePlotter(default_figsize = (12, 4), logger=logger)
    # Define the variables to be plotted
    data = {"P_ekf_scaled": P_ekf_scaled,}
    # Dynamically create variable_groups as a dictionary
    variable_groups = {"P_ekf_scaled": numeric_columns}
    # Calculate the x_limits
    x_limits = [(start_time, end_time) for _ in range(len(variable_groups))]
    # Dynamically generate a colors dictionary with lists of colors
    keys = list(data.keys())
    colormap = plt.cm.get_cmap('tab10', len(keys))
    #colors = {"P_ekf_scaled": [colormap(i)] for i, key in enumerate(keys)}
    # Use the FlexiblePlotter to create the plot
    plotter.plot_multiple_variables_across_dfs_with_labels(
        data=data,
        variables=variable_groups,
        x_axis_format="datetime",
        x_limit=x_limits[0],
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,
        grid=True,
        show_plot=False,
        y_limit=(0, 0.3),
        #colors=colors,
    )
    # Save the figure
    plt.savefig(os.path.join(output_path, f"{modelname}_{date_string}_{'&'.join(data.keys())}_cov_ALT.png"), bbox_inches='tight', dpi=300)
    # Log the output path
    logger.info(f"Plot saved to {output_path}")

    #------------------------------------------------------------------------------------------------------------# 
    # PLOT VFA
    if 'y_df_prior' in locals() and 'y_df_opt' in locals():
        prior = y_df_prior
        posterior = y_df_opt
    # Define the data and variables for plotting
    data = {
        'sim_ref_past': df_model_ref_past,
        'sim_past': df_model_past,
        'ekf_past': df_ekf_past,
        'data': df_data_past_off,
        'sim_ref': y_df_ref,
        'sim': y_df,
        'ekf': all_df_ekf,
        'previous_pred': y_df_previouspred,
        'z_star_past': z_star_past_df
    }
    variables = [
        [('sim_ref_past', 'S2'), ('sim_past', 'S2'), ('ekf_past', 'S2'), ('data', 'VFA (mol-basis)'), ('previous_pred', 'S2'), ('z_star_past', 'S2')],
        [('sim_ref', 'S2')],
        [('sim', 'S2')],
        [('ekf', 'S2'), ('ekf', 'lb_S2'), ('ekf', 'ub_S2')]
    ]
    if 'y_df_prior' in locals() and 'y_df_opt' in locals():
        data.update({'prior': prior, 'posterior': posterior})
        variables.append([('prior', 'S2')])
        variables.append([('posterior', 'S2')])
    data.update({'prediction': x_star_df, 
                'z_star_nom': z_star_nom_df})
    variables.append([('prediction', 'S2'),
                    ('z_star_nom', 'S2')])
    if ancillary.formulation == 'ancillaryonline':
        data.update({'z_star_anc': z_star_df})
        variables.append([('z_star_anc', 'S2')])
    # Instantiate plotter
    plotter = FlexiblePlotter(default_figsize = (10, 4), logger=logger)
    # Generate plots
    plotter.generate_plots_with_incremental_variables(
        data=data,
        variables=variables,
        output_folder = os.path.join(directory, testname, 'plots'),
        run_index=index_control_run,
        title = 'S2',
        y_label=r"Total Volatile Fatty Acids (mmol L$^{-1}$)",
        y_limit=(0,40),
        y_ticks = 5,
        x_limit=(model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght']), ancillary.integrator_parameters['end_timestamp']),
        colors= {'sim_ref_past': ['grey'], 'sim_past': ['grey'], 'ekf_past': ['grey'], 'data': ['red'], 
                'sim_ref':['grey'],'sim':['black'], 'ekf': ['blue',"lightblue", "lightblue"], 'prediction': ['black'], 'z_star_nom': ['grey'], 'z_star_anc': ['grey'], 'z_star_past': ['pink']}, #[plt.cm.tab10(np.linspace(0, 1, 3))[i] for i in range(3)]
        alphas={'sim_ref_past': [0.5], 'sim_past': [0.5], 'ekf_past': [0.5], 'data': [1], 
                'sim_ref':[1],'sim':[1], 'ekf': [1,0.8,0.8], 'prediction': [0.8], 'prior': [0.2], 'posterior': [0.5], 'z_star_nom': [0.2], 'z_star_anc': [0.8], 'z_star_past': [1]},
        linestyles={'sim_ref_past': ['--'], 'sim_past': ['--'], 'ekf_past': ['-'], 'sim_ref':['--'],'sim':['--'], 'ekf': ['-','--', '--'], 'z_star_past': ['--']},
        plot_types={'sim_ref_past': ['line'], 'sim_past': ['line'], 'ekf_past': ['line'], 'data': ['scatter'], 
                    'sim_ref':['line'],'sim':['line'], 'ekf': ['line'], 'prediction': ['line'], 'prior': ['line'], 
                    'posterior': ['line'], 'previous_pred':['scatter'], 'z_star_nom': ['scatter'], 'z_star_anc': ['line'], 'z_star_past': ['line']},
        x_axis_format="datetime",
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,
        figsize=(10, 4),
        grid=True,
        vertical_lines=[model.integrator_parameters['end_timestamp']],
        horizontal_lines=[opti_param_values['bnds_x'][6][0], opti_param_values['bnds_x'][6][1],
                          opti_param_values['bnds_z'][6][0], opti_param_values['bnds_z'][6][1]],
        legend_param = {
            "loc": "upper left",                                                # Note: position the legend outside the plot
            "bbox_to_anchor": (1.05, 1),                                        # Note: place the legend outside the axes
            "borderaxespad": 0.0,                                               # Note: reduce padding between the axes and the legend
            "frameon": False                                                    # Note: remove the legend box frame
        },
        layout_param = {
            "pad": 1.0,                                                         # Note: adjust padding to make space for the legend
            "rect": [0, 0, 0.85, 1]                                             # Note: shrink the plot area to fit the legend
    }
    )
    #------------------------------------------------------------------------------------------------------------# 
    # PLOT THE OTHER STATES
    bounds_x = np.concatenate((opti_param_values['bnds_x'][0:6], opti_param_values['bnds_x'][7:]), axis=0)
    bounds_z = np.concatenate((opti_param_values['bnds_z'][0:6], opti_param_values['bnds_z'][7:]), axis=0)
    for name, yaxislimits, bnds_x, bnds_z in zip(['Xh1', 'Xh2', 'Xh3', 'X1', 'X2', 'S1', 'csi', 'xM_gb'],[[0,6],[0,3],[0,1],[0,6],[0,6],[0,0.2],[-15,10],[0.4,0.6]],
                                                bounds_x, bounds_z):
        variables = [
            [('sim_ref_past', name), ('sim_past', name), ('ekf_past', name), ('previous_pred', name), ('z_star_past', name)],
            [('sim_ref', name)],
            [('sim', name)],
            [('ekf', name), ('ekf', f'lb_{name}'), ('ekf', f'ub_{name}')]
        ]
        if 'y_df_prior' in locals() and 'y_df_opt' in locals():
            data.update({'prior': prior, 'posterior': posterior})
            variables.append([('prior', name)])
            variables.append([('posterior', name)])
        # Define the data and variables for plotting
        data.update({'prediction': x_star_df, 
                    'z_star_nom': z_star_nom_df})
        variables.append([('prediction', name),
                        ('z_star_nom', name)])
        if ancillary.formulation == 'ancillaryonline':
            data.update({'z_star_anc': z_star_df})
            variables.append([('z_star_anc', name)])
        # Instantiate plotter
        plotter = FlexiblePlotter(default_figsize = (10, 4), logger=logger)
        # Generate plots
        plotter.generate_plots_with_incremental_variables(
            data=data,
            variables=variables,
            output_folder = os.path.join(directory, testname, 'plots'),
            run_index=index_control_run,
            title = name,
            y_label=f"{name} default U.o.M.",
            y_limit=(yaxislimits[0],yaxislimits[1]),
            y_ticks = min(10, max(2, int((yaxislimits[1] - yaxislimits[0]) // ((yaxislimits[1] - yaxislimits[0]) / 5 if (yaxislimits[1] - yaxislimits[0]) > 0 else 1)))),
            x_limit=(model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght']), ancillary.integrator_parameters['end_timestamp']),
            colors= {'sim_ref_past': ['grey'], 'sim_past': ['grey'], 'ekf_past': ['grey'], 'data': ['red'], 
                    'sim_ref':['grey'],'sim':['black'], 'ekf': ['blue',"lightblue", "lightblue"], 'prediction': ['black'], 'z_star_nom': ['grey'], 'z_star_anc': ['grey'], 'z_star_past': ['pink']}, #[plt.cm.tab10(np.linspace(0, 1, 3))[i] for i in range(3)]
            alphas={'sim_ref_past': [0.5], 'sim_past': [0.5], 'ekf_past': [0.5], 'data': [1], 
                    'sim_ref':[1],'sim':[1], 'ekf': [1,0.8, 0.8], 'prediction': [0.8], 'prior': [0.2], 'posterior': [0.5], 'z_star_nom': [0.2], 'z_star_anc': [0.8], 'z_star_past': [1]},
            linestyles={'sim_ref_past': ['--'], 'sim_past': ['--'], 'ekf_past': ['-'], 'sim_ref':['--'],'sim':['--'], 'ekf': ['-','dashed', 'dashed'], 'z_star_past': ['--']},
            plot_types={'sim_ref_past': ['line'], 'sim_past': ['line'], 'ekf_past': ['line'], 'data': ['scatter'], 
                        'sim_ref':['line'],'sim':['line'], 'ekf': ['line'], 'prediction': ['line'], 'prior': ['line'], 
                        'posterior': ['line'], 'previous_pred':['scatter'], 'z_star_nom': ['scatter'], 'z_star_anc': ['line'], 'z_star_past': ['line']},
            x_axis_format="datetime",
            x_axis_ticks_format="%d-%m %H",
            tick_rotation=90,
            interval=12,
            figsize=(10, 4),
            grid=True,
            vertical_lines=[model.integrator_parameters['end_timestamp']],
            horizontal_lines=[bnds_x[0], bnds_x[1],
                            bnds_z[0], bnds_z[1],],
            legend_param = {
                "loc": "upper left",                                            # Note: position the legend outside the plot
                "bbox_to_anchor": (1.05, 1),                                    # Note: place the legend outside the axes
                "borderaxespad": 0.0,                                           # Note: reduce padding between the axes and the legend
                "frameon": False                                                # Note: remove the legend box frame
            },
            layout_param = {
                "pad": 1.0,                                                     # Note: adjust padding to make space for the legend
                "rect": [0, 0, 0.85, 1]                                         # Note: shrink the plot area to fit the legend
        }
        )
    #------------------------------------------------------------------------------------------------------------#
    # PLOT METHANE FLOWRATE
    if 'y_df_prior' in locals() and 'y_df_opt' in locals():
        prior = y_df_prior
        posterior = y_df_opt
    # Define the data and variables for plotting
    data = {
        'setpoint': y_ref,
        'sim_ref_past': df_model_ref_past,
        'sim_past': df_model_past,
        'ekf_past': df_ekf_past,
        'data': df_data_past,
        'sim_ref': y_df_ref,
        'sim': y_df,
        'ekf': all_df_ekf,
        'previous_pred': y_df_previouspred,
        'y_star_nom_past': y_star_nom_past_df
    }
    variables = [
        [('setpoint', 'Qch4_ref'), ('sim_ref_past', 'qM_gb'), ('sim_past', 'qM_gb'), ('ekf_past', 'qM_gb'), ('data', 'ch4_rate'), ('previous_pred', 'qM_gb'), ('y_star_nom_past', 'qM_gb')],
        [('sim_ref', 'qM_gb')],
        [('sim', 'qM_gb')],
        [('ekf', 'qM_gb')]
    ]
    if 'y_df_prior' in locals() and 'y_df_opt' in locals():
        data.update({'prior': prior, 'posterior': posterior})
        variables.append([('prior', 'qM_gb')])
        variables.append([('posterior', 'qM_gb')])
    data.update({'prediction': y_star_df, 
                'y_star_nom': y_star_nom_df})
    variables.append([('prediction', 'qM_gb'),
                    ('y_star_nom', 'qM_gb')])
    if ancillary.formulation == 'ancillaryonline':
        data.update({'y_star_anc': y_star_anc_df})
        variables.append([('y_star_anc', 'qM_gb')])
    # Generate plots
    plotter.generate_plots_with_incremental_variables(
        data=data,
        variables=variables,
        output_folder = os.path.join(directory, testname, 'plots'),
        run_index=index_control_run,
        title = 'qM_gb',
        y_label=r"CH4 flowrate (mmol L$^{-1}$ day$^{-1}$)",
        y_limit=(10,70),
        y_ticks = 5,
        x_limit=(model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght']), ancillary.integrator_parameters['end_timestamp']),
        colors= {'sim_ref_past': ['grey'], 'sim_past': ['grey'], 'ekf_past': ['grey'], 'data': ['red'], 
                'sim_ref':['grey'],'sim':['black'], 'ekf': ['blue'], 'prediction': ['black'], 'setpoint': ['violet'], 'y_star_nom': ['grey'], 'y_star_anc': ['grey'], 'y_star_nom_past': ['pink']}, #[plt.cm.tab10(np.linspace(0, 1, 3))[i] for i in range(3)]
        alphas={'sim_ref_past': [0.5], 'sim_past': [0.5], 'ekf_past': [0.5], 'data': [0.8], 
                'sim_ref':[1],'sim':[1], 'ekf': [1], 'prediction': [0.8], 'prior': [0.2], 'posterior': [0.5], 'y_star_nom': [0.2], 'y_star_anc': [0.8], 'y_star_nom_past': [1]},
        linestyles={'sim_ref_past': ['--'], 'sim_past': ['--'], 'ekf_past': ['-'], 'sim_ref':['--'],'sim':['--'], 'ekf': ['-'], 'setpoint': ['--'], 'y_star_nom_past': ['--']},
        plot_types={'sim_ref_past': ['line'], 'sim_past': ['line'], 'ekf_past': ['line'], 'data': ['line'], 
                    'sim_ref':['line'],'sim':['line'], 'ekf': ['line'], 'prediction': ['line'], 'prior': ['line'], 'posterior': ['line'], 'previous_pred':['scatter'],
                    'setpoint': ['line'], 'y_star_nom': ['scatter'], 'y_star_anc': ['line'], 'y_star_nom_past': ['line']},
        x_axis_format="datetime",
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,
        figsize=(10, 4),
        grid=True,
        vertical_lines=[model.integrator_parameters['end_timestamp']],
        legend_param = {
            "loc": "upper left",                                                # Note: position the legend outside the plot
            "bbox_to_anchor": (1.05, 1),                                        # Note: place the legend outside the axes
            "borderaxespad": 0.0,                                               # Note: reduce padding between the axes and the legend
            "frameon": False                                                    # Note: remove the legend box frame
        },
        layout_param = {
            "pad": 1.0,                                                         # Note: adjust padding to make space for the legend
            "rect": [0, 0, 0.85, 1]                                             # Note: shrink the plot area to fit the legend
    }
    )
    #------------------------------------------------------------------------------------------------------------# 
    # PLOT CO2CH4 BIOGAS RATIO
    if 'y_df_prior' in locals() and 'y_df_opt' in locals():
        prior = y_df_prior
        posterior = y_df_opt
    # Define the data and variables for plotting
    data = {
        'setpoint': y_ref,
        'sim_ref_past': df_model_ref_past,
        'sim_past': df_model_past,
        'ekf_past': df_ekf_past,
        'data': df_data_past,
        'sim_ref': y_df_ref,
        'sim': y_df,
        'ekf': all_df_ekf,
        'previous_pred': y_df_previouspred,
        'y_star_nom_past': y_star_nom_past_df
    }
    variables = [
        [('setpoint', 'co2ch4_ref'), ('sim_ref_past', 'co2ch4_gb'), ('sim_past', 'co2ch4_gb'), ('ekf_past', 'co2ch4_gb'), ('data', 'co2ch4'), ('previous_pred', 'co2ch4_gb'), ('y_star_nom_past', 'co2ch4_gb')],
        [('sim_ref', 'co2ch4_gb')],
        [('sim', 'co2ch4_gb')],
        [('ekf', 'co2ch4_gb')]
    ]
    if 'y_df_prior' in locals() and 'y_df_opt' in locals():
        data.update({'prior': prior, 'posterior': posterior})
        variables.append([('prior', 'co2ch4_gb')])
        variables.append([('posterior', 'co2ch4_gb')])
    data.update({'prediction': y_star_df, 
                'y_star_nom': y_star_nom_df})
    variables.append([('prediction', 'co2ch4_gb'),
                    ('y_star_nom', 'co2ch4_gb')])
    if ancillary.formulation == 'ancillaryonline':
        data.update({'y_star_anc': y_star_anc_df})
        variables.append([('y_star_anc', 'co2ch4_gb')])
    # Generate plots
    plotter.generate_plots_with_incremental_variables(
        data=data,
        variables=variables,
        output_folder = os.path.join(directory, testname, 'plots'),
        run_index=index_control_run,
        title = 'co2ch4_gb',
        y_label="Biogas composition ratio (-)",
        y_limit=(0.5,1.5),
        y_ticks = 5,
        x_limit=(model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght']), ancillary.integrator_parameters['end_timestamp']),
        colors= {'sim_ref_past': ['grey'], 'sim_past': ['grey'], 'ekf_past': ['grey'], 'data': ['red'], 
                'sim_ref':['grey'],'sim':['black'], 'ekf': ['blue'], 'prediction': ['black'], 'setpoint': ['violet'], 'y_star_nom': ['grey'], 'y_star_anc': ['grey'], 'y_star_nom_past': ['pink']}, #[plt.cm.tab10(np.linspace(0, 1, 3))[i] for i in range(3)]
        alphas={'sim_ref_past': [0.5], 'sim_past': [0.5], 'ekf_past': [0.5], 'data': [1], 
                'sim_ref':[1],'sim':[1], 'ekf': [1], 'prediction': [0.8], 'prior': [0.2], 'posterior': [0.5], 'y_star_nom': [0.2], 'y_star_anc': [0.8], 'y_star_nom_past': [1]},
        linestyles={'sim_ref_past': ['--'], 'sim_past': ['--'], 'ekf_past': ['-'], 'sim_ref':['--'],'sim':['--'], 'ekf': ['-'], 'setpoint': ['--'], 'y_star_nom_past': ['--']},
        plot_types={'sim_ref_past': ['line'], 'sim_past': ['line'], 'ekf_past': ['line'], 'data': ['line'], 
                    'sim_ref':['line'],'sim':['line'], 'ekf': ['line'], 'prediction': ['line'], 'prior': ['line'], 'posterior': ['line'], 'previous_pred':['scatter'],
                    'setpoint': ['line'], 'y_star_nom': ['scatter'], 'y_star_anc': ['line'], 'y_star_nom_past': ['line']},
        x_axis_format="datetime",
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,
        figsize=(10, 4),
        grid=True,
        vertical_lines=[model.integrator_parameters['end_timestamp']],
        legend_param = {
            "loc": "upper left",                                                # Note: position the legend outside the plot
            "bbox_to_anchor": (1.05, 1),                                        # Note: place the legend outside the axes
            "borderaxespad": 0.0,                                               # Note: reduce padding between the axes and the legend
            "frameon": False                                                    # Note: remove the legend box frame
        },
        layout_param = {
            "pad": 1.0,                                                         # Note: adjust padding to make space for the legend
            "rect": [0, 0, 0.85, 1]                                             # Note: shrink the plot area to fit the legend
    }
    )
    #------------------------------------------------------------------------------------------------------------# 
    # PLOT CONTROL ACTION
    # Define the data and variables for plotting
    data = {
        'past': df_u_past,
        'previous_optimal': Uk_guess_df,
        'v_star_past': v_star_past_df,
        'optimal': u_star_df,
        'v_star_nom': v_star_nom_df,
    }
    variables = [
        [('past', 'Uk'), ('previous_optimal', 'Uk'), ('v_star_past', 'Uk')],
        [('optimal', 'Uk')],
        [('v_star_nom', 'Uk')]
    ]
    # Generate plots
    plotter.generate_plots_with_incremental_variables(
        data=data,
        variables=variables,
        output_folder = os.path.join(directory, testname, 'plots'),
        run_index=index_control_run,
        title = 'u',
        y_label=r"Tomato flowrate (mL day$^{-1}$)",
        y_limit=(0,400),
        y_ticks = 5,
        x_limit=(model.integrator_parameters['end_timestamp'] - timedelta(days=model.integrator_parameters['update_window_lenght']), ancillary.integrator_parameters['end_timestamp']),
        #x_limit=[None],
        colors= {'past': ['grey'], 'optimal': ['black'], 'previous_optimal': ['grey'], 'v_star_past': ['pink'], 'v_star_nom': ['grey']}, #[plt.cm.tab10(np.linspace(0, 1, 3))[i] for i in range(3)]
        alphas={'past': [1], 'optimal': [0.8], 'previous_optimal': [0.5], 'v_star_past': [1], 'v_star_nom': [0.2]},
        linestyles={'past': ['-'], 'v_star_past': ['--']},
        markers={'v_star_past': ['o']},
        plot_types={'past': ['step'], 'optimal': ['scatter'], 'previous_optimal': ['scatter'], 'v_star_past': ['scatter'], 'v_star_nom': ['scatter']},
        x_axis_format="datetime",
        x_axis_ticks_format="%d-%m %H",
        tick_rotation=90,
        interval=12,
        figsize=(10, 4),
        grid=True,
        vertical_lines=[model.integrator_parameters['end_timestamp']],
        horizontal_lines=[opti_param_values['bnds_u'][0], opti_param_values['bnds_u'][1]],
        legend_param = {
            "loc": "upper left",                                                # Note: position the legend outside the plot
            "bbox_to_anchor": (1.05, 1),                                        # Note: place the legend outside the axes
            "borderaxespad": 0.0,                                               # Note: reduce padding between the axes and the legend
            "frameon": False                                                    # Note: remove the legend box frame
        },
        layout_param = {
            "pad": 1.0,                                                         # Note: adjust padding to make space for the legend
            "rect": [0, 0, 0.85, 1]                                             # Note: shrink the plot area to fit the legend
    }
    )
    
    #------------------------------------------------------------------------------------------------------------# 
    # # In[15]: CONVERT CONTROL ACTION TO ON/OFF TIMES FOR THE RASPBERRY PI RELAY
    # u_current = u_star_df['Uk'].iloc[0]
    # u_max = integrator_parameters['nmpc']['opti_parameters']['bnds_u'][1]
    # dt = integrator_parameters['nmpc']['control_interval']*3600
    # pump_dose_per_minute = 20                                                     # Note: mL/min (add it to integrator_parameters?)
    # period = 20                                                                   # Note: secc
    # tuple_seconds_ini, on_periods_tot, conversion_error = convert_u_to_pwm(u_current, u_max, dt, pump_dose_per_minute, period)
    # rounded_dosage = period*on_periods_tot*pump_dose_per_minute/60*24/dt*3600     # Note: mL/day
    # pwm_star_dict = dict(zip(['on_sec', 'off_sec', 'tot_on_periods', 'rounded_dosage'],
    #                 np.array([tuple_seconds_ini[0], tuple_seconds_ini[1], on_periods_tot, rounded_dosage])))
    # pwm_star_df = pd.DataFrame([pwm_star_dict])                                   # Note: duplicate last element
    # pwm_star_df.insert(0, 'Timestamp', integrator_parameters['nmpc']['start_timestamp'])
    # save_df_with_check(pwm_star_df, os.path.join(directory, testname, "Input", f'{ancillary.modelname}_pwm_actual_{ancillary.problemtype}.csv'), log=True)
    
    #------------------------------------------------------------------------------------------------------------# 
    # In[16]: UPDATE START/END TIMESTAMPS AND SAVE TO INTEGRATOR.JSON 
    integrator_parameters['model']['start_timestamp'] = integrator_parameters['model']['end_timestamp']
    integrator_parameters['model']['end_timestamp'] = integrator_parameters['model']['start_timestamp'] + timedelta(hours=ancillary.integrator_parameters['control_interval'])
    integrator_parameters['ekf']['start_timestamp'] = integrator_parameters['ekf']['end_timestamp']
    integrator_parameters['ekf']['end_timestamp'] = integrator_parameters['ekf']['start_timestamp'] + timedelta(hours=ancillary.integrator_parameters['control_interval'])
    integrator_parameters['nmpc']['start_timestamp'] = integrator_parameters['model']['end_timestamp']
    integrator_parameters['nmpc']['end_timestamp'] = integrator_parameters['model']['end_timestamp'] + timedelta(hours=ancillary.integrator_parameters['control_interval']*ancillary.integrator_parameters['Hp'])
    update_parameters('integrator.json', integrator_parameters['model'],'model', log=True)
    update_parameters('integrator.json', integrator_parameters['ekf'],'ekf', log=True)
    update_parameters('integrator.json', integrator_parameters['nmpc'],'nmpc', log=True)
    
    #------------------------------------------------------------------------------------------------------------# 
    # In[17]: UPDATE INDEX CONTROL RUN AND SAVE TO TXT
    index_control_run += 1
    index_control_run_df['IndexRun'] = index_control_run
    index_control_run_df.to_csv("index_control_run.txt", index=False)
    
    logger.info('############################################### END TUBE MAIN ###############################################')
    
    #------------------------------------------------------------------------------------------------------------#
# In[18]: MAIN FUNCTION
if __name__ == "__main__":
    # Allow the script to be executed from the terminal
    if len(sys.argv) != 2:
        print("Usage: python test_closedloop_tube_main.py <modelname>")
        sys.exit(1)

    modelname = sys.argv[1]
    main(modelname)

# # Set back to the original directory                                              # Uncomment these lines if you want to run this code as stand alone
# os.chdir(os.path.dirname(__file__))
# logger.info(f"The current working directory is now set to: {os.getcwd()}")