{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dcf378e6",
   "metadata": {},
   "source": [
    "# Optimisation de charge de véhicules électriques par MOPSO assisté par IA\n",
    "### Par Amine AIT MOUSSA, Siham DAANOUNI, et Clément DELAMOTTE\n",
    "\n",
    "---\n",
    "\n",
    "**Objectif :** Minimiser le coût de charge, l'insatisfaction utilisateur et la tension sur le réseau.\n",
    "\n",
    "**Méthode :** Utilisation d'un algorithme MOPSO (Multi-Objective Particle Swarm Optimization), couplé à un modèle de régression (Surrogate Model) pour accélérer l'évaluation des particules."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea73e020",
   "metadata": {},
   "source": [
    "## 1. Imports and plot function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5cf6f9f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "from mopso import MOPSO \n",
    "from surrogate_handler import SurrogateHandler\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "54c5e14f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D # Nécessaire pour la 3D\n",
    "\n",
    "def plot_pareto_3d(archive, model_type:str):\n",
    "    fig = plt.figure(figsize=(12, 8))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "    # Extraction des scores depuis l'archive\n",
    "    # f_best[0] = Coût, f_best[1] = Insatisfaction, f_best[2] = Stress Réseau\n",
    "    f1 = [p.f_best[0] for p in archive] \n",
    "    f2 = [p.f_best[1] for p in archive] \n",
    "    f3 = [p.f_best[2] for p in archive] \n",
    "\n",
    "    # Création du nuage de points\n",
    "    img = ax.scatter(f1, f2, f3, c=f3, cmap='viridis', s=60, edgecolors='black')\n",
    "    \n",
    "    ax.set_xlabel('Coût (€)')\n",
    "    ax.set_ylabel('Insatisfaction (SoC manquant)')\n",
    "    ax.set_zlabel('Pic Réseau (kW)')\n",
    "    ax.set_title(f'Front de Pareto des Solutions Non-Dominées ({model_type})')\n",
    "    \n",
    "    # Barre de couleur\n",
    "    cbar = fig.colorbar(img, ax=ax, pad=0.1)\n",
    "    cbar.set_label('Intensité du Pic Réseau (kW)')\n",
    "    \n",
    "    # Sauvegarde et affichage\n",
    "    filename = f\"{model_type}_pareto_3d.png\"\n",
    "    plt.savefig(filename)\n",
    "    print(f\"Graphique sauvegardé sous : {filename}\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee8b04e4",
   "metadata": {},
   "source": [
    "## 2. SmartMOPSO class\n",
    "We created a new class that heritate the classic MOPSO, but we added :\n",
    "- in the initialization the surrogate type\n",
    "- a new implementation to the iteration function, to use the surrogate for predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f783ba32",
   "metadata": {},
   "outputs": [],
   "source": [
    "class SmartMOPSO(MOPSO):\n",
    "    def __init__(self, model_type=None, **kwargs):\n",
    "        super().__init__(**kwargs)\n",
    "        \n",
    "        # Initialize Surrogate Handler if model_type is provided\n",
    "        self.use_surrogate = (model_type is not None)\n",
    "        if self.use_surrogate:\n",
    "            self.surrogate_handler = SurrogateHandler(model_type)\n",
    "\n",
    "            # Pre-fill with initial particle data\n",
    "            for p in self.particles:\n",
    "                self.surrogate_handler.add_data(p.x, p.f_current[1])\n",
    "\n",
    "    def iterate(self):\n",
    "        train_freq = 10 # Retrain every 10 iterations\n",
    "        \n",
    "        # Check if retraining is needed\n",
    "        if self.use_surrogate and (self.t % train_freq == 0):\n",
    "            self.surrogate_handler.train()\n",
    "\n",
    "        # Determine if AI prediction should be used\n",
    "        use_ai = (self.use_surrogate and \n",
    "                  self.surrogate_handler.is_trained and \n",
    "                  self.t % train_freq != 0)\n",
    "\n",
    "        # Main loop (overriding original logic to manage control flow)\n",
    "        for t in range(self.t): \n",
    "            self.select_leader()\n",
    "            \n",
    "            for i in range(self.n):\n",
    "                # Movement \n",
    "                self.particles[i].update_velocity(self.leader.x, self.c1, self.c2, self.w)\n",
    "                self.particles[i].update_position()\n",
    "                self.particles[i].keep_boudaries(self.A_max)\n",
    "\n",
    "                if use_ai:\n",
    "                    # Fast exact calculation (f1, f3)\n",
    "                    f1 = self.particles[i].f1(self.prices)\n",
    "                    f3 = self.particles[i].f3()\n",
    "                    \n",
    "                    # Slow prediction (f2) by using Surrogate\n",
    "                    f2_pred = self.surrogate_handler.predict(self.particles[i].x)\n",
    "                    \n",
    "                    # Inject scores without running the expensive 'updating_socs'\n",
    "                    self.particles[i].f_current = [f1, f2_pred, f3]\n",
    "                    \n",
    "                else:\n",
    "                    # Standard Calculation (Slow and Exact)\n",
    "                    self.particles[i].updating_socs(self.socs, self.capacities)\n",
    "                    self.particles[i].evaluate(self.prices, self.socs, self.socs_req, self.times)\n",
    "                    \n",
    "                    # Capture data for AI training\n",
    "                    if self.use_surrogate:\n",
    "                        self.surrogate_handler.add_data(self.particles[i].x, self.particles[i].f_current[1])\n",
    "\n",
    "                self.particles[i].update_best()\n",
    "            \n",
    "            self.update_archive()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83075d35",
   "metadata": {},
   "source": [
    "## 3. Utility functions\n",
    "Manage CSV files reading and calculate physical constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a7789c8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_elec_prices(csv_file:str, sep:str=';'):\n",
    "    elec_df = pd.read_csv(filepath_or_buffer=csv_file, sep=sep, skipinitialspace=True)\n",
    "\n",
    "    # Mean of Winter and Summer of 2025 electric prices (Euros/MWh)\n",
    "    elec_mean = (elec_df['Winter 2025'].mean() + elec_df['Summer 2025'].mean())/2\n",
    "\n",
    "    # Standard variation of Winter and Summer of 2025 electric prices (Euros/MWh)\n",
    "    elec_std = (elec_df['Winter 2025'].std() + elec_df['Summer 2025'].std())/2\n",
    "\n",
    "    elec_mean = elec_mean / 1000\n",
    "    elec_std = elec_std / 1000\n",
    "\n",
    "\n",
    "    print(f'Electricity prices:\\n - Mean: ${elec_mean}€/Mwh\\n - Std:  ${elec_std}€/Mwh')\n",
    "    return elec_mean, elec_std\n",
    "\n",
    "def generate_capacities(csv_file:str, nb_vehicles:int, seed:int=42, sep:str=';'):\n",
    "    cap_df = pd.read_csv(filepath_or_buffer=csv_file, sep=sep)\n",
    "\n",
    "    # Getting back all kind of battery capacities with unique values\n",
    "    all_capacities = cap_df['Battery Capacity kwh'].dropna().unique()\n",
    "\n",
    "    # Extracting random values for generating the array of capacities\n",
    "    capacities = pd.Series(all_capacities).sample(n=nb_vehicles, random_state=seed)\n",
    "\n",
    "    print(f'Capacities of vehicles (kwh): ${capacities}')\n",
    "    return capacities.tolist()\n",
    "\n",
    "def get_power_constants(nb_vehicles:int, nb_consumers:int=67000000):\n",
    "    mean_consumption = (87028 + 46847 + 52374 + 29819)/4 # Mean of consumption in France in 2025 (estimate according to data/grid_capacity.txt)\n",
    "    sim_ratio = nb_vehicles / nb_consumers # Ratio to reduce A_max of simulation to realistic restrictions\n",
    "\n",
    "    a_max = sim_ratio * mean_consumption\n",
    "    x_max = a_max / nb_vehicles # For init, uniform charging/discharging for every vehicle\n",
    "    x_min = -x_max\n",
    "    return a_max, x_max, x_min"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b510ed",
   "metadata": {},
   "source": [
    "## 4. Execution Function\n",
    "Execute the pipeline on a scenario with simulation parameters by using the SmartMOPSO class for the classic MOPSO or for MOPSO and Surrogate. It does the calculation for the functions f1 and f3, and calculate or predict by AI the f2 function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a8797755",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_scenario(scenario_name, capacities:list, price_mean:float, price_std:float, model_type=None, n:int=20, t:int=30, w:float=0.4, c1:float=0.3, c2:float=0.2, archive_size:int=10, nb_vehicles:int=10, delta_t:int=60, nb_of_ticks:int=48):\n",
    "    A_MAX, X_MAX, X_MIN = get_power_constants(nb_vehicles=nb_vehicles)\n",
    "\n",
    "\n",
    "    print(f\"\\n--- Launching Scenario: {scenario_name} ---\")\n",
    "    # Simulation parameters\n",
    "    params = {\n",
    "        'A_max': A_MAX, 'price_mean': price_mean, 'price_std': price_std,\n",
    "        'capacities': capacities, 'n': n, 't': t, \n",
    "        'w': w, 'c1': c1, 'c2': c2,\n",
    "        'nb_vehicles': nb_vehicles, 'delta_t': delta_t, 'nb_of_ticks': nb_of_ticks,\n",
    "        'x_min':X_MIN, 'x_max':X_MAX\n",
    "    }\n",
    "    \n",
    "    # Instantiate extended class\n",
    "    optimizer = SmartMOPSO(model_type=model_type, **params)\n",
    "    \n",
    "    start_time = time.time()\n",
    "    \n",
    "    # Run simulation\n",
    "    optimizer.iterate()\n",
    "        \n",
    "    end_time = time.time()\n",
    "    duration = end_time - start_time\n",
    "    \n",
    "    # Retrieve best f2 (e.g. from archive)\n",
    "    best_f2 = min([p.f_best[1] for p in optimizer.archive]) if optimizer.archive else 0\n",
    "    \n",
    "    print(f\"Finished in {duration:.2f} seconds.\")\n",
    "    print(f\"Best f2 found: {best_f2:.4f}\")\n",
    "    \n",
    "    return duration, best_f2, optimizer.archive"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdbd5f06",
   "metadata": {},
   "source": [
    "## 5. Scenario comparison\n",
    "\n",
    "We use here 3 configurations to evaluate our surrogate models :\n",
    "1.  **No AI (Baseline) :** Exact calculation, without any approximation.\n",
    "2.  **MLP (Multi-Layer Perceptron) :** Using a neural network to approximate $f_2$.\n",
    "3.  **RF (Random Forest) :** Using random forest to approximate $f_2$.\n",
    "\n",
    "The goal is to compare the **execution time** and **solutions quality**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6cdd4953",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Electricity prices:\n",
      " - Mean: $0.08469423999999999€/Mwh\n",
      " - Std:  $0.043112482841363334€/Mwh\n",
      "Capacities of vehicles (kwh): $33     83.0\n",
      "36     20.0\n",
      "4      72.8\n",
      "13     35.8\n",
      "30     77.4\n",
      "26     17.3\n",
      "6      65.0\n",
      "27     11.4\n",
      "24     82.0\n",
      "15     87.0\n",
      "17     93.4\n",
      "8      65.5\n",
      "16     79.2\n",
      "12     58.0\n",
      "19     26.0\n",
      "9      39.2\n",
      "32     90.0\n",
      "0      95.0\n",
      "25    200.0\n",
      "5      75.0\n",
      "dtype: float64\n",
      "\n",
      "--- Launching Scenario: Only MOPSO ---\n",
      "Finished in 0.40 seconds.\n",
      "Best f2 found: 5.3493\n",
      "\n",
      "--- Launching Scenario: With MLP ---\n",
      "Finished in 0.68 seconds.\n",
      "Best f2 found: 6.3994\n",
      "\n",
      "--- Launching Scenario: With Random Forest ---\n",
      "Finished in 0.52 seconds.\n",
      "Best f2 found: 5.2397\n",
      "\n",
      "--- Launching Scenario: Only MOPSO ---\n",
      "Finished in 0.96 seconds.\n",
      "Best f2 found: 4.5495\n",
      "\n",
      "--- Launching Scenario: With MLP ---\n",
      "Finished in 1.30 seconds.\n",
      "Best f2 found: 5.3998\n",
      "\n",
      "--- Launching Scenario: With Random Forest ---\n",
      "Finished in 1.12 seconds.\n",
      "Best f2 found: 4.4896\n",
      "\n",
      "--- Launching Scenario: Only MOPSO ---\n",
      "Finished in 9.63 seconds.\n",
      "Best f2 found: 5.4597\n",
      "\n",
      "--- Launching Scenario: With MLP ---\n",
      "Finished in 10.31 seconds.\n",
      "Best f2 found: 6.9999\n",
      "\n",
      "--- Launching Scenario: With Random Forest ---\n",
      "Finished in 9.77 seconds.\n",
      "Best f2 found: 5.7496\n",
      "\n",
      "--- Launching Scenario: Only MOPSO ---\n",
      "Finished in 18.79 seconds.\n",
      "Best f2 found: 6.6090\n",
      "\n",
      "--- Launching Scenario: With MLP ---\n",
      "Finished in 20.60 seconds.\n",
      "Best f2 found: 5.8495\n",
      "\n",
      "--- Launching Scenario: With Random Forest ---\n",
      "Finished in 22.21 seconds.\n",
      "Best f2 found: 3.4493\n",
      "\n",
      "--- Launching Scenario: Only MOPSO ---\n",
      "Finished in 196.01 seconds.\n",
      "Best f2 found: 5.0795\n",
      "\n",
      "--- Launching Scenario: With MLP ---\n",
      "Finished in 206.84 seconds.\n",
      "Best f2 found: 5.2493\n",
      "\n",
      "--- Launching Scenario: With Random Forest ---\n",
      "Finished in 202.25 seconds.\n",
      "Best f2 found: 5.4393\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# CSV files\n",
    "elec_price_csv = 'data/elec_prices.csv'\n",
    "capacity_csv = 'data/vehicle_capacity.csv'\n",
    "\n",
    "# Global Simulation parameters\n",
    "T = 30 # Number of iterations (for the particles)\n",
    "W = 0.4 # Inertia (for exploration)\n",
    "C1 = 0.3 # Individual trust\n",
    "C2 = 0.2 # Social trust\n",
    "ARC_SIZE = 10 # Archive size\n",
    "nb_vehicle = 20\n",
    "\n",
    "P_MEAN, P_STD = calculate_elec_prices(elec_price_csv)\n",
    "CAPACITIES = generate_capacities(capacity_csv, nb_vehicles=nb_vehicle)\n",
    "\n",
    "NB_TICKS = 48\n",
    "DELTA = 60\n",
    "\n",
    "results = {\n",
    "    'MOPSO':[],\n",
    "    'MLP': [],\n",
    "    'RF': []\n",
    "}\n",
    "nb_particles = [20,50,500,1000,10000]\n",
    "\n",
    "for k in range(len(nb_particles)):\n",
    "    # 1. Without Surrogate (Baseline)\n",
    "    d1, f1_score, _ = run_scenario(\n",
    "        \"Only MOPSO\", \n",
    "        capacities=CAPACITIES, \n",
    "        price_mean=P_MEAN, \n",
    "        price_std=P_STD, \n",
    "        nb_vehicles=nb_vehicle,  # Important pour la cohérence\n",
    "        model_type=None,\n",
    "        n=nb_particles[k]\n",
    "    )\n",
    "    results['MOPSO'].append((d1, f1_score))\n",
    "\n",
    "    # 2. With MLP\n",
    "    d2, f2_score, _ = run_scenario(\n",
    "        \"With MLP\", \n",
    "        capacities=CAPACITIES, \n",
    "        price_mean=P_MEAN, \n",
    "        price_std=P_STD, \n",
    "        nb_vehicles=nb_vehicle,\n",
    "        model_type='mlp',\n",
    "        n=nb_particles[k]\n",
    "    )\n",
    "    results['MLP'].append((d2, f2_score))\n",
    "\n",
    "    # 3. With Random Forest\n",
    "    d3, f3_score, _ = run_scenario(\n",
    "        \"With Random Forest\", \n",
    "        capacities=CAPACITIES, \n",
    "        price_mean=P_MEAN, \n",
    "        price_std=P_STD, \n",
    "        nb_vehicles=nb_vehicle,\n",
    "        model_type='rf',\n",
    "        n=nb_particles[k]\n",
    "    )\n",
    "    results['RF'].append((d3, f3_score))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "165b48cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== SUMMARY ===\n",
      "Mode            | Time (s)   | Best f2   \n",
      "---------------------------------------------\n",
      "MOPSO          _20              | 0.40       | 5.3493    \n",
      "MOPSO          _50              | 0.96       | 4.5495    \n",
      "MOPSO          _500             | 9.63       | 5.4597    \n",
      "MOPSO          _1000            | 18.79      | 6.6090    \n",
      "MOPSO          _10000           | 196.01     | 5.0795    \n",
      "MLP            _20              | 0.68       | 6.3994    \n",
      "MLP            _50              | 1.30       | 5.3998    \n",
      "MLP            _500             | 10.31      | 6.9999    \n",
      "MLP            _1000            | 20.60      | 5.8495    \n",
      "MLP            _10000           | 206.84     | 5.2493    \n",
      "RF             _20              | 0.52       | 5.2397    \n",
      "RF             _50              | 1.12       | 4.4896    \n",
      "RF             _500             | 9.77       | 5.7496    \n",
      "RF             _1000            | 22.21      | 3.4493    \n",
      "RF             _10000           | 202.25     | 5.4393    \n"
     ]
    }
   ],
   "source": [
    "# --- DISPLAY RESULTS ---\n",
    "print(\"\\n=== SUMMARY ===\")\n",
    "print(f\"{'Mode':<15} | {'Time (s)':<10} | {'Best f2':<10}\")\n",
    "print(\"-\" * 45)\n",
    "for k, v in results.items():\n",
    "    for i in range(len(nb_particles)):\n",
    "        print(f\"{k:<15}_{nb_particles[i]:<15} | {v[i][0]:<10.2f} | {v[i][1]:<10.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "314f8d5a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 100140 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def plot_time_benchmark(nb_particles_list, results_dict):\n",
    "    \n",
    "    t_mopso = [item[0] for item in results_dict['MOPSO']]\n",
    "    t_mlp   = [item[0] for item in results_dict['MLP']]\n",
    "    t_rf    = [item[0] for item in results_dict['RF']]\n",
    "\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    \n",
    "    plt.plot(nb_particles_list, t_mopso, 'o-', label='Sans IA (MOPSO)', color='#1f77b4', linewidth=2)\n",
    "    plt.plot(nb_particles_list, t_mlp,   's--', label='Avec MLP', color='#ff7f0e', linewidth=2)\n",
    "    plt.plot(nb_particles_list, t_rf,    '^-.', label='Avec Random Forest', color='#2ca02c', linewidth=2)\n",
    "    \n",
    "    plt.title(\"Temps d'exécution selon le nombre de particules\", fontsize=14, fontweight='bold')\n",
    "    plt.xlabel(\"Nombre de Particules\", fontsize=12)\n",
    "    plt.ylabel(\"Temps (s)\", fontsize=12)\n",
    "    plt.grid(True, linestyle=':', alpha=0.7)\n",
    "    plt.legend(fontsize=11)\n",
    "    \n",
    "\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "    \n",
    "plot_time_benchmark(nb_particles, results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1e3fa2ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 118980 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_f2_benchmark(nb_particles_list, results_dict):\n",
    "    s_mopso = [item[1] for item in results_dict['MOPSO']]\n",
    "    s_mlp   = [item[1] for item in results_dict['MLP']]\n",
    "    s_rf    = [item[1] for item in results_dict['RF']]\n",
    "\n",
    "\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    \n",
    "    plt.plot(nb_particles_list, s_mopso, 'o-', label='Sans IA (MOPSO)', color='#1f77b4', linewidth=2)\n",
    "    plt.plot(nb_particles_list, s_mlp,   's--', label='Avec MLP', color='#ff7f0e', linewidth=2)\n",
    "    plt.plot(nb_particles_list, s_rf,    '^-.', label='Avec Random Forest', color='#2ca02c', linewidth=2)\n",
    "    \n",
    "    plt.title(\"Meilleur Score F2 (Convergence) selon le nombre de particules\", fontsize=14, fontweight='bold')\n",
    "    plt.xlabel(\"Nombre de Particules (log scale)\", fontsize=12)\n",
    "    plt.ylabel(\"Meilleur F2 Score\", fontsize=12)\n",
    "    plt.grid(True, linestyle=':', alpha=0.7)\n",
    "    plt.legend(fontsize=11)\n",
    "    \n",
    "    plt.xscale('log') \n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_f2_benchmark(nb_particles, results)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Optim_Metaheuristique (3.11.14)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
