Line Modeling Simulations
Originally Contributed by: Rodrigo Henriquez-Auba and José Daniel Lara
Introduction
This tutorial will introduce an example of considering dynamic lines in PowerSimulationsDynamics.
This tutorial presents a simulation of a three-bus system, with an infinite bus (represented as a voltage source behind an impedance) at bus 1, a one d- one q- machine on bus 2 and an inverter of 19 states, as a virtual synchronous machine at bus 3. The perturbation will be the trip of two of the three circuits (triplicating its resistance and impedance) of the line that connects bus 1 and bus 3. This case also consider a dynamic line model for connection between buses 2 and 3. We will compare it against a system without dynamic lines.
It is recommended to check the OMIB tutorial first, since that includes more details and explanations on all definitions and functions.
Step 1: Package Initialization
julia> using PowerSimulationsDynamicsjulia> using PowerSystemsjulia> using PowerNetworkMatricesjulia> using PowerSystemCaseBuilderjulia> using Sundialsjulia> using Plots
PowerSystemCaseBuilder.jl is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function build_system. For more details visit PowerSystemCaseBuilder Documentation
Step 2: Data creation
Load the system using PowerSystemCaseBuilder.jl:
julia> threebus_sys = build_system(PSIDSystems, "Three Bus Dynamic data Example System")System ┌───────────────────┬─────────────┐ │ Property │ Value │ ├───────────────────┼─────────────┤ │ Name │ │ │ Description │ │ │ System Units Base │ SYSTEM_BASE │ │ Base Power │ 100.0 │ │ Base Frequency │ 60.0 │ │ Num Components │ 19 │ └───────────────────┴─────────────┘ Static Components ┌─────────────────┬───────┬────────────────────────┬───────────────┐ │ Type │ Count │ Has Static Time Series │ Has Forecasts │ ├─────────────────┼───────┼────────────────────────┼───────────────┤ │ ACBus │ 3 │ false │ false │ │ Arc │ 3 │ false │ false │ │ Area │ 1 │ false │ false │ │ Line │ 3 │ false │ false │ │ LoadZone │ 1 │ false │ false │ │ Source │ 1 │ false │ false │ │ StandardLoad │ 3 │ false │ false │ │ ThermalStandard │ 2 │ false │ false │ └─────────────────┴───────┴────────────────────────┴───────────────┘ Dynamic Components ┌──────────────────┬───────┐ │ Type │ Count │ ├──────────────────┼───────┤ │ DynamicGenerator │ 1 │ │ DynamicInverter │ 1 │ └──────────────────┴───────┘
In addition, we will create a new copy of the system on which we will simulate the same case, but will consider dynamic lines:
julia> threebus_sys_dyn = deepcopy(threebus_sys);
Step 3: Create the fault and simulation on the Static Lines system
First, we construct the perturbation, by properly computing the new Ybus on the system:
julia> #Make a copy of the original system sys2 = deepcopy(threebus_sys) #Triplicates the impedance of the line named "BUS 1-BUS 3-i_1"System ┌───────────────────┬─────────────┐ │ Property │ Value │ ├───────────────────┼─────────────┤ │ Name │ │ │ Description │ │ │ System Units Base │ SYSTEM_BASE │ │ Base Power │ 100.0 │ │ Base Frequency │ 60.0 │ │ Num Components │ 19 │ └───────────────────┴─────────────┘ Static Components ┌─────────────────┬───────┬────────────────────────┬───────────────┐ │ Type │ Count │ Has Static Time Series │ Has Forecasts │ ├─────────────────┼───────┼────────────────────────┼───────────────┤ │ ACBus │ 3 │ false │ false │ │ Arc │ 3 │ false │ false │ │ Area │ 1 │ false │ false │ │ Line │ 3 │ false │ false │ │ LoadZone │ 1 │ false │ false │ │ Source │ 1 │ false │ false │ │ StandardLoad │ 3 │ false │ false │ │ ThermalStandard │ 2 │ false │ false │ └─────────────────┴───────┴────────────────────────┴───────────────┘ Dynamic Components ┌──────────────────┬───────┐ │ Type │ Count │ ├──────────────────┼───────┤ │ DynamicGenerator │ 1 │ │ DynamicInverter │ 1 │ └──────────────────┴───────┘julia> fault_branches = get_components(ACBranch, sys2)ACBranch Counts: Line: 3julia> for br in fault_branches if get_name(br) == "BUS 1-BUS 3-i_1" br.r = 3 * br.r br.x = 3 * br.x b_new = (from = br.b.from / 3, to = br.b.to / 3) br.b = b_new end end #Obtain the new Ybusjulia> Ybus_fault = Ybus(sys2).data #Define Fault: Change of YBus3×3 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 9 stored entries: 0.91954-10.9011im -0.689655+8.27586im -0.229885+2.75862im -0.689655+8.27586im 0.714334-8.78642im -0.0246792+1.11056im -0.229885+2.75862im -0.0246792+1.11056im 0.254564-3.33585imjulia> Ybus_change = NetworkSwitch( 1.0, #change at t = 1.0 Ybus_fault, #New YBus );
Now, we construct the simulation:
julia> #Time span of our simulation tspan = (0.0, 30.0) #Define Simulation(0.0, 30.0)julia> sim = Simulation( ResidualModel, #Type of model used threebus_sys, #system pwd(), #folder to output results tspan, #time span Ybus_change, #Type of perturbation )Simulation Summary ┌─────────────────────────┬────────────────┐ │ Property │ Value │ ├─────────────────────────┼────────────────┤ │ Status │ BUILT │ │ Simulation Type │ Residual Model │ │ Initialized? │ Yes │ │ Multimachine system? │ No │ │ Time Span │ (0.0, 30.0) │ │ Number of States │ 33 │ │ Number of Perturbations │ 1 │ └─────────────────────────┴────────────────┘
We can obtain the initial conditions as:
julia> #Print the initial states. It also give the symbols used to describe those states. show_states_initial_value(sim)Voltage Variables ==================== BUS 1 ==================== Vm 1.02 θ 0.0 ==================== BUS 2 ==================== Vm 1.0142 θ -0.0247 ==================== BUS 3 ==================== Vm 1.0059 θ 0.05 ==================== ==================== Differential States generator-102-1 ==================== eq_p 0.6478 ed_p 0.6672 δ 0.9386 ω 1.0 Vf 1.0781 Vr1 0.0333 Vr2 -0.1941 Vm 1.0142 ==================== Differential States generator-103-1 ==================== θ_oc 0.4573 ω_oc 1.0 q_oc -0.4453 ξd_ic 0.0013 ξq_ic 0.0004 γd_ic 0.0615 γq_ic -0.0138 ϕd_ic 0.8765 ϕq_ic -0.1978 vd_pll 0.8986 vq_pll 0.0 ε_pll 0.0 θ_pll 0.2354 ir_cnv 0.7462 ii_cnv 0.757 vr_filter 0.8738 vi_filter 0.2095 ir_filter 0.7617 ii_filter 0.6923 ====================
Step 4: Run the simulation of the Static Lines System
julia> #Run the simulation execute!( sim, #simulation structure IDA(), #Sundials DAE Solver dtmax = 0.02, #Maximum step size )SIMULATION_FINALIZED::BUILD_STATUS = 6
Step 5: Store the solution
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 │ 1607 │ │ Number of States │ 33 │ │ Total solve time │ 2.425772422 │ └────────────────────────────┴─────────────┘julia> series2 = get_voltage_magnitude_series(results, 102)([0.0, 0.001, 0.002, 0.004, 0.008, 0.016, 0.032, 0.052000000000000005, 0.07200000000000001, 0.09200000000000001 … 29.83954153464911, 29.85954153464911, 29.87954153464911, 29.89954153464911, 29.91954153464911, 29.93954153464911, 29.95954153464911, 29.97954153464911, 29.999541534649108, 30.0], [1.0142000000000002, 1.0142000000001268, 1.0142000000001326, 1.0142000000001927, 1.01420000000014, 1.014200000000125, 1.0142000000001052, 1.0142000000000875, 1.014200000000071, 1.0142000000000564 … 1.0140816770731866, 1.0140823112693007, 1.014082942782033, 1.0140835715610201, 1.014084197556437, 1.0140848207190878, 1.0140854410003357, 1.0140860583521254, 1.0140866727269917, 1.0140866867751974])julia> zoom = [ (series2[1][ix], series2[2][ix]) for (ix, s) in enumerate(series2[1]) if (s > 0.90 && s < 1.6) ];
Step 3.1: Create the fault and simulation on the Dynamic Lines system
An important aspect to consider is that DynamicLines must not be considered in the computation of the Ybus. First we construct the Dynamic Line, by finding the Line named "BUS 2-BUS 3-i_1", and then adding it to the system.
julia> # get component return the Branch on threebus_sys_dyn named "BUS 2-BUS 3-i_1" dyn_branch = DynamicBranch(get_component(Line, threebus_sys_dyn,"BUS 2-BUS 3-i_1")) # Adding a dynamic line will immediately remove the static line from the system.BUS 2-BUS 3-i_1 (DynamicBranch): branch: BUS 2-BUS 3-i_1 (Line) n_states: 2 states: [:Il_R, :Il_I] internal: InfrastructureSystems.InfrastructureSystemsInternaljulia> add_component!(threebus_sys_dyn, dyn_branch)
Similarly, we construct the Ybus fault by creating a copy of the original system, but removing the Line "BUS 2-BUS 3-i_1" to avoid considering it in the Ybus:
julia> #Make a copy of the original system sys3 = deepcopy(threebus_sys); #Remove Line "BUS 2-BUS 3-i_1"julia> remove_component!(Line, sys3, "BUS 2-BUS 3-i_1") #Triplicates the impedance of the line named "BUS 1-BUS 2-i_1"julia> fault_branches2 = get_components(Line, sys3)Line Counts: Line: 2julia> for br in fault_branches2 if get_name(br) == "BUS 1-BUS 3-i_1" br.r = 3 * br.r br.x = 3 * br.x b_new = (from = br.b.from / 3, to = br.b.to / 3) br.b = b_new end end #Obtain the new Ybusjulia> Ybus_fault_dyn = Ybus(sys3).data #Define Fault: Change of YBus3×3 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 7 stored entries: 0.91954-10.9011im -0.689655+8.27586im -0.229885+2.75862im -0.689655+8.27586im 0.689655-8.17586im ⋅ -0.229885+2.75862im ⋅ 0.229885-2.72529imjulia> Ybus_change_dyn = NetworkSwitch( 1.0, #change at t = 1.0 Ybus_fault_dyn, #New YBus )NetworkSwitch(1.0, sparse([1, 2, 3, 4, 5, 6, 1, 2, 4, 5 … 5, 6, 1, 2, 4, 5, 1, 3, 4, 6], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2 … 4, 4, 5, 5, 5, 5, 6, 6, 6, 6], [0.9195402298850576, -0.6896551724137931, -0.2298850574712644, 10.901149425287358, -8.275862068965518, -2.7586206896551726, -0.6896551724137931, 0.6896551724137931, -8.275862068965518, 8.175862068965518 … -0.6896551724137931, -0.2298850574712644, 8.275862068965518, -8.175862068965518, -0.6896551724137931, 0.6896551724137931, 2.7586206896551726, -2.7252873563218394, -0.2298850574712644, 0.2298850574712644], 6, 6))
Step 4.1: Run the simulation of the Dynamic Lines System
Now, we construct the simulation:
julia> # Define Simulation sim_dyn = Simulation( ResidualModel, #Type of model used threebus_sys_dyn, #system pwd(), #folder to output results (0.0, 30.0), #time span Ybus_change_dyn, #Type of perturbation )Simulation Summary ┌─────────────────────────┬────────────────┐ │ Property │ Value │ ├─────────────────────────┼────────────────┤ │ Status │ BUILT │ │ Simulation Type │ Residual Model │ │ Initialized? │ Yes │ │ Multimachine system? │ No │ │ Time Span │ (0.0, 30.0) │ │ Number of States │ 35 │ │ Number of Perturbations │ 1 │ └─────────────────────────┴────────────────┘
julia> # Run the simulation execute!( sim_dyn, #simulation structure IDA(), #Sundials DAE Solver dtmax = 0.02, #Maximum step size )SIMULATION_FINALIZED::BUILD_STATUS = 6
We can obtain the initial conditions as:
julia> #Print the initial states. It also give the symbols used to describe those states. show_states_initial_value(sim_dyn)Voltage Variables ==================== BUS 1 ==================== Vm 1.02 θ 0.0 ==================== BUS 2 ==================== Vm 1.0141 θ -0.0155 ==================== BUS 3 ==================== Vm 0.9686 θ 0.1354 ==================== ==================== Differential States generator-102-1 ==================== eq_p 0.681 ed_p 0.6477 δ 0.907 ω 1.0 Vf 1.1115 Vr1 0.0355 Vr2 -0.2001 Vm 1.0141 ==================== Differential States generator-103-1 ==================== θ_oc 0.5658 ω_oc 1.0 q_oc -0.3799 ξd_ic 0.0014 ξq_ic 0.0003 γd_ic 0.0595 γq_ic -0.014 ϕd_ic 0.8484 ϕq_ic -0.201 vd_pll 0.8719 vq_pll 0.0 ε_pll -0.0 θ_pll 0.3332 ir_cnv 0.715 ii_cnv 0.7768 vr_filter 0.824 vi_filter 0.2851 ir_filter 0.7361 ii_filter 0.7158 ==================== ==================== Line Current States ==================== Line BUS 2-BUS 3-i_1 Il_R -0.16141 Il_I -0.06381 ====================
Step 5.1: Store the solution
julia> results_dyn = read_results(sim_dyn)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 │ 1721 │ │ Number of States │ 35 │ │ Total solve time │ 2.298442545 │ └────────────────────────────┴─────────────┘julia> series2_dyn = get_voltage_magnitude_series(results_dyn, 102);julia> zoom_dyn = [ (series2_dyn[1][ix], series2_dyn[2][ix]) for (ix, s) in enumerate(series2_dyn[1]) if (s > 0.90 && s < 1.6) ];
Step 6.1: Compare the solutions:
We can observe the effect of Dynamic Lines
julia> plot(series2_dyn, label = "V_gen_dyn");julia> plot!(series2, label = "V_gen_st", xlabel = "Time [s]", ylabel = "Voltage [pu]");
that looks quite similar. The differences can be observed in the zoom plot:
julia> plot(zoom_dyn, label = "V_gen_dyn");julia> plot!(zoom, label = "V_gen_st", xlabel = "Time [s]", ylabel = "Voltage [pu]");
