Integrazione codebase RSE

files/energy/input/DSM_optimizer/func/func.py

@@ -0,0 +1,459 @@
+import numpy as np
+import random
+import yaml
+# Sezione di print: 
+import pandas as pd
+import plotly.graph_objects as go
+
+###################################################################################################################################################
+
+def print_result(perc_econd,n_devices,total_energy_shared,total_incentive,total_energy_consumption,total_energy_cost,perc_error,*param):
+    print("\nGA OPT\n")
+    print("Percentuale condivisione: {:.2f} %".format(float(perc_econd)*100)) 
+    print("Percentuale condivisione per singolo dispositivo: {:.2f} %".format(float(perc_econd/n_devices)*100))
+    print(f'Total energy sherred: {total_energy_shared}')
+    print(f'Incentivo totale: {total_incentive}')
+    print(f'Total energy consumption: {total_energy_consumption}')
+    print(f'Costo totale: {total_energy_cost}')                
+    print("Optimization in percentuale su consumo: {:.2f}%".format(float((perc_error)*100)))
+    print('popolazione {}, generazione: {}, mutazione: {}'.format(*param))
+
+###################################################################################################################################################
+
+def stampa_result(dataframe, immission_profile, num_intervals, day, show_plot = True):
+    
+    # # print per tutti i dispositiv: 
+    df = pd.DataFrame({
+        'immission_profile': immission_profile,
+        
+    }, index=np.arange(num_intervals))
+   
+    # # Create figure:
+    x = pd.date_range(start='2018-01-01 00:00+00:00', freq='15min', inclusive='left', periods=97)   
+    fig = go.Figure()
+
+    # tutte queste sono cumulate, quindi stacked area
+    for column in dataframe.columns:
+        fig.add_trace(go.Scatter(
+            x = x, y = dataframe[column], mode='lines',
+                    name = column,
+                    line_shape = 'hv',
+                    textposition='top center',stackgroup='one'))
+    
+    # Inserisco il profilo di immissione: 
+    fig.add_trace(go.Scatter(
+            x = x, y = df['immission_profile'], mode='lines',
+                    name = 'immission_profile',
+                    line_shape = 'hv',
+                    textposition='top center'))
+
+    fig.update_yaxes(title_text = "Power (kW)")          # verificare se è effettivamente la potenza: 
+    fig.update_xaxes(title_text = 'Tempo (HH:MM)', dtick = 1000*60*30, range = [x[0], x[96]])
+    fig.update_layout(title= f'Day {day}', xaxis_tickformat = '%H:%M', plot_bgcolor='white', legend = dict(orientation='h',yanchor="top", y = -0.4, xanchor="left", x=0.01)) #barmode='stack') # si setta la modalità di rappresentazione del grafico a barre, in questo caso avranno una visualizzazione di tipo "stack", ovvero impilate le une sulle altre
+    
+    if show_plot:
+        fig.show()
+
+    #########################################################################################################
+
+    config = yaml.safe_load(open("config.yml", 'r'))
+    path = config['forldername_graphs_DSM_optimizer']
+
+    title = "DSM optimization day " + str(day)
+
+    fig.write_html(path + title + ".html")
+    fig.write_image(path + title + ".png", width=1000, height=1100/13.2*5, scale = 4)
+
+    #########################################################################################################
+
+    # print delle tipologie d'utente:
+    # fig = go.Figure()
+    
+    # Definizione dei raggruppamenti basati sulle stringhe specifiche
+    # user_types = ['single', 'coppia', 'famiglia_3p', 'famiglia_4p', 'famiglia_mag4p']
+
+    # Dizionario per memorizzare i risultati delle somme per ogni gruppo
+    # group_sums = {}
+
+    # Iterazione sui tipi di utenti per raggruppare e sommare le colonne
+    # for user_type in user_types:
+    #     # Selezione delle colonne che contengono la stringa specifica
+    #     columns_to_group = [col for col in dataframe.columns if user_type in col]
+        
+    #     # Somma delle colonne selezionate
+    #     df_grouped_sum = dataframe[columns_to_group].sum(axis=1)
+        
+    #     # Memorizzazione del risultato nel dizionario
+    #     group_sums[user_type] = df_grouped_sum
+
+    # Creazione di un nuovo DataFrame dai risultati delle somme
+    # df_grouped_sums = pd.DataFrame(group_sums)
+    
+    # # tutte queste sono cumulate, quindi stacked area
+    # for column in df_grouped_sums.columns:
+    #     fig.add_trace(go.Scatter(
+    #         x = x, y = df_grouped_sums[column], mode='lines',
+    #                 name = column,
+    #                 line_shape = 'hv', #'spline',
+    #                 textposition='top center',stackgroup='one'))
+    
+    # # Inserisco il profilo di immissione: 
+    # fig.add_trace(go.Scatter(
+    #         x = x, y = df['immission_profile'], mode='lines',
+    #                 name = 'immission_profile',
+    #                 line_shape = 'hv', #'spline',
+    #                 textposition='top center'))
+
+    # fig.update_yaxes(title_text = "Power (kW)")          # verificare se è effettivamente la potenza: 
+    # fig.update_xaxes(title_text = 'Tempo (HH:MM)', dtick = 1000*60*30, range = [x[0], x[96]])
+    # fig.update_layout(title= 'GA Optimization 2', xaxis_tickformat = '%H:%M', plot_bgcolor='white', legend = dict(orientation='h',yanchor="top", y = -0.4, xanchor="left", x=0.01)) #barmode='stack') # si setta la modalità di rappresentazione del grafico a barre, in questo caso avranno una visualizzazione di tipo "stack", ovvero impilate le une sulle altre
+    # #fig.show()
+
+    # Print con aggregati: 
+    # fig = go.Figure()
+     
+    # df = pd.DataFrame({
+    #     'autonsumption_profile': autoconsumption_profile,
+    #     'pv_profile': pv_profile
+    # }, index=np.arange(num_intervals))
+    # # Inserisco i profili di produzione: 
+    # fig.add_trace(go.Scatter(
+    #         x = x, y = df['pv_profile'], mode='lines',
+    #                 name = 'pv_profile',
+    #                 line_shape = 'hv', #'spline',
+    #                 textposition='top center'))
+    
+    # fig.add_trace(go.Scatter(
+    #     x = x, y = df['autonsumption_profile'], mode='lines',
+    #             name = 'autonsumption_profile',
+    #             line_shape = 'hv', #'spline',
+    #             textposition='top center',stackgroup='one'))
+
+    # # Inserisco il profilo dei dispositivi: 
+    # for column in df_grouped_sums.columns:
+    #     fig.add_trace(go.Scatter(
+    #         x = x, y = df_grouped_sums[column], mode='lines',
+    #                 name = column,
+    #                 line_shape = 'hv', #'spline',
+    #                 textposition='top center',stackgroup='one'))
+    
+    # fig.update_yaxes(title_text = "Power (kW)")          # verificare se è effettivamente la potenza: 
+    # fig.update_xaxes(title_text = 'Tempo (HH:MM)', dtick = 1000*60*30, range = [x[0], x[96]])
+    # fig.update_layout(title= 'GA Optimization 2', xaxis_tickformat = '%H:%M', plot_bgcolor='white', legend = dict(orientation='h',yanchor="top", y = -0.4, xanchor="left", x=0.01)) #barmode='stack') # si setta la modalità di rappresentazione del grafico a barre, in questo caso avranno una visualizzazione di tipo "stack", ovvero impilate le une sulle altre
+    # fig.update_layout(autosize=False,width=1400,height=800)
+    # html_str = fig.to_html(full_html=False, include_plotlyjs='cdn')
+    
+    # #fig.show()
+
+###################################################################################################################################################
+
+def stampa_best_solution_gen(fig, best_solution_gen,params, n_device):
+    # create the dataframe:
+    
+    df = pd.DataFrame({
+        'fitness_values': best_solution_gen,
+        
+    }, index=np.arange(len(best_solution_gen)))
+    
+    x = df.index
+    
+    fig.add_trace(go.Scatter(
+            x = x, y = df['fitness_values'], mode='lines',
+                    name = 'ft_p_{}_g_{}_m_{}'.format(*params),
+                    line_shape = 'hv',
+                    textposition='top center'))
+
+    fig.update_yaxes(title_text = "fitness_values")
+    fig.update_xaxes(title_text = 'Number of Gen')
+    fig.update_layout(title= f'GA Optimization Device for {n_device} devices', plot_bgcolor='white', legend = dict(orientation='h',yanchor="top", y = -0.4, xanchor="left", x=0.01)) #barmode='stack') # si setta la modalità di rappresentazione del grafico a barre, in questo caso avranno una visualizzazione di tipo "stack", ovvero impilate le une sulle altre
+
+###################################################################################################################################################
+
+def intervallo_funzionamento(start,end,n):
+    giorno = np.zeros(n,int)
+    for i in range(start,end):
+        giorno[i]=1 
+    return giorno
+
+###################################################################################################################################################
+
+# Funzione per generare un profilo di consumo casuale di lunghezza specificata
+def generate_random_consumption_profile(length):
+    return np.random.rand(length)
+
+###################################################################################################################################################
+
+def first_nonzero_index(arr):
+    """funzione per trovare l'indice del primo elemento non nullo di un array
+
+    Args:
+        arr (array): array di numpy
+
+    Returns:
+        int: indice del primo elemento non nullo, 0 se l'array è vuoto
+    """
+    nonzero_indices = np.nonzero(arr)[0]
+    return nonzero_indices[0] if nonzero_indices.size > 0 else 0
+
+###################################################################################################################################################
+
+# Function to calculate total energy cost for a device on/off schedule
+# Precedente
+# def calculate_device_cost(device_schedule, immission_profile, max_power_contract, energy_cost_per_hour):
+    
+    # total_energy_cost = 0
+    # total_energy_consumption = 0
+    # total_energy_shared      = 0
+    # incentivo = 0.11
+    # total_cost = 0
+    
+    # num_intervals = len(device_schedule[0])
+    # for t in range(num_intervals):
+    #     total_power = sum(device[t] for device in device_schedule)
+        
+                
+    #     for device in device_schedule:
+    #         # Verifica del vincolo di potenza contrattuale massima
+    #         if device[t] > max_power_contract:
+    #             return float('inf')  # Penalità per violazione del vincolo
+        
+    #     total_energy_cost = total_power*energy_cost_per_hour[t]
+    #     total_energy_shared = min(total_power, immission_profile[t])*incentivo
+            
+    #     # # Se supero energia immessa:
+    #     # if total_power > immission_profile[t]:
+    #     #     total_energy_cost += total_power
+
+    #     # Calcolo del costo considerando l'immissione nel profilo
+    #     total_cost += total_energy_cost - total_energy_shared
+    # return total_cost
+
+###################################################################################################################################################
+
+# Post rev. NO check
+# def calculate_device_cost(device_schedule, immission_profile, energy_cost_per_hour):
+#     device_schedule = np.array(device_schedule)
+#     total_power = device_schedule.sum(axis=0)
+
+#     total_energy_cost = total_power * np.array(energy_cost_per_hour)
+#     total_energy_shared = np.minimum(total_power, immission_profile) * 0.11
+#     total_cost = total_energy_cost.sum() - total_energy_shared.sum()   
+#     return total_cost
+
+###################################################################################################################################################
+
+def calculate_device_cost(device_schedule, immission_profile, energy_cost_per_hour):
+    
+    # device_schedule = np.asarray(device_schedule)
+    # total_power = device_schedule.sum(axis=0)
+    total_power = np.sum(device_schedule, axis=0)
+    total_energy_cost = np.dot(total_power, energy_cost_per_hour)
+    
+    total_energy_shared = np.dot(np.minimum(total_power, immission_profile), 0.11)
+    
+    total_cost = total_energy_cost.sum() - total_energy_shared.sum()
+    return total_cost
+
+# Objective function to calculate total energy cost for all solutions
+# def objective_function(solutions):
+#     total_costs = []
+#     for solution in solutions:
+#         total_cost = np.sum(calculate_device_cost(solution, consumption_profiles))
+#         total_costs.append(total_cost)
+#     return total_costs
+
+###################################################################################################################################################
+
+# Funzione per generare una soluzione casuale rispettando il profilo di consumo di ciascun dispositivo
+# PRE rev. 
+def generate_random_solution(num_intervals, consumption_profiles, allowed_intervals, possible_starts):
+    solution = []
+    for profile in consumption_profiles:
+        # Seleziona casualmente un punto di inizio nel profilo di consumo
+        # Integro ora il vincolo sugli istanti di tempo permessi:            
+        start_index = random.choice(possible_starts)
+        # start_index = np.random.randint(0, num_intervals - len(profile))
+        # Costruisci il programma di utilizzo del dispositivo utilizzando il profilo di consumo
+        device_schedule = [0] * num_intervals
+        device_schedule[start_index:start_index + len(profile)] = profile
+        solution.append(device_schedule)
+    return solution
+
+###################################################################################################################################################
+
+# def generate_random_solution(num_intervals, consumption_profiles, allowed_intervals, possible_starts):
+#     solution = np.zeros((len(consumption_profiles), num_intervals))
+#     profile_lengths = np.array([len(profile) for profile in consumption_profiles])
+    
+#     for i, profile in enumerate(consumption_profiles):
+#         start_index = np.random.choice(possible_starts)
+#         solution[i, start_index:start_index + profile_lengths[i]] = profile
+    
+#     return solution
+
+###################################################################################################################################################
+
+# #Mutation function to randomly modify an existing solution
+def mutate_solution(solution, mutation_rate, consumption_profiles, allowed_intervals, num_intervals, possible_starts):
+    for i, device_schedule in enumerate(solution):
+        if np.random.rand() < mutation_rate:
+            # Seleziona casualmente un nuovo punto di inizio per il profilo di consumo
+            # Integro ora il vincolo sugli istanti di tempo permessi:            
+            profile = consumption_profiles[i]            
+            new_start_index = random.choice(possible_starts)
+            #new_start_index = np.random.randint(0, len(device_schedule) - len(profile))
+            # Aggiorna il programma di utilizzo del dispositivo
+            new_device_schedule = [0] * len(device_schedule)
+            new_device_schedule[new_start_index:new_start_index + len(profile)] = profile
+            solution[i] = new_device_schedule
+    return solution
+
+###################################################################################################################################################
+
+# def mutate_solution(solution, mutation_rate, consumption_profiles, allowed_intervals, num_intervals, possible_starts):
+#     profile_lengths = np.array([len(profile) for profile in consumption_profiles])
+    
+#     for i, device_schedule in enumerate(solution):
+#         if np.random.rand() < mutation_rate:
+#             profile = consumption_profiles[i]
+#             new_start_index = np.random.choice(possible_starts)
+#             solution[i, :] = 0
+#             solution[i, new_start_index:new_start_index + profile_lengths[i]] = profile
+    
+#     return solution
+
+###################################################################################################################################################
+
+# PRE 
+#Algoritmo genetico
+def genetic_algorithm(num_intervals, population_size, generations, mutation_rate, consumption_profiles,immission_profile, max_power_contract, energy_cost_per_hour, allowed_intervals):
+    
+    # calcolo possibili istanti attivazione:
+    for profile in consumption_profiles:
+        possible_starts = [i for i in allowed_intervals if i + len(profile)-1 <= allowed_intervals[-1]] #ultimo valore di allowed intervals
+        possible_starts = sorted(possible_starts, key=lambda x: -sum(immission_profile[x:x + len(profile)]))
+        
+    # Inizializza la popolazione iniziale di soluzioni casuali
+    population = [generate_random_solution(num_intervals, consumption_profiles, allowed_intervals, possible_starts) for _ in range(population_size)]
+    
+    # per analisi, eliminare poi: 
+    best_solution_gen_value =[]
+    
+    # Esegui le iterazioni per un numero fisso di generazioni
+    for generation in range(generations):
+        # Calcola il valore di fitness per ogni soluzione nella popolazione
+        # fitness_values = objective_function(population, consumption_profiles)
+        fitness_values = [np.sum(calculate_device_cost(solution, immission_profile, energy_cost_per_hour)) for solution in population]
+        
+        # Seleziona i genitori in base al loro valore di fitness
+        parents = [population[i] for i in np.argsort(fitness_values)[:int(population_size/2)]]  
+        
+        # Genera nuovi individui attraverso incrocio (crossover)
+        children = []
+        for _ in range(population_size - len(parents)):
+            parent1, parent2 = random.sample(parents, 2)
+            crossover_point = np.random.randint(num_intervals)
+            child = parent1[:crossover_point] + parent2[crossover_point:]
+            #child = np.concatenate((parent1[:, :crossover_point], parent2[:, crossover_point:]), axis=1)
+            children.append(child)
+
+        # Aggiorna la popolazione con nuovi individui generati attraverso crossover e mutazione
+        population = parents + children
+        population = [mutate_solution(solution, mutation_rate, consumption_profiles, allowed_intervals, num_intervals, possible_starts) for solution in population]
+        
+        # Elitismo: conserva la migliore soluzione
+        best_solution = min(population, key=lambda sol: calculate_device_cost(sol, immission_profile, energy_cost_per_hour))
+        population[0] = best_solution
+        
+        # Introduci diversità ogni 10 generazioni
+        if generation % 5 == 0:
+            for i in range(len(population)):
+                if np.random.rand() < 0.1:
+                    population[i] = generate_random_solution(num_intervals, consumption_profiles,allowed_intervals, possible_starts)    
+        
+        
+        # Per analisi ---> da eliminare poi
+        best_solution_gen_value.append(fitness_values[np.argmin(fitness_values)])
+        
+    # Calcola il valore di fitness per ogni soluzione nella popolazione finale
+    fitness_values = [np.sum(calculate_device_cost(solution, immission_profile, energy_cost_per_hour)) for solution in population]
+    
+    # Trova e restituisci la migliore soluzione
+    best_solution_index = np.argmin(fitness_values)
+    best_solution = population[best_solution_index]
+    best_cost = fitness_values[best_solution_index]
+    
+    return best_solution, best_cost, best_solution_gen_value
+
+###################################################################################################################################################
+
+# POST using Concurrency: too slow
+
+# from concurrent.futures import ProcessPoolExecutor
+# def calculate_fitness(solution, immission_profile, energy_cost_per_hour):
+#     return calculate_device_cost(solution, immission_profile, energy_cost_per_hour)
+
+###################################################################################################################################################
+
+# def fitness_wrapper(args):
+#     solution, immission_profile, energy_cost_per_hour = args
+#     return calculate_fitness(solution, immission_profile, energy_cost_per_hour)
+
+###################################################################################################################################################
+
+# # # # POST rev No Check
+# def genetic_algorithm(num_intervals, population_size, generations, mutation_rate, consumption_profiles, immission_profile, max_power_contract, energy_cost_per_hour, allowed_intervals):
+#     population = [generate_random_solution(num_intervals, consumption_profiles, allowed_intervals, immission_profile) for _ in range(population_size)]
+    
+#     best_solution_gen_value = []
+    
+#     with ProcessPoolExecutor() as executor:
+#         for generation in range(generations):
+#             fitness_values = list(executor.map(fitness_wrapper, [(sol, immission_profile, energy_cost_per_hour) for sol in population]))
+            
+#             parents = [population[i] for i in np.argsort(fitness_values)[:int(population_size / 2)]]
+            
+#             children = []
+#             for _ in range(population_size - len(parents)):
+#                 parent1, parent2 = random.sample(parents, 2)
+#                 crossover_point = np.random.randint(num_intervals)
+#                 if len(parent1) != len(parent2):
+#                     parent1, parent2 = np.broadcast_to(parent1, (len(parent2), len(parent2[0]))), np.broadcast_to(parent2, (len(parent1), len(parent1[0])))
+#                 child = np.concatenate((parent1[:crossover_point], parent2[crossover_point:]), axis=0)
+#                 children.append(child)
+            
+#             population = parents + children
+#             population = parents + children
+#             population = [mutate_solution(solution, mutation_rate, consumption_profiles, allowed_intervals, num_intervals, immission_profile) for solution in population]
+            
+#             best_solution = min(population, key=lambda sol: fitness_wrapper((sol, immission_profile, energy_cost_per_hour)))
+#             population[0] = best_solution
+            
+#             if generation % 5 == 0:
+#                 for i in range(len(population)):
+#                     if np.random.rand() < 0.1:
+#                         population[i] = generate_random_solution(num_intervals, consumption_profiles, allowed_intervals, immission_profile)
+            
+#             best_solution_gen_value.append(min(fitness_values))
+    
+#     fitness_values = [fitness_wrapper((sol, immission_profile, energy_cost_per_hour)) for sol in population]
+    
+#     best_solution_index = np.argmin(fitness_values)
+#     best_solution = population[best_solution_index]
+#     best_cost = fitness_values[best_solution_index]
+    
+#     return best_solution, best_cost, best_solution_gen_value
+
+###################################################################################################################################################
+
+
+
+###################################################################################################################################################
+
+
+
+###################################################################################################################################################
+
+
+
+###################################################################################################################################################

files/energy/input/DSM_optimizer/func/ott.py

@@ -0,0 +1,260 @@
+import numpy as np
+import pandas as pd
+import plotly.graph_objects as go
+from files.energy.input.DSM_optimizer.func.func import genetic_algorithm, stampa_result
+
+###################################################################################################################################################
+
+def ottimizzazione(immission_profile, consumption_profiles,num_intervals, day, create_plot = True, show_plot = True):
+    
+    ##### Trovare valori permessi immissione nei quali andare a inserire elettrodomestici:
+    # Trovare gli indici dei valori maggiori di zero
+       
+    allowed_intervals = np.where(np.array(immission_profile) > 0)[0].tolist()
+    max_power_contract = 3
+
+    # Grid search parameters
+    n_devices = len(consumption_profiles)
+    
+    # consumo totale:
+    #consumption_total = sum([sum(consumption_profile) for consumption_profile in consumption_profiles])
+    consumption_total = sum(map(np.sum, consumption_profiles))
+    # Define energy cost (e.g., cost per kWh)
+    # energy_cost_per_hour = [0.3]* num_intervals #- intervallo_funzionamento(35,60,96)*0.2 # Electricity price per kWh
+    energy_cost_per_hour = np.full(num_intervals, 0.3)
+    # Generation profile
+    # # Trova il valore massimo
+    # max_value = max(pv_profile)
+
+    # # Normalizza la lista
+    # normalized_profile = [x / max_value for x in pv_profile]
+
+    # # scalo in funzione della taglia: 
+    # taglia_FV = max(2,1*(math.ceil(n_devices/4)))
+    # #taglia_FV = consumption_total/(sum(normalized_profile)*0.04*n_devices)   # la taglia deve rendere il rapporto tra cosnumo singolo dispositivo e totale energia immessa = 0.6
+        
+    # generation_profile = [x * taglia_FV for x in normalized_profile]   
+    
+    # Immission profile
+    # """ Questa parte è da modificare quando avremo la produzione fotovoltaica"""
+    # immission_profile = [
+    #     max(0, generation_profile[i] - immission_profile[i])
+    #     for i in range(num_intervals)
+    # ]
+    
+    # Calcoliamo il rapporto tra il consumo e l'immissionee, per valutare percentuale possibili condivisione sul totale: 
+
+    # immission_total = sum(immission_profile)
+    immission_total = np.sum(immission_profile)
+    
+    perc_econd= min(consumption_total, immission_total)/immission_total  # varia tra 0 e 1 
+    perc_error = 0
+    perc_error_1 = 0
+    
+    # Calcoliamo il rapporto tra il consumo e l'immissionee, per valutare percentuale possibili condivisione sul totale: 
+    
+    if perc_econd/n_devices < 0.04 and perc_econd < 0.65:  
+        # Grid search parameters
+        population_size = min(50,max(10,n_devices*2))
+        generations = 20
+        mutation_rate = 0.05
+        
+    elif perc_econd/n_devices < 0.06 and perc_econd < 0.55:
+        # Grid search parameters
+        population_size = min(50,n_devices*3)
+        generations = max(20,min(50,n_devices*2))
+        mutation_rate = 0.05
+    
+    elif perc_econd/n_devices < 0.01 and perc_econd > 0.65:
+        # Grid search parameters
+        """"da mnodificare---->"""
+        population_size = 20
+        generations = 20
+        mutation_rate = 0.05
+    
+    else:
+        # Grid search parameters
+        population_size = min(50,n_devices*3)
+        generations = max(20,min(50,n_devices*6))
+        mutation_rate = 0.05
+
+    # print(f"\nNUMBER DEVICE(S): {n_devices}")
+    
+    #################### SOLUTIONS:####################
+    # Creazione del dizionario di input alla funzione.
+
+    params = {
+            "num_intervals": num_intervals,
+            "population_size": population_size,
+            "generations": generations,
+            "mutation_rate": mutation_rate,
+            "consumption_profiles": consumption_profiles,
+            "immission_profile": immission_profile,
+            "max_power_contract": max_power_contract,
+            "energy_cost_per_hour": energy_cost_per_hour,
+            "allowed_intervals": allowed_intervals
+        }
+
+    # vera e propria funzione di ottimizzazione:
+    best_solution = genetic_algorithm(**params)
+    param = [population_size, generations, mutation_rate]
+    
+    # funzione per il plot
+    #################### PLOT:####################    
+    # create the dataframe devices:
+    # Possiamo differenziarlo in funzione della tipologia di utenza e del nome del dispositivo: ad esempio user_n_single_washing_machine ecc..
+    # magari creo una funzione ad hoc per questo, vediamo dopo  ? 
+    
+    ##############
+    # Creazione delle colonne in modo compatto
+    u = 0
+    #columns = [f'user_{u}_{key}' for key, count in sched_dev.items() for _ in range(1, count + 1)]
+    columns = []    
+    for count in range(int(n_devices/2)):
+        for key in ['washing_machine','dish_washer']:
+            columns.append(f'user_{u}_{key}')
+        u = u +1
+
+    # Creazione del DataFrame
+    df_plot = pd.DataFrame(data=np.transpose(best_solution[0]),
+                        index=np.arange(num_intervals),
+                        columns=columns)
+    
+    # df_plot = pd.DataFrame(data=np.transpose(best_solution[0]),
+    #                     index=np.arange(num_intervals),
+    #                     columns=[f'user_{u=u+1}_{key}' for u in range(1, n_devices + 1) for key in sched_dev.keys()])
+    
+    ##############
+
+    # Crea il DataFrame immissione:
+    df = pd.DataFrame({
+        'immission_profile': immission_profile,
+        
+    }, index=np.arange(num_intervals))
+
+    # Voglio calcolarmi l'energia totale condivisa. Con ritorno economico dato da incentivo totale: 
+    # max(immissione[t], sum(prelievo utenze domestiche[t])*incentivo for 
+    # Calcolare la somma degli elementi su ogni riga per il primo DataFrame
+    sum_df1_1 = df_plot.sum(axis=1)
+    # Sottrarre il valore corrispondente del secondo DataFrame
+    energy_shared_1 = np.minimum(sum_df1_1.values, df['immission_profile'].values)
+    
+    # Creare un nuovo DataFrame con una colonna chiamata 'energy_shared'
+    #df_energy_shared_1 = pd.DataFrame({'energy_shared': energy_shared_1})
+    #print(f'Vector of energy sherred: {df_energy_shared}')
+    #total_energy_shared_1 = df_energy_shared_1['energy_shared'].sum()
+    total_energy_shared_1 = np.sum(energy_shared_1)
+    total_incentive_1 = total_energy_shared_1*0.11
+    total_energy_consumption_1 = np.sum(sum_df1_1)
+    total_energy_cost_1 = total_energy_consumption_1*0.3    
+    #total_immission_energy_1 = df['immission_profile'].sum()
+    total_immission_energy_1 = np.sum(immission_profile)
+    perc_error_1 = total_energy_shared_1/min(total_immission_energy_1,total_energy_consumption_1)
+    
+    if  perc_error_1<=1:  
+                 
+        e_cond_i = 0
+        e_cond_f =  0 
+        e_cons_profile = np.zeros(96)
+        
+        # Crea il DataFrame dispositivi in surplus:
+        df_user_over = pd.DataFrame({                        
+                                    }, 
+                                    index=np.arange(num_intervals))
+        # Ciclo
+        for column in df_plot.columns:
+            e_cons_profile += df_plot[column]
+            #e_cond_f = sum([min(immission_profile[t], e_cons_profile[t])for t in range(len(e_cons_profile))])
+            e_cond_f = np.sum(np.minimum(immission_profile, e_cons_profile))
+            
+            if e_cond_f == e_cond_i:
+                df_user_over[column] = df_plot[column]
+                e_cons_profile -= df_plot[column]
+            
+            #if e_cond_f < sum(e_cons_profile):
+            if e_cond_f < np.sum(e_cons_profile):
+                
+                df_user_over[column] = df_plot[column]
+                e_cons_profile -= df_plot[column]
+                e_cond_f = e_cond_i
+            
+            e_cond_i = e_cond_f
+
+        
+        if not df_user_over.dropna().empty:
+            df_total_new =df_plot.drop(columns= df_user_over.columns)
+                       
+            df_total_new_sum = df_total_new.sum(axis=1)
+            
+            # Trasformiamo per vettorializzare:
+            # Calcolo del nuovo profilo di immissione utilizzando Pandas
+            #immission_profile_series = pd.Series(immission_profile)
+            # Calcolo del nuovo profilo di immissione
+            #immission_profile_new = (immission_profile_series - df_total_new_sum).clip(lower=0).tolist()
+            immission_profile_new = np.clip(immission_profile - df_total_new_sum.values, 0, None)
+            
+            ##### Trovare valori permessi immissione nei quali andare a inserire elettrodomestici:
+            # Trovare valori permessi di immissione
+            allowed_intervals = np.where(np.array(immission_profile) > 0)[0].tolist()
+
+            # Trasposizione e filtraggio dei profili di consumo
+            consumption_profiles = [df_user_over.transpose().values[t][df_user_over.transpose().values[t] > 0] for t in range(len(df_user_over.columns))]
+
+            
+            population_size = min(50,len(consumption_profiles)*20)
+            generations     = min(50,len(consumption_profiles)*20)
+            
+            params = {
+                "num_intervals": num_intervals,
+                "population_size": population_size,
+                "generations": generations,
+                "mutation_rate": mutation_rate,
+                "consumption_profiles": consumption_profiles,
+                "immission_profile": immission_profile_new,
+                "max_power_contract": max_power_contract,
+                "energy_cost_per_hour": energy_cost_per_hour,
+                "allowed_intervals": allowed_intervals
+            }
+             
+            # second GA
+            best_solution_2 = genetic_algorithm(**params)
+            param = [population_size, generations, mutation_rate]
+            # create the dataframe devices:
+            df_user_opt_2 = pd.DataFrame(data=np.transpose(best_solution_2[0]),
+                        index=np.arange(num_intervals),
+                        columns= df_user_over.columns)
+   
+            df_total_new = pd.concat([df_total_new, df_user_opt_2], axis=1)
+            sum_df1 = df_total_new.sum(axis=1)
+            # Sottrarre il valore corrispondente del secondo DataFrame
+            energy_shared = np.minimum(sum_df1.values, df['immission_profile'].values)
+            # Creare un nuovo DataFrame con una colonna chiamata 'energy_shared'
+            # df_energy_shared = pd.DataFrame({'energy_shared': energy_shared})
+            # total_energy_shared = df_energy_shared['energy_shared'].sum()
+            total_energy_shared = np.sum(energy_shared)
+            total_incentive = total_energy_shared*0.11
+            #total_energy_consumption = sum_df1.sum()
+            total_energy_consumption = np.sum(sum_df1)
+            total_energy_cost = total_energy_consumption*0.3
+            #total_immission_energy = df['immission_profile'].sum()
+            #total_immission_energy = total_immission_energy_1.copy()
+            perc_error = total_energy_shared/min(total_immission_energy_1,total_energy_consumption)
+            
+           
+    # teniamo solo quello con percentuale condivisione migliore, quindi: 
+    if perc_error < perc_error_1:
+        #print("Return first")
+        if create_plot:
+            stampa_result(df_plot, immission_profile, num_intervals, day, show_plot)
+        #print_result(perc_econd,n_devices,total_energy_shared_1,total_incentive_1,total_energy_consumption_1,total_energy_cost_1,perc_error_1,*param)
+        #perc_condivisione_utilizzo = total_energy_shared_1/total_energy_consumption_1
+        return df_plot
+    
+    else:
+        #print("Return second")
+        if create_plot:
+            stampa_result(df_total_new, immission_profile, num_intervals, day, show_plot)
+        #print_result(perc_econd,n_devices,total_energy_shared,total_incentive,total_energy_consumption,total_energy_cost,perc_error,*param)
+        #perc_condivisione_utilizzo = total_energy_shared/total_energy_consumption # type: ignore
+        return df_total_new
+