Open
Description
Hello!
I have a power grid model I am experimenting with.
Example at the bottom of this post (Note that the model is generated automatically:) )
Currently, when structural_simplify is used the model fails with the following message:
julia> PowerGridsReal__System99Simulate()
ERROR: UndefVarError: L3_port_omegaRef not defined
Stacktrace:
[1] macro expansion
@ C:\Users\johti17\.julia\packages\SymbolicUtils\qulQp\src\code.jl:394 [inlined]
[2] macro expansion
@ C:\Users\johti17\.julia\packages\Symbolics\FGTCH\src\build_function.jl:519 [inlined]
[3] macro expansion
@ C:\Users\johti17\.julia\packages\SymbolicUtils\qulQp\src\code.jl:351 [inlined]
[4] macro expansion
@ C:\Users\johti17\.julia\packages\RuntimeGeneratedFunctions\KrkGo\src\RuntimeGeneratedFunctions.jl:129 [inlined]
[5] macro expansion
@ .\none:0 [inlined]
[6] generated_callfunc
@ .\none:0 [inlined]
[7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x7dad26c4, 0x36373ff1, 0x91d469e3, 0x44dcb6f8, 0xe2027209)})(::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, ::Vector{Float64}, ::Vector{Float64}, ::Float64)
@ RuntimeGeneratedFunctions C:\Users\johti17\.julia\packages\RuntimeGeneratedFunctions\KrkGo\src\RuntimeGeneratedFunctions.jl:117
However, it works fine if dae_index_lowering is used instead:
#reducedSystem = structural_simplify(
# firstOrderSystem;
# simplify = false,
# allow_parameter = true,
#)
reducedSystem = dae_index_lowering(
firstOrderSystem
)The reason for the issue seems to be that structurally_simplify removes the
variables involved in the events (L3_port_omegaRef and L2_port_omegaRef):
events = [
[-L3_port_omegaRef - 0.0 ~ 0] => [ifCond1 ~ true],
[-L2_port_omegaRef - 0.0 ~ 0] => [ifCond2 ~ true],
[-L2_port_omegaRef - 0.0 ~ 0] => [ifCond3 ~ true],
[-L3_port_omegaRef - 0.0 ~ 0] => [ifCond4 ~ true],
]
Is there a way to trick structurally_simplify to not remove a set of variables specified by the user?
I suppose one issue is that the simplification is done before the events are added to the system
Cheers,
John
Complete model:
using ModelingToolkit
using DifferentialEquations
import ModelingToolkit.IfElse
function PowerGridsReal__System99CallbackSet(aux)
return CallbackSet()
end
function PowerGridsReal__System99Model(
tspan = (0.0, 1.0))
ModelingToolkit.@variables t
D = ModelingToolkit.Differential(t)
parameters = ModelingToolkit.@parameters(begin
G1_V
G1_Ta
G1_droop
G1_p
T1_B
T2_B
end)
begin
variableConstructors = Function[]
begin
function generateStateVariables1()
(
:t,
:(ifCond1(t)),
:(ifCond2(t)),
:(ifCond3(t)),
:(ifCond4(t)),
:(T1_B_act(t)),
:(T1_closed(t)),
:(T1_open(t)),
:(T1_close(t)),
:(T2_B_act(t)),
:(T2_closed(t)),
:(T2_open(t)),
:(T2_close(t)),
:(G1_theta(t)),
:(G1_omega(t)),
)
end
push!(variableConstructors, generateStateVariables1)
end
begin
function generateAlgebraicVariables1()
(
:t,
:(G1_port_v_re(t)),
:(G1_port_v_im(t)),
:(G1_port_i_re(t)),
:(G1_port_i_im(t)),
:(G1_port_omegaRef(t)),
:(G1_Ps(t)),
:(G1_Pc(t)),
:(G1_Pe(t)),
:(L1_port_v_re(t)),
:(L1_port_v_im(t)),
:(L1_port_i_re(t)),
:(L1_port_i_im(t)),
:(L1_port_omegaRef(t)),
:(L1_P(t)),
:(L1_Q(t)),
:(L2_port_v_re(t)),
:(L2_port_v_im(t)),
:(L2_port_i_re(t)),
:(L2_port_i_im(t)),
:(L2_port_omegaRef(t)),
:(L2_P(t)),
:(L2_Q(t)),
:(L3_port_v_re(t)),
:(L3_port_v_im(t)),
:(L3_port_i_re(t)),
:(L3_port_i_im(t)),
:(L3_port_omegaRef(t)),
:(L3_P(t)),
:(L3_Q(t)),
:(T1_Pab(t)),
:(T1_port_a_v_re(t)),
:(T1_port_a_v_im(t)),
:(T1_port_a_i_re(t)),
:(T1_port_a_i_im(t)),
:(T1_port_a_omegaRef(t)),
:(T1_port_b_v_re(t)),
:(T1_port_b_v_im(t)),
:(T1_port_b_i_re(t)),
:(T1_port_b_i_im(t)),
:(T1_port_b_omegaRef(t)),
:(T2_Pab(t)),
:(T2_port_a_v_re(t)),
:(T2_port_a_v_im(t)),
:(T2_port_a_i_re(t)),
:(T2_port_a_i_im(t)),
:(T2_port_a_omegaRef(t)),
:(T2_port_b_v_re(t)),
:(T2_port_b_v_im(t)),
:(T2_port_b_i_re(t)),
:(T2_port_b_i_im(t)),
)
end
push!(variableConstructors, generateAlgebraicVariables1)
end
begin
function generateAlgebraicVariables2()
(
:t,
:(T2_port_b_omegaRef(t)),
:(ifEq_tmp33(t)),
:(ifEq_tmp30(t)),
:(ifEq_tmp31(t)),
:(ifEq_tmp32(t)),
)
end
push!(variableConstructors, generateAlgebraicVariables2)
end
end
componentVars = []
for constructor in variableConstructors
res = eval(
ModelingToolkit.Symbolics._parse_vars(
"CustomCall",
Real,
constructor(),
),
)
push!(componentVars, res[2:end])
end
vars = collect(Iterators.flatten(componentVars))
pars = Dict(
G1_V => float(1.0),
G1_Ta => float(10.0),
G1_droop => float(0.05),
G1_p => float(0),
T1_B => float(-5.0),
T2_B => float(-5.0),
)
startEquationComponents = []
begin
startEquationConstructors = Function[]
begin
function generateStartEquations0()
[
ifCond1 => false,
ifCond2 => false,
ifCond3 => false,
ifCond4 => false,
T1_closed => true,
T1_B_act => -5.0,
T2_closed => true,
T2_B_act => -5.0,
T1_closed => 0.0,
T1_open => 0.0,
T1_close => 0.0,
T2_closed => 0.0,
T2_open => 0.0,
T2_close => 0.0,
G1_port_v_re => 1.0,
G1_port_v_re => 0.0,
G1_port_v_im => 0.0,
G1_port_i_re => 0.0,
G1_port_i_im => 0.0,
G1_port_omegaRef => 0.0,
G1_Ps => 0.0,
G1_Pc => 0.0,
G1_Pe => 0.0,
L1_port_v_re => 1.0,
L1_port_v_re => 0.0,
L1_port_v_im => 0.0,
L1_port_i_re => 0.0,
L1_port_i_im => 0.0,
L1_port_omegaRef => 0.0,
L1_P => 0.0,
L1_Q => 0.0,
L2_port_v_re => 1.0,
L2_port_v_re => 0.0,
L2_port_v_im => 0.0,
L2_port_i_re => 0.0,
L2_port_i_im => 0.0,
L2_port_omegaRef => 0.0,
L2_P => 0.0,
L2_Q => 0.0,
L3_port_v_re => 1.0,
L3_port_v_re => 0.0,
L3_port_v_im => 0.0,
L3_port_i_re => 0.0,
L3_port_i_im => 0.0,
L3_port_omegaRef => 0.0,
L3_P => 0.0,
L3_Q => 0.0,
T1_Pab => 0.0,
T1_port_a_v_re => 1.0,
T1_port_a_v_re => 0.0,
]
end
push!(startEquationConstructors, generateStartEquations0)
end
begin
function generateStartEquations1()
[
T1_port_a_v_im => 0.0,
T1_port_a_i_re => 0.0,
T1_port_a_i_im => 0.0,
T1_port_a_omegaRef => 0.0,
T1_port_b_v_re => 1.0,
T1_port_b_v_re => 0.0,
T1_port_b_v_im => 0.0,
T1_port_b_i_re => 0.0,
T1_port_b_i_im => 0.0,
T1_port_b_omegaRef => 0.0,
T2_Pab => 0.0,
T2_port_a_v_re => 1.0,
T2_port_a_v_re => 0.0,
T2_port_a_v_im => 0.0,
T2_port_a_i_re => 0.0,
T2_port_a_i_im => 0.0,
T2_port_a_omegaRef => 0.0,
T2_port_b_v_re => 1.0,
T2_port_b_v_re => 0.0,
T2_port_b_v_im => 0.0,
T2_port_b_i_re => 0.0,
T2_port_b_i_im => 0.0,
T2_port_b_omegaRef => 0.0,
ifEq_tmp33 => 0.0,
ifEq_tmp30 => 0.0,
ifEq_tmp31 => 0.0,
ifEq_tmp32 => 0.0,
G1_theta => 0.0,
G1_omega => 1.0,
T1_closed => true,
T1_B_act => -5.0,
T2_closed => true,
T2_B_act => -5.0,
]
end
push!(startEquationConstructors, generateStartEquations1)
end
end
for constructor in startEquationConstructors
push!(startEquationComponents, constructor())
end
initialValues = collect(Iterators.flatten(startEquationComponents))
equationComponents = []
begin
equationConstructors = Function[]
begin
G1_V = float(1.0)
G1_Ta = float(10.0)
G1_droop = float(0.05)
G1_p = float(0)
T1_B = float(-5.0)
T2_B = float(-5.0)
function generateEquations0()
[
0 ~ L1_port_omegaRef - T1_port_a_omegaRef,
0 ~ L1_port_omegaRef - G1_port_omegaRef,
0 ~ G1_port_v_im - L1_port_v_im,
0 ~ G1_port_v_im - T1_port_a_v_im,
0 ~ L1_port_v_re - G1_port_v_re,
0 ~ L1_port_v_re - T1_port_a_v_re,
0 ~ T2_port_a_omegaRef - T1_port_b_omegaRef,
0 ~ T2_port_a_omegaRef - L2_port_omegaRef,
0 ~ T2_port_a_v_im - T1_port_b_v_im,
0 ~ T2_port_a_v_im - L2_port_v_im,
0 ~ T2_port_a_v_re - L2_port_v_re,
0 ~ T2_port_a_v_re - T1_port_b_v_re,
0 ~ T2_port_b_omegaRef - L3_port_omegaRef,
0 ~ L3_port_v_im - T2_port_b_v_im,
0 ~ L3_port_v_re - T2_port_b_v_re,
0 ~ G1_port_i_re + L1_port_i_re + T1_port_a_i_re,
0 ~ G1_port_i_im + L1_port_i_im + T1_port_a_i_im,
0 ~ L2_port_i_re + T1_port_b_i_re + T2_port_a_i_re,
0 ~ L2_port_i_im + T1_port_b_i_im + T2_port_a_i_im,
0 ~ L3_port_i_re + T2_port_b_i_re,
0 ~ L3_port_i_im + T2_port_b_i_im,
D(G1_theta) ~
314.1592653589793G1_omega -
314.1592653589793G1_port_omegaRef,
D(G1_omega) ~ -(((G1_Pe - G1_Pc) - G1_Ps) / (G1_Ta * G1_omega)),
0 ~ G1_port_v_re - G1_V * cos(G1_theta),
0 ~ G1_port_v_im - G1_V * sin(G1_theta),
0 ~
G1_Pe +
G1_port_i_im * G1_port_v_im +
G1_port_i_re * G1_port_v_re,
0 ~ G1_Pc + (G1_omega - 1.0) / G1_droop,
0 ~ G1_port_omegaRef - G1_omega,
0 ~
(
L1_port_i_im * L1_port_v_im +
L1_port_i_re * L1_port_v_re
) - L1_P,
0 ~
(L1_port_i_re * L1_port_v_im - L1_Q) -
L1_port_i_im * L1_port_v_re,
0 ~ ifEq_tmp30,
0 ~ ifEq_tmp31,
0 ~ ifEq_tmp32,
0 ~ ifEq_tmp33,
0 ~ T1_port_a_omegaRef - T1_port_b_omegaRef,
0 ~ T1_port_a_i_re + T1_port_b_i_re,
0 ~ T1_port_a_i_im + T1_port_b_i_im,
0 ~
T1_port_a_i_re +
T1_B_act * (T1_port_a_v_im - T1_port_b_v_im),
0 ~
T1_port_a_i_im -
T1_B_act * (T1_port_a_v_re - T1_port_b_v_re),
0 ~
(T1_Pab - T1_port_a_i_im * T1_port_a_v_im) -
T1_port_a_i_re * T1_port_a_v_re,
0 ~ T2_port_a_omegaRef - T2_port_b_omegaRef,
0 ~ T2_port_a_i_re + T2_port_b_i_re,
0 ~ T2_port_a_i_im + T2_port_b_i_im,
0 ~
T2_port_a_i_re +
T2_B_act * (T2_port_a_v_im - T2_port_b_v_im),
0 ~
T2_port_a_i_im -
T2_B_act * (T2_port_a_v_re - T2_port_b_v_re),
0 ~
(T2_Pab - T2_port_a_i_im * T2_port_a_v_im) -
T2_port_a_i_re * T2_port_a_v_re,
-0.0 ~ G1_Ps - 1.0,
-0.0 ~ L1_P - 0.8,
0 ~ L1_Q,
-0.0 ~ L2_P - 0.1,
]
end
push!(equationConstructors, generateEquations0)
end
begin
G1_V = float(1.0)
G1_Ta = float(10.0)
G1_droop = float(0.05)
G1_p = float(0)
T1_B = float(-5.0)
T2_B = float(-5.0)
function generateEquations1()
[
0 ~ L2_Q,
-0.0 ~ L3_P - 0.1,
0 ~ L3_Q,
D(T1_B_act) ~ 0,
D(T1_closed) ~ 0,
D(T1_open) ~ 0,
D(T1_close) ~ 0,
D(T2_B_act) ~ 0,
D(T2_closed) ~ 0,
D(T2_open) ~ 0,
D(T2_close) ~ 0,
D(ifCond1) ~ 0,
D(ifCond2) ~ 0,
D(ifCond3) ~ 0,
D(ifCond4) ~ 0,
0 ~
ifelse(
ifCond1 == true,
(L3_port_i_re * L3_port_v_im - L3_Q) -
L3_port_i_im * L3_port_v_re,
L3_port_i_im,
) - ifEq_tmp33,
0 ~
ifelse(
ifCond2 == true,
(
L2_port_i_im * L2_port_v_im +
L2_port_i_re * L2_port_v_re
) - L2_P,
L2_port_i_re,
) - ifEq_tmp30,
0 ~
ifelse(
ifCond3 == true,
(L2_port_i_re * L2_port_v_im - L2_Q) -
L2_port_i_im * L2_port_v_re,
L2_port_i_im,
) - ifEq_tmp31,
0 ~
ifelse(
ifCond4 == true,
(
L3_port_i_im * L3_port_v_im +
L3_port_i_re * L3_port_v_re
) - L3_P,
L3_port_i_re,
) - ifEq_tmp32,
]
end
push!(equationConstructors, generateEquations1)
end
end
for constructor in equationConstructors
push!(equationComponents, constructor())
end
eqs = collect(Iterators.flatten(equationComponents))
events = [
[-L3_port_omegaRef - 0.0 ~ 0] => [ifCond1 ~ true],
[-L2_port_omegaRef - 0.0 ~ 0] => [ifCond2 ~ true],
[-L2_port_omegaRef - 0.0 ~ 0] => [ifCond3 ~ true],
[-L3_port_omegaRef - 0.0 ~ 0] => [ifCond4 ~ true],
]
nonLinearSystem = ODESystem(
eqs,
t,
vars,
parameters;
name = :(
$(Symbol("PowerGridsReal.System99"))
),
continuous_events = events,
)
firstOrderSystem = nonLinearSystem
reducedSystem = structural_simplify(
firstOrderSystem;
simplify = false,
allow_parameter = true,
)
# reducedSystem = dae_index_lowering(
# firstOrderSystem
# )
local event_p = [1.0, 10.0, 0.05, G1_p, -5.0, -5.0]
local discreteVars = collect(
values(
ModelingToolkit.OrderedDict(
T1_closed => true,
T1_B_act => -5.0,
T2_closed => true,
T2_B_act => -5.0,
T1_closed => 0.0,
T1_open => 0.0,
T1_close => 0.0,
T2_closed => 0.0,
T2_open => 0.0,
T2_close => 0.0,
),
),
)
event_p = vcat(event_p, discreteVars)
local aux = Vector{Any}(undef, 3)
aux[1] = event_p
aux[2] = Float64[]
aux[3] = reducedSystem
callbacks =
PowerGridsReal__System99CallbackSet(aux)
problem = ModelingToolkit.ODEProblem(
reducedSystem,
initialValues,
tspan,
pars,
callback = callbacks,
)
return (problem, callbacks, initialValues, reducedSystem, tspan, pars, vars)
end
end
function PowerGridsReal__System99Simulate(
tspan = (0.0, 1.0);
solver = Rosenbrock23(),
)
(
PowerGridsReal__System99Model_problem,
callbacks,
ivs,
PowerGridsReal__System99Model_ReducedSystem,
tspan,
pars,
vars,
) = PowerGridsReal__System99Model(tspan)
solve(
PowerGridsReal__System99Model_problem,
solver,
)
end
PowerGridsReal__System99Simulate()