Linearization regression

[baggepinnen]

The linearization workflow below took

julia> @time batch_ss(model, :u, :y, ops);
  0.108045 seconds (731.06 k allocations: 60.930 MiB, 0.00% compilation time)

on MTKv8 while on master it takes

3.217503 seconds (2.82 M allocations: 1.931 GiB, 1.94% gc time, 0.00% compilation time)

Most of this time appears to be in this call.

It also appears as if the result is wrong in MTK v9

I have other examples where the regression is over 100x

using ModelingToolkit, ModelingToolkitStandardLibrary, ControlSystemsMTK, ControlSystemsBase
using ModelingToolkitStandardLibrary.Blocks

@parameters t;
D = Differential(t);
rc = 0.25 # Reference concentration 

@mtkmodel MixingTank begin
    @parameters begin
        c0  = 0.8,   [description = "Nominal concentration"]
        T0  = 308.5, [description = "Nominal temperature"]
        a1  = 0.2674
        a21 = 1.815
        a22 = 0.4682
        b   = 1.5476
        k0  = 1.05e14
        ϵ   = 34.2894
    end

    @variables begin
        gamma(t),     [description = "Reaction speed"]
        xc(t) = c0,   [description = "Concentration"]
        xT(t) = T0,   [description = "Temperature"]
        xT_c(t) = T0, [description = "Cooling temperature"]
    end

    @components begin
        T_c = RealInput()
        c = RealOutput()
        T = RealOutput()
    end

    begin
        τ0   = 60
        wk0  = k0/c0
        wϵ   = ϵ*T0
        wa11 = a1/τ0
        wa12 = c0/τ0
        wa13 = c0*a1/τ0
        wa21 = a21/τ0
        wa22 = a22*T0/τ0
        wa23 = T0*(a21 - b)/τ0
        wb   = b/τ0
    end
    @equations begin
        gamma ~ xc*wk0*exp( -wϵ/xT)
        D(xc) ~ -wa11*xc - wa12*gamma + wa13
        D(xT) ~ -wa21*xT + wa22*gamma + wa23 + wb*xT_c

        xc ~ c.u
        xT ~ T.u
        xT_c ~ T_c.u
    end
end

Ftf = tf(1, [(100), 1])^3
Fss = ss(Ftf)

"Compute initial state that yields y0 as output"
function init_filter(y0)
    (; A,B,C,D) = Fss
    Fx0 = -A\B*y0
    @assert C*Fx0 ≈ [y0] "C*Fx0*y0 ≈ y0 failed, got $(C*Fx0*y0) ≈ $(y0)]" 
    Fx0
end
RefFilter(; y0, name) = ODESystem(Fss; name, x0=init_filter(y0))

@mtkmodel InverseControlledTank begin
    begin
        c0 = 0.8    #  "Nominal concentration
        T0 = 308.5  #  "Nominal temperature
        x10 = 0.42  
        x20 = 0.01 
        u0 = -0.0224 

        c_start = c0*(1-x10)        # Initial concentration
        T_start = T0*(1+x20)        # Initial temperature
        c_high_start = c0*(1-0.72)  # Reference concentration
        T_c_start = T0*(1+u0)       # Initial cooling temperature
    end
    @components begin
        ref             = Constant(k=0.25) # Concentration reference
        ff_gain         = Gain(k=1) # To allow turning ff off
        controller      = PI(gainPI.k=10, T=500)
        tank            = MixingTank(xc=c_start, xT = T_start, c0=c0, T0=T0)
        inverse_tank    = MixingTank(xc=c_start, xT = T_start, c0=c0, T0=T0)
        feedback        = Feedback()
        add             = Add()
        filter          = RefFilter(y0=c_start) # Initialize filter states to the initial concentration
        noise_filter    = FirstOrder(k=1, T=1, x=T_start)
        limiter         = Limiter(y_max=370, y_min=250) # Saturate the control input 
    end
    @equations begin
        connect(ref.output, :r, filter.input)
        connect(filter.output, inverse_tank.c)
        connect(inverse_tank.T_c, ff_gain.input)
        connect(ff_gain.output, :uff, limiter.input)
        connect(limiter.output, add.input1)
        connect(controller.ctr_output, :u, add.input2)
        connect(add.output, :u_tot, tank.T_c)
        connect(inverse_tank.T, feedback.input1)
        connect(tank.T, :y, noise_filter.input)
        connect(noise_filter.output, feedback.input2)
        connect(feedback.output, :e, controller.err_input)
    end
end;
@named model = InverseControlledTank()
cm = complete(model)

op = Dict(
    D(cm.inverse_tank.xT) => 1,
    cm.tank.xc => 0.65
)

@time linearize(model, :u, :y; op)
using BenchmarkTools

ops = fill(op, 100)
@time batch_ss(model, :u, :y, ops)

[baggepinnen]

[baggepinnen]

In MTK v8, the call to init is very fast

julia> @btime init(prob, OrdinaryDiffEq.Rodas5P())
  61.636 μs (361 allocations: 74.64 KiB)

and produces a zero residual (last two equations)

julia> prob.f(res.u, prob.p, 0)
9-element Vector{Float64}:
 -0.0
  0.00042497264388805433
 -0.0013812304089120954
  0.11162675106913678
 -0.0
 -0.0033437500000000004
  0.0
 -6.622879952788985e-9
  3.034388078286968e-8

The new init produces huge residuals