ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

Linearization regression

Open baggepinnen opened this issue 1 month ago • 2 comments

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 avatar Jun 03 '24 16:06 baggepinnen