Skip to content

Code generation contains Differential when equations contain variable time parameter #3480

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bradcarman opened this issue Mar 20, 2025 · 13 comments · May be fixed by #3493
Open

Code generation contains Differential when equations contain variable time parameter #3480

bradcarman opened this issue Mar 20, 2025 · 13 comments · May be fixed by #3493
Assignees
Labels
bug Something isn't working

Comments

@bradcarman
Copy link
Contributor

Take the following

@component function FilteredInput(; name, x0=0, T=0.1)
    params = @parameters begin
      k(t) = x0
      T = T
    end
    vars = @variables begin
      x(t) = k
      dx(t) = 0
      ddx(t)
    end
    systems = []
    eqs = [
      D(x) ~ dx
      D(dx) ~ ddx
      dx ~ (k - x)/T
    ]
    return ODESystem(eqs, t, vars, params; systems, name)
end

@mtkbuild sys = FilteredInput()
ModelingToolkit.observed(sys)

gives...

julia> ModelingToolkit.observed(sys)
3-element Vector{Equation}:
 dx(t) ~ (k(t) - x(t)) / T
 dxˍt(t) ~ (-dx(t) + Differential(t)(k(t))) / T
 ddx(t) ~ dxˍt(t)

Note the Differential(t)(k(t))) which cannot be computed.

Manifest:

[961ee093] ModelingToolkit v9.66.0
@bradcarman bradcarman added the bug Something isn't working label Mar 20, 2025
@baggepinnen
Copy link
Contributor

baggepinnen commented Mar 21, 2025

One problem here is that if k is a function of time, then the named variable ddx will depend on the time derivative of k, D(k), and there is no way for MTK to know what the derivative of this time-dependent parameter is. The problem goes away if you don't name ddx, and instead write the model as

julia> @component function FilteredInput(; name, x0=0, T=0.1)
           params = @parameters begin
             k(t) = x0
             T = T
           end
           vars = @variables begin
             x(t) = k
             dx(t) = 0
           end
           systems = []
           eqs = [
             D(x) ~ dx
             D(dx) ~ (k - x)/T
           ]
           return ODESystem(eqs, t, vars, params; systems, name)
       end

julia> @mtkbuild sys = FilteredInput()

julia> ModelingToolkit.observed(sys)
Equation[]

since then MTK has no reason to make ddx available as an observed, and the detivative of k never appears in any expression, it's only being integrated.

Furthermore, if you intend for k(t) to be piecewise constant, the derivative of k is infinite or ill defined at the time points when k changes.

@AayushSabharwal
Copy link
Member

Making parameter derivatives zero should handle this?

@ChrisRackauckas
Copy link
Member

We do need to make sure that when a parameter is a function, say a datainterpolation, that we still handle that case correctly.

@baggepinnen
Copy link
Contributor

Making parameter derivatives zero should handle this?

"handle"... if the derivative is zero then there's no need to have the parameter depend on time. If the parameter truly does vary with time it's wrong to just assume that the derivative is zero.

@ChrisRackauckas
Copy link
Member

Well the issue here is that discrete variables and parameter functions need to do this separately.

@baggepinnen
Copy link
Contributor

baggepinnen commented Mar 24, 2025

In both cases is it wrong to assume that the derivative is zero. A discrete-time variable that is piecewise constant has an ill-defined derivative (in the non-distribution sense at least) at the time point when it changes value, if the model is such that it depends on this derivative, the model is probably wrong. If the parameter is a function of time and MTK knows the derivative of this function it is fine, otherwise erroring is the only reasonable option.

@hersle
Copy link
Contributor

hersle commented Mar 25, 2025

See also #2823 (from #2823 (comment) and below). A good solution to this issue would probably also solve that issue.

@AayushSabharwal
Copy link
Member

We do need to handle parameter derivatives in the explicit case at least.

A discrete-time variable that is piecewise constant has an ill-defined derivative (in the non-distribution sense at least) at the time point when it changes value

Isn't that the point of it being discrete-time, though? The callback handles the instantaneous change in value, and the assumption of the derivative being zero holds at all other time points. Is there a case where it would give an incorrect result?

At the very least, I think we can allow explicitly specifying parameter derivatives and erroring when they're required but not provided.

@baggepinnen
Copy link
Contributor

Is there a case where it would give an incorrect result?

It hides a likely model error, MTK is not supposed to have to differentiate such a discrete-time parameter. Take the model in the OP as an example, ddx contains the time derivative of k which is an impulse, this is not mechanically feasible and the k should have been added to the acceleration rather than the velocity, i.e., a model error. One cannot perturb a velocity other than through a force, and forces affect the acceleration, they only affect velocity through integration of the acceleration.

Change the domain from mechanical to electrical or any other and you get the same problem, MTK should error if you try to impulsively disturb a variable in any other way than to explicitly change the variable in a callback. Note, it's k being modified by the callback here, not ddx, which would be implicitly disturbed impulsively.

@AayushSabharwal
Copy link
Member

Interesting. So time-dependent parameter derivatives should error unless the user explicitly provides a derivative for them?

@baggepinnen
Copy link
Contributor

If the derivative is of a parameter that changes discretely, then yes, that's my opinion at least. In the FMU case, the output is probably a continuous variable with a well-defined time derivative, and it's our problem that we call that a parameter 😅

The modelica crowd has a number of recent papers on how to do the right thing with Dirac impulses, but I don't think any tool implements this research yet.

@AayushSabharwal
Copy link
Member

In the FMU case, the output is probably a continuous variable with a well-defined time derivative, and it's our problem that we call that a parameter

Well yes and no. CoSimulation FMUs by nature are meant to be handled via callbacks, and their outputs/states are discrete variables. Discrete variables are stored in the parameter struct. If a model uses the output of such an FMU in an equation that MTK differentiates, then this is still valid. The FMU allows obtaining the output derivative, we just don't have a way to hook into it yet. The question then arises what to do if the FMU doesn't support calculating a high enough order derivative. Error?

@ChrisRackauckas
Copy link
Member

In both cases is it wrong to assume that the derivative is zero. A discrete-time variable that is piecewise constant has an ill-defined derivative (in the non-distribution sense at least) at the time point when it changes value, if the model is such that it depends on this derivative, the model is probably wrong. If the parameter is a function of time and MTK knows the derivative of this function it is fine, otherwise erroring is the only reasonable option.

This is missing context. It's like saying the derivative of the PDE equation cannot be directly calculated because there's a boundary condition which the derivative equation may not satisfy. The extra context is the domain. The definition of the eqs for the continuous variables are that those are defined on the domain of t where no clocks are being triggered. So yes, the derivative of a discrete variable x(t) is not zero. But for any continuous equation in eqs, which is defined as being valid only on the open interval time spans for which the discrete behavior is not triggered, the derivative is zero. So if we want to mince words, sure the derivative of x(t) is not zero, but it is correct to substitute D(x) => 0 for discretes in the context of the eqs array.

It hides a likely model error, MTK is not supposed to have to differentiate such a discrete-time parameter. Take the model in the OP as an example, ddx contains the time derivative of k which is an impulse, this is not mechanically feasible and the k should have been added to the acceleration rather than the velocity, i.e., a model error. One cannot perturb a velocity other than through a force, and forces affect the acceleration, they only affect velocity through integration of the acceleration.

There are mathematically reasonable equations to be considered for this though. I mean, you can have a pendulum where the length L is a discrete value that changes periodically. Mathematically that's a perfectly reasonable set of equations, and D(L) => 0 would be required in order for index reduction to work. That is precisely the case for which this rule is needed. That said, it's fine to want the ability to switch this off because I agree there could be something wonky in the model. But in general it's mathematically allowable and correct and so we should at least do the math right.

MTK should error if you try to impulsively disturb a variable in any other way than to explicitly change the variable in a callback.

You can change L only in callbacks and that's precisely the case where you need D(L) => 0 for index reduction. So yes, I agree MTK should error if there's an impulse anywhere except a callback, but that's orthogonal to this.

CoSimulation FMUs by nature are meant to be handled via callbacks, and their outputs/states are discrete variables. Discrete variables are stored in the parameter struct.

Yes, Cosimulation FMUs are clocked, so they have discrete inputs and outputs. https://www.claytex.com/tech-blog/fmi-basics-co-simulation-vs-model-exchange/ makes this very clear with a nice image.

A nice way to do this would be to add metadeta for dummy derivative substitutes, which when not supplied is 0. Then when going through the index reduction process, you check if this is supplied, if so you put that expression in (which may just be a different output of an FMU for example), and if not it's 0.

Even then the derivative substitution to zero is not incorrect, as within the mathematical with a finite dt clocked FMU it's correct that at any place where f is called the derivative is zero. Though it is interesting to think about if you had an FMU for L(t) in the pendulum which was continuously changing, the derivative corresponds to different implicit constraints. The implicit constraint with the zero derivative ensures that the energy drift of the actual discretized simulation is zero. The implicit constraint with the non-zero derivative means that the energy drift against the true solution is zero, which will bias the discretized simulation in a way that could be good. So you want to use that information if you can, but it's not entirely clear it's better in this context since it depends on the definition of "the correct energy" here...

The FMU allows obtaining the output derivative, we just don't have a way to hook into it yet. The question then arises what to do if the FMU doesn't support calculating a high enough order derivative. Error?

Yes, when that comes up that's definitely an error.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants