Add control systems examples to docs

[zekeriyasari]

Include some examples from the control theory to the documentation of Jusdl.

[zekeriyasari]

#21

[ValentinKaisermayer]

I propose the following example:

PI-Controller for a System with Actuator Saturation

The following example visualizes the wind-up effect.

using Jusdl, Plots

"""
System.
"""
@def_ode_system mutable struct Sys{RH,RO,ST,IP,OP} <: AbstractODESystem 
    righthandside::RH = function ode(dx, x, u, t)
        dx[1] = x[2]
        dx[2] = -1. * x[1] - 0.5 * x[2] + u[1](t)
    end    
    readout::RO = (x, u, t) -> 0.9 * x[1] + x[2]
    state::ST = [0., 0.]
    input::IP = Inport() 
    output::OP = Outport() 
end

"""
Time continuos PI controller.
"""
@def_ode_system mutable struct PIController{RH,RO,ST,IP,OP} <: AbstractODESystem 
    kp::Float64 = 1.
    Tn::Float64 = 0.1
    righthandside::RH = (dx, x, u, t, kp = kp, Tn = Tn) -> (dx[1] = kp / Tn * u[1](t))    
    readout::RO = (x, u, t, kp = kp) -> kp * u[1](t) + x[1]
    state::ST = [0.]
    input::IP = Inport() 
    output::OP = Outport() 
end

"""
Time continuos PI controller with actuator saturation.
"""
@def_ode_system mutable struct PIControllerWithSat{RH,RO,ST,IP,OP} <: AbstractODESystem 
    kp::Float64 = 1.
    Tn::Float64 = 0.1
    lb::Float64 = -1.
    ub::Float64 = 1.
    righthandside::RH = (dx, x, u, t, kp = kp, Tn = Tn, lb = lb, ub = ub) -> (dx[1] = kp / Tn * u[1](t))    
    readout::RO = (x, u, t, kp = kp, ub = ub, lb = lb) -> max(min(kp * u[1](t) + x[1], ub), lb)
    state::ST = [0.]
    input::IP = Inport() 
    output::OP = Outport() 
end

"""
Time continuos PI controller actuator saturation and anti-windup.
"""
@def_ode_system mutable struct PIControllerWithSatAndAntiWindup{RH,RO,ST,IP,OP} <: AbstractODESystem 
    kp::Float64 = 1.
    Tn::Float64 = 0.1
    lb::Float64 = -1.
    ub::Float64 = 1.
    Ta::Float64 = 0.1
    righthandside::RH = function ode(dx, x, u, t, kp=kp, Tn=Tn, Ta=Ta, lb=lb, ub=ub)
        dx[1] = kp / Tn * u[1](t) + 1 / Ta * (-(x[1] + kp * u[1](t)) + max(min(kp * u[1](t) + x[1], ub), lb))
    end
    readout::RO = (x, u, t, kp = kp, lb = lb, ub = ub) -> max(min(kp * u[1](t) + x[1], ub), lb)
    state::ST = [0.]
    input::IP = Inport() 
    output::OP = Outport() 
end

# Construct the model
@defmodel model1 begin
    @nodes begin
        step = StepGenerator(delay=1.) 
        sys = Sys() 
        sum = Adder(signs=(+, -))
        pi_contr = PIController(kp=3., Tn=0.5)
        writer = Writer(input=Inport(3))
        mem = Memory(delay=0.01)
    end
    @branches begin
        step[1] => sum[1]
        sys[1] => mem[1]
        mem[1] => sum[2]
        sum[1] => pi_contr[1]
        pi_contr[1] => sys[1]
        step[1] => writer[1]
        sys[1] => writer[2]
        pi_contr[1] => writer[3]
    end
end

# Construct the model without anti-windup measure
@defmodel model2 begin
    @nodes begin
        step = StepGenerator(delay=1.)  
        sys = Sys() 
        sum = Adder(signs=(+, -))
        pi_contr = PIControllerWithSat(kp=3., Tn=0.5, lb=-1.5, ub=1.5)
        writer = Writer(input=Inport(3))
        mem = Memory(delay=0.01)
    end
    @branches begin
        step[1] => sum[1]
        sys[1] => mem[1]
        mem[1] => sum[2]
        sum[1] => pi_contr[1]
        pi_contr[1] => sys[1]
        step[1] => writer[1]
        sys[1] => writer[2]
        pi_contr[1] => writer[3]
    end
end

# Construct the model with anti-windup measure
@defmodel model3 begin
    @nodes begin
        step = StepGenerator(delay=1.) 
        sys = Sys() 
        sum = Adder(signs=(+, -))
        pi_contr = PIControllerWithSatAndAntiWindup(kp=3., Tn=0.5, lb=-1.5, ub=1.5, Ta=0.1)
        writer = Writer(input=Inport(3))
        mem = Memory(delay=0.01)
    end
    @branches begin
        step[1] => sum[1]
        sys[1] => mem[1]
        mem[1] => sum[2]
        sum[1] => pi_contr[1]
        pi_contr[1] => sys[1]
        step[1] => writer[1]
        sys[1] => writer[2]
        pi_contr[1] => writer[3]
    end
end

# call simulate
sim1 = simulate!(model1, 0, 0.01, 6.)
sim2 = simulate!(model2, 0, 0.01, 6.)
sim3 = simulate!(model3, 0, 0.01, 6.)

# Read and plot data 
t, x = read(getnode(model1, :writer).component)
p1 = plot(t, x[:, 1], label="r(t)", xlabel="t", legend=:bottomright)
plot!(p1, t, x[:, 2], label="y(t) without actuator saturation", xlabel="t")
p2 = plot(t, x[:, 3], label="u(t) without actuator saturation", xlabel="t")

t, x = read(getnode(model2, :writer).component)
plot!(p1, t, x[:, 2], label="y(t) without anti-windup", xlabel="t")
plot!(p2, t, x[:, 3], label="u(t) without anti-windup", xlabel="t")

t, x = read(getnode(model3, :writer).component)
plot!(p1, t, x[:, 2], label="y(t) with anti-windup", xlabel="t")
plot!(p2, t, x[:, 3], label="u(t) with anti-windup", xlabel="t")
plot(p1, p2, layout=(2, 1))

[zekeriyasari]

Looks great :) I wonder the details of the subject (i.e. the wind-up effect). Would you mind providing me some docs or related sources?

[ValentinKaisermayer]

Sure!
The Matlab page is not bad. And also the Wikipedia article.

A short description from me:
When designing control loops, it should always be noted that in practical applications the manipulated variable is limited in terms of amount, i.e. it cannot take on any value. If the output of the PI controller is in this limit, the control loop is essentially in an "open-loop". Since, no change of the measured variable influences the manipulated variable. The output of the PI-controller is limited but the value of the integrator of the I-part integrates on and on. Now, if the control error changes sign it can take a long time to "down-integrate" the value of the integrator. Hence, large over- and undershoots are typical for the windup effect. There are two main methods in order to mitigate this effect: back-calculation or clamping of the integrator.
In the example above the system can not follow the step change and due to the actuator limit the windup effect occurs. If one would not demand a step change or at least a smaller step change the windup effect may not occur.

Hence, this is a very relevant topic for practical implementations and should therefore be always included in simulation studies when designing controllers.

[zekeriyasari]

Thank you for the devotion of your time and explanation.

[ValentinKaisermayer]

You are welcome! I really like this package. I think it has great potential.

[zekeriyasari]

Thank you for your kind words. I just want to let you know that the package is under renaming, but I could not come up with a good name yet. I would appreciate any suggestions. In case you are interested, you may join in the discussion in #27