{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Line Modeling simulation with [PowerSimulationsDynamics.jl](https://github.com/NREL-SIIP/PowerSimulationsDynamics.jl)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Originally Contributed by**: José Daniel Lara"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Introduction"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This tutorial will introduce the modeling of an inverter with Virtual Inertia in a multi-machine\n",
    "model of the system. We will load the data directly from PSS/e dynamic files"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The tutorial uses a modified 14-bus system on which all the synchronous machines have been\n",
    "substituted by generators with ESAC1A AVR's and no Turbine Governors."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In the first portion of the tutorial we will simulate the system with the original data and\n",
    "cause a line trip between Buses 2 and 4. In the second part of the simulation, we will switch\n",
    "generator 6 with a battery using an inverter and perform the same fault."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Load the packages"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SIIPExamples # Only needed for the tutorial, comment if you want to run\n",
    "import DisplayAs # Only needed for the tutorial\n",
    "using PowerSimulationsDynamics\n",
    "using PowerSystems\n",
    "using Logging\n",
    "using Sundials\n",
    "using Plots\n",
    "gr()\n",
    "PSD = PowerSimulationsDynamics"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Create the system"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "file_dir = joinpath(\n",
    "    dirname(dirname(pathof(SIIPExamples))),\n",
    "    \"script\",\n",
    "    \"4_PowerSimulationsDynamics_examples\",\n",
    "    \"Data\",\n",
    ")\n",
    "\n",
    "sys = System(joinpath(file_dir, \"14bus.raw\"), joinpath(file_dir, \"dyn_data.dyr\"))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define Simulation Problem with a 20 second simulation period and the branch trip at t = 1.0"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sim = PSD.Simulation(\n",
    "    ResidualModel, #Type of model used\n",
    "    sys,         #system\n",
    "    file_dir,       #path for the simulation output\n",
    "    (0.0, 20.0), #time span\n",
    "    BranchTrip(1.0, Line, \"BUS 02-BUS 04-i_4\");\n",
    "    console_level = Logging.Info,\n",
    ")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now that the system is initialized, we can verify the system states for potential issues."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "show_states_initial_value(sim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We execute the simulation with an additional tolerance for the solver set at 1e-8."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "PSD.execute!(sim, IDA(); abstol = 1e-8)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Using `PowerSimulationsDynamics` tools for exploring the results, we can plot all the voltage\n",
    "results for the buses"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "result = read_results(sim)\n",
    "p = plot()\n",
    "for b in get_components(Bus, sys)\n",
    "    voltage_series = get_voltage_magnitude_series(result, get_number(b))\n",
    "    plot!(\n",
    "        p,\n",
    "        voltage_series;\n",
    "        xlabel = \"Time\",\n",
    "        ylabel = \"Voltage Magnitude [pu]\",\n",
    "        label = \"Bus - $(get_name(b))\",\n",
    "    )\n",
    "end\n",
    "img = DisplayAs.PNG(p) # This line is only needed because of literate use display(p) when running locally"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also explore the frequency of the different generators"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "p2 = plot()\n",
    "for g in get_components(ThermalStandard, sys)\n",
    "    state_series = get_state_series(result, (get_name(g), :ω))\n",
    "    plot!(\n",
    "        p2,\n",
    "        state_series;\n",
    "        xlabel = \"Time\",\n",
    "        ylabel = \"Speed [pu]\",\n",
    "        label = \"$(get_name(g)) - ω\",\n",
    "    )\n",
    "end\n",
    "img = DisplayAs.PNG(p2) # This line is only needed because of literate use display(p2) when running locally"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "It is also possible to explore the small signal stability of this system we created. However,\n",
    "Since a simulation has already taken place, we need to reset the model."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "res = small_signal_analysis(sim; reset_simulation = true)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The eigenvalues can be explored visually"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "scatter(res.eigenvalues; legend = false)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Modifying the system and adding storage"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Reload the system for this example"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sys = System(joinpath(file_dir, \"14bus.raw\"), joinpath(file_dir, \"dyn_data.dyr\"))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We want to remove the generator 6 and the dynamic component attached to it."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "thermal_gen = get_component(ThermalStandard, sys, \"generator-6-1\")\n",
    "remove_component!(sys, get_dynamic_injector(thermal_gen))\n",
    "remove_component!(sys, thermal_gen)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can now define our storage device and add it to the system"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "storage = GenericBattery(\n",
    "    name = \"Battery\",\n",
    "    bus = get_component(Bus, sys, \"BUS 06\"),\n",
    "    available = true,\n",
    "    prime_mover = PrimeMovers.BA,\n",
    "    active_power = 0.6,\n",
    "    reactive_power = 0.16,\n",
    "    rating = 1.1,\n",
    "    base_power = 25.0,\n",
    "    initial_energy = 50.0,\n",
    "    state_of_charge_limits = (min = 5.0, max = 100.0),\n",
    "    input_active_power_limits = (min = 0.0, max = 1.0),\n",
    "    output_active_power_limits = (min = 0.0, max = 1.0),\n",
    "    reactive_power_limits = (min = -1.0, max = 1.0),\n",
    "    efficiency = (in = 0.80, out = 0.90),\n",
    ")\n",
    "\n",
    "add_component!(sys, storage)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "A good sanity check it running a power flow on the system to make sure all the components\n",
    "are properly scaled and that the system is properly balanced. We can use `PowerSystems` to\n",
    "perform this check. We can get the results back and perform a sanity check"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "res = solve_powerflow(sys)\n",
    "res[\"bus_results\"]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "After verifying that the system works, we can define our inverter dynamics and add it to the\n",
    "battery that has already been stored in the system."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "inverter = DynamicInverter(\n",
    "    name = get_name(storage),\n",
    "    ω_ref = 1.0, # ω_ref,\n",
    "    converter = AverageConverter(rated_voltage = 138.0, rated_current = 100.0),\n",
    "    outer_control = OuterControl(\n",
    "        VirtualInertia(Ta = 2.0, kd = 400.0, kω = 20.0),\n",
    "        ReactivePowerDroop(kq = 0.2, ωf = 1000.0),\n",
    "    ),\n",
    "    inner_control = CurrentControl(\n",
    "        kpv = 0.59,     #Voltage controller proportional gain\n",
    "        kiv = 736.0,    #Voltage controller integral gain\n",
    "        kffv = 0.0,     #Binary variable enabling the voltage feed-forward in output of current controllers\n",
    "        rv = 0.0,       #Virtual resistance in pu\n",
    "        lv = 0.2,       #Virtual inductance in pu\n",
    "        kpc = 1.27,     #Current controller proportional gain\n",
    "        kic = 14.3,     #Current controller integral gain\n",
    "        kffi = 0.0,     #Binary variable enabling the current feed-forward in output of current controllers\n",
    "        ωad = 50.0,     #Active damping low pass filter cut-off frequency\n",
    "        kad = 0.2,\n",
    "    ),\n",
    "    dc_source = FixedDCSource(voltage = 600.0),\n",
    "    freq_estimator = KauraPLL(\n",
    "        ω_lp = 500.0, #Cut-off frequency for LowPass filter of PLL filter.\n",
    "        kp_pll = 0.084,  #PLL proportional gain\n",
    "        ki_pll = 4.69,   #PLL integral gain\n",
    "    ),\n",
    "    filter = LCLFilter(lf = 0.08, rf = 0.003, cf = 0.074, lg = 0.2, rg = 0.01),\n",
    ")\n",
    "add_component!(sys, inverter, storage)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "These are the current system components:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sys"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define Simulation problem using the same parameters:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sim = PSD.Simulation(\n",
    "    ResidualModel, #Type of model used\n",
    "    sys,         #system\n",
    "    file_dir,       #path for the simulation output\n",
    "    (0.0, 20.0), #time span\n",
    "    BranchTrip(1.0, Line, \"BUS 02-BUS 04-i_4\");\n",
    "    console_level = Logging.Info,\n",
    ")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can verify the small signal stability of the system before running the simulation:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "res = small_signal_analysis(sim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "scatter(res.eigenvalues)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We execute the simulation"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "PSD.execute!(sim, IDA(); abstol = 1e-8)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Using `PowerSimulationsDynamics` tools for exploring the results, we can plot all the voltage\n",
    "results for the buses"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "result = read_results(sim)\n",
    "p = plot()\n",
    "for b in get_components(Bus, sys)\n",
    "    voltage_series = get_voltage_magnitude_series(result, get_number(b))\n",
    "    plot!(\n",
    "        p,\n",
    "        voltage_series;\n",
    "        xlabel = \"Time\",\n",
    "        ylabel = \"Voltage Magnitude [pu]\",\n",
    "        label = \"Bus - $(get_name(b))\",\n",
    "    )\n",
    "end\n",
    "img = DisplayAs.PNG(p) # This line is only needed because of literate use display(p) when running locally"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also explore the frequency of the different static generators and storage"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "p2 = plot()\n",
    "for g in get_components(ThermalStandard, sys)\n",
    "    state_series = get_state_series(result, (get_name(g), :ω))\n",
    "    plot!(\n",
    "        p2,\n",
    "        state_series;\n",
    "        xlabel = \"Time\",\n",
    "        ylabel = \"Speed [pu]\",\n",
    "        label = \"$(get_name(g)) - ω\",\n",
    "    )\n",
    "end\n",
    "state_series = get_state_series(result, (\"Battery\", :ω_oc))\n",
    "plot!(p2, state_series; xlabel = \"Time\", ylabel = \"Speed [pu]\", label = \"Battery - ω\")\n",
    "img = DisplayAs.PNG(p2) # This line is only needed because of literate use display(p2) when running locally"
   ],
   "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
}