GG causes internal StackOverflow with ModelingToolkit

[jd-lara]

MWE causes error below.

using ModelingToolkit
const MTK = ModelingToolkit

params = MTK.@parameters begin
        t
        # AC side quantities
        ωb      # Base Frequency
        # Grid impadance
        lg      # Grid reactance
        rg      # Grid resistance
        #Reference set-point input
        pʳ      # Active Power Setpoint
        qʳ      # Reactive Power Setpoint
        vʳ      # Voltage Setpoint
        ωʳ      # Frequency Setpoint
        # Load at rated voltage
        vl      # Load rated voltage
        pl      # Load rated power
        # Filter parameters
        lf      # Filter reactance
        cf      # Filter capacitance
        rf      # Filter Resistance
        ωz      # Filtering frequency
        # Transformer Parameters
        rt      # Transformer resistance
        lt      # Transformer reactance
        # OuterControl Loops
        Dp      # Active Power Droop
        Dq      # Reactive Power Droop
        # SRF Current Control
        kip     # Current control propotional gain
        kii     # Current control integral gain
        kffi    # Current control differential gain
        # SRF Voltage Control
        kvp     # Voltage control propotional gain
        kvi     # Voltage control integral gain
        kffv    # Voltage control differential gain
        # Virtual Impedance
        rv
        lv
        # DC Source Parameters
        leq     # Equivalent inductance (i.e. battery inductance and DC/DC converter inductance)
        req     # Equivalent resistance
        vb      # Battery Voltage
        cdc     # Dc-side capacitance
        # DC/DC converter controller parameters
        vdcʳ    # DC Voltage reference
        kpvb    # DC/DC Voltage control integral gain
        kivb    # DC/DC Voltage control propotional gain
        kpib    # DC/DC Current control propotional gain
        kiib    # DC/DC Current control Integral gain
        a1      # First coefficient of Pade approximation
        a2      # Second co-efficient of Pade approxmiation
        Ts      # DC/DC controller time delay
    end

    variables = MTK.@variables begin
        eg_d #d-axis capacitor filter voltage
        eg_q #q-axis capacitor filter voltage
        is_d #d-axis current flowing into filter
        is_q #q-axis current flowing into filter
        ig_d #d-axis current flowing into grid
        ig_q #q-axis current flowing into grid
        pf   #Filtered active power measurement
        qf   #Filtered reactive power measurement
        ξ_d  #d-axis integrator term for outer AC/DC PI controller
        ξ_q  #q-axis integrator term for outer AC/DC PI controller
        γ_d  #d-axis integrator term for inner AC/DC PI controller
        γ_q  #d-axis integrator term for inner AC/DC PI controller
        vdc  #DC Voltage
        ibat #Battery Current
        η    #Integrator term for outer DC/DC PI controller
        κ    #Integrator term for inner DC/DC PI controller
        # TODO: Verify in the nomenclature equation is the appropiate for each term of the Pade approximation
        M    # First term for Pade approx
        L    # Second term for Pade approx
    end

    # Expressions
    pm = eg_d * ig_d + eg_q * ig_q  # AC Active Power Calculation
    qm = -eg_d * ig_q + eg_q * ig_d # AC Reactive Power Calculation
    ω_a = ωʳ + Dp * (pʳ - pf)  # Active Power Drop
    # TODO: Original model had pf here. Verify
    v_hat = vʳ + Dq * (qʳ - qf) # Reactive Power Drop
    v_iref_d = v_hat - rv * ig_d + ω_a * lv * ig_q # d-axis virtual impedance equation
    v_iref_q = -rv * ig_q - ω_a * lv * ig_d # q-axis virtual impedance equation
    i_hat_d = kvp * (v_iref_d - eg_d) + kvi * ξ_d - ω_a * cf * eg_q # Inner voltage controller d PI
    i_hat_q = kvp * (v_iref_q - eg_q) + kvi * ξ_q + ω_a * cf * eg_d # Inner voltage controller q PI
    v_md = kip * (i_hat_d - is_d) + kii * γ_d - ω_a * lf * is_q # Inner current controller d PI
    v_mq = kip * (i_hat_q - is_q) + kii * γ_q + ω_a * lf * is_d # Inner current controller q PI
    p_inv = v_md * is_d + v_mq * is_q # Active power drawn from inverter
    q_inv = -v_md * is_q + v_mq * is_d # Reactive power drawn from inverter
    v_gd = (vl^2 / pl) * ig_d
    v_gq = (vl^2 / pl) * ig_q
    i_ref = kpvb * (vdcʳ - vdc) + kivb * η
    i_in = (vb * ibat - ibat^2 * req) / vdc
    d_dc = (-a2 / Ts) * M + kpib * (i_ref - i_in) + kiib * κ

    model_rhs = [
        ### Grid forming equations
        #𝜕eg_d/𝜕t
        (ωb / cf) * (is_d - ig_d) + ω_a * ωb * eg_q
        #𝜕eg_q/𝜕t
        (ωb / cf) * (is_q - ig_q) - ω_a * ωb * eg_d
        #𝜕is_d/𝜕t
        (ωb / lf) * (v_md - eg_d) - (rf * ωb / lf) * is_d + ωb * ω_a * is_q
        #𝜕is_q/𝜕t
        (ωb / lf) * (v_mq - eg_q) - (rf * ωb / lf) * is_q - ωb * ω_a * is_d
        #𝜕ig_d/𝜕t
        (ωb / lt) * (eg_d - v_gd) - (rt * ωb / lt) * ig_d + ωb * ω_a * ig_q
        #𝜕ig_q/𝜕t
        (ωb / lt) * (eg_q - v_gq) - (rt * ωb / lt) * ig_q - ωb * ω_a * ig_d
        #𝜕pf/𝜕t
        ωz * (pm - pf)
        #𝜕qf/𝜕t
        ωz * (qm - qf)
        #𝜕ξ_d/𝜕t
        v_iref_d - eg_d
        #𝜕ξ_q/𝜕t
        v_iref_q - eg_q
        #𝜕γ_d/𝜕t
        i_hat_d - is_d
        #𝜕γ_q/𝜕t
        i_hat_q - is_q
        ### DC-side equations
        #∂vdc/∂t
        ωb * ((i_in - p_inv / (2 * vdc)) / (cdc))
        #∂ib/∂t
        (ωb / leq) * (vb - req * ibat - (1 - d_dc) * vdc)
        #∂η/∂t
        vdcʳ - vdc # Integrator for DC/DC outer PI controller
        #∂κ/dt
        i_ref - i_in # Integrator for DC/DC inner PI controller
        # ∂M/dt
        (-a1 / Ts) * M + (-a2 / Ts^2) * L + kpib * (i_ref - i_in) + kiib * ibat # First term in Pade approximation
        # ∂M/dt
        M # Second term in Pade approx.
    ]

_eqs = zeros(length(model_rhs)) .~ model_rhs
_nl_system =  MTK.NonlinearSystem(_eqs, [variables...], [params...][2:end])
jac = MTK.generate_jacobian(_nl_system, expression = Val{false})[2]

u0 = [
    0.9808903498431045
 -0.09574546685726393
  0.48560976954113544
 -0.022652021638989283
  0.4785215581979559
 -0.09526912125100888
  0.47849876511951955
  0.04763229168221858
  0.0006501651605950486
 -0.00012944174083017514
  0.06869560693368727
 -0.006700239365187476
  1.01
  0.33065335095454984
  0.059307890471113116
  0.027324656609971815
  1.706771821458912e-29
  2.7823532298747112e-8
]
param_eval = (out, params) -> jac(out, u0, params)
n= length(u0)
J = zeros(n, n)
parameter_values =
[
    376.99111843077515
   0.2
   0.01
   0.5
   0.0
   1.0
   1.0
   1.0
   0.5
   0.08
   0.074
   0.003
  37.69911184307752
   0.01
   0.2
   0.02
   0.001
   1.27
  14.3
   0.0
   0.59
 736.0
   1.0
   0.0
   0.2
   0.7918317967460096
   0.0
   0.7246376811594203
   0.08974273574244604
   1.01
   0.6
   4.0
   0.3863
  10.34
   6.0
  12.0
   0.0003125
]
param_eval(J, parameter_values)

This is the error

Internal error: encountered unexpected error in runtime:
StackOverflowError()
intersect at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:0
intersect_invariant at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:2742
intersect at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:3002
intersect_tuple at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:2642 [inlined]
intersect at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:2958
intersect_unionall_ at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:2491
intersect_unionall at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:2539
intersect at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:0
intersect_all at /Users/julia/buildbot/worker/package_macos64/build/src/subtype.c:3047
jl_type_intersection_env_s at 

I had to cut the bulk of the StackOverFlow stack trace since Github as a limit in the number of Characters.

cc @ChrisRackauckas

[thautwarm]

Hello, would you mind using files to share the full trace with me? I have to detect the GG related issues and bugs and report them to the language devs.

[jd-lara]

@thautwarm https://gist.github.com/jd-lara/e5bff07a22cf224e374d9f3f2d73a702

[thautwarm]

Thanks @jd-lara , and I'm sorry that you did suffer a lot from the extremely long trace...
Currently according to the trace I'd think it's because the type is too big for analyze and makes the compiler crash. GG encodes julia expressions into julia types, and the latter can be really large when your generated code is large.
Anyway I'll try to make experiments and report an issue to julialang dev.
Thanks again for your report, and appologize for your bad experience with GG.

[thautwarm]

I made some optimizations to avoid drastic AST -> Type conversions.
Maybe it works for this issue, but I don't have time now for testing.
@YingboMa Maybe you could give me a hand if you're not busy?

[thautwarm]

still unsolved yet..

[YingboMa]

This particular case can be solved by just not generating all the 0s. We should be able to solve this downstream.

julia> symjac = MTK.calculate_jacobian(_nl_system);

julia> J = sparse(Float64.(isequal.(symjac, 0)));

julia> jac = MTK.generate_jacobian(_nl_system, expression = Val{false}, sparse = true)[2]

julia> param_eval = (out, params) -> jac(out, u0, params);

julia> param_eval(J, parameter_values)

julia> J
18×18 SparseMatrixCSC{Float64,Int64} with 228 stored entries:
  [1 ,  1]  =  -377.153
  [6 ,  1]  =  -8243.38
  [10,  1]  =  443.061
  [14,  1]  =  1884.96
  [15,  1]  =  18.0398
  [16,  1]  =  3.59156
  [17,  1]  =  -1.0
...

[jd-lara]

@YingboMa yeah that makes sense, it will force us to allocate the Jacobian matrix again to calculate the eigenvalues though.... but I guess it isn't that bad of a solution.