Incorrect codegen for SDE system

[baggepinnen]

The following example generates a SDE system where the g function fails with

julia> sol = solve(prob, SRIW1())
ERROR: UndefVarError: `wind₊alt` not defined

looking at the generated code, there are still several symbolic variables in there

julia> prob.g.g_iip
RuntimeGeneratedFunction(#=in ModelingToolkit=#, #=using ModelingToolkit=#, :((ˍ₋out, ˍ₋arg1, ˍ₋arg2, t)->begin
    ˍ₋out[1] = (/)((*)((*)((*)((*)(0.32808400000000004, ˍ₋arg2[2]), (^)((+)(0.177, (*)(0.0027001313199999997, wind₊alt(t))), 0.7999999999999999)), NaNMath.sqrt((/)((*)((*)(2.0886476139744556, wind₊alt(t)), wind₊V(t)), (^)((+)(0.177, (*)(0.0027001313199999997, wind₊alt(t))), 1.2)))), wind₊V(t)), (*)(3.28084, wind₊alt(t)))

using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq, LinearAlgebra
using ControlSystemsBase, StochasticDiffEq
using ModelingToolkitStandardLibrary.Blocks
import ModelingToolkit: t_nounits as t, D_nounits as D
function DrydenWind2(; name)
    @parameters begin
        W20 = 4.572, [description = "Mean wind at 20 feet"]
        b = 10, [description = "Wing span"]
    end
    @variables begin
        u(t) = 0, [description = "Wind velocity in x direction"]
        v(t) = 0, [description = "Wind velocity in x direction"]
        w(t) = 0, [description = "Wind velocity in x direction"]
        Vx(t)[1:2] = zeros(2), [description = "State for wind speed in y direction"]
        Wx(t)[1:2] = zeros(2), [description = "State for wind speed in z direction"]
        p(t) = 0, [description = "Angular rate roll"]
        q(t) = 0, [description = "Angular rate pitch"]
        r(t) = 0, [description = "Angular rate yaw"]
        Q(t)[1:3] = zeros(3), [description = "States for angular rate pitch"]
        R(t)[1:3] = zeros(3), [description = "States for angular rate yaw"]
        σu(t) = 0, [description = "Turbulence intensity parameters"]
        σv(t) = 0, [description = "Turbulence intensity parameters"]
        σw(t) = 0, [description = "Turbulence intensity parameters"]
        Lu(t) = 0, [description = "Scaling parameters"]
        Lv(t) = 0, [description = "Scaling parameters"]
        Lw(t) = 0, [description = "Scaling parameters"]
        V(t), [description = "speed"]
        alt(t) = 0, [description = "Altitude"]
    end
    @brownian η
    Vx, Wx, Q, R = collect.((Vx, Wx, Q, R))
    systems = @named begin
        altitudeIn = RealInput()
        velocityIn = RealInput()
    end
    D = Differential(t)
    alt_ft = alt * 3.28084
    W20_ft = W20 * 3.28084
    σw = 0.1 * W20_ft
    σu = σw / (0.177 + 0.000823 * alt_ft)^0.4
    σv = σu
    Lu = alt_ft / (0.177 + 0.000823 * alt_ft)^1.2
    Lv = Lu
    Lw = alt_ft

    # Logitudinal
    Au = (σu * (sqrt(2 * Lu / pi * V)))

    mAp = (σw * (sqrt(0.8 / V)))
    nAp = (pi / (4 * b))^(1 / 6)
    dAp1 = Lw^(1 / 3)
    dAp2 = 4 * b / pi * V
    ssP = ControlSystemsBase.ss(mAp * ControlSystemsBase.tf([nAp], [dAp1 * dAp2, dAp1]))
    Ap = vec(ssP.A)
    Bp = vec(ssP.B)

    # Lateral
    mAv = (σv * (sqrt(Lv / pi * V)))
    nAv = sqrt(3) * Lv / V
    dAv = Lv / V
    ssV = ControlSystemsBase.ss(
        mAv * ControlSystemsBase.tf([nAv, 1], [dAv^2, 2 * dAv, 1]),
        balance = false)
    Av = ssV.A
    Bv = ssV.B
    Cv = ssV.C

    nAr = 1 / V
    dAr = 3 * b / pi * V
    tfr = ControlSystemsBase.tf([-nAr, 0], [dAr, 1])
    ssR = ControlSystemsBase.ss(
        tfr * mAv * ControlSystemsBase.tf([nAv, 1], [dAv^2, 2 * dAv, 1]),
        balance = false)
    Ar = ssR.A
    Br = ssR.B
    Cr = ssR.C

    # Vertical
    mAw = (σw * (sqrt(Lw / pi * V)))
    nAw = sqrt(3) * Lw / V
    dAw = Lw / V
    ssW = ControlSystemsBase.ss(
        mAw * ControlSystemsBase.tf([nAw, 1], [dAw^2, 2 * dAw, 1]),
        balance = false)
    Aw = ssW.A
    Bw = ssW.B
    Cw = ssW.C

    nAq = 1 / V
    dAq = 4 * b / pi * V
    tfq = ControlSystemsBase.tf([nAq, 0], [dAq, 1])
    ssQ = ControlSystemsBase.ss(
        tfq * mAw * ControlSystemsBase.tf([nAw, 1], [dAw^2, 2 * dAw, 1]),
        balance = false)
    Aq = ssQ.A
    Bq = ssQ.B
    Cq = ssQ.C

    eqs = [
        V ~ velocityIn.u
        alt ~ altitudeIn.u

        D(u) ~ (V / Lu) * (Au * η - u)
        D(p) ~ dot(Ap, p) + dot(Bp, η)
        D.(Vx) .~ vec(Av * Vx + Bv * η)
        D.(R) .~ vec((Ar * R) + Br * η)
        D.(Wx) .~ vec(Aw * Wx + Bw * η)
        D.(Q) .~ vec(Aq * Q + Bq * η)
        v ~ (Cv * Vx)[1]
        w ~ (Cw * Wx)[1]
        q ~ (Cq * Q)[1]
        r ~ (Cr * R)[1]
    ]
    System(eqs, t; systems, name)
end

@mtkmodel WindWithInput begin
	@components begin
		altitude = Constant(k=6)
        velocity = Constant(k=1)
		wind = DrydenWind2()
	end
	@equations begin
		connect(altitude.output, :ua, wind.altitudeIn)
        connect(velocity.output, :uv, wind.velocityIn)
	end
end;

tspan = (0.0,100.0)
@named model = WindWithInput()
cmodel = complete(model)
ssys = structural_simplify(model)
prob = SDEProblem(ssys, [], tspan)
sol = solve(prob, SRIW1())

cc @shobhitvoleti