Operating point passed to `linearization_function` closure has no effect

[baggepinnen]

I get the same result from linearization no matter what operating point I'm sending in. The first equality below incorrectly returns true, while the second correctly returns false. I also do not get S1 == S1f like I would have expected

using ModelingToolkit, ControlSystemsBase, Test, ControlSystemsMTK
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

begin
	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

	# Create an MTK-compatible constructor 
	RefFilter(; y0, name) = ODESystem(Fss; name, x0=init_filter(y0))
end
@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 = Gain(k=1)
		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)

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

op2 = Dict(
    D(cm.inverse_tank.xT) => 1,
    cm.tank.xc => 0.45
)

output = :y
lin_fun, ssys = Blocks.get_sensitivity_function(model, output)
matrices1 = linearize(ssys, lin_fun, op=op1)
matrices2 = linearize(ssys, lin_fun, op=op2)
S1f = ss(matrices1...)
S2f = ss(matrices2...)
@test S1f != S2f


matrices1, ssys = Blocks.get_sensitivity(model, output; op=op1)
matrices2, ssys = Blocks.get_sensitivity(model, output; op=op2)
S1 = ss(matrices1...)
S2 = ss(matrices2...)
@test S1 != S2

@test S1 == S1f
@test S2 == S2f

[ChrisRackauckas]

Can you double check this with the latest Symbolics and SymbolicUtils?

[AayushSabharwal]

This doesn't seem to be related to hashconsing

[AayushSabharwal]

The issue is multi-fold:

• get_u0 here doesn't recognize D(cm.inverse_tank.xT) for some reason. This is easily fixable by using the new process_SciMLProblem
• Whatever u0 you get out of it doesn't matter because linearization solves an initialization problem if the algebraic conditions aren't satisfied, that solve fails but instead of erroring it just uses the initial conditions from the failed solve
• The initialization fails because it's built using the defaults of the system which contain inconsistent values. Some of these values are meant to be overridden by op, but since that op isn't provided to linearization_function it isn't respected. E.g. model.tank.xc has a default of 0.4 which becomes an equation tank.xc ~ 0.4 in initialization, and the tank.xc => 0.45 in op isn't respected. Note that this is fixed by passing initialize = false to get_sensitivity_function (which forwards it to linearization_function), but that also avoids checking that algebraic conditions are satisfied.

There are multiple things we should do to solve the initialization issue. Firstly, this should use the existing OverrideInit infrastructure in SciMLBase instead of hand-rolling it. Next, it should error on failing initialization and use CheckInit if initialize = false. Lastly, we should enable passing var => nothing to op to remove defaults in the initialization system.