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

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.

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())

Unfortunately there is currently no delay component in the standard library

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?

Yes, you can program a component that implements the Padé approximant - Wikipedia …

See also: RATIONAL APPROXIMATION OF TIME DELAY

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

@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

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

1 Like

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. Unfortunately, ControlSystems cannot model both time delays and static nonlinearities at the same time, only one or the other :confused:

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.

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

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?

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…

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

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).

Here’s something to get you started

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)

3 Likes

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

Thanks a lot!

1 Like