Class III - An introduction to PowerModels.jl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class III - An introduction to PowerModels.jl\n",
    "\n",
    "*Los Alamos National Laboratory Grid Science Winter School, 2019*\n",
    "\n",
    "Welcome! This tutorial will introduce you to the basics of the [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl) package. If you haven't yet, work through [Class I - An introduction  to Julia](Class%20I%20-%20An%20introduction%20to%20Julia.ipynb) and [Class II - An introduction  to JuMP](Class%20II%20-%20An%20introduction%20to%20JuMP.ipynb) first.\n",
    "\n",
    "As in Class II, run the following magic sauce to check we're good to go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Pkg\n",
    "Pkg.activate(@__DIR__)\n",
    "Pkg.instantiate()\n",
    "println(\"Excellent! Everything is good to go!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background \n",
    "\n",
    "This [presentation](https://github.com/lanl-ansi/tutorial-grid-science-2019/blob/master/assets/infrastructure_optimization_in_julia.pdf) provides motivation and context for this notebook.\n",
    "\n",
    "Some additional informaiton is also available at,\n",
    "* [PowerModels Documentation](https://lanl-ansi.github.io/PowerModels.jl/stable/)\n",
    "* [PowerModels Network Model Update](https://youtu.be/j7r4onyiNRQ)\n",
    "* [PSCC 2018 Presentation](https://youtu.be/AEEzt3IjLaM)\n",
    "* [JuMP Dev 2017 Presentation](https://youtu.be/W4LOKR7B4ts)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with the Network Model\n",
    "\n",
    "A simple 5-bus network model is provided with this tutorial.  The matpower file can be viewed [here](../edit/data/pglib_opf_case5_pjm.m).\n",
    "\n",
    "The `parse_file` function is used to load a text files into the PowerModels data model, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PowerModels\n",
    "\n",
    "data = PowerModels.parse_file(\"data/pglib_opf_case5_pjm.m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case the file parser generated a few informational of messages.  \n",
    "\n",
    "The first `info` message indicates that the file contains a data block named `area` that is not part of the Matpower data specification.  Still the parser will load this data using some standards for [Matpower data extentions](https://lanl-ansi.github.io/PowerModels.jl/stable/network-data/#Working-with-Matpower-Data-Files-1).\n",
    "\n",
    "The second block of `warn` messages indicates that PowerModels has detected that the generator cost model could be simplified.  In this case a linear function specified as a quadratic function was transformed into a linear function.  PowerModels makes a best effort to clean data while parsing, but the best practice is to clean datasets until they do not produce any warning messages in PowerModels.\n",
    "\n",
    "Let us use the silence command to ignore similar warnings for the rest of this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PowerModels.silence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets look at some of the network data.  The data for the first bus, load and branch can be viewed as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"bus\"][\"1\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"load\"][\"1\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the \"pd\" value is 3.0, while in the Matpower file it was 300.  By default all data in InfrastructureModels is in per-unit.  The functions `make_mixed_units` and `make_per_unit` can be used for switching back and forth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"branch\"][\"1\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's great, but looking at components one-by-one can get boring fast.\n",
    "\n",
    "All InfrastructureModels packages provide a `print_summary` function that prints a table-like summary of the network data to the terminal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PowerModels.print_summary(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dictionary comprehensions also provide an easy way to take slices of the data.\n",
    "\n",
    "For example lets inspect the generation profile and output capabilities using comprehensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[gen[\"pg\"] for (i,gen) in data[\"gen\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[gen[\"pg\"] + gen[\"qg\"]im for (i,gen) in data[\"gen\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[(gen[\"qmin\"], gen[\"qmax\"]) for (i,gen) in data[\"gen\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving Optimal Power Flow (OPF)\n",
    "\n",
    "Before we can solve a power model, we need a solver.  Lets start with the simple linear programming solver HiGHS and try to solve a DC OPF model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using HiGHS\n",
    "lp_solver = HiGHS.Optimizer\n",
    "\n",
    "result = solve_dc_opf(data, lp_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result object contains a variety of useful information about the optimization problem solved, including the objective value and wall clock runtime.  Detailed documentation of the result dictionary is available [here](https://lanl-ansi.github.io/PowerModels.jl/stable/result-data/), however most interesting point is the `solution` data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result[\"solution\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not every useful, lets try `print_summary`,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PowerModels.print_summary(result[\"solution\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks great!  But `vm` and `qg` are missing...  On to AC OPF!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = solve_ac_opf(data, lp_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Quiz Question\n",
    "That was ugly... What happened?!?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets try using an NLP solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Ipopt\n",
    "nlp_solver = optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0)\n",
    "\n",
    "result = solve_ac_opf(data, nlp_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It worked!  Lets look at the solution,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PowerModels.print_summary(result[\"solution\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution make sense intuitively.  The voltage magnitude on bus 3 is pushed to the upper limit, presumably to minimize line losses.  Total generation costs are similar to DC-OPF which seems reasonable.\n",
    "\n",
    "The solution seems reasonable, but is it optimal?  To test that we need a convex relaxation!  But first a brief digression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A little secret, the functions `solve_dc_opf` and `solve_ac_opf` are short hands for a much more general `solve_opf` function which supports many variants of optimal power flow problem.  Specifically,\n",
    "```\n",
    "solve_dc_opf(data, solver) => solve_opf(data, DCPPowerModel, solver)\n",
    "solve_dc_opf(data, solver) => solve_opf(data, ACPPowerModel, solver)\n",
    "```\n",
    "\n",
    "Try the following code with different PowerModels.  Suggested models include `ACPPowerModel`,`ACRPowerModel`,`DCPPowerModel`,`DCPLLPowerModel`,`LPACCPowerModel`, `SOCWRPowerModel`, `QCWRPowerModel`.  A fairly complete list of implemented models is available [here](https://lanl-ansi.github.io/PowerModels.jl/stable/formulation-details/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = solve_opf(data, ACPPowerModel, nlp_solver)\n",
    "PowerModels.print_summary(result[\"solution\"])\n",
    "println()\n",
    "println(\"Objective: $(result[\"objective\"])\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point you might be curious what is PowerModels doing under the hood.  In short, it builds and solves JuMP models.  The [PowerModelsAnnex](https://github.com/lanl-ansi/PowerModelsAnnex.jl/tree/master/src/model) includes some examples of the kind of JuMP models that PowerModels when you call `run_opf`.\n",
    "\n",
    "Another way to see what PowerModels is doing is to pass your own JuMP model, then you can inspect it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuMP\n",
    "\n",
    "m = Model()\n",
    "result = solve_opf(data, DCPPowerModel, nlp_solver, jump_model=m)\n",
    "\n",
    "println(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing Optimality Gaps\n",
    "\n",
    "Getting back to the question of AC-OPF solution quality.  An easy way to certify AC-OPF solution quality is to use a convex-relaxation to compute a lower bound.  If the \"gap\" between the solution and the lower bound is small we are happy.\n",
    "\n",
    "Let build a helper function for computing optimality gaps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function soc_gap(data_file, nlp_solver)\n",
    "    data = PowerModels.parse_file(data_file)\n",
    "    ac_result = solve_opf(data, ACPPowerModel, nlp_solver)\n",
    "    soc_result = solve_opf(data, SOCWRPowerModel, nlp_solver)\n",
    "    \n",
    "    # it is important to check the solver terminated nominally!\n",
    "    @assert  ac_result[\"termination_status\"] == OPTIMAL ||  ac_result[\"termination_status\"] == LOCALLY_SOLVED\n",
    "    @assert soc_result[\"termination_status\"] == OPTIMAL || soc_result[\"termination_status\"] == LOCALLY_SOLVED\n",
    "    \n",
    "    # by convention optimality gaps are given as a percentage\n",
    "    gap = 100*(ac_result[\"objective\"] - soc_result[\"objective\"])/ac_result[\"objective\"]\n",
    "    \n",
    "    return (gap=gap, ub_time=ac_result[\"solve_time\"], lb_time=soc_result[\"solve_time\"]) \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soc_gap(\"data/pglib_opf_case5_pjm.m\", nlp_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A 14% gap, that's HUGE!  It is not clear if the AC-OPF solution is optimal or not.  Clearly this case requires more study.\n",
    "\n",
    "In the mean time, try the `soc_gap` function on a few more of the network cases included with this tutorial.  Available datasets are `pglib_opf_case14_ieee.m`, `pglib_opf_case30_ieee.m`, `pglib_opf_case73_ieee_rts.m`, `pglib_opf_case118_ieee.m`, `pglib_opf_case200_tamu.m`, `pglib_opf_case500_tamu.m`, `pglib_opf_case1354_pegase.m`, `pglib_opf_case2383wp_k.m` (larger may take a few minutes to run).  The complete PGLib-OPF benchmark library with over 100 cases is available [here](https://github.com/power-grid-lib/pglib-opf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soc_gap(\"data/pglib_opf_case14_ieee.m\", nlp_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discussion Point\n",
    "\n",
    "Can you draw any conclusions from computing optimality gaps on a variety of cases?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Simple Contingency Study\n",
    "\n",
    "Component `status` parameters can be useful for exploring a variety of component failure scenarios.  Lets use it to do a simple generator N-1 study on the 5-bus example.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function single_generator_loss(data_file, nlp_solver)\n",
    "    data = PowerModels.parse_file(data_file)\n",
    "    \n",
    "    result = solve_opf(data, ACPPowerModel, nlp_solver)\n",
    "\n",
    "    # it is important to check the solver terminated nominally!\n",
    "    @assert  result[\"termination_status\"] == OPTIMAL || result[\"termination_status\"] == LOCALLY_SOLVED\n",
    "\n",
    "    println(\"base cost: $(result[\"objective\"])\")\n",
    "    \n",
    "    for (i,gen) in sort(data[\"gen\"])\n",
    "        # remove the generator\n",
    "        gen[\"gen_status\"] = 0\n",
    "        \n",
    "        result = solve_opf(data, ACPPowerModel, nlp_solver)\n",
    "        \n",
    "        println(\"gen $(i) out: $(result[\"objective\"]), $(result[\"termination_status\"])\")\n",
    "\n",
    "        # restore the generator\n",
    "        gen[\"gen_status\"] = 1\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_generator_loss(\"data/pglib_opf_case5_pjm.m\", nlp_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simple study seems to suggest that this network is not N-1 secure and cannot support the loss of generator 3 or 5.  Let's compute how much these two generators contribute to the total generation capacity in the data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pg_capacity = sum(gen[\"pmax\"] for (i,gen) in data[\"gen\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "100*data[\"gen\"][\"3\"][\"pmax\"]/pg_capacity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "100*data[\"gen\"][\"5\"][\"pmax\"]/pg_capacity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems reasonable that losing over 30% of the total generating capacity would be problematic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving Optimal Transmission Switching (OTS)\n",
    "\n",
    "Because InfrastructureModels is built on JuMP, is it fairly easy to encode challenging mixed integer nonlinear optimization problems.  Solving these problems however presents additional challenges.\n",
    "\n",
    "Lets consider the OTS problem.  An extension of OPF where branches can be removed from the network (i.e. a binary decion).  Lets try solving OPF and OTS on the 5-bus example.  Because the OTS problem has discrete variables we will need to use a MIP solver for that one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = PowerModels.parse_file(\"data/pglib_opf_case5_pjm.m\")\n",
    "\n",
    "opf_result = solve_opf(data, DCPPowerModel, lp_solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using HiGHS\n",
    "\n",
    "mip_solver = optimizer_with_attributes(HiGHS.Optimizer, \"output_flag\" => false)\n",
    "\n",
    "ots_result = solve_ots(data, DCPPowerModel, mip_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at those massive cost savings!!!  Lets look at the solution to see what happened."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PowerModels.print_summary(ots_result[\"solution\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that removing branch 5 from the network will remove congestion caused by Ohm's Law and yield the observed savings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### OTS Home Work\n",
    "\n",
    "Explore the AC-OTS problem on the 5-bus case.  Build the model with commands like `run_ots(data, ACPPowerModel, <solver>)` and `run_ots(data, SOCWRPowerModel, <solver>)`.  Use solvers like Juniper/Bonmin as heuristics to find upper bound solutions and solvers like Juniper/Pavito to find lower bounds.  Compare the gaps similar to how we did for the OPF probem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Topics\n",
    "\n",
    "PowerModels includes a variety of useful features that are outside the scope of this brief tutorial.  These features include,\n",
    "* Extracting Duals/LMPs from linear models using output options\n",
    "* Network models with HVDC lines or storage devices (see https://github.com/lanl-ansi/energy-storage-example)\n",
    "* Optimization Based Bound Tightening algorithms\n",
    "* Using mulitnetworks for time series analysis\n",
    "* Scalable SCOPF with PTDF cuts\n",
    "* Modeling multi-phase networks (e.g. see [PowerModelsDistribution](https://github.com/lanl-ansi/PowerModelsDistribution.jl))\n",
    "* A variety of experimental features in [PowerModelsAnnex](https://github.com/lanl-ansi/PowerModelsAnnex.jl)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.3",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}