add JuliaC example (JuliaControl#13)

* add JuliaC example

* add c test program

* make wrappers `2ccallable` again

* juliac demo working!

Co-authored-by: Cody Tapscott <topolarity@tapscott.me>

* update juliac readme instrucitons

* add note about julia from julia

---------

Co-authored-by: Cody Tapscott <topolarity@tapscott.me>

Author: baggepinnen

Project.toml

@@ -13,7 +13,8 @@ julia = "1.7"
 AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AllocCheck", "Test", "ControlSystemsBase", "FixedPointNumbers"]
+test = ["AllocCheck", "Test", "ControlSystemsBase", "FixedPointNumbers", "JET"]

README.md

@@ -195,6 +195,36 @@ plot([res, res_fp], plotu=true, lab=["Float64" "" string(T) ""]); ylabel!("u + d
 
 The fixed-point controller behaves roughly the same in this case, but artifacts are clearly visible. If the number of bits used for the fractional part is decreased, the controller will start to misbehave.
 
+## Compilation using JuliaC
+The file `examples/juliac/juliac_pid.jl` contains a JuliaC-compatible interface that can be compiled into a C-callable shared library using JuliaC. To compile the file, run the following from the `examples/juliac` folder:
+```bash
+julia +nightly --project <PATH_TO_JULIA_REPO>/julia/contrib/juliac.jl --output-lib juliac_pid --experimental --trim=unsafe-warn --compile-ccallable juliac_pid.jl
+```
+where `<PATH_TO_JULIA_REPO>` should be replaced with the path to the Julia repository on your system. The command will generate a shared library `juliac_pid` that can be called from C. The file `examples/juliac/juliac_pid.h` contains the C-compatible interface to the shared library. The C program may be compiled with a command like
+```bash
+export LD_LIBRARY_PATH=<PATH_TO_JULIA_REPO>/julia/usr/lib:$LD_LIBRARY_PATH
+gcc -o pid_program test_juliac_pid.c -I <PATH_TO_JULIA_REPO>/julia/usr/include/julia -L<PATH_TO_JULIA_REPO>/julia/usr/lib -ljulia -ldl
+```
+and then run by
+```bash
+./pid_program
+```
+which should produce the output
+```
+DiscretePIDs/examples/juliac> ./pid_program 
+Loading juliac_pid.so
+Loaded juliac_pid.so
+Finding symbols
+Found all symbols!
+calculate_control! returned: 1.000000
+calculate_control! returned: 2.000000
+calculate_control! returned: 3.000000
+calculate_control! returned: 3.000000
+calculate_control! returned: 3.000000
+```
+At the time of writing, this requires a nightly version of julia
+
+
 ## See also
 - [TrajectoryLimiters.jl](https://github.com/baggepinnen/TrajectoryLimiters.jl) To generate dynamically feasible reference trajectories with bounded velocity and acceleration given an instantaneous reference $r(t)$ which may change abruptly.
 - [SymbolicControlSystems.jl](https://github.com/JuliaControl/SymbolicControlSystems.jl) For C-code generation of LTI systems.

examples/juliac/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DiscretePIDs = "c1363496-6848-4723-8758-079b737f6baf"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

examples/juliac/juliac_pid.jl

@@ -0,0 +1,50 @@
+module JuliacPID
+import DiscretePIDs
+import DiscretePIDs: DiscretePID
+import Base.@ccallable
+
+const T = Float64 # The numeric type used by the controller
+
+# Set the initial PID parameters here
+const pid = DiscretePIDs.DiscretePID(; K = T(1), Ti = 1, Td = false, Ts = 1)
+
+
+@ccallable function calculate_control!(r::T, y::T, uff::T)::T
+    DiscretePIDs.calculate_control!(pid, r, y, uff)::T
+end
+
+@ccallable function set_K!(K::T, r::T, y::T)::Cvoid
+    DiscretePIDs.set_K!(pid, K, r, y)
+    nothing
+end
+
+@ccallable function set_Ti!(Ti::T)::Cvoid
+    DiscretePIDs.set_Ti!(pid, Ti)
+    nothing
+end
+
+@ccallable function set_Td!(Td::T)::Cvoid
+    DiscretePIDs.set_Td!(pid, Td)
+    nothing
+end
+
+@ccallable function reset_state!()::Cvoid
+    DiscretePIDs.reset_state!(pid)
+    nothing
+end
+
+# @ccallable function main()::Cint
+#     println(Core.stdout, "I'm alive and well")
+#     u = calculate_control!(0.0, 0.0, 0.0)
+#     println(Core.stdout, u)
+
+#     Cint(0)
+# end
+
+
+end
+
+# compile using something like this, modified to suit your local paths
+# cd(@__DIR__)
+# run(`/home/fredrikb/repos/julia/julia --project --experimental /home/fredrikb/repos/julia/contrib/juliac.jl --output-lib juliac_pid --experimental --trim=unsafe-warn --compile-ccallable juliac_pid.jl`)
+# run(`ls -ltrh`)

examples/juliac/test_juliac_pid.c

@@ -0,0 +1,74 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <dlfcn.h>
+// #include <julia.h>
+
+// Path to julia binary folder
+#define JULIA_PATH "/home/fredrikb/repos/julia/usr/bin/" // NOTE: modify this path
+
+// Path to juliac compiled shared object file
+#define LIB_PATH "/home/fredrikb/.julia/dev/DiscretePIDs/examples/juliac/juliac_pid.so" // NOTE: modify this path
+
+
+
+// Define the types of the julia @ccallable functions
+typedef void (*jl_init_with_image_t)(const char *bindir, const char *sysimage);
+typedef double (*calculate_control_t)(double r, double y, double uff);
+typedef void (*set_K_t)(double K, double r, double y);
+typedef void (*set_Ti_t)(double Ti);
+typedef void (*set_Td_t)(double Td);
+typedef void (*reset_state_t)();
+
+
+int main() {
+
+    // Load the shared library
+    printf("Loading juliac_pid.so\n");
+    void *lib_handle = dlopen(LIB_PATH, RTLD_LAZY);
+    if (!lib_handle) {
+        fprintf(stderr, "Error: Unable to load library %s\n", dlerror());
+        exit(EXIT_FAILURE);
+    }
+    printf("Loaded juliac_pid.so\n");
+
+    // Locate the julia functions function
+    printf("Finding symbols\n");
+    jl_init_with_image_t jl_init_with_image = (jl_init_with_image_t)dlsym(lib_handle, "jl_init_with_image");
+
+    calculate_control_t calculate_control = (calculate_control_t)dlsym(lib_handle, "calculate_control!");
+    set_K_t set_K = (set_K_t)dlsym(lib_handle, "set_K!");
+    set_Ti_t set_Ti = (set_Ti_t)dlsym(lib_handle, "set_Ti!");
+    set_Td_t set_Td = (set_Td_t)dlsym(lib_handle, "set_Td!");
+    reset_state_t reset_state = (reset_state_t)dlsym(lib_handle, "reset_state!");
+
+
+    if (jl_init_with_image == NULL || calculate_control == NULL) {
+        char *error = dlerror();
+        fprintf(stderr, "Error: Unable to find symbol: %s\n", error);
+        exit(EXIT_FAILURE);
+    }
+    printf("Found all symbols!\n");
+
+    // Init julia
+    jl_init_with_image(JULIA_PATH, LIB_PATH);
+
+    // Trivial test program that computes a few control outputs and modifies K
+    double r = 1.0, y = 0.0, uff = 0.0;
+    double result = calculate_control(r, y, uff);
+    printf("calculate_control! returned: %f\n", result);
+    result = calculate_control(r, y, uff);
+    printf("calculate_control! returned: %f\n", result);
+    set_K(0.0, r, y);
+    for (int i = 0; i < 3; ++i) {
+        result = calculate_control(r, y, uff);
+        printf("calculate_control! returned: %f\n", result);
+    }
+
+    // jl_atexit_hook(0);
+    return 0;
+}
+
+
+// Compile this C program using a command like the one above, modified to suit your paths
+// export LD_LIBRARY_PATH=/home/fredrikb/repos/julia/usr/lib:$LD_LIBRARY_PATH
+// gcc -o pid_program test_juliac_pid.c -I /home/fredrikb/repos/julia/usr/include/julia -L/home/fredrikb/repos/julia/usr/lib -ljulia -ldl

examples/juliac/test_juliac_pid.jl

@@ -0,0 +1,34 @@
+# NOTE: it is currently not possible to call a julia-produced shared-library from julia.
+# To test the compiled shared library, see test_juliac_pid.c instead.
+cd(@__DIR__)
+
+const T = Float64
+@info("Loading juliac_pid.so")
+lib = Libc.Libdl.dlopen("/home/fredrikb/.julia/dev/DiscretePIDs/examples/juliac/juliac_pid.so")
+@info("Loaded juliac_pid.so, finding calculate_control!")
+const calc = Libc.Libdl.dlsym(lib, :calculate_control!)
+@info("Found calculate_control!")
+
+function pid(r::T, y::T, uff::T)
+    ccall(calc, T, (T, T, T), r, y, uff)
+end
+
+pid(0.0, 0.0, 0.0) # test
+
+using ControlSystemsBase, Plots
+Tf = 15   # Simulation time
+Ts = 0.01 # sample time
+
+P   = c2d(ss(tf(1, [1, 1])), Ts) # Process to be controlled, discretized using zero-order hold
+
+ctrl = function(x,t)
+    y = (P.C*x)[] # measurement
+    d = 1         # disturbance
+    r = 0         # reference
+    u = pid(T(r), T(y), T(0))
+    u + d # Plant input is control signal + disturbance
+end
+
+res = lsim(P, ctrl, Tf)
+
+plot(res, plotu=true); ylabel!("u + d", sp=2)

src/DiscretePIDs.jl

@@ -71,7 +71,7 @@ u = calculate_control!(pid, r, y, uff) # Equivalent to the above
 - `D`: Derivative part
 - `yold`: Last measurement signal
 
-See also [`calculate_control!`](@ref), [`set_K!`](@ref), [`set_Ti!`](@ref), [`set_Td!`](@ref)
+See also [`calculate_control!`](@ref), [`set_K!`](@ref), [`set_Ti!`](@ref), [`set_Td!`](@ref), [`reset_state!`](@ref).
 """
 function DiscretePID(;
     K::T  = 1f0,
@@ -115,6 +115,8 @@ end
     set_K!(pid::DiscretePID, K, r, y)
 
 Update `K` in the PID controller. This function takes the current reference and measurement as well in order to provide bumpless transfer. This is realized by updating the internal state `I`.
+
+Note: Due to the bumpless transfer, setting ``K = 0`` does not imply that the controller output will be 0 if the integral state is non zero. To reset the controller state, call `reset_state!(pid)`.
 """
 function set_K!(pid::DiscretePID, K, r, y)
     Kold = pid.K
@@ -124,6 +126,7 @@ function set_K!(pid::DiscretePID, K, r, y)
         pid.bi = K * pid.Ts / pid.Ti
         pid.I = pid.I + Kold*(pid.b*r - y) - K*(pid.b*r - y)
     end
+    nothing
 end
 
 """
@@ -139,6 +142,7 @@ function set_Ti!(pid::DiscretePID{T}, Ti) where T
     else
         pid.bi = zero(T)
     end
+    nothing
 end
 
 """
@@ -151,6 +155,7 @@ function set_Td!(pid::DiscretePID, Td)
     pid.Td = Td
     pid.ad = Td / (Td + pid.N * pid.Ts)
     pid.bd = pid.K * pid.N * pid.ad
+    nothing
 end
 
 
@@ -193,6 +198,7 @@ function reset_state!(pid::DiscretePID)
     pid.I = zero(pid.I)
     pid.D = zero(pid.D)
     pid.yold = zero(pid.yold)
+    nothing
 end
 
 end

test/runtests.jl

@@ -2,6 +2,7 @@ using DiscretePIDs
 using Test
 using ControlSystemsBase
 using AllocCheck
+using JET
 
 @testset "DiscretePIDs.jl" begin
 
@@ -95,6 +96,10 @@ reset_state!(pid)
 res3 = lsim(P, ctrl, Tf)
 @test res3.y == res2.y
 
+@test_opt pid(1.0, 1.0)
+@test_opt pid(1.0, 1.0, 1.0)
+# @report_call pid(1.0, 1.0)
+
 ## Test with FixedPointNumbers
 using FixedPointNumbers
 T = Fixed{Int16, 10} # 16-bit signed fixed-point with 10 bits for the fractional part