IDA solver error (possibly) for a DAE problem with complex number

New to Julia
1

Hi there. I have a problem with this code here. I’m trying to simulate a three-phase short circuit fault in a power system according to the instructions in Example 2 here:

https://github.com/NREL-Sienna/PowerSimulationsDynamics.jl/blob/main/docs/src/perturbations.md

The power system is also taken from the SIIPExamples repository. The main thing to do for simulating this scenario is using the NetworkSwitch function, which accepts complex numbers as input, and I don’t want to define new method for this function because the numbers are meant to be complex. Anyway, in this scenario, we introduce two perturbations and feed it into a simulation function. My code executes up until the last line. When I attempt to execute the simulation using IDA solver, I get this error:

Error: Execution failed
│   exception =
│    BoundsError: attempt to access 28×28 SparseMatrixCSC{Float64, Int64} with 257 stored entries at index [785]
│    Stacktrace:
│      [1] throw_boundserror(A::SparseMatrixCSC{Float64, Int64}, I::Tuple{Int64})
│        @ Base ./abstractarray.jl:737
│      [2] checkbounds
│        @ ./abstractarray.jl:702 [inlined]
│      [3] _setindex!
│        @ ./abstractarray.jl:1430 [inlined]
│      [4] setindex!
│        @ ./abstractarray.jl:1396 [inlined]
│      [5] (::PowerSimulationsDynamics.var"#27#28"{NetworkSwitch})

which is strange because the apply_fault and clear_fault are both 28*28 matrices, which mean that there should only be 784 cells. I double checked all the matrices and made sure that they all have proper dimensions and the code is still not running.

What was confusing to me was that IDA solver works with non-complex values. But I have to have complex values. I have tried several other solvers but could not find anything that works. I’m a newbie to Julia so I suspect that the problem is with something different and I have no clue. Can you help me?

This is the code:

using SIIPExamples
using PowerSimulationsDynamics
import PowerSimulationsDynamics as PSD
using PowerSystems
using Sundials
using SparseArrays
using OrdinaryDiffEq


# Load the power system model
file_dir = joinpath(dirname(dirname(pathof(SIIPExamples))), "script", "4_PowerSimulationsDynamics_examples", "Data")
sys = System(joinpath(file_dir, "14bus.raw"), joinpath(file_dir, "dyn_data.dyr"))

# Retrieve the Y-bus matrix and bus mapping
(ybus_matrix, bus_mapping) = PowerSimulationsDynamics._get_ybus(sys)

# Transform ybus_matrix into Complex type
ybus_matrix = SparseMatrixCSC{ComplexF64, Int64}(ybus_matrix)

# Store the original Y-Bus Matrix
original_ybus_matrix = copy(ybus_matrix)

# The original Y-Bus Matrix is already in Complex type
typeof(original_ybus_matrix)

# Modify Y-Bus Matrix for Fault Condition
bus_index = bus_mapping[6]
large_value = 1.0e4 + 0.0im
ybus_matrix[bus_index, bus_index] += large_value

# Create Sparse Matrices for Fault Condition and Clearing Fault
fault_admittance = sparse(ybus_matrix)
clear_fault_matrix = sparse(original_ybus_matrix)

# Define NetworkSwitch Events
time_of_fault = 1.0  # Time when the fault occurs
time_of_clear = 1.1  # Time when the fault is cleared

# fault_admittance = SparseMatrixCSC{Float64, Int64}(fault_admittance)
apply_fault = NetworkSwitch(time_of_fault, fault_admittance)
clear_fault = NetworkSwitch(time_of_clear, clear_fault_matrix)

three_phase_fault = [apply_fault, clear_fault]


# Define and Run Simulation
sim = PSD.Simulation(
    ResidualModel,  # Type of model used
    sys,            # System
    file_dir,       # Path for the simulation output
    (0.0, 10.0),    # Time span
    three_phase_fault,  # Events
)


# Execute the simulation with IDA
PSD.execute!(sim, IDA(); abstol = 1e-8)
2

Just split the input to real and complex parts

3

Inputs to which function? Can you elaborate?

4

The initial condition, I. Just write it with reals and combine to complex within f

5

Well they should eventually be combined for NetworkSwitch function, as this function only supports complex. But I have no problems with the data before this function.

6

Yes cool so then if you just pack it it’s solved in 3 lines of code, so it’s at least a solution for now

7

Sorry but I don’t understand what you’re saying. Can you explain what you mean?

8 
function f(u)
  _f(ComplexF64.(u[1,:], u[2,:]))
end

Makes f take in real numbers and reformat it for a function _f that wants complex numbers. It’s just a reinterpretation.

Of course you can make that even faster with reinterpret(Vector{ComplexF64}, u), interpreting it differently though

9

Thank you, but I still don’t see the point. I have a simulation scenario with 2 perturbations with complex values. What is the point of separating the two parts? Do you mean I should solve the simulation twice (once for real values and the other time for imaginary)?

Can you tell me what should I do with this f? I would be very happy if you elaborate on your suggestion as I’m confused.

10

That makes it real valued

11

What makes what real valued? Can you please use longer sentences and be more specific?

12

Complex numbers are just two real numbers. Write a function that takes in 2N real numbers, converts them into Complex numbers, does the computation in Complex, and then converts the result back to real.

2 Likes