{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "One Machine against Infinite Bus (OMIB) simulation with [PowerSimulationsDynamics.jl](https://github.com/NREL-SIIP/PowerSimulationsDynamics.jl)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Originally Contributed by**: Rodrigo Henriquez and José Daniel Lara"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Introduction"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This tutorial will introduce you to the functionality of `PowerSimulationsDynamics`\n",
    "for running power system dynamic simulations."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This tutorial presents a simulation of a two-bus system with an infinite bus\n",
    "(represented as a voltage source behind an impedance) at bus 1, and a classic\n",
    "machine on bus 2. The perturbation will be the trip of one of the two circuits\n",
    "(doubling its resistance and impedance) of the line that connects both buses."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Dependencies"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SIIPExamples #hide\n",
    "using PowerSimulationsDynamics\n",
    "using PowerSystems\n",
    "using Sundials\n",
    "using Plots\n",
    "gr()\n",
    "PSD = PowerSimulationsDynamics"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "`PowerSystems` (abbreviated with `PSY`) is used to properly define the data structure and establish an equilibrium\n",
    "point initial condition with a power flow routine, while `Sundials` is\n",
    "used to solve the problem defined in `PowerSimulationsDynamics`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Load the system\n",
    "_The following command requires that you have executed the\n",
    "[dynamic systems data example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/2_PowerSystems_examples/09_loading_dynamic_systems_data.ipynb)\n",
    "previously to generate the json file._"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "file_dir = joinpath(\n",
    "    dirname(dirname(pathof(SIIPExamples))),\n",
    "    \"script\",\n",
    "    \"4_PowerSimulationsDynamics_examples\",\n",
    "    \"Data\",\n",
    ")\n",
    "omib_sys = System(joinpath(file_dir, \"omib_sys.json\"))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Build the simulation and initialize the problem"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The next step is to create the simulation structure. This will create the indexing\n",
    "of our system that will be used to formulate the differential-algebraic system of\n",
    "equations. To do so, it is required to specify the perturbation that will occur in\n",
    "the system. `PowerSimulationsDynamics` supports three types of perturbations:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- Network Switch: Change in the Y-bus values.\n",
    "- Branch Trip: Disconnects a line from the system.\n",
    "- Change in Reference Parameter"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we will use a Branch Trip perturbation, that is modeled by modifying the\n",
    "specifying which line we want to trip. In this case we disconnect one of the lines\n",
    "that connects BUS 1 and BUS 2, named \"BUS 1-BUS 2-i_1\"."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "With this, we are ready to create our simulation structure:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "time_span = (0.0, 30.0)\n",
    "perturbation_trip = BranchTrip(1.0, Line, \"BUS 1-BUS 2-i_1\")\n",
    "sim = PSD.Simulation(ResidualModel, omib_sys, pwd(), time_span, perturbation_trip)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This will automatically initialize the system by running a power flow\n",
    "and update `V_ref`, `P_ref` and hence `eq_p` (the internal voltage) to match the\n",
    "solution of the power flow. It will also initialize the states in the equilibrium,\n",
    "which can be printed with:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "show_states_initial_value(sim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To examine the calculated initial conditions, we can export them into a dictionary:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "x0_init = PSD.get_initial_conditions(sim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Run the Simulation"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Finally, to run the simulation we simply use:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "PSD.execute!(\n",
    "    sim, #simulation structure\n",
    "    IDA(), #Sundials DAE Solver\n",
    "    dtmax = 0.02,\n",
    "); #Arguments: Maximum timestep allowed"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In some cases, the dynamic time step used for the simulation may fail. In such case, the\n",
    "keyword argument `dtmax` can be used to limit the maximum time step allowed for the simulation."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Exploring the solution"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "First, we need to load the simulation results into memory"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "results = read_results(sim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "`PowerSimulationsDynamics` has two functions to obtain different\n",
    "states of the solution:\n",
    " - `get_state_series(results, (\"generator-102-1\", :δ))`: can be used to obtain the solution as\n",
    "a tuple of time and the required state. In this case, we are obtaining the rotor angle `:δ`\n",
    "of the generator named `\"generator-102-1\"`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "angle = get_state_series(results, (\"generator-102-1\", :δ));\n",
    "plot(angle, xlabel = \"time\", ylabel = \"rotor angle [rad]\", label = \"rotor angle\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `get_voltage_magnitude_series(results, 102)`: can be used to obtain the voltage magnitude as a\n",
    "tuple of time and voltage. In this case, we are obtaining the voltage magnitude at bus 102\n",
    "(where the generator is located)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "volt = get_voltage_magnitude_series(results, 102);\n",
    "plot(volt, xlabel = \"time\", ylabel = \"Voltage [pu]\", label = \"V_2\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Optional: Small Signal Analysis"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "`PowerSimulationsDynamics` uses automatic differentiation to compute the reduced Jacobian\n",
    "of the system for the differential states. This can be used to analyze the local stability\n",
    "of the linearized system. We need to re-initialize our simulation:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sim2 = PSD.Simulation(ResidualModel, omib_sys, pwd(), time_span)\n",
    "\n",
    "small_sig = small_signal_analysis(sim2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The `small_sig` result can report the reduced jacobian for $\\delta$ and $\\omega$,"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "small_sig.reduced_jacobian"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "and can also be used to report the eigenvalues of the reduced linearized system:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "small_sig.eigenvalues"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  },
  "kernelspec": {
   "name": "julia-1.7",
   "display_name": "Julia 1.7.2",
   "language": "julia"
  }
 },
 "nbformat": 4
}