structurally_simplify removes algebraic variables involved in events

[JKRT]

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()

[YingboMa]

Is there a way to trick structurally_simplify to not remove a set of variables specified by the user?

You can mark the variable as irreducible.

julia> using ModelingToolkit

julia> @variables t x(t);

julia> @named sys = ODESystem([x ~ 0]); equations(structural_simplify(sys))
Equation[]

julia> @variables t x(t) [irreducible = true];

julia> @named sys = ODESystem([x ~ 0]); equations(structural_simplify(sys))
1-element Vector{Equation}:
 0 ~ -x(t)

[YingboMa]

Here's a better way to generate the variables.

julia> t = Symbolics.variable(:t, T = Real);

julia> var_names = [:x, :y]
2-element Vector{Symbol}:
 :x
 :y

julia> vars = map(n -> Symbolics.variable(n, T = Symbolics.FnType{Tuple{Real}, Real})(t), var_names)
2-element Vector{Num}:
 x(t)
 y(t)