One Machine against Infinite Bus (OMIB) Simulation

Originally Contributed by: Rodrigo Henriquez-Auba and José Daniel Lara

Introduction

This tutorial will introduce you to the functionality of PowerSimulationsDynamics for running power system dynamic simulations.

This tutorial presents a simulation of a two-bus system with an infinite bus (represented as a voltage source behind an impedance) at bus 1, and a classic machine on bus 2. The perturbation will be the trip of one of the two circuits (doubling its resistance and impedance) of the line that connects both buses.

Dependencies

julia> using PowerSimulationsDynamics
julia> PSID = PowerSimulationsDynamicsPowerSimulationsDynamics
julia> using PowerSystemCaseBuilder
julia> using PowerSystems
julia> const PSY = PowerSystemsPowerSystems
julia> using Sundials
julia> using Plots
julia> gr()Plots.GRBackend()

PowerSystems (abbreviated with PSY) is used to properly define the data structure and establish an equilibrium point initial condition with a power flow routine, while Sundials is used to solve the problem defined in PowerSimulationsDynamics.

Load the system

We load the system using PowerSystemCaseBuilder.jl:

julia> omib_sys = build_system(PSIDSystems, "OMIB System")System
┌───────────────────┬─────────────┐
│ Property          │ Value       │
├───────────────────┼─────────────┤
│ Name              │             │
│ Description       │             │
│ System Units Base │ SYSTEM_BASE │
│ Base Power        │ 100.0       │
│ Base Frequency    │ 60.0        │
│ Num Components    │ 10          │
└───────────────────┴─────────────┘

Static Components
┌─────────────────┬───────┬────────────────────────┬───────────────┐
│ Type            │ Count │ Has Static Time Series │ Has Forecasts │
├─────────────────┼───────┼────────────────────────┼───────────────┤
│ ACBus           │ 2     │ false                  │ false         │
│ Arc             │ 1     │ false                  │ false         │
│ Area            │ 1     │ false                  │ false         │
│ Line            │ 2     │ false                  │ false         │
│ LoadZone        │ 1     │ false                  │ false         │
│ Source          │ 1     │ false                  │ false         │
│ ThermalStandard │ 1     │ false                  │ false         │
└─────────────────┴───────┴────────────────────────┴───────────────┘

Dynamic Components
┌──────────────────┬───────┐
│ Type             │ Count │
├──────────────────┼───────┤
│ DynamicGenerator │ 1     │
└──────────────────┴───────┘

Build the simulation and initialize the problem

The next step is to create the simulation structure. This will create the indexing of our system that will be used to formulate the differential-algebraic system of equations. To do so, it is required to specify the perturbation that will occur in the system. PowerSimulationsDynamics supports multiple types of perturbations. See Perturbations

Here, we will use a Branch Trip perturbation, that is modeled by modifying the specifying which line we want to trip. In this case we disconnect one of the lines that connects BUS 1 and BUS 2, named "BUS 1-BUS 2-i_1".

With this, we are ready to create our simulation structure:

julia> time_span = (0.0, 30.0)(0.0, 30.0)
julia> perturbation_trip = BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1")BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1")
julia> sim = PSID.Simulation(
           ResidualModel, # Type of formulation
           omib_sys, # System
           mktempdir(), # Output directory
           time_span,
           perturbation_trip)Simulation Summary
┌─────────────────────────┬────────────────┐
│ Property                │ Value          │
├─────────────────────────┼────────────────┤
│ Status                  │ BUILT          │
│ Simulation Type         │ Residual Model │
│ Initialized?            │ Yes            │
│ Multimachine system?    │ No             │
│ Time Span               │ (0.0, 30.0)    │
│ Number of States        │ 6              │
│ Number of Perturbations │ 1              │
└─────────────────────────┴────────────────┘

This will automatically initialize the system by running a power flow and update V_ref, P_ref and hence eq_p (the internal voltage) to match the solution of the power flow. It will also initialize the states in the equilibrium, which can be printed with:

julia> show_states_initial_value(sim)Voltage Variables
====================
BUS 1
====================
Vm 1.05
θ 0.0
====================
BUS 2
====================
Vm 1.04
θ 0.0229
====================
====================
Differential States
generator-102-1
====================
δ 0.1685
ω 1.0
====================

Run the Simulation

Finally, to run the simulation we simply use:

julia> execute!(
           sim, #simulation structure
           IDA(), #Sundials DAE Solver
           dtmax = 0.02, #Arguments: Maximum timestep allowed
       );

In some cases, the dynamic time step used for the simulation may fail. In such case, the keyword argument dtmax can be used to limit the maximum time step allowed for the simulation.

Exploring the solution

First, we need to load the simulation results into memory:

julia> results = read_results(sim)Simulation Results Summary
┌────────────────────────────┬─────────────┐
│ Property                   │ Value       │
├────────────────────────────┼─────────────┤
│ System Base Power [MVA]    │ 100.0       │
│ System Base Frequency [Hz] │ 60.0        │
│ Time Span                  │ (0.0, 30.0) │
│ Total Time Steps           │ 1512        │
│ Number of States           │ 6           │
│ Total solve time           │ 1.525229078 │
└────────────────────────────┴─────────────┘

PowerSimulationsDynamics has two functions to obtain different states of the solution:

• get_state_series(results, ("generator-102-1", :δ)): can be used to obtain the solution as a tuple of time and the required state. In this case, we are obtaining the rotor angle :δ of the generator named "generator-102-1"`.

julia> angle = get_state_series(results, ("generator-102-1", :δ));
julia> plot(angle, xlabel = "time", ylabel = "rotor angle [rad]", label = "rotor angle")Plot{Plots.GRBackend() n=1}

• get_voltage_magnitude_series(results, 102): can be used to obtain the voltage magnitude as a tuple of time and voltage. In this case, we are obtaining the voltage magnitude at bus 102 (where the generator is located).

julia> volt = get_voltage_magnitude_series(results, 102);
julia> plot(volt, xlabel = "time", ylabel = "Voltage [pu]", label = "V_2")Plot{Plots.GRBackend() n=1}

Optional: Small Signal Analysis

PowerSimulationsDynamics uses automatic differentiation to compute the reduced Jacobian of the system for the differential states. This can be used to analyze the local stability of the linearized system. We need to re-initialize our simulation:

julia> sim2 = Simulation(ResidualModel, omib_sys, mktempdir(), time_span)Simulation Summary
┌─────────────────────────┬────────────────┐
│ Property                │ Value          │
├─────────────────────────┼────────────────┤
│ Status                  │ BUILT          │
│ Simulation Type         │ Residual Model │
│ Initialized?            │ Yes            │
│ Multimachine system?    │ No             │
│ Time Span               │ (0.0, 30.0)    │
│ Number of States        │ 6              │
│ Number of Perturbations │ 0              │
└─────────────────────────┴────────────────┘
julia> small_sig = small_signal_analysis(sim2)The system is small signal stable

The small_sig result can report the reduced jacobian for $\delta$ and $\omega$,

julia> small_sig.reduced_jacobian2×2 Matrix{Float64}:
  0.0       376.991
 -0.466763   -0.317662

and can also be used to report the eigenvalues of the reduced linearized system:

julia> small_sig.eigenvalues2-element Vector{ComplexF64}:
 -0.15883100381194407 - 13.264252847627304im
 -0.15883100381194407 + 13.264252847627304im