don't error on implicit parameter equation

[hexaeder]

I use MTK to generate "open loop" models. In closed loop loop, it is possible that the input of a system directly depends on its output, i.e. $u_a = f(y_a)$.

 ┌────────┐y_a = u_b┌────────┐  
 │System A├──────→──┤System B│
 └────────┘         └────────┘
     └─────────←────────┘
     u_a       =      y_b

This means you may implicitly force $y_a$ by constraining $u_a$. Under those conditions, it makes sense to define implicit equations on inputs (which are handled as parameters internally). MTK@10 does not allow this anymore, because it sees an equation containing only parameters and requires it to be an explicit parameter equation.

Instead of throwing, I think it should just leave the equation as it is.
Alternatively, one could check if any of the parameters is an input and throws otherwise.

Example:

    @mtkmodel ImplicitForcing begin
        @variables begin
            u(t), [description = "Input Variable", input=true]
            y(t), [description = "fully implicit output", output=true]
        end
        @equations begin
            0 ~ sqrt(u) # implicitly forces output y because u=f(y) in  closed loop
        end
    end
    @named implicit = ImplicitForcing()
    mtkcompile(implicit; inputs = unbound_inputs(implicit)) # errors

ERROR: ArgumentError: LHS of parameter dependency equation must be a single parameter. Found 0.

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

[AayushSabharwal]

The MWE doesn't error for me on MTK#master?

julia> using ModelingToolkit
[ Info: Auto-importing using ModelingToolkit: t_nounits as t, D_nounits as D, ModelingToolkit as MTK, observed

julia> @mtkmodel ImplicitForcing begin
           @variables begin
               u(t), [description = "Input Variable", input=true]
               y(t), [description = "fully implicit output", output=true]
           end
           @equations begin
               0 ~ sqrt(u) # implicitly forces output y because u=f(y) in  closed loop
           end
       end
ModelingToolkit.Model{typeof(__ImplicitForcing__), Dict{Symbol, Any}}(__ImplicitForcing__, Dict{Symbol, Any}(:variables => Dict{Symbol, Dict{Symbol, Any}}(:y => Dict(:output => true, :type => Real, :description => "fully implicit output"), :u => Dict(:type => Real, :description => "Input Variable", :input => true)), :kwargs => Dict{Symbol, Dict}(:y => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real), :u => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real)), :independent_variable => :t, :equations => Any["0 ~ sqrt(u)"]), false)

julia> @named implicit = ImplicitForcing()
Model implicit:
Equations (1):
  1 standard: see equations(implicit)
Unknowns (2): see unknowns(implicit)
  u(t): Input Variable
  y(t): fully implicit output

julia> mtkcompile(implicit; inputs = ModelingToolkit.unbound_inputs(implicit)) # errors
Model implicit:
Parameters (1): see parameters(implicit)
  u(t): Input Variable

[hexaeder]

I'm confused. On d1246c0 i get the error I reported. On master from 1h ago (847c41f) I reproduce your behavior and the equation vanishes completely. Thats also strange behavior because I think it cannot be resolved for u?

Now MTK seems to just ignore the equation and is happy to solve with a parameter which does not fullfil the constraint. See below MWE where I default u to 1 which is incompatible with the equations so it should throw init failure.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D, unbound_inputs
using OrdinaryDiffEqTsit5
@mtkmodel ImplicitForcing begin
    @variables begin
        u(t)=1, [description = "Input Variable", input=true]
        y(t), [description = "fully implicit output", output=true]
    end
    @equations begin
        0 ~ sqrt(u)
    end
end
@named implicit = ImplicitForcing()
sys = mtkcompile(implicit; inputs = unbound_inputs(implicit)) # test for no error

prob = ODEProblem(sys, [], (0.0, 1.0))
sol = solve(prob, Tsit5()) # happy to solve even though constrsaint is not fulfilled

[AayushSabharwal]

It won't automatically throw an initial failure. The equation needs to be moved to initialization_eqs.

[hexaeder]

I mean i don't think the equation should be dropped, and why is it? I would expect it to be a a system with one state y(t) and one equation 0~sqrt(u) and the DAE solver would throw init failure because it can't find a y which satisfies the equation.

[AayushSabharwal]

I mean i don't think the equation should be dropped, and why is it?

Because it doesn't involve any variables.

would expect it to be a a system with one state y(t) and one equation 0~sqrt(u) and the DAE solver would throw init failure because it can't find a y which satisfies the equation.

Exactly. It should throw an initialization failure. So the equation should be an initialization equation. It's not just "one state" and "one equation". The state has nothing to do with that particular equation.

[hexaeder]

But isn't the fact that an input variable is transformed into a parameter in structural simplify more of an implementation detail? My whole point is that I think equations containing inputs cannot be handled the same way as equations only containing fixed parameters.

[AayushSabharwal]

Hmmm. Alright, so now the question is where this equation disappears to.

[hexaeder]

I guess the recent change which drops the equations rather than keeping it was introduced in #3681 (or any other PR in the last 2 days but i think this one is the most promient as it changes alot in how tearing is handled)

[AayushSabharwal]

Okay so the way this worked previously is arguably by accident.

First, the system is technically overdetermined since it has one equation for zero unknowns (y doesn't count because it isn't involved in any equations and is removed by structural_simplify/mtkcompile). However, because the equation is a pure parameter equation it passes the consistency check since the check only considers equations that involve some unknowns. This behavior seems intentional.

Now, tearing for systems goes through tearing_reassemble for both fully- and under/over-determined systems. In v9, tearing_reassemble sorts all equations topologically and codegens them. To handle underdetermined systems, tearing_reassemble is passed extra information (full_var_eq_matching) specifically in the under/over-determined case. This extra information acts as a flag telling tearing that "hey the system isn't fully determined". However, the exact same information is available in the fully determined case, it just isn't passed. #3681 passes all information in both cases because it is necessary for BLT sorting, and differentiates the two cases by passing the fully_determined keyword to tearing_reassemble.

Thus in the MWE of this PR, the input equation is an "extra" equation that should be handled by the fully_determined = false logic. However, it escapes the consistency check and lands in the fully_determined = true case which ignores it since it isn't part of the BLT sorting and fully-determined systems don't have equations that aren't part of the BLT sorting.

I believe the fix here is to make the consistency check more strict and simplify this MWE with fully_determined = false.
This will then also require the change to process_parameter_equations in this PR. However, I think the check there should ensure that implicit parameter-only equations contain inputs.

[hexaeder]

Thanks for the clarification, this explains why it just worked in v9.

Regarding the bigger picture: the actual problem seems to be the handling of fully implicit outputs.
In this case, 0 ~ sqrt(u) is in something like the limit M*D(y) ~ sqrt(u) for M->0. For input-output systems it kinda makes sens to have "fully implicit" outputs as they appear naturally for M->0. In this sense, the system is also not overdetermined, because it has a free output variable and one equation. However I cannot tell MTK to use outputs = [y] because y does not show up in the equations anymore.

Previously I hacked around this problem using

_to_zero(x) = 0
@register_symbolic _to_zero(x)
eqs = _tozero(y) ~ sqrt(u)

which "tricks" MTK into doing the right thing but is not very obvious to users. This still seems to work on master:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D, unbound_inputs, unbound_outputs

_tozero(x) = 0
@register_symbolic _tozero(x)

@mtkmodel ImplicitForcing begin
    @variables begin
        u(t), [description = "Input Variable", input=true]
        y(t), [description = "fully implicit output", output=true]
    end
    @equations begin
        _tozero(y) ~ sqrt(u) # implicitly forces output y because u=f(y) in  closed loop
    end
end
@named implicit = ImplicitForcing()
mtkcompile(implicit; inputs = unbound_inputs(implicit), outputs=unbound_outputs(implicit)) # works

Maybe the better solution for the problem would be to somehow allow to specify fully implicit outputs on the mtkcompile level, rather than dropping the output as i did before.