ERROR: CheckInit specified but initialization not satisfied

[bradcarman]

In the below example, I'm getting an error that appears to be about initialization. However, I can see that initialization works just fine. However, when I attempt to solve the problem I get an error, even when specifying the initial states.

MWE

using ModelingToolkit
using OrdinaryDiffEq
using Plots

using ModelingToolkit: D_nounits as D, t_nounits as t

L = 0.001

@mtkmodel System begin
    @parameters begin
        g = 9.807
        rho = 1000
        area = π*1^2

        # pump 
        c_1 = 1*g*rho
        c_2 = -10
        power(t) = 0
    end
    @variables begin
        p_1(t)
        p_2(t)
        h(t) = 10
        flow_in(t), [guess=0]
        flow_out(t)
    end
    @equations begin
        # pump
        # p_2 - p_1 ~ c_1*power + c_2*flow_in
        (p_2 - p_1)*flow_in ~ power

        # tank
        p_2 ~ rho*g*h
        D(h)*area ~ flow_in - flow_out

        # inputs
        p_1 ~ 0
        flow_out ~ 100 * L
    end
    @continuous_events begin
        [h ~ 5] => [power ~ 10]
        [h ~ 10] => [power ~ 0]
    end
end

@mtkbuild sys = System()

initsys = ModelingToolkit.generate_initializesystem(sys)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0])
initsol = solve(initprob)

sts = unknowns(sys)
u0 = sts .=> initsol[sts]
prob = ODEProblem(sys, u0, (0, 3600))
sol = solve(prob) #ERROR: CheckInit specified but initialization not satisfied. normresid = 0.004765618901311116 > abstol = 1.0e-6

[bradcarman]

OK, I see a hint, if I move the pump equation to a simple expression, then the problem solves fine. See below. Why should I need to do that?

@mtkmodel System begin
    @parameters begin
        g = 9.807
        rho = 1000
        area = π*1^2

        # pump 
        c_1 = 1*g*rho
        c_2 = -10
        power(t) = 0
    end
    @variables begin
        p_1(t)
        p_2(t)
        h(t) = 10
        # flow_in(t), [guess=0]
        flow_out(t)
    end
    begin
        # pump equation
        flow_in = power/(p_2 - p_1)
    end
    @equations begin
        # pump
        #(p_2 - p_1)*flow_in ~ power

        # tank
        p_2 ~ rho*g*h
        D(h)*area ~ flow_in - flow_out

        # inputs
        p_1 ~ 0
        flow_out ~ 20 * L
    end
    @continuous_events begin
        [h ~ 5] => [power ~ 2500]
        [h ~ 10] => [power ~ 0]
    end
end

[aml5600]

@bradcarman I have hit this problem frequently now too. Mine were related to the continuous events that I had, the check was happening after an event had occurred (although my events were terminating, I would still hit it).

It was difficult to address because abstol is used in both the solving tolerance, and init checking, so could not directly change that...needed to drop reltol pretty low to avoid this issue.

[bradcarman]

For reference, this is solvable by Modelica without issue...

model MTK
	parameter Real g = 9.807;
	parameter Real rho = 1000;
	parameter Real area = 3.14;
	parameter Real eta = 0.65;
	Real power(start=0);
	Real p_1;
	Real p_2;
	Real h(start = 10);
	Real flow_out;
	Real flow_in;
equation
	flow_in*(p_2 - p_1) = power*eta;
	p_2 = g*rho*h;
	area*der(h) = -flow_out + flow_in;
	p_1 = 0;
	flow_out = 0.02;

        der(power) = 0;

	when h <= 5 then
	  reinit(power, 5000);
	end when;
	
	when h >= 10 then
	  reinit(power, 0);
	end when;
end MTK;

[ChrisRackauckas]

Okay, it's not an issue, but the error message is not informative enough. So let me explain a little bit. @aml5600 it's the opposite of your original problem.

@mtkmodel System begin
    @parameters begin
        g = 9.807
        rho = 1000
        area = π*1^2

        # pump 
        c_1 = 1*g*rho
        c_2 = -10
        power(t) = 0
    end
    @variables begin
        p_1(t)
        p_2(t)
        h(t) = 10
        flow_in(t), [guess=0]
        flow_out(t)
    end
    @equations begin
        # pump
        # p_2 - p_1 ~ c_1*power + c_2*flow_in
        (p_2 - p_1)*flow_in ~ power

        # tank
        p_2 ~ rho*g*h
        D(h)*area ~ flow_in - flow_out

        # inputs
        p_1 ~ 0
        flow_out ~ 100 * L
    end
    @continuous_events begin
        [h ~ 5] => [power ~ 10, flow_in ~ power / (p_2 - p_1)]
        [h ~ 10] => [power ~ 0]
    end
end

That's what you're effectively doing. Modelica is changing your flow_in behind the scenes without declaring it will: it's allowed to change any algebraic variable in reinit. We have seen many cases where this can be rather dangerous though. One case is where a user changes an algebraic variable in a callback, the trivial behavior is then it just undoes your change and 👍 that is absolutely correct under algebraic variable initialization.

So what we have done is the following, and note we're half way through this process. We have first changed all of the callbacks to default to CheckInit(). What this means is, by default we assume that you give us a consistent system after a callback. You did not, flow_in ~ power / (p_2 - p_1) fails, so done. We should put a better error message on top of the standard one to contextualize it, but that's always the safest default. "I wanted to change this one thing", "sorry that's impossible, here's why".

We have next made it possible to change the callback_initializealg, that's what #3197 is all about. So if you do want it to reinitialize by changing all of the algebraic variables, you can pass callback_initializealg = BrownBasicInit() to get the behavior of Modelica, and for full initialization it's ShampineInit(). We just improved the messaging around that about an hour ago SciML/SciMLBase.jl#886 but we need to contextualize that a bit more in the context when it's a callback.

Now the real thing is to also fix the Modelica problem of changing algebraic variables. This is also a part of this rollout, let me explain. When you define a symbolic callback, what you're really doing is defining a larger implicit discrete system. Let's take the example of manually moving a pendulum. You might say D(x) = integrator.u[1] + 1 in code, but what you actually mean is the following:

x ~ x_old + 1
x^2 + y^2 ~ 1

Right now our symbolic callback would do this via x ~ x + 1, but it's very underdetermined what that actually means. What we really want to say here is, update x using the previous value to be one more, and also solve for y to get a consistent system. When I write it out like this, it should be obvious why Modelica has a bad behavior here in the case where y is the chosen variable: then it would say "explicitly update x ~ x + 1, now solve for x^2 + y^2 ~ 1, and if you do that then xgoes right back to where it is (this is whyreinit` in Modelica has the restriction that it only allows for changes to differential variables).

So how do we handle this sanely? Well what we can do is we can make a new primitive, the ImplicitDiscreteSystem. #2077 this is just a very useful primitive anyways, things like a handrolled implicit euler would be very straightforward with such an object so it's good for embedded codegen things. If we have this in place, then what we can do is change the symbolic callback interface to be an ImplicitDiscreteSystem, where you give us x ~ prev(x) + 1 (where we need to define the shifting operator in this context to make it simpler). What it would then do is always append the algebraic equations and make an explicit discrete system:

x ~ x_old + 1
x^2 + y^2 ~ 1

where the prev are basically just parameters. It would append the algebraic equations and by default it can append that the unknowns to solve for includes the changed values and it always appends the algebraic equations. So in the case where x is the differential variable, this would be an implicit system for x and y. The tearing of the nonlinear system would then make the first thing explicit, and you'd get effectively the same behavior as Modelica in this case. In the case where y is the differential variable, you'd have to declare that it's solvable in this system and solving this would still be a correct way to do it.

Now notice that if you do that nonlinear solve in the callback, then the callback is always guaranteed to be a consistent state, and therefore CheckInit() will pass. As such, you can see that is then consistent with in this direction of "assume the callback gives you something consistent, check that it is, and give an option for a different basic initializations" because now any symbolic callback will want that same CheckInit() initial behavior. But then you're always guaranteed correctness and it has well-defined more general semantics and debugging tools, since you'll be able to then investigate the true meaning behind your callback (because of what it means implicitly to the algebraic variables).

But again, this is in-progress and the MTK v10 and we don't have it done yet. We have already added in all of the extra correctness tests to make sure no callbacks continue with incorrect states, and this error that you're seeing, but now we are fixing up the error messages and moving the default handling to implicit discrete.