# ModelingToolkit -- can I use \`min\` and \`max\` functions?

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934>\
**Category:** Modelling & Simulations\
**Created:** [April 10, 2025, 2:53pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934 "2025-04-10T14:53:23Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 2:53pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/1 "2025-04-10T14:53:24Z")

</div>

I have written a PID controller in ModelingToolkit (yes, I know it probably exists…). I need to clip the control output `u` to lie in the range [u\_min, u\_max]. So I try to do:

```julia
# Constraining the unconstrained control signal u_u to the valid range [u_min, u_max]

u ~ max(u_min, min(u_max,u_u))

```

Problem is: when `u_u` reaches `u_max`, the simulation crashes.

A corresponding Modelica code doesn’t have any problems with doing a similar statement.

**Questions** :

1. Is there a problem with using `min` and `max`? (note: the code works while `u_u` lies in the interval `(u_min, u_max)`, so there must be some kind of allowance for these functions.
2. Or is the problem related to using the solver I use? … I have used `Rodas5()`.

At some stage, I also want to include anti-windup in the form of locking the integrator state while `u_u` lies outside of `[u_min, u_max]`.

I include my PID controller below.

```julia
# PID controller
# ===============
@mtkmodel control_PID begin
    # Model parameters
    @parameters begin
        # Controller parameters
        K_p = 1, [description = "Proportional gain"]
        T_i = 1, [description = "Integral time"]
        T_d = 1, [description = "Derivative time"]
        eta = 0.1, [description = "Ratio, error filter time constant/derivative time, -"]
        tau_f = T_d/eta, [description = "Filter time constant"]
        #
        kappa_p = K_p*(T_i+T_d)/T_i, [description = "Parallel P-gain"]
        kappa_i = K_p/T_i, [description = "Parallel I-gain"]
        kappa_d = K_p*T_d, [description = "Parallel D-gain"]
        #
        u_0 = 0.5, [description = "Assumed steady control signal"]
        u_min = 0.0, [description = "Minimum valid control signal"]
        u_max = 1.0, [description = "Maximum valid control signal"]
    end
    #
    # Model variables, with initial values needed
    @variables begin
        xi_i(t)=u_0/kappa_i,[description = "Integral state"]
        xi_f(t)=0, [description = "Filter state"]
        e(t), [input=true,description = "Control error"]
        u_u(t), [description = "Unconstrained control signal"]
        u(t), [output=true,description = "Control signal"]
    end
    #
    # Model equations
    @equations begin
            Dt(xi_i) ~ e 
            Dt(xi_f) ~ -(xi_f - e)/tau_f 
            u_u ~ kappa_p*e + kappa_i*xi_i + kappa_d*(e - xi_f/tau_f)
     # Constraining the unconstrained control signal u_u to the valid range [u_min, u_max]
        u ~ max(u_min, min(u_max,u_u))
    end
end
;

```

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [April 10, 2025, 3:46pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/2 "2025-04-10T15:46:44Z")

</div>

I have used both `clamp` and `max` in my MTK models, and they work fine. I am using the `FBDF` solver which is meant for stiff systems. Which solver are you using in Modelica? Maybe try using a similar solver in MTK. Also, using a [ContinuousCallback](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/) would maybe help here.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 4:37pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/3 "2025-04-10T16:37:29Z")

</div>

Here is the response with the `Rodas5()` solver:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/b/8b288a3cb5b9258438d737fb3ed350bf3cb6bee5.png)

So apparently the clipping works.

The `FBDF()` solver does slightly better, i.e., simulates the system a fraction of a second longer. [The abscissa ticks indicate time in _minutes_.]

Modelica uses the default solver, which is `dassl` with a tolerance of 1e-6.

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [April 10, 2025, 4:53pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/4 "2025-04-10T16:53:22Z")

</div>

Try plotting the differential variables `xi_i` or `xi_f`. I think the clamp might not be the reason the model is failing.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 5:00pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/5 "2025-04-10T17:00:51Z")

</div>

OK – I asked MS Co-pilot about which is the best DAE solver for Julia… Here is the answer:

* * *

For solving stiff differential-algebraic equations (DAEs) in Julia, several solvers stand out depending on your requirements:

1. **IDA from Sundials.jl** : Ideal for high-performance solutions with Float64 problems.
2. **DFBDF from OrdinaryDiffEq.jl** : A Julia-native option, great for adaptive-order and adaptive-time BDF methods.
3. **RadauIIA5** : Excellent for high-accuracy solutions.
4. **DASSL.jl** : Implements the DASSL algorithm for stiff DAEs using variable step-size BDF.

Each solver has its strengths, so the best choice depends on your problem’s specifics, such as mass matrix form, accuracy requirements, and computational constraints.

* * *

For almost all of these (don’t recall the exception), I am told that the solver is incompatible with the problem… (if one of them was compatible, it failed…).

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 5:11pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/6 "2025-04-10T17:11:43Z")

</div>

OK - I switched to `clamp` – the syntax is more elegant than the combination of `min` and `max`.

Still, the simulation crash with too large a disturbance (driving the control signal into saturation).

Btw… is there a trick to lock the integrator when `u_u` (the “un-clamped” control signal) is saturated? (Anti-windup).

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [April 10, 2025, 5:14pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/7 "2025-04-10T17:14:14Z")

</div>

> [@BLI](#):
>
> OK – I asked MS Co-pilot about which is the best DAE solver for Julia… Here is the answer:

Use this instead:

> **[SciMLBenchmarks.jl: Benchmarks for Scientific Machine Learning (SciML) and...](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/)**

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 5:19pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/8 "2025-04-10T17:19:28Z")

</div>

OK… I tried to _remove_ the `clamp` statement, thus leaving the control signal unconstrained.

Now the simulation crashes at an earlier time than when I constrain the signal.

I guess this indicates that the reason for the crashing is _not_ the constraint on the control signal, but something else.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 5:23pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/9 "2025-04-10T17:23:04Z")

</div>

Note: I have specified `tstops` at the locations where the changes in inputs (a disturbance) happens. But, of course, the crash is slightly after the `tstops` value.

Also, the system is rather oscillatoric.

Typical control signal response when I start with disturbance at 85% of max value, jumps up to 91%, jumps back down to 85%, then down to 50% and finally down to 20%.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/6/76328fed544b7ad94f1ca8fbb04407b0f65a892d.png)

If I start by jumping up from 85% to 92%, the system crashes. [But OpenModelica handles this case.]

* * *

I have tried to change the tolerance. If I increase the `abstol` to `1e-8`, the initialization of the problem fails.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 10, 2025, 7:21pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/10 "2025-04-10T19:21:52Z")

</div>

You might want to apply `IfLifting`, which is a custom pass that changes `min` and `max` to `if` statements that then get lowered to callbacks. That can make it more robust.

```julia
structural_simplify(sys, additional_passes = [IfLifting])

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 8:31pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/11 "2025-04-10T20:31:47Z")

</div>

Currently, I do:

```julia
@mtkbuild mod_L = Model_Large()

```

How can I insert `additional_passes = [IfLifting]` into this?

(`Large` does not mean that the model is large, but that I consider large disturbances…)

* * *

Btw, I now tend to suspect that it is model accuracy that is the problem… the error message I get is:

```julia
Warning: At t=138.6494837605311, dt was forced below floating point epsilon 2.842170943040401e-14, and step error estimate = 0.13043281293703687. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64)

```

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [April 10, 2025, 8:43pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/12 "2025-04-10T20:43:53Z")

</div>

> **[ModelingToolkit Language: Modeling with @mtkmodel, @connectors and @mtkbuild...](https://docs.sciml.ai/ModelingToolkit/stable/basics/MTKLanguage/#Build-structurally-simplified-models:)**
>
> Documentation for ModelingToolkit.jl.

But if removing the `clamp` didn’t solve the issue, this probably won’t help. Are you sure all of your formulas are correct?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 10, 2025, 9:08pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/13 "2025-04-10T21:08:18Z")

</div>

Yes, I am pretty sure the model is correct. With jumping up to 91%, the model output is virtually identical to what OpenModelica gives. I can go up to 91.5%, the result is still the same as in OpenModelica. But MTK/Julia crashes if I go up to 92%, while OpenModelica still works. But increasing the disturbance to more than 92% of the max disturbance, makes OpenModelica also crash.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 10, 2025, 9:30pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/14 "2025-04-10T21:30:18Z")

</div>

> [@BLI](#):
>
> ```julia
> @mtkbuild mod_L = Model_Large()
> 
> ```
> 
> How can I insert `additional_passes = [IfLifting]` into this?

Remove the `@mtkbuild` and do the structural simplification as a separate step.

> [@BLI](#):
>
> Yes, I am pretty sure the model is correct. With jumping up to 91%, the model output is virtually identical to what OpenModelica gives. I can go up to 91.5%, the result is still the same as in OpenModelica. But MTK/Julia crashes if I go up to 92%, while OpenModelica still works. But increasing the disturbance to more than 92% of the max disturbance, makes OpenModelica also crash.

This is because it if-lifts by default. We don’t yet because the pass is pretty new and experimental, but we’re getting there.

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [April 10, 2025, 11:09pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/15 "2025-04-10T23:09:14Z")

</div>

I think your system is slightly abusive of modeling systems. Your `u_u` has a direct feedthrough input from `e` plus integrated input in `xi_i` and `xi_f`. So the integrator has to deal with two+ sources of danger.

My guess is that MTK correctly detects the first rail. But once saturated it is still doing hard work to deal with fluctuations in the input. Adaptive solvers try to decrease step size to satisfy tolerances, except that feedthrough is unaffected by step size, hence the tiny time steps which error out. In the days before MTK/Modelica/Simulink, this would be a nightmare to code by hand, and stiff solvers wouldn’t help because they want to deal with stiff dynamics not stiff feedthrough.

You may be causing yourself issues by implementing `u_u` and `u` open-loop, where technically they’re not going through any dynamics. It might actually help to use them in closed-loop, because that at least should make the system causal (no feed-through). Then the system is merely stiff and not also DAE.

> [@BLI](#):
>
> Btw… is there a trick to lock the integrator when `u_u` (the “un-clamped” control signal) is saturated? (Anti-windup).

See this previous thread on anti-windup, with similar questions about why Modelica was able to handle it well. I had speculated that it has some clever bits to allow this:

> It appears Modelica can do things several ways. (1) Naively using adaptive ODE solver, where corners are tight-stepped. (2) Aware of corners and degree of smoothness which may help with polynomial order and stepping. (3) High-level symbolically aware when a signal is at limit which may help avoid if statements or discrete callbacks.

My impression is that MTK has adaptive time steps and stiff solvers, and can convert `min` to callbacks (if-lifting), but I don’t think it has Modelica’s degree of symbolic awareness.

As a workaround, you could implement back-calculation manually. It applies a feedback term that should slow (reduce stiffness in) the dynamics when saturated. But I would try to get rid of the feed-through first.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 11, 2025, 8:58am UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/16 "2025-04-11T08:58:37Z")

</div>

> [@apo383](#):
>
> I think your system is slightly abusive of modeling systems. Your `u_u` has a direct feedthrough input from `e` plus integrated input in `xi_i` and `xi_f`. So the integrator has to deal with two+ sources of danger.

What I included is a fairly standard implementation of a parallel form PID controller; no abuse. `xi_f` is just the state in a low-pass filter needed to make the derivative action realizable. I just included the controller – it is connected with the real system that I try to control (not shown, somewhat more complex) in a feedback structure.

* * *

Hm… I’m testing a possible anti-windup scheme. It seems to work, at least with the `Rodas5()` integrator; `FBDF()` gives an error message:

```julia
# PID controller
# ===============
@mtkmodel control_PID begin
    # Model parameters
    # -- structure parameters
    @structural_parameters begin
        # anti_windup = true: lock integrator state during clipping
        # anti-windup = false: no locking of integrator state
        anti_windup = false
    end
    #
    # -- parameters
    @parameters begin
        # Controller parameters
        K_p = 1, [description = "Proportional gain"]
        T_i = 1, [description = "Integral time"]
        T_d = 1, [description = "Derivative time"]
        eta = 0.1, [description = "Ratio, error filter time constant/derivative time, -"]
        tau_f = T_d/eta, [description = "Filter time constant"]
        #
        kappa_p = K_p*(T_i+T_d)/T_i, [description = "Parallel P-gain"]
        kappa_i = K_p/T_i, [description = "Parallel I-gain"]
        kappa_d = K_p*T_d, [description = "Parallel D-gain"]
        #
        u_0 = 0.5, [description = "Assumed steady control signal"]
        u_min = 0, [description = "Minimum valid control signal"]
        u_max = 1, [description = "Maximum valid control signal"]
    end
    #
    # Model variables, with initial values needed
    @variables begin
        xi_i(t)=u_0/kappa_i,[description = "Integral state"]
        xi_f(t)=0, [description = "Filter state"]
        e(t), [input=true,description = "Control error"]
        u_u(t), [description = "Unconstrained control signal"]
        u(t), [output=true,description = "Control signal"]
    end
    #
    # Model equations
    @equations begin
        # Parallel implmenentation
        if !anti_windup
            # No anti-windup
            Dt(xi_i) ~ e
        else
            # Anti-windup of integrator state
            Dt(xi_i) ~ e*(1-sign(abs(u_u - clamp(u_u,u_min,u_max))))
        end
        Dt(xi_f) ~ -(xi_f - e)/tau_f 
        u_u ~ kappa_p*e + kappa_i*xi_i + kappa_d*(e - xi_f/tau_f)
        # Constraining the unconstrained control signal u_u to the valid range [u_min, u_max]
        u ~ clamp(u_u, u_min, u_max) #max(u_min, min(u_max,u_u))
    end
end

```

I have put a severe constraint on the control signal, i.e., `u_max = 0.78` (the normal would be `u_max=1.0` – I do this just to test the algorithm):

CONTROL signal  
Without anti-windup:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/5/85a663fc2d9f0237d47789ead764a0430812a5b2.png)

With anti-windup:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/7/e730076e94d60c5bcd858475d21886508ef6a77a.png)

INTEGRATOR state  
Without anti-windup:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/2/82bfe1e4e635b4c0d1db0d3104ddcaf65af80843.png)

With anti-windup - note how the integrator state is frozen from ca. 2 min to 2.8 min or so; that is when this controller is clipping.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/c/7c92ec3def565e845c1bb9b7c1370c9fb724dda8.png)

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [April 11, 2025, 6:00pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/17 "2025-04-11T18:00:29Z")

</div>

> [@BLI](#):
>
> fairly standard implementation of a parallel form PID controller; no abuse

That’s fine, my problem is with `u_u` and `u` as algebraic outputs. They both depend on `e`, a non-differentiable input. DAE solvers usually differentiate algebraic equations to make them differential, impossible with `e`. Perhaps MTK is smart enough to eliminate `u_u` and `u` and not represent them since they’ll be integrated, but IMO only in closed-loop. In open loop it thinks you want `u` as an output, which needs to be represented somehow. That’s what I mean by feed-through being problematic.

In any case, your anti-windup example is very promising. ~~~Back-calculation’s feedback decreases the stiffness~~~ [EDIT: Noticed you’re clamping `xi_i`], making it easier for the integrator to deal with bumpy `e` [and avoiding a large integral].

However that raises the question why your no-anti-windup case works. Didn’t you complain about the integrator failing, or was that only in open loop? Perhaps your main issue was merely slow (not failing) integrator in closed-loop, and you used open loop for debugging, which failed? In which case the open loop is not a good example because it algebraically depends on `e`, whereas the closed-loop system is merely stiff. And MTK was smart enough to eliminate or deal with `u_u` and `u` since they get integrated.

My brief scan of OpenModelica code suggests they have a couple ways to deal with this, including some semantic/symbolic awareness. Their anti-windup code leverages this and has option to enforce a polynomial smoothness, or lock the integrator, or back-calculate. I believe Simulink has the latter two. I also believe MTK does not (yet) go as far as Modelica.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 11, 2025, 6:22pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/18 "2025-04-11T18:22:11Z")

</div>

> [@apo383](#):
>
> However that raises the question why your no-anti-windup case works.

There is no mystery with this. PID controllers without anti-windup may work well even when the control signal clips. However, they then tend to integrate up an unnecessarily large (absolute value) integrator state which takes a long time to reverse/“empty”, and thus leads to slow approach to steady state after lengthy periods of clipping. Also, this build-up may extend the period of clipping.

The main idea of anti-windup is to lock the integrator state while the control signal clips, to avoid this unnecessary build-up of integrator state.

* * *

My initial problem was how to implement (i) clipping, and (ii) anti-windup. I think I have now solved both of these, although solver `FBDF` fails when I use anti-windup. Possibly, I can remedy that problem by Chris’ proposal to use a stage which insert if-sentences in the underlying code. I haven’t had time to look into that.

* * *

Then I had a secondary problem in that the MTK/Julia code can not handle as large disturbance as the similar (Open-)Modelica code can. The difference is small, though. I observed that the MTK/Julia code can handle up to 91.5% of the maximum disturbance for my model (while doing so, the control signal stays below the maximum value, thus there is no clipping), while the Modelica code can handle up to 92%, and in doing so actually clips the control signal.

The reason for this could be either (i) the default solver in OpenModelica is superior to what I have found for MTK/Julia, or (ii) some slight differences in the Modelica code. I will check for such differences one more time – it is somewhat complicated on my 13 inch laptop screen to have both codes up at the same time… If there are differences, these will be tiny, because up to 91.5% disturbance, the response is virtually indistinguishable.

* * *

Anyways, I have learned from these discussions! Thanks to all you guys who have contributed with suggestions.

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [April 11, 2025, 8:14pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/19 "2025-04-11T20:14:45Z")

</div>

I agree with what you say. The closed-loop system is causal and better shows its issue is indeed with wind-up. I was suspicious of the open-loop example because it’s acausal, perhaps leading to the crash without clipping. I think open- and closed-loop crashed at different times and slightly different reasons. (Also I should have noticed you’re clamping in closed loop rather than back-calculating.)

[This previous thread](https://discourse.julialang.org/t/handling-of-saturated-integrator-with-anti-windeup/124819) discusses the differences with OpenModelica. They have nonlinear blocks that are saturation (limiter) aware:

> <https://github.com/modelica/ModelicaStandardLibrary/blob/daa10221dfd96caf4f4309190733715109ee5a3d/Modelica/Blocks/Nonlinear.mo#L33C2-L38C18>

This code allows several ways to saturate. It’s implemented explicitly with `if u > uMax...`. The `smooth(0,` is option to smooth with a given polynomial order, currently disabled with 0. The other branch with `homotopy` is for reinitializing integrator states after hitting a discrete change, difficult when algebraic equations contain `if`. They can iteratively solve starting with simplified knowledge of the rail.

I would expect MTK if-lifting to be great for identifying discrete changes. Not sure about initialization, but it’s possible to furnish MTK with guesses for solving algebraic equations, where the first guess is that system is stuck on rail.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 13, 2025, 2:44pm UTC](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934/20 "2025-04-13T14:44:38Z")

</div>

> [@apo383](#):
>
> I was suspicious of the open-loop example because it’s acausal, perhaps leading to the crash without clipping. I think open- and closed-loop crashed at different times and slightly different reasons

Please note: I have not shown any results in open-loop.

[Next page](https://discourse.julialang.org/t/modelingtoolkit-can-i-use-min-and-max-functions/127934.md?page=2)