Fixes a bug when calling Jacobian for a system with delays (#382)

* bug fix * add test
NREL-Sienna · Jul 2, 2024 · 3e1fe3e · 3e1fe3e
1 parent a001d10
commit 3e1fe3e
Show file tree

Hide file tree

Showing 3 changed files with 48 additions and 1 deletion.
diff --git a/src/base/jacobian.jl b/src/base/jacobian.jl
@@ -184,10 +184,11 @@ function JacobianFunctionWrapper(
         end
     jac = zeros(n, n)
     if sparse_retrieve_loop > 0
+        h(p, t; idxs = nothing) = typeof(idxs) <: Number ? x0[idxs] : x0
         for _ in 1:sparse_retrieve_loop
             temp = zeros(n, n)
             rng_state = copy(Random.default_rng())
-            Jf(temp, x0 + Random.randn(n))
+            Jf(temp, x0 + Random.randn(n), h, 0.0)
             copy!(Random.default_rng(), rng_state)
             jac .+= abs.(temp)
         end

diff --git a/test/test_base.jl b/test/test_base.jl
@@ -691,6 +691,24 @@ end
     end
 end
 
+@testset "Jacobian with Delays" begin
+    for model_type in [MassMatrixModel]
+        model_type = MassMatrixModel
+        omib_sys_copy = deepcopy(omib_sys)
+        add_degov_to_omib!(omib_sys_copy)
+        @test_throws IS.ConflictingInputsError get_jacobian(model_type, omib_sys_copy, -1)
+
+        jac = get_jacobian(model_type, omib_sys_copy)
+        @test isa(jac, Function)
+        @test isa(jac.Jv, PSID.SparseArrays.SparseMatrixCSC{Float64, Int64})
+
+        jac_dense = get_jacobian(model_type, omib_sys_copy, 0)
+        @test isa(jac_dense, Function)
+        @test !isa(jac_dense.Jv, PSID.SparseArrays.SparseMatrixCSC{Float64, Int64})
+        @test isa(jac_dense.Jv, Matrix{Float64})
+    end
+end
+
 @testset "Test 01 OMIB ResidualModel" begin
     path = mktempdir()
     try

diff --git a/test/utils/data_utils.jl b/test/utils/data_utils.jl
@@ -50,3 +50,31 @@ function get_ybus_fault_threebus_sys(sys)
     Ybus_fault = PNM.Ybus(fault_branch, sorted_buses)[:, :]
     return Ybus_fault
 end
+
+function add_degov_to_omib!(omib_sys)
+    gen = get_component(ThermalStandard, omib_sys, "generator-102-1")
+    dyn_gen = get_component(DynamicGenerator, omib_sys, "generator-102-1")
+    new_gov = PSY.DEGOV(;
+        T1 = 0.0,
+        T2 = 0.0,
+        T3 = 0.0,
+        K = 18.0,
+        T4 = 12.0,
+        T5 = 5.0,
+        T6 = 0.2,
+        Td = 0.5,
+        P_ref = 0.0,
+    )
+    dyn_gen_new = DynamicGenerator(;
+        name = get_name(dyn_gen),
+        ω_ref = get_ω_ref(dyn_gen),
+        machine = get_machine(dyn_gen),
+        shaft = get_shaft(dyn_gen),
+        avr = get_avr(dyn_gen),
+        prime_mover = new_gov,
+        pss = get_pss(dyn_gen),
+        base_power = get_base_power(dyn_gen),
+    )
+    remove_component!(omib_sys, dyn_gen)
+    add_component!(omib_sys, dyn_gen_new, gen)
+end