examples/Microgrid_sizing_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b14f854f-51e9-4253-8946-d7ca39628ed0",
   "metadata": {},
   "source": [
    "# Microgrid sizing optimization\n",
    "\n",
    "Following the first example about [Microgrid performance simulation](Microgrid_Wind-Solar.ipynb)\n",
    "using [Microgrids.jl](https://github.com/Microgrids-X/Microgrids.jl),\n",
    "this notebook delves into **optimal microgrid sizing**.\n",
    "This means finding the optimal size (rated power or energy) of each component by solving an optimization problem where the objective function is composed of Key Performance Indicators (KPIs) related to the economic and technical performance of the project.\n",
    "\n",
    "Like in the previous example, the sizing is done for a microgrid project with *wind* and *solar* sources, plus a *battery* and a *dispatchable generator*.\n",
    "\n",
    "<img alt='schematic of a microgrid with a symbolic representation of the notion of sizing' src='./images/microgrid_sizing.png' style='height:18em'>\n",
    "\n",
    "The main steps of the optimal sizing process are:\n",
    "1. set up the problem (in particular loading project data like prices and time series)\n",
    "2. run the optimization\n",
    "3. analyze the results\n",
    "\n",
    "## About the optimization formulation and process\n",
    "\n",
    "### Optimization formulation\n",
    "\n",
    "In this illustrative sizing example, the optimization objective (criterion) is:\n",
    "- minimize the Levelized Costs of Electricity `lcoe` (€/kWh)\n",
    "- while keeping the load shedding rate below some threshold (ex.: shedding rate ≤ 1%)\n",
    "\n",
    "but other indicators could be used as well: rate of renewables above some threshold...\n",
    "\n",
    "The above problem is a constrained optimization problem. Since we'll use unconstrained optimization solvers, we will reformulate it as: objective = lcoe + penalty for excess of shedding.\n",
    "\n",
    "### Optimization process\n",
    "\n",
    "This is a “blackbox” type of optimization problem: the optimization algorithm calls iteratively the Microgrid simulator with sizings to be tested until convergence (or reaching the maximum allowed number of iterations):\n",
    "\n",
    "<img alt='schematic of microgrid sizing optimization process: iterations call of the simulator by optimizer' src='./images/MG_optim_diagram.png' style='height:18em'>\n",
    "\n",
    "The code is meant for using the [NLopt](https://nlopt.readthedocs.io/en/latest/) library of optimization solvers (through its [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl/) wrapper), but it could be adapted to other optimizers with minor changes (see remark below about the modularity of the code)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c2f04250-290c-4c38-93e6-074702ec0700",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Microgrids\n",
    "using NLopt # optimization solvers\n",
    "using Printf # pretty print results\n",
    "using BenchmarkTools # for timing performance, optional"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e20bad8-81a3-4a6b-8c1e-6848a20c0d55",
   "metadata": {},
   "source": [
    "## Load Microgrid project data\n",
    "\n",
    "Loading parameters and time series for a Microgrid project with *wind* and *solar* sources, plus a *battery* and a *dispatchable generator*. \n",
    "Values gathered from the [Microgrid_Wind-Solar.ipynb](Microgrid_Wind-Solar.ipynb) notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a078bb9d-37a9-45f0-bb00-4404a44dcb30",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data definition for Microgrid with wind, solar, storage and generator...\n"
     ]
    }
   ],
   "source": [
    "include(\"./Microgrid_Wind-Solar_data.jl\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9431f37f-4a44-42f7-a017-4935125731f2",
   "metadata": {},
   "source": [
    "## Setting up the cost function (criterion) to be optimized\n",
    "\n",
    "The key coding objective is to **encapsulate** the microgrid simulator (`simulate` function of `Microgrids.jl` package) into an objective function that can be called by the optimization algorithm, that is which respects its expected interface (here NLopt).\n",
    "\n",
    "To increase the modularity which facilitates using optimization solvers others that NLopt's we implement the encapsulation by **3 nested functions**:\n",
    "\n",
    "1. Simulation of Microgrid project described by a sizing vector `x` (vector input) → returns all simulation statistics\n",
    "2. Extract KPIs of interest to build a multi-objective criterion: here lcoe and shedding rate\n",
    "3. Combine these KPIs as one mono-objective criterion: here LCOE + penalty if shedding rate > shed_max\n",
    "   - and match the interface expected by NLopt's optimizers\n",
    "\n",
    "but if one cares more about compactness, this could be assembled into one single function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adb0b6bb-54f7-4849-b1b0-b236681be8d9",
   "metadata": {},
   "source": [
    "### Wrapper of the Microgrid simulator\n",
    "\n",
    "accept a sizing vector `x`, then create all `Microgrids.jl` components and call `simulate`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b1b119b9-0e9c-447d-b708-ef64e11b0d0d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "simulate_microgrid"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Simulate the performance of a Microgrid project of size `x`\n",
    "with x=[power_rated_gen, energy_rated_sto, power_rated_pv, power_rated_wind]\n",
    "\n",
    "Returns stats, costs\n",
    "\"\"\"\n",
    "    function simulate_microgrid(x)\n",
    "    project = Project(lifetime, discount_rate, timestep, \"€\")\n",
    "    # Split decision variables (converted MW → kW):\n",
    "    power_rated_gen = x[1]*1000\n",
    "    energy_rated_sto = x[2]*1000\n",
    "    power_rated_pv = x[3]*1000\n",
    "    power_rated_wind = x[4]*1000\n",
    "\n",
    "    # Create components\n",
    "    gen = DispatchableGenerator(power_rated_gen,\n",
    "        fuel_intercept, fuel_slope, fuel_price,\n",
    "        investment_price_gen, om_price_gen, lifetime_gen,\n",
    "        load_ratio_min,\n",
    "        replacement_price_ratio, salvage_price_ratio, fuel_unit)\n",
    "    batt = Battery(energy_rated_sto,\n",
    "        investment_price_sto, om_price_sto, lifetime_sto, lifetime_cycles,\n",
    "        charge_rate, discharge_rate, loss_factor_sto, SoC_min, SoC_ini,\n",
    "        replacement_price_ratio, salvage_price_ratio)\n",
    "    pv = Photovoltaic(power_rated_pv, irradiance,\n",
    "        investment_price_pv, om_price_pv,\n",
    "        lifetime_pv, derating_factor_pv,\n",
    "        replacement_price_ratio, salvage_price_ratio)\n",
    "    windgen = WindPower(power_rated_wind, cf_wind,\n",
    "        investment_price_wind, om_price_wind,\n",
    "        lifetime_wind,\n",
    "        replacement_price_ratio, salvage_price_ratio)\n",
    "    mg = Microgrid(project, Pload, gen, batt, [pv, windgen])\n",
    "\n",
    "    # Launch simulation:\n",
    "    traj, stats, costs = simulate(mg)\n",
    "\n",
    "    return stats, costs\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5321132-a5bf-467b-839b-99e36a2e781c",
   "metadata": {},
   "source": [
    "Test of the simulator wrapper (on a baseline sizing):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d765dc6a-eeda-4727-b259-479176670dc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1.8, 5.0, 3.0, 0.9], 0.22924812869928668, 21.89002772908652)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Baseline sizing: same as in Microgrid_Wind-Solar.ipynb notebook\n",
    "x_base = [power_rated_gen, energy_rated_sto, power_rated_pv, power_rated_wind]/1000.\n",
    "# run simulation:\n",
    "stats, costs = simulate_microgrid(x_base)\n",
    "x_base, costs.lcoe, costs.npc/1e6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e51a19f-0f68-44e4-ae09-ad11806f5747",
   "metadata": {},
   "source": [
    "### Build the objective functions (criteria)\n",
    "\n",
    "- first bi-objective function x ↦ (lcoe, shedding rate)(x)\n",
    "- then wrapped into a mono objective x ↦ J(x) by using a penalty for the excess of shedding rate\n",
    "  - and match the interface expected by NLopt's optimizers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "470f5630-322c-497f-bac9-beb87f014f71",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "obj_multi"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"Multi-objective criterion for microgrid performance: lcoe, shedding rate\"\n",
    "function obj_multi(x)\n",
    "    stats, costs = simulate_microgrid(x)\n",
    "    # Extract KPIs of interest\n",
    "    lcoe = costs.lcoe # $/kWh\n",
    "    shed_rate = stats.shed_rate # in [0,1]\n",
    "    return lcoe, shed_rate\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6b3eb2d2-1513-424b-baa0-a49ca5247bb4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "obj"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Mono-objective criterion: LCOE + penalty if shedding rate > `shed_max`\n",
    "\n",
    "with signature adapted to NLopt with `grad` as 2nd argument\n",
    "\n",
    "load shedding penalty threshold `shed_max` should be in [0,1[\n",
    "\"\"\"\n",
    "function obj(x, grad, shed_max, w_shed_max=1e5)\n",
    "    lcoe, shed_rate = obj_multi(x)\n",
    "    over_shed = shed_rate - shed_max\n",
    "    if over_shed > 0.0\n",
    "        penalty = w_shed_max*over_shed\n",
    "    else\n",
    "        penalty = 0.0\n",
    "    end\n",
    "    J = lcoe + penalty\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c57f7af9-ea2f-40a8-b502-ba98b472410c",
   "metadata": {},
   "source": [
    "### Tests the objective functions\n",
    "\n",
    "Sizing being tested:\n",
    "- baseline sizing from the simulation notebook: perfect quality of service (QoS) with zero load shedding\n",
    "- baseline modified with a halved generator sizing: very good QoS with a bit of load shedding → not penalized\n",
    "- small PV and small wind generators only: low LCOE (i.e. the production-only LCOE of these sources) but but extremely bad QoS → huge penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8de2153d-559f-4a60-880e-e8ea5fa3b810",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Base. multi: (0.22924812869928668, 0.0) mono: 0.22924812869928668\n",
      "Shed. multi: (0.20671000160924613, 0.009602858175478355) mono: 0.20671000160924613\n",
      "PV.   multi: (0.10149685980966677, 0.923547868561659) mono: 91354.8883530257\n",
      "Wind. multi: (0.10040264224635914, 0.74395737102815) mono: 73395.83750545724\n"
     ]
    }
   ],
   "source": [
    "# Test:\n",
    "shed_max = 0.01 # in [0,1]\n",
    "\n",
    "x_shed = [power_rated_gen/2, energy_rated_sto, power_rated_pv, power_rated_wind]/1000.\n",
    "x_pv   = [0. 0. 500.   0.]/1000.\n",
    "x_wind = [0. 0.   0. 500.]/1000.\n",
    "\n",
    "println(\"Base. multi: \", obj_multi(x_base), \" mono: \", obj(x_base,[], shed_max))\n",
    "println(\"Shed. multi: \", obj_multi(x_shed), \" mono: \", obj(x_shed,[], shed_max))\n",
    "println(\"PV.   multi: \", obj_multi(x_pv), \" mono: \", obj(x_pv,[], shed_max))\n",
    "println(\"Wind. multi: \", obj_multi(x_wind), \" mono: \", obj(x_wind,[], shed_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfbf1ccf-a7df-474a-991c-9437bd83d8ab",
   "metadata": {},
   "source": [
    "## Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15efd32-c2ab-4d9a-98cf-5918cdf0b9b1",
   "metadata": {},
   "source": [
    "### Setting up the optimization problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87ccea54-0d83-4b01-8f7e-de37c8dfa0d9",
   "metadata": {},
   "source": [
    "bounds of the design space and starting point: derived from maximal load power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5ac92d7e-5de1-4ddf-a050-6897dbc8af2d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  2.0484\n",
       " 17.07\n",
       " 17.07\n",
       "  8.535"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pload_max = maximum(Pload)\n",
    "\n",
    "xmin = [0., 0., 1e-3, 0.] # 1e-3 instead of 0.0, because LCOE is NaN if ther is exactly zero generation\n",
    "x0 = [1.0, 3.0, 3.0, 2.0] * (Pload_max/1000)\n",
    "xmax = [1.2, 10.0, 10.0, 5.0] * (Pload_max/1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c2bc984-bae8-403f-8130-d50efc8f5717",
   "metadata": {},
   "source": [
    "Optionally ban some choices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9f65e094-4c58-4e9b-b4ee-b626f6504f6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solar power forbidden (optional)\n",
    "#x0[3] = 1e-3\n",
    "#xmax[3] = 1e-3\n",
    "# Wind power forbidden (optional)\n",
    "#x0[4] = 0.\n",
    "#xmax[4] = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1bd17b9-c02d-4280-8715-c46efe6c03f4",
   "metadata": {},
   "source": [
    "Check cost function on `xmin` and `xmax`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6d5ef67e-5d40-4bb0-a13d-65dc3d57bceb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((0.10149685980963595, 0.9998470957371233), (0.8229416738277807, 0.0))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obj_multi(xmin), obj_multi(xmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3daceed1-2315-4eca-b7bd-b3a1e31cff09",
   "metadata": {},
   "source": [
    "### Wrapper of the optimization process\n",
    "\n",
    "This is an optional step, but recommended to explore easily the impact of the many parameters taken by optimization algorithms.\n",
    "\n",
    "See [optimization termination conditions](https://nlopt.readthedocs.io/en/latest/NLopt_Introduction/#termination-conditions) in NLopt documention for the meaning of `xtol_rel`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "08491363-e1ce-47d3-a425-4bc83888be12",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "optim_mg"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Optimize sizing of microgrid based on the `obj` function\n",
    "\n",
    "Parameters:\n",
    "- `x0`: initial sizing (for the algorithms which need them)\n",
    "- `shed_max`: load shedding penalty threshold (same as in `obj`)\n",
    "- `algo` could be one of LN_SBPLX, GN_DIRECT, GN_ESCH...\n",
    "- `maxeval`: maximum allowed number of calls to the objective function,\n",
    "  that is to the microgrid simulation\n",
    "- `xtol_rel`: termination condition based on relative change of sizing, see NLopt doc.\n",
    "- `srand`: random number generation seed (for algorithms which use some stochastic search)\n",
    "\n",
    "Problem bounds are taken as the global variables `xmin`, `xmax`,\n",
    "but could be added to the parameters as well.\n",
    "\"\"\"\n",
    "function optim_mg(x0, shed_max, algo=:GN_CRS2_LM, maxeval=1000, xtol_rel=1e-4, srand=1)\n",
    "    nx = length(x0) # number of optim variables\n",
    "    opt = Opt(algo, nx)\n",
    "    NLopt.srand(srand)\n",
    "    \n",
    "    opt.lower_bounds = xmin\n",
    "    opt.upper_bounds = xmax\n",
    "    opt.min_objective = (x, grad) -> obj(x, grad, shed_max)\n",
    "    opt.xtol_rel = xtol_rel\n",
    "    opt.maxeval = maxeval\n",
    "    \n",
    "    (fopt, xopt, ret) = optimize(opt, x0)\n",
    "    return xopt, ret, opt.numevals\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11264431-19b1-4349-bc90-19d88bba1b87",
   "metadata": {},
   "source": [
    "### Run optimization & analyze results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d32fcbf4-42fc-42aa-a4ad-1b8ca24a984d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GN_CRS2_LM algo: MAXEVAL_REACHED after 1000 iterations. \n",
      "x*=[889.0, 2013.0, 1570.9, 1354.0]\n",
      "LCOE*: 0.19072074165935016\n",
      "shed*: 0.009961277810994507\n"
     ]
    }
   ],
   "source": [
    "algo = :GN_CRS2_LM # could be one of GN_CRS2_LM, GN_DIRECT, GN_ESCH, LN_SBPLX...\n",
    "shed_max = 0.01 # in [0,1]\n",
    "maxeval=1000\n",
    "xopt, ret, numevals = optim_mg(x0, shed_max, algo, maxeval)\n",
    "\n",
    "@printf(\"%s algo: %s after %d iterations. \\nx*=\", algo, ret, numevals)\n",
    "println(round.(xopt*1000; digits=1)) # kW\n",
    "lcoe_opt, shed_rate_opt = obj_multi(xopt)\n",
    "println(\"LCOE*: \", lcoe_opt)\n",
    "println(\"shed*: \", shed_rate_opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58772966-74b1-431c-9f84-e8efe7352130",
   "metadata": {},
   "source": [
    "optional local \"polishing\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "78b3e3f8-0de8-47f3-a8a4-40358feef172",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LN_SBPLX polish: XTOL_REACHED after 260 iterations. \n",
      "x*=[888.5, 1881.3, 1570.9, 1354.0]\n",
      "LCOE*: 0.1906609830904905\n",
      "shed*: 0.009999998875359207\n"
     ]
    }
   ],
   "source": [
    "algo_polish = :LN_SBPLX\n",
    "xopt_polish, ret, numevals = optim_mg(xopt, shed_max, algo_polish)\n",
    "\n",
    "@printf(\"%s polish: %s after %d iterations. \\nx*=\", algo_polish, ret, numevals)\n",
    "println(round.(xopt_polish*1000; digits=1)) # kW\n",
    "lcoe_opt, shed_rate_opt = obj_multi(xopt_polish)\n",
    "println(\"LCOE*: \", lcoe_opt)\n",
    "println(\"shed*: \", shed_rate_opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f543939c-c26d-40fb-aefa-129dff57ae1d",
   "metadata": {},
   "source": [
    "Retrieve all performance statistics of the optimized sizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1e046d9e-21ec-4834-9bf3-5bc1c9be1f0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NPC*: 18.03 M$ (compared to 21.89 M$ in baseline)\n",
      "rate of renewables: 75.3 %\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "OperationStats with fields:\n",
       "- served_energy: 6.7075e6 kWh\n",
       "- shed_energy: 67487.0 kWh\n",
       "- shed_max: 680.56 kW\n",
       "- shed_hours: 370.0 h\n",
       "- shed_duration_max: 21.0 h\n",
       "- shed_rate: 0.0099613 in [0,1]\n",
       "- gen_energy: 1.6592e6 kWh\n",
       "- gen_hours: 3462.0 h\n",
       "- gen_fuel: 398200.0 L\n",
       "- storage_cycles: 190.11 \n",
       "- storage_char_energy: 401830.0 kWh\n",
       "- storage_dis_energy: 363560.0 kWh\n",
       "- storage_loss_energy: 38269.0 kWh\n",
       "- spilled_energy: 1.2495e6 kWh\n",
       "- spilled_max: 1941.1 kW\n",
       "- spilled_rate: 0.1972 in [0,1]\n",
       "- renew_potential_energy: 6.3361e6 kWh\n",
       "- renew_energy: 5.0866e6 kWh\n",
       "- renew_rate: 0.75264 in [0,1]\n"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats_opt, costs_opt = simulate_microgrid(xopt)\n",
    "@printf(\"NPC*: %.2f M\\$ (compared to %.2f M\\$ in baseline)\\n\", costs_opt.npc/1e6, costs.npc/1e6)\n",
    "@printf(\"rate of renewables: %.1f %%\\n\", stats_opt.renew_rate*100)\n",
    "# Display all operation statistics\n",
    "stats_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "464f74e6-2939-4bed-ac6f-87ec8ef99660",
   "metadata": {},
   "source": [
    "#### Summary of some results for 1% max shedding rate\n",
    "\n",
    "Global search algorithms:\n",
    "- **CRS2**: 0.190 720 \\$/kWh (with shedding ~ threshold).\n",
    "  x* = [889.0, 2013.0, 1570.9, 1354.0] → *best global algorithm choice so far, although effect of seed not studied*\n",
    "- **DIRECT**: 0.191 297 \\$/kWh (with shedding < threshold) in 1000 iter.\n",
    "  x* = [910.4, 1837.8, 1543.3, 1375.7] \n",
    "    - \"polishing\" with SBPLX: 0.190 650 \\$/kWh (with shedding ~ threshold) in 207 extra iter.\n",
    "      x* = [887.1, 1992.4, 1543.3, 1375.6]\n",
    "    - increasing DIRECT maxeval to 10k doesn't bring much improvement (0.191 258 \\$/kWh)\n",
    "      compared to SBPLX polishing\n",
    "- **ESCH**: 0.194 884 \\$/kWh with shedding OK (with shedding < threshold) in 1000 iter.\n",
    "  x* = [1028.1, 1769.6, 1420.5, 1365.1] → *much slower to converge than DIRECT or CRS*\n",
    "\n",
    "Local search algorithms:\n",
    "- **SUBPLX**: 0.190 630 \\$/kWh with (with shedding ~ threshold) **in 479 iter** (<maxeval).\n",
    "  x* = [888.6, 1873.5, 1558.5, 1355.9]\n",
    "  - very competitive, but too sensitive to initialization `x0`?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fec0632-2c4b-4683-8d53-c8062fd238d7",
   "metadata": {},
   "source": [
    "#### Case with perfect quality of service (QoS)\n",
    "\n",
    "Another optimization case with perfect quality of service (QoS), that is zero load shedding  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8f4a08f2-96d5-47f1-a9fb-41a419051be8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GN_CRS2_LM algo: MAXEVAL_REACHED after 1002 iterations. \n",
      "x*=[1561.8, 2821.2, 1758.3, 1464.5]\n",
      "LCOE*: 0.2074272180878332\n",
      "shed*: 0.0\n"
     ]
    }
   ],
   "source": [
    "algo = :GN_CRS2_LM # could be one of LN_SBPLX, GN_DIRECT, GN_CRS2_LM, GN_ESCH...\n",
    "shed_max = 0.00 # in [0,1]\n",
    "maxeval=1000\n",
    "xopt, ret, numevals = optim_mg(x0, shed_max, algo, maxeval)\n",
    "\n",
    "@printf(\"%s algo: %s after %d iterations. \\nx*=\", algo, ret, numevals)\n",
    "println(round.(xopt*1000; digits=1)) # kW\n",
    "lcoe_opt, shed_rate_opt = obj_multi(xopt)\n",
    "println(\"LCOE*: \", lcoe_opt)\n",
    "println(\"shed*: \", shed_rate_opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d4a8f44-eebd-49cd-bcf3-2e44fd3372c7",
   "metadata": {},
   "source": [
    "#### Performance (duration) of the optimization\n",
    "\n",
    "About 2 s for 10k iterations, i.e. 0.2 ms/iter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9dc290cb-ac52-44db-8ee3-c54dce0fe4be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.021157 seconds (1.50 M allocations: 7.882 GiB, 3.09% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([1.5931999999999993, 2.586127114769105, 1.7946595138202435, 1.422494440730462], :MAXEVAL_REACHED, 10000)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time xopt, ret, numevals = optim_mg(x0, shed_max, :GN_DIRECT, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce32248d-39de-448a-829c-376e1eaaa550",
   "metadata": {},
   "source": [
    "Fact: optimization time relates to simulation time with Toptim ~ Tsim × neval\n",
    "\n",
    "[Benchmarktools](https://juliaci.github.io/BenchmarkTools.jl/)'s `@btime` is need for accurate timing of the simulator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "16330858-ca64-45e6-b641-d197924e401d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  168.679 μs (147 allocations: 826.34 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.31411996962098027, 0.0)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime obj_multi(x0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df6cf8f-a878-4718-b8f8-a90a62aa9aee",
   "metadata": {},
   "source": [
    "## Optimization with Automatic Differentiation *(TO BE WRITTEN)*\n",
    "\n",
    "As demonstrated in our conference article:\n",
    "\n",
    "de Godoy Antunes, P. Haessig, C. Wang, and R. Chouhy Leborgne, “Optimal Microgrid Sizing using Gradient-based Algorithms with Automatic Differentiation,” ISGT Europe 2022, Novi Sad, Serbia, 2022. https://dx.doi.org/10.1109/ISGT-Europe54678.2022.9960498, Archived at https://hal.archives-ouvertes.fr/hal-03370004\n",
    "\n",
    "Gradient of optimization criteria computed with [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).\n",
    "This allows using **gradient-based optimization algorithms** like LBFGS or SLSQP (the latter supports constraints) which converge in less iterations.\n",
    "\n",
    "This requires using the `smoothing` parameter (ε in the article is `Smoothing.transition`) of the Microgrid simulator to smooth out some discontinuities.\n",
    "\n",
    "DOCUMENTATION TO BE WRITTEN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5d974d16-f0fc-4f9a-a48c-1f7b2bc48786",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "search: \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1mm\u001b[22m\u001b[0m\u001b[1mu\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1me\u001b[22m \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1mm\u001b[22m\u001b[0m\u001b[1mu\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1me\u001b[22m_microgrid\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "simulate(mg::Microgrid, smoothing::Smoothing=NoSmoothing)\n",
       "\\end{verbatim}\n",
       "Simulate the technical and economic performance of a Microgrid \\texttt{mg}.\n",
       "\n",
       "Discontinuous computations can optionally be smoothed (relaxed) using the \\texttt{Smoothing} parameters \\texttt{transition} and \\texttt{gain} (see their doc). Smoothing is recommended when using gradient-based optimization.\n",
       "\n",
       "Returns:\n",
       "\n",
       "\\begin{itemize}\n",
       "\\item Operational trajectories from \\texttt{sim\\_operation} (should be optional in future version)\n",
       "\n",
       "\n",
       "\\item Operational statistics from \\texttt{sim\\_operation}\n",
       "\n",
       "\n",
       "\\item Microgrid project costs from \\texttt{sim\\_economics}\n",
       "\n",
       "\\end{itemize}\n"
      ],
      "text/markdown": [
       "```\n",
       "simulate(mg::Microgrid, smoothing::Smoothing=NoSmoothing)\n",
       "```\n",
       "\n",
       "Simulate the technical and economic performance of a Microgrid `mg`.\n",
       "\n",
       "Discontinuous computations can optionally be smoothed (relaxed) using the `Smoothing` parameters `transition` and `gain` (see their doc). Smoothing is recommended when using gradient-based optimization.\n",
       "\n",
       "Returns:\n",
       "\n",
       "  * Operational trajectories from `sim_operation` (should be optional in future version)\n",
       "  * Operational statistics from `sim_operation`\n",
       "  * Microgrid project costs from `sim_economics`\n"
      ],
      "text/plain": [
       "\u001b[36m  simulate(mg::Microgrid, smoothing::Smoothing=NoSmoothing)\u001b[39m\n",
       "\n",
       "  Simulate the technical and economic performance of a Microgrid \u001b[36mmg\u001b[39m.\n",
       "\n",
       "  Discontinuous computations can optionally be smoothed (relaxed) using the\n",
       "  \u001b[36mSmoothing\u001b[39m parameters \u001b[36mtransition\u001b[39m and \u001b[36mgain\u001b[39m (see their doc). Smoothing is\n",
       "  recommended when using gradient-based optimization.\n",
       "\n",
       "  Returns:\n",
       "\n",
       "    •  Operational trajectories from \u001b[36msim_operation\u001b[39m (should be optional in\n",
       "       future version)\n",
       "\n",
       "    •  Operational statistics from \u001b[36msim_operation\u001b[39m\n",
       "\n",
       "    •  Microgrid project costs from \u001b[36msim_economics\u001b[39m"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?simulate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b788d2-a5c7-4444-9f45-fefa8929d6fe",
   "metadata": {},
   "source": [
    "---\n",
    "## *How we imagine the optimization interface in a future version...*\n",
    "\n",
    "Desirable feature: *dynamically select* which variables are optimized (e.g. using `getproperty` to dynamically access structure fields...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2a51ac3b-4493-472b-8b0f-ac924edf5de7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Tuple{Vector{Symbol}, Float64, Float64, Float64}}:\n",
       " ([:generator, :power_rated], 1000.0, 0.0, 2000.0)\n",
       " ([:battery, :energy_max], 3000.0, 0.0, 6000.0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#mg0 = baseline microgrid project, with given price and time series data\n",
    "opti_params = [ # list of :component/variable name, initial, min, max values\n",
    "    ([:generator,:power_rated], 1000.0, 0., 2000.), \n",
    "    ([:battery,:energy_max],    3000.0, 0., 6000.)\n",
    "]\n",
    "#optim_mg(mg0, opti_params)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.4",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}