# Model a first-order-plus-delay system in ModelingToolkit

**URL:** <https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [February 4, 2024, 10:21am UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695 "2024-02-04T10:21:43Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![VPBML](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vpbml/32/212617_2.png) [@VPBML](https://discourse.julialang.org/u/VPBML)\
**Post date:** [February 4, 2024, 10:21am UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/1 "2024-02-04T10:21:44Z")

</div>

Hi,  
I wonder if there is any straightforward way to model a system consisting of a first-order-plus-delay plant T(s) = K.exp(-delay.s)/(tau.s + 1) and a PI controller using components from ModelingToolkitStandardLibrary. I built a system with first-order plant T(s) = K/(tau.s + 1) and a PI controller, but do not know how to introduce the delay part. I would really appreciate it if you guys can give me some guidance.

```julia
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using OrdinaryDiffEq

@parameters t

@named P = FirstOrder(k = 1, T = 1, x = 2.0)
@named C = LimPI(k = 1.0, T = 1.0, Ta = 1.0, int__x = 10.0, u_max = 100, u_min = 0)
@named ref = Blocks.Constant(k = 2.0)
@named feedback = Blocks.Feedback()

connections = [connect(ref.output, feedback.input1)
            connect(P.output, feedback.input2)
            connect(feedback.output, C.err_input)
            connect(C.ctr_output, :plant_input, P.input)]

closed_loop = ODESystem(connections, t, systems = [P,C,ref,feedback], name = :closed_loop)
sys = structural_simplify(closed_loop)
prob = ODEProblem(sys, Pair[], (0.0, 20.0))
sol = solve(prob, Rodas5())

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2024, 10:39am UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/2 "2024-02-04T10:39:41Z")

</div>

Unfortunately there is currently no delay component in the standard library

---

<div class="post-metadata">

**Author:** ![VPBML](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vpbml/32/212617_2.png) [@VPBML](https://discourse.julialang.org/u/VPBML)\
**Post date:** [February 4, 2024, 2:41pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/3 "2024-02-04T14:41:27Z")

</div>

> [@VPBML](#):
>
> model a system consisting of a first-order-plus-delay plant T(s) = K.exp(-delay.s)/(tau.s + 1) and a PI controller

Thanks. Is there any way I can model a system consisting of a first-order-plus-delay plant T(s) = K.exp(-delay.s)/(tau.s + 1) and a PI controller in ModelingToolkit?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [February 4, 2024, 2:48pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/4 "2024-02-04T14:48:18Z")

</div>

Yes, you can program a component that implements the [Padé approximant - Wikipedia](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) …

See also: [RATIONAL APPROXIMATION OF TIME DELAY](https://www2.humusoft.cz/www/papers/tcp09/035_hanta.pdf)

> **[Custom Components · ModelingToolkitStandardLibrary.jl](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/custom_component/)**
>
> Documentation for ModelingToolkitStandardLibrary.jl.

Sometimes it is easier just to write the differential equations in ModelingToolkit notation without using any components… A PI controller could look like this:

```julia
@variables t Pgc(t) e_ω(t)
@parameters P=1 I=1

D = Differential(t)

# PI speed controller, e_ω is the input and Pgc the output
D(Pgc) ~ P * D(e_ω) + I * e_ω,

```

s-function for Padé delay: [https://ris.utwente.nl/ws/portalfiles/portal/134422804/Some\_remarks\_on\_Pade-approximations.pdf](https://ris.utwente.nl/ws/portalfiles/portal/134422804/Some_remarks_on_Pade-approximations.pdf)

If you know the s-function of your approximated delay, you still need to transfer it to a differential equation, see: [Single Diff Eq → Transfer Function](https://lpsa.swarthmore.edu/Representations/SysRepTransformations/TF2SDE.html)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2024, 4:23pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/5 "2024-02-04T16:23:04Z")

</div>

I’m not sure if there is an interface to the delay solvers exposed from MTK. Since you are this far talking about a linear system and a linear controller you could model and simulate this using [ControlSystems.jl](https://juliacontrol.github.io/ControlSystems.jl/stable/examples/smith_predictor/). Unfortunately, ControlSystems cannot model _both_ time delays _and_ static nonlinearities at the same time, only one or the other 😕

Padé approximants introduce artifacts that are usually undesirable for time-domain simulations, so only use this if you are okay with the non-minimum phase zeros they introduce.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [February 4, 2024, 4:31pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/6 "2024-02-04T16:31:56Z")

</div>

> [@baggepinnen](#):
>
> non-minimum phase zeros they introduce.

Does this approximation introduces a non-minimum phase zero?

\frac{1}{1+sT}

or this one?

\frac{6-2sT}{6 + 4sT + (sT)²}

See: [https://ris.utwente.nl/ws/portalfiles/portal/134422804/Some\_remarks\_on\_Pade-approximations.pdf](https://ris.utwente.nl/ws/portalfiles/portal/134422804/Some_remarks_on_Pade-approximations.pdf)

---

<div class="post-metadata">

**Author:** ![VPBML](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vpbml/32/212617_2.png) [@VPBML](https://discourse.julialang.org/u/VPBML)\
**Post date:** [February 4, 2024, 5:22pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/7 "2024-02-04T17:22:57Z")

</div>

Thanks for the answer and suggestion. I took a quick look and ControlSystems.jl seems to be a better choice. Do you know if we can simulate the behavior of the above system (FOPD plant with PI controller) with the presence of output noise somehow?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [February 4, 2024, 5:58pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/8 "2024-02-04T17:58:51Z")

</div>

If your system is linear, use ControlSystems.jl, otherwise ModelingToolkit… You can also derive a linearized model from ModelingToolkit and then use that in ControlSystems…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2024, 6:10pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/9 "2024-02-04T18:10:59Z")

</div>

Yes you can, but this has to be done in discrete time. You need to

1. Discretize the plant model using the function `c2d`
2. Compute the closed-loop transfer function from output noise n to plant output z, this is given by (assuming measured output y = z + n)  
G\_n = \dfrac{-PC}{1 + PC}

which can be computed by calling `Gn = feedback(P*C) # Ignoring the - sign` where `P` is your discretized plant and `C` is a discrete-time PI controller. You can create this using the function `pid` with the keyword argument `Ts` to set the sample rate to get a discrete-time controller.

Once you have the transfer function `Gn`, you simulate this with random noise as input

```julia
N = 1000 # The number of input samples
n = randn(1, N) # Random noise input
result = lsim(Gn, n)
using Plots
plot(result)

```

If you further have a reference input, you can add r to n since the transfer function is the same from reference to output as from measurement noise to output (except for the - sign from noise, but this can be ignored since the noise is random anyways).

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2024, 6:17pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/10 "2024-02-04T18:17:21Z")

</div>

Here’s something to get you started

```julia
using ControlSystemsBase, Plots
τ = 1 # Delay time
Ts = 0.1 # Sampling time

P = tf(1, [1, 1]) * delay(τ)
Pd = c2d(P, Ts) # Discretized plant using ZoH discretization with sample time Ts
C = pid(1, 1; Ts) # A very well-tuned PI controller
Gn = feedback(Pd*C) # The closed-loop transfer function from r and n to plant output

N = 1000 # The number of input samples
n = 0.1 * randn(1, N) # Random noise input with std 0.1
r = sign.(sin.(0.2 * (0:Ts:(N-1)*Ts)))' # Square wave reference input

result = lsim(Gn, r - n) # - n in case n is not actually random
plot(result)
plot!(result.t, r', label="Reference", color=:black, linestyle=:dash, linewidth=2)

```

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

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [February 4, 2024, 7:07pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/11 "2024-02-04T19:07:09Z")

</div>

The latter yes – it has a NMP zero at 3/T.

---

<div class="post-metadata">

**Author:** ![VPBML](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vpbml/32/212617_2.png) [@VPBML](https://discourse.julialang.org/u/VPBML)\
**Post date:** [February 5, 2024, 5:34pm UTC](https://discourse.julialang.org/t/model-a-first-order-plus-delay-system-in-modelingtoolkit/109695/12 "2024-02-05T17:34:35Z")

</div>

Thanks a lot!