Merge pull request #77 from baggepinnen/highpass

allow `FirstOrder` to operate as highpass instead of lowpass

Author: ChrisRackauckas

src/Blocks/continuous.jl

@@ -61,15 +61,22 @@ function Derivative(; name, k=1, T, x_start=0.0)
 end
 
 """
-    FirstOrder(; name, k=1, T, x_start=0.0)
+    FirstOrder(; name, k=1, T, x_start=0.0, lowpass=true)
 
-A first-order filter with a single real pole in `s = -T` and gain `k`. The transfer function
+A first-order filter with a single real pole in `s = -T` and gain `k`. If `lowpass=true` (default), the transfer function
 is given by `Y(s)/U(s) = `
 ```
    k   
 ───────
 sT + 1
 ```
+and if `lowpass=false`, by
+```
+sT + 1 - k   
+──────────
+  sT + 1
+```
+
 
 # Parameters:
 - `k`: Gain
@@ -82,15 +89,15 @@ sT + 1
 
 See also [`SecondOrder`](@ref)
 """
-function FirstOrder(; name, k=1, T, x_start=0.0)
+function FirstOrder(; name, k=1, T, x_start=0.0, lowpass=true)
     T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive"))
     @named siso = SISO()
     @unpack u, y = siso
     sts = @variables x(t)=x_start
     pars = @parameters T=T k=k
     eqs = [
         D(x) ~ (k*u - x) / T
-        y ~ x
+        lowpass ? y ~ x : y ~ k*u - x
     ]
     extend(ODESystem(eqs, t, sts, pars; name=name), siso)
 end

test/Blocks/continuous.jl

@@ -43,7 +43,7 @@ end
 @testset "PT1" begin
     pt1_func(t, k, T) = k * (1 - exp(-t / T)) # Known solution to first-order system
 
-    k, T = 1.0, 0.1
+    k, T = 1.2, 0.1
     @named c = Constant(; k=1)
     @named pt1 = FirstOrder(; k=k, T=T)
     @named iosys = ODESystem(connect(c.output, pt1.input), t, systems=[pt1, c])
@@ -52,6 +52,15 @@ end
     sol = solve(prob, Rodas4())
     @test sol.retcode == :Success
     @test sol[pt1.output.u] ≈ pt1_func.(sol.t, k, T) atol=1e-3
+
+    # Test highpass feature
+    @named pt1 = FirstOrder(; k=k, T=T, lowpass=false)
+    @named iosys = ODESystem(connect(c.output, pt1.input), t, systems=[pt1, c])
+    sys = structural_simplify(iosys)
+    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
+    sol = solve(prob, Rodas4())
+    @test sol.retcode == :Success
+    @test sol[pt1.output.u] ≈ k .- pt1_func.(sol.t, k, T) atol=1e-3
 end
 
 @testset "PT2" begin