Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add parameters to PID controllers #336

Merged
merged 5 commits into from
Oct 16, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 45 additions & 23 deletions src/Blocks/continuous.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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

"""
Expand Down Expand Up @@ -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),
Expand All @@ -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

"""
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

"""
Expand Down Expand Up @@ -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"]
Expand Down
2 changes: 1 addition & 1 deletion src/Blocks/nonlinear.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
10 changes: 5 additions & 5 deletions test/Blocks/continuous.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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])
Expand Down
2 changes: 1 addition & 1 deletion test/Mechanical/translational_modelica.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading