optimisation_batterie_EOD.py

# -*- coding: utf-8 -*-
"""
Created on Mon Nov 25 11:54:31 2024

@author: Despouysoli
"""


import pandas as pd
import highspy
from linopy import Model
import matplotlib.pyplot as plt
import numpy as np
import os





def preparer_modele (dico_site, df_chroniques, delta_t):
    
    # Constantes
    # =========================================================================
    
    # Fraction d'heure correspondant à un pas de temps
    tau = delta_t / 60
    # Puissance consommée par les auxiliaires de la batterie
    p_aux = dico_site['conso_veille'] * dico_site['e_max'] * 1e-3
    assert p_aux < dico_site['p_max_batt'] and -p_aux > dico_site['p_min_batt']
    # Recherche de la plus grande plage de limitations (s'il en existe)
    fenetre_maximale = 1
    for limitation in dico_site['Limitations spécifiques']:
        if limitation["fenetre"] > fenetre_maximale:
            fenetre_maximale = limitation["fenetre"]
    # big M pour les puissances de batteries, puis pour les energies
    big_m_p_batt = abs(dico_site['p_max_batt']) + p_aux + 1
    big_m_energies = big_m_p_batt * tau
    
    # Création du modèle et des dimensions
    # =========================================================================

    # Définition des indices
    timesteps = pd.Index(range(-fenetre_maximale, df_chroniques.shape[0]), name="timestep")
    # Creation du modèle
    m = Model()
    
    # Variables
    # =========================================================================
    
    # ---------- Puissances ----------
    # Puissances injectées / soutirées par le site (MW)
    p_site = m.add_variables(name='p_site', coords=[timesteps], dims=['timestep'])
    # Puissances injectées / soutirées par la batterie (MW)
    p_batt = m.add_variables(name='p_batterie', coords=[timesteps], dims=['timestep'])
    # Puissances spilled par le site (MW)
    p_spill = m.add_variables(name='p_spilled', coords=[timesteps], dims=['timestep'],
                              lower=0,
                              upper=df_chroniques.p_inj.clip(lower=0).reindex(timesteps).fillna(0))

    # ---------- Energies ----------
    # Energie ***vue de la batterie *** pendant un timestep où elle charge (MWh)
    e_charge = m.add_variables(name='e_charge', coords=[timesteps], dims=['timestep'],
                              lower=0, )
                              # upper=-tau * df_chroniques.p_min.reindex(timesteps).fillna(0) * dico_site['eta_c'] / 100)
    # Energie ***vue de la batterie *** pendant un timestep où elle se décharge (MWh)
    e_decharge = m.add_variables(name='e_decharge', coords=[timesteps], dims=['timestep'],
                                 lower=0, )
                                 # upper=tau * df_chroniques.p_max.reindex(timesteps).fillna(0) / dico_site['eta_d'] * 100)
    # Energie de la batterie ***à la fin de l'heure courante*** (MWh)
    e = m.add_variables(name='energie', coords=[timesteps], dims=['timestep'],
                        lower=dico_site['soc_min'] * dico_site['e_max'] / 100, 
                        upper=dico_site['soc_max'] * dico_site['e_max'] / 100)
    # # Etat de la batterie : en charge / non
    en_charge = m.add_variables(name='en_charge', coords=[timesteps], dims=['timestep'],
                                binary=True)
    
    # Contraintes
    # =========================================================================
    
    # ---------- Puissance de la batterie ----------
    # Puissance de la batterie avant l'optimisation
    m.add_constraints(p_batt.where(timesteps < 0) == 0,
                      name='contrainte p_batt pendant les timesteps négatifs')
    # Puissance minimale de la batterie et des auxiliaires pendant l'optimisation
    m.add_constraints(p_batt.where(timesteps >= 0) >= dico_site['p_min_batt'] + p_aux,
                      name='borne inf pour p_batt+p_aux')
    # Puissance maximale de la batterie et des auxiliaires pendant l'optimisation
    m.add_constraints(p_batt.where(timesteps >= 0) <= dico_site['p_max_batt'] + p_aux,
                      name='borne sup pour p_batt+p_aux')
    
    # ---------- Puissance spilled ----------
    # ==> intégré dans la déclaration de la bariable p_spill
    
    # ---------- Puissance site ----------
    # Calcul / définition de p_site
    lhs = p_site - p_batt + p_spill - df_chroniques.p_inj.reindex(timesteps).fillna(0) + p_aux
    m.add_constraints(lhs == 0,
                      mask = lhs.data.timestep >= 0,
                      name='définition de p_site')
    # Contrainte de borne inf sur p_site
    m.add_constraints(p_site >= df_chroniques.p_min.clip(lower=-dico_site['p_racc']).reindex(timesteps).fillna(-dico_site['p_racc']),
                      name='définition de borne inf sur p_site')
    # Contrainte de borne sup sur p_site
    m.add_constraints(p_site <= df_chroniques.p_max.clip(upper=dico_site['p_racc']).reindex(timesteps).fillna(dico_site['p_racc']),
                      name='définition de borne sup sur p_site')
    
    # ---------- Energie chargée / déchargée à chaque pas de temps ----------
    # Contraintes sur les puissances de charge en prenant en compte le rendement
    m.add_constraints(e_charge  >= -p_batt * dico_site['eta_c'] / 100 * tau, 
                      name='contrainte charge effective')    
    # Contraintes sur les puissances de décharge en prenant en compte le rendement
    m.add_constraints(e_decharge  >= p_batt / dico_site['eta_d'] * 100 * tau, 
                      name='contrainte decharge effective')    
    
    # ---------- Non-simultaneite entre charge et decharge ----------
    m.add_constraints(e_decharge * dico_site['eta_d'] / 100 - e_charge / dico_site['eta_c'] * 100 == p_batt * tau, 
                      name='contrainte liant charge et decharge (non-simultaneite -- puissance)')
    # Contraintes pour forcer la valeur de "en_charge"
    m.add_constraints(en_charge * big_m_p_batt  + p_batt >= 0, 
                      name='p_batt<0 ==> en_charge==1')
    m.add_constraints(-en_charge * big_m_p_batt  -p_batt + big_m_p_batt >= 0, 
                      name='p_batt>0 ==> en_charge==0')
    # Si en_charge = 1, alors e_decharge = 0, et si en_charge = 0, alors e_charge = 0
    m.add_constraints( e_decharge <= big_m_energies*(1-en_charge), 
                      name='en_charge = 1 ==> e_decharge = 0')
    m.add_constraints( e_charge <= big_m_energies*en_charge, 
                      name='en_charge = 0 ==> e_charge = 0')

    # ---------- Couplage entre énergies ----------
    lhs = e - e.shift(timestep=+1) - e_charge + e_decharge
    mask = lhs.data.timestep >= 0
    m.add_constraints(lhs=lhs, sign='=', rhs=0,  mask=mask,
                      name='contrainte etat de charge')

    # ---------- Contraintes pour quelques pas de temps en énergies ----------
    # Energie périodique (hors pas de temps négatifs et dernier pas de temps)
    if dico_site['periode_soc'] is not None and dico_site['soc_periodique'] is not None:
        mask = (e.data.timestep % dico_site['periode_soc'] == 0) & (e.data.timestep > 0) &\
            (e.data.timestep != e.data.timestep[-1])
        m.add_constraints(e == dico_site['soc_periodique'] * dico_site['e_max'] / 100,  
                          mask=mask, name='contrainte etat de charge périodique')
    # Energie initiale
    if dico_site['soc_initiale'] is not None:
        m.add_constraints(e == dico_site['soc_initiale'] * dico_site['e_max'] / 100, 
                          mask=e.data.timestep < 0,
                          name='contrainte etat de charge initial')
    # Energie finale
    if dico_site['soc_finale'] is not None:
        m.add_constraints(e == dico_site['soc_finale'] * dico_site['e_max'] / 100, 
                          mask=e.data.timestep == e.data.timestep[-1],
                          name='contrainte etat de charge finale')

    # ---------- Limitations spécifiques ----------
    for limitation in dico_site['Limitations spécifiques']:
        if limitation['type'] == 'charge':
            m.add_constraints(e_charge.rolling(timestep=limitation['fenetre']).sum() <= limitation['valeur'], 
                              name=f"Charge cumulee sur une fenetre glissante de {limitation['fenetre']} pas de temps inf. a {limitation['valeur']}")
        elif limitation['type'] == 'décharge':
            m.add_constraints(e_decharge.rolling(timestep=limitation['fenetre']).sum() <= limitation['valeur'], 
                              name=f"Decharge cumulee sur une fenetre glissante de {limitation['fenetre']} pas de temps inf. a {limitation['valeur']}")
        else:
            raise Exception(f"Les limitations spéciales de type {limitation['type']} ne sont pas gérées")

    return m, timesteps, tau



def generer_dataframe_solution(m, dico_site, df_chroniques):
    # Mise en forme des résultats et ajout des données d'entrée pour constituer
    # un DataFrame de synthèse
    df_synthese = m.solution.to_dataframe()
    df_synthese['p_inj'] = df_chroniques.p_inj
    # labels_chg_dchg = ['energie_charge_batterie', 'energie_decharge_batterie']
    # labels_chg_dchg_cumules = ['energie_charge_cumulee', 'energie_decharge_cumulee']
    # df_synthese[labels_chg_dchg_cumules] = df_synthese.rolling(window=window, min_periods=window).sum().loc[:, labels_chg_dchg]
    df_synthese = df_synthese.loc[ df_synthese.index >= 0 ]
    df_synthese.loc[:, 'p_min'] = df_chroniques.p_min
    df_synthese.loc[:, 'p_max'] = df_chroniques.p_max
    df_synthese.loc[:, 'prix'] = df_chroniques.prix
    df_synthese.loc[:, 'gains'] = (df_synthese.loc[:, 'prix'] * df_synthese.loc[:, 'p_site']).round(2)
    df_synthese.loc[:, 'gains_cumules'] = df_synthese.loc[:, 'gains'].cumsum()
    df_synthese = df_synthese.reindex(range(-1, df_synthese.shape[0]))
    if dico_site['soc_initiale'] is not None:
        df_synthese.loc[-1, 'energie'] = dico_site['soc_initiale'] * dico_site['e_max'] / 100
    # Arrondi à la 2ie décimale
    ordered_cols = ['p_batterie', 'p_spilled', 'p_inj', 'p_site', 'p_min', 'p_max', 
                    'e_charge', 'e_decharge', 'energie', 'en_charge', 
                    # 'e_charge', 'e_decharge', 'energie',  
                    'prix', 'gains', 'gains_cumules']
    df_synthese = df_synthese[ordered_cols].round(2)
    df_synthese.en_charge = df_synthese.en_charge.astype(bool)
    return df_synthese.sort_index()



def affichage_resultats (df_synthese, dico_site):
    fig, axs = plt.subplots(4, 1, sharex=True)
    
    axs[0].plot(df_synthese.index, df_synthese.prix, label="Prix électricité (€)")
    axs[0].grid(True)
    axs[0].legend(bbox_to_anchor=(1, 1), loc=1, borderaxespad=0)
    
    axs[1].plot(df_synthese.index, df_synthese.gains_cumules, label="Gains cumulés (€)")
    axs[1].grid(True)
    axs[1].legend(bbox_to_anchor=(1, 1), loc=1, borderaxespad=0)

    axs[2].fill_between(df_synthese.index, df_synthese.p_min, df_synthese.p_max, 
                        color='blue', alpha=0.2, 
                        label="Puissance site admissible")
    axs[2].plot(df_synthese.index, [dico_site['p_racc'] for _ in df_synthese.index], 
                color='red', linestyle="dashed", label="+Pracc")
    axs[2].plot(df_synthese.index, [-dico_site['p_racc'] for _ in df_synthese.index], 
                color='red', linestyle="dashed", label="-Pracc")

    axs[2].plot(df_synthese.index, df_synthese.p_site, linestyle="solid", label="Puissance site")
    axs[2].plot(df_synthese.index, df_synthese.p_batterie, linestyle="solid", label="Puissance batterie")
    axs[2].plot(df_synthese.index, df_synthese.p_inj, linestyle="solid", label="Prod/conso intra-site")
    axs[2].plot(df_synthese.index, df_synthese.p_spilled, linestyle="solid", label="Puissance écrêtée")
    axs[2].grid(True)
    axs[2].legend(bbox_to_anchor=(1, 1), loc=1, borderaxespad=0)
    
    # vect_e_min = pd.Series(index=df_synthese.index, data=[e_min for _ in df_synthese.index])
    # vect_e_max = pd.Series(index=df_synthese.index, data=[e_max for _ in df_synthese.index])
    # axs[3].fill_between(df_synthese.index, vect_e_min, vect_e_max, color='red', alpha=0.2, label="Plage admissible pour l'énergie")
    axs[3].plot(df_synthese.index, df_synthese.energie, linestyle="solid", color='red', label="Energie")
    axs[3].grid(True)
    axs[3].legend(bbox_to_anchor=(1, 1), loc=1, borderaxespad=0)

    plt.show()








def maximisation_gains_puis_minimisation_ene(dico_site, df_chroniques, delta_t, log_to_console=True):
    # Préparation du modèle avec ses variables et ses contraintes
    m, timesteps, tau = preparer_modele(dico_site, df_chroniques, delta_t)
    
    # Fonction objectif
    objectif_gains = - (m.variables['p_site'] * df_chroniques.prix.reindex(timesteps).fillna(0) ).sum()
    objectif_p_spilled = m.variables['p_spilled'].sum()
    alpha = 0.9 * df_chroniques.prix[df_chroniques.prix>0].min()
    m.add_objective(objectif_gains + alpha * objectif_p_spilled)
    
    # C'est parti !
    status = m.solve(log_to_console=log_to_console)
    
    # On trouve quelquechose, ou pas ?
    if status != ('warning', 'infeasible'): # Solution trouvée
        return generer_dataframe_solution(m, dico_site, df_chroniques)
    else: # Problème infaisable
        return None
    

def minimisation_ene_puis_maximisation_gains(dico_site, df_chroniques, delta_t, log_to_console=True):
    # Préparation du modèle avec ses variables et ses contraintes
    m, timesteps, tau = preparer_modele(dico_site, df_chroniques, delta_t)
    
    # Fonction objectif
    objectif_gains = - (m.variables['p_site'] * df_chroniques.prix.reindex(timesteps).fillna(0) ).sum()
    objectif_p_spilled = m.variables['p_spilled'].sum()

    facteur_marge = 1.5
    alpha = facteur_marge * 2 * dico_site['p_racc']  * df_chroniques.prix.abs().max()
    m.add_objective(objectif_gains + alpha * objectif_p_spilled)
    
    # C'est parti !
    status = m.solve(log_to_console=log_to_console)
    
    # On trouve quelquechose, ou pas ?
    if status != ('warning', 'infeasible'): # Solution trouvée
        return generer_dataframe_solution(m, dico_site, df_chroniques)
    else: # Problème infaisable
        return None
    

def minimisation_ene_puis_report_injections(dico_site, df_chroniques, delta_t, log_to_console=True):
    # Préparation du modèle avec ses variables et ses contraintes
    m, timesteps, tau = preparer_modele(dico_site, df_chroniques, delta_t)
    
    # Fonction objectif
    objectif_p_spilled = m.variables['p_spilled'].sum()
    objectif_puissance_spilled_normalisee = (60/delta_t)/dico_site['p_max_batt']*m.variables['energie'].sum() 
    alpha = 1e-3
    m.add_objective(objectif_p_spilled + alpha*objectif_puissance_spilled_normalisee)
    
    # C'est parti !
    status = m.solve(log_to_console=log_to_console)
    
    # On trouve quelquechose, ou pas ?
    if status != ('warning', 'infeasible'): # Solution trouvée
        return generer_dataframe_solution(m, dico_site, df_chroniques)
    else: # Problème infaisable
        return None