add parameters to PID controllers (#336)

* add parameters to PID controllers

* workaround for scoping issue

* remove a_end parameter

* test: fix depwarn

* format

---------

Co-authored-by: Aayush Sabharwal <[email protected]>

Author: baggepinnen

src/Blocks/continuous.jl

@@ -258,13 +258,21 @@ See also [`LimPID`](@ref)
     with_D = !isequal(Td, false)
     @named err_input = RealInput() # control error
     @named ctr_output = RealOutput() # control signal
-    !isequal(Ti, false) &&
-        (Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0")))
-    !isequal(Td, false) &&
-        (Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0")))
-    Nd > 0 || throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0"))
+    @symcheck Ti ≥ 0 ||
+              throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0"))
+    @symcheck Td ≥ 0 ||
+              throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0"))
+    @symcheck Nd > 0 ||
+              throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0"))
+
+    pars = @parameters begin
+        k = k, [description = "Proportional gain"]
+        Ti = Ti, [description = "Integrator time constant"]
+        Td = Td, [description = "Derivative time constant"]
+        Nd = Nd, [description = "Derivative limit"]
+    end
 
-    @named gainPID = Gain(k)
+    @named gainPID = Gain(; k)
     @named addPID = Add3()
     if with_I
         @named int = Integrator(k = 1 / Ti, x = int__x)
@@ -304,7 +312,7 @@ See also [`LimPID`](@ref)
     else
         push!(eqs, connect(Dzero.output, addPID.input3))
     end
-    ODESystem(eqs, t, [], []; name = name, systems = sys)
+    ODESystem(eqs, t, [], pars; name = name, systems = sys)
 end
 
 """
@@ -334,15 +342,22 @@ The simplified expression above is given without the anti-windup protection.
               throw(ArgumentError("Time constant `Ta` has to be strictly positive"))
     @symcheck T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive"))
     @symcheck u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max"))
+    pars = @parameters begin
+        k = k, [description = "Proportional gain"]
+        T = T, [description = "Integrator time constant"]
+        Ta = Ta, [description = "Tracking time constant"]
+        u_max = u_max, [description = "Upper saturation limit"]
+        u_min = u_min, [description = "Lower saturation limit"]
+    end
     @named err_input = RealInput() # control error
     @named ctr_output = RealOutput() # control signal
-    @named gainPI = Gain(k)
+    @named gainPI = Gain(; k)
     @named addPI = Add()
     @named addTrack = Add()
     @named int = Integrator(k = 1 / T, x = int__x)
     @named limiter = Limiter(y_max = u_max, y_min = u_min)
     @named addSat = Add(k1 = 1, k2 = -1)
-    @named gainTrack = Gain(1 / Ta)
+    @named gainTrack = Gain(k = 1 / Ta)
     sys = [err_input, ctr_output, gainPI, addPI, int, addTrack, limiter, addSat, gainTrack]
     eqs = [
         connect(err_input, addPI.input1),
@@ -357,7 +372,7 @@ The simplified expression above is given without the anti-windup protection.
         connect(addTrack.output, int.input),
         connect(int.output, addPI.input2)
     ]
-    ODESystem(eqs, t, [], []; name = name, systems = sys)
+    ODESystem(eqs, t, [], pars; name = name, systems = sys)
 end
 
 """
@@ -408,18 +423,25 @@ where the transfer function for the derivative includes additional filtering, se
         Ti = k / Ti
         Td = Td / k
     end
-    0 ≤ wp ≤ 1 ||
-        throw(ArgumentError("wp out of bounds, got $(wp) but expected wp ∈ [0, 1]"))
-    0 ≤ wd ≤ 1 ||
-        throw(ArgumentError("wd out of bounds, got $(wd) but expected wd ∈ [0, 1]"))
-    !isequal(Ti, false) &&
-        (Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0")))
-    !isequal(Td, false) &&
-        (Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0")))
+    @symcheck Ti ≥ 0 ||
+              throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0"))
+    @symcheck Td ≥ 0 ||
+              throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0"))
     @symcheck u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max"))
     @symcheck Nd > 0 ||
               throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0"))
 
+    pars = @parameters begin
+        k = k, [description = "Proportional gain"]
+        Ti = Ti, [description = "Integrator time constant"]
+        Td = Td, [description = "Derivative time constant"]
+        wp = wp, [description = "Set-point weighting in the proportional part"]
+        wd = wd, [description = "Set-point weighting in the derivative part"]
+        Ni = Ni, [description = "Anti-windup tracking gain"]
+        Nd = Nd, [description = "Derivative limit"]
+        u_max = u_max, [description = "Upper saturation limit"]
+        u_min = u_min, [description = "Lower saturation limit"]
+    end
     @named reference = RealInput()
     @named measurement = RealInput()
     @named ctr_output = RealOutput() # control signal
@@ -431,7 +453,7 @@ where the transfer function for the derivative includes additional filtering, se
         if with_AWM
             @named addI = Add3(k1 = 1, k2 = -1, k3 = 1)
             @named addSat = Add(k1 = 1, k2 = -1)
-            @named gainTrack = Gain(1 / (k * Ni))
+            @named gainTrack = Gain(k = 1 / (k * Ni))
         else
             @named addI = Add(k1 = 1, k2 = -1)
         end
@@ -492,7 +514,7 @@ where the transfer function for the derivative includes additional filtering, se
         push!(eqs, connect(Dzero.output, addPID.input2))
     end
 
-    ODESystem(eqs, t, [], []; name = name, systems = sys)
+    ODESystem(eqs, t, [], pars; name = name, systems = sys)
 end
 
 """
@@ -611,13 +633,13 @@ See also [`StateSpace`](@ref) which handles MIMO systems, as well as [ControlSys
             description = "Denominator coefficients of transfer function (e.g., `s² + 2ωs + ω^2` is specified as [1, 2ω, ω^2])"
         ]
         bb[1:(nbb + nb)] = [zeros(nbb); b]
-        d = bb[1] / a[1], [description = "Direct feedthrough gain"]
     end
+    d = bb[1] / a[1]# , [description = "Direct feedthrough gain"]
 
     a = collect(a)
-    @parameters a_end = ifelse(a[end] > 100 * symbolic_eps(sqrt(a' * a)), a[end], 1.0)
+    a_end = ifelse(a[end] > 100 * symbolic_eps(sqrt(a' * a)), a[end], 1.0)
 
-    pars = [collect(b); a; collect(bb); d; a_end]
+    pars = [collect(b); a; collect(bb)]
     @variables begin
         x(t)[1:nx] = zeros(nx),
         [description = "State of transfer function on controller canonical form"]

src/Blocks/nonlinear.jl

@@ -19,7 +19,7 @@ Limit the range of a signal.
 @component function Limiter(; name, y_max, y_min = y_max > 0 ? -y_max : -Inf)
     @symcheck y_max ≥ y_min || throw(ArgumentError("`y_min` must be smaller than `y_max`"))
     m = (y_max + y_min) / 2
-    @named siso = SISO(u_start = m, y_start = m) # Default signals to center of saturation to minimize risk of saturation while linearizing etc.
+    siso = SISO(u_start = m, y_start = m, name = :siso) # Default signals to center of saturation to minimize risk of saturation while linearizing etc.
     @unpack u, y = siso
     pars = @parameters y_max=y_max [description = "Maximum allowed output of Limiter $name"] y_min=y_min [
         description = "Minimum allowed output of Limiter $name"

test/Blocks/continuous.jl

@@ -430,7 +430,7 @@ end
         @named pt1 = TransferFunction(b = [1.2], a = [3.14, 1])
         @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
         sys = structural_simplify(iosys)
-        prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+        prob = ODEProblem(sys, Pair[], (0.0, 100.0))
         sol = solve(prob, Rodas4())
         @test sol.retcode == Success
         @test sol[pt1.output.u]≈pt1_func.(sol.t, 1.2, 3.14) atol=1e-3
@@ -440,7 +440,7 @@ end
         @named pt1 = TransferFunction(b = [1.2], a = [3.14, 0])
         @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
         sys = structural_simplify(iosys)
-        prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+        prob = ODEProblem(sys, Pair[], (0.0, 100.0))
         sol = solve(prob, Rodas4())
         @test sol.retcode == Success
         @test sol[pt1.output.u] ≈ sol.t .* (1.2 / 3.14)
@@ -464,7 +464,7 @@ end
         @named pt1 = TransferFunction(b = [w^2], a = [1, 2d * w, w^2])
         @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
         sys = structural_simplify(iosys)
-        prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+        prob = ODEProblem(sys, Pair[], (0.0, 100.0))
         sol = solve(prob, Rodas4())
         @test sol.retcode == Success
         @test sol[pt1.output.u]≈pt2_func.(sol.t, k, w, d) atol=1e-3
@@ -474,7 +474,7 @@ end
         @named pt1 = TransferFunction(b = [1, 0], a = [1, 1])
         @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
         sys = structural_simplify(iosys)
-        prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+        prob = ODEProblem(sys, Pair[], (0.0, 100.0))
         sol = solve(prob, Rodas4())
         @test sol.retcode == Success
         @test sol[pt1.output.u]≈1 .- pt1_func.(sol.t, 1, 1) atol=1e-3
@@ -484,7 +484,7 @@ end
         @named pt1 = TransferFunction(b = [2.7], a = [pi])
         @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
         sys = structural_simplify(iosys)
-        prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+        prob = ODEProblem(sys, Pair[], (0.0, 100.0))
         sol = solve(prob, Rodas4())
         @test sol.retcode == Success
         @test all(==(2.7 / pi), sol[pt1.output.u])

test/Mechanical/translational_modelica.jl

@@ -105,5 +105,5 @@ end
     prob = ODEProblem(sys, [], (0, 2pi))
     sol = solve(prob, Rodas4())
     tv = 0:0.1:(2pi)
-    @test sol(tv, idxs = sys.mass.s)≈@.(2sin(2pi * tv * 3)) atol=1e-2
+    @test sol(tv, idxs = sys.mass.s).u≈@.(2sin(2pi * tv * 3)) atol=1e-2
 end