Winch model in Matlab: examples, tests and documentation

.gitignore

@@ -46,3 +46,4 @@ install_version_1.13.txt
 install_version_1.14.txt
 bin/kps-image-1.13.so
 *.bak
+examples/matlab/Manifest*.toml

docs/make.jl

@@ -15,6 +15,7 @@ makedocs(
         "Flight Path Controller" => "flight_path_controller.md",
         "Examples"  => "examples.md",
         "Projects and Turbulence" => "projects.md",
+        "Winch Model in Matlab" => "matlab.md",
         "Learning Control"  => "learning.md",
     ],
     checkdocs = :exports,

docs/src/examples.md

@@ -39,6 +39,11 @@ The primary example demonstrating the complete control stack.
 You can select one of the four, pre-defined projects, run the simulation, and on the left you have menu options to 
 print the statistics or show one of the many 2D plots, visualizing the results.
 
+## [autopilot_matlab.jl](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/autopilot_matlab.jl) — Autopilot with the Winch Model in Matlab
+
+Runs autopilot.jl with the winch model implemented in Matlab instead of the winch model of the
+selected project. Needs a Matlab installation, see [Winch Model in Matlab](@ref).
+
 ## [parking_4p.jl](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/parking_4p.jl) — Parking Controller
 
 Demonstrates how to bring a four-line kite to a stable parking position at

docs/src/matlab.md

@@ -0,0 +1,188 @@
+# Winch Model in Matlab
+
+The model of the winch can be implemented in Matlab instead of Julia. The rest of the simulation
+stays in Julia: the kite and tether model (KiteModels.jl), the winch controller
+(WinchControllers.jl), the flight path controller and the state machine `SystemStateControl`.
+The Matlab winch is driven by `SystemStateControl` through the Julia winch controller for the full
+power cycle: parking, power production, reel-in and depower.
+
+The Matlab function
+[`examples/matlab/async_machine_step.m`](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/matlab/async_machine_step.m)
+implements the same model as the `AsyncMachine` of
+[WinchModels.jl](https://github.com/aenarete/WinchModels.jl): an asynchronous machine with gearbox,
+drum, friction and brake. It is a starting point for your own Matlab winch models.
+
+## Requirements
+
+- Matlab R2020b or newer, with a license that allows starting the Matlab engine. Matlab Coder is
+  not needed. Octave is not supported.
+- On Linux, `csh` (the Matlab engine is started with it).
+- WinchModels.jl and KiteModels.jl with support for the `MatlabMachine`. Until they are released,
+  use the development versions (branch `matlab-machine`).
+
+## Setup
+
+The Matlab examples have their own environment, `examples/matlab`, because
+[MATLAB.jl](https://github.com/JuliaInterop/MATLAB.jl) cannot be installed without Matlab. In a
+Julia session started in the root folder of KiteControllers.jl:
+
+```julia
+using Pkg
+Pkg.activate("examples/matlab")
+ENV["MATLAB_ROOT"] = "/usr/local/MATLAB/R2021a"   # only if the wrong Matlab is found
+Pkg.instantiate()
+Pkg.build("MATLAB")
+```
+
+Packages in the system image cannot be replaced by other versions. To work on development versions
+of WinchModels.jl or KiteModels.jl, start Julia without the system image with
+`bin/run_julia --nosysimage`.
+
+## Usage
+
+### Autopilot with the Matlab winch
+
+Choose `autopilot_matlab` in the menu of the examples (`menu()`), or run:
+
+```julia
+include("examples/autopilot_matlab.jl")
+```
+
+This runs
+[autopilot.jl](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/autopilot.jl)
+with the Matlab winch model instead of the winch model of the selected project. If the active
+environment is not `examples/matlab`, for example when started from the menu, the script starts a
+new Julia process for it, without system image, and waits until the window is closed.
+
+### Selecting the Matlab winch in the settings
+
+In your own scripts, set the winch model in the settings file of the project:
+
+```yaml
+    winch_model: "MatlabMachine" # or AsyncMachine or TorqueControlledMachine
+```
+
+and load MATLAB.jl and add the folder of the Matlab function to the Matlab path before the first
+time step:
+
+```julia
+using MATLAB, KiteModels
+mat"addpath('examples/matlab')"
+```
+
+`KPS3` and `KPS4` then create a `MatlabMachine`. An unknown value of `winch_model` stops with an
+error.
+
+## How it works
+
+The acceleration of the winch is calculated in the residual of the DAE solver, which evaluates it
+many times per time step. A call of Matlab costs about a millisecond, so Matlab is called only once
+per time step, in `update_winch_state!` of WinchModels.jl:
+
+- Julia passes the reel-out speed and the tether force at the start of the time step, the set speed
+  of the winch controller and the state of the brake and of the rate limiter of the set speed.
+- `async_machine_step.m` updates the brake and the rate limited set speed and returns the
+  acceleration and its partial derivatives with respect to speed and force.
+- Within the time step, the residual uses the first order expansion
+  `acc0 + dacc_dv*(speed - speed0) + dacc_dF*(force - force0)`, evaluated in Julia.
+
+### Linearization within the time step
+
+Within one time step, the set speed, the state of the brake and the parameters are constant, so
+the acceleration of the winch is a function of the reel-out speed `v` and the tether force `F`
+only:
+
+```text
+acc = f(v, F)
+```
+
+At the start of each time step, Matlab evaluates this function once at the current point
+`(v0, F0)` and returns its value and its two slopes:
+
+```text
+acc0    = f(v0, F0)
+dacc_dv = ∂f/∂v (v0, F0)
+dacc_dF = ∂f/∂F (v0, F0)
+```
+
+During the time step, the DAE solver needs the acceleration many times, at slightly different
+values of `v` and `F`. Instead of calling Matlab, it uses the tangent plane through this point,
+the first two terms of the Taylor series (first order expansion):
+
+```text
+acc ≈ acc0 + dacc_dv * (v - v0) + dacc_dF * (F - F0)
+```
+
+The error of this approximation depends on how curved `f` is:
+
+- In the force there is no error: the force enters the torque balance linearly, so `dacc_dF` is
+  constant.
+- In the speed there is an error: the motor torque depends on the slip `δ = ω_sync - ω` as
+  `δ / (R² + L²δ²)`, and the Coulomb friction is a smooth sign function of the speed. The error grows
+  approximately with `½ ∂²f/∂v² (v - v0)²`, quadratically with the change of the speed within the
+  time step and with the curvature of the torque-slip curve.
+- With the brake engaged there is no error: the acceleration is `brake_acc * v`, which is linear.
+
+The error is therefore largest when the speed changes a lot within one time step in a strongly
+curved part of the torque-slip curve, for example at the transitions between reel-out and reel-in.
+Over full power cycles the effect is small, see [Accuracy](@ref). If more accuracy were needed,
+Matlab could also return the second derivative `∂²f/∂v²`, and Julia would add the quadratic term;
+this still needs only one call of Matlab per time step.
+
+### Calling Matlab
+
+The parameters (gear ratio, drum radius, inertia, friction, maximal acceleration, motor data) are
+sent to Matlab once, before the first time step. The Matlab engine is called by the package
+extension `WinchModelsMATLABExt`, which is loaded by `using MATLAB`: one evaluation per time step
+with the inputs as text, without the exception check of MATLAB.jl, and one read of the result
+vector. Errors in Matlab are reported as Julia errors.
+
+## Accuracy
+
+[`examples/matlab/compare_full_cycle.jl`](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/matlab/compare_full_cycle.jl)
+flies the power cycles of autopilot.jl, without viewer and turbulence, with the Julia
+`AsyncMachine` and with the `MatlabMachine`, for 460 s:
+
+| project | energy, Julia | energy, Matlab | difference |
+| --- | --- | --- | --- |
+| hydra10_951 | 654.2 Wh | 654.3 Wh | 0.014 % |
+| hydra20_426 | 752.3 Wh | 752.4 Wh | 0.018 % |
+| hydra20_600 | 1759.4 Wh | 1760.2 Wh | 0.044 % |
+| hydra20_920 | 1555.1 Wh | 1557.8 Wh | 0.177 % |
+
+Mean and maximal tether force differ by at most 0.13 %. The difference comes from the
+linearization within the time step.
+[`examples/matlab/test_async_machine_step.jl`](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/matlab/test_async_machine_step.jl)
+checks that the Matlab function itself agrees with the Julia model to 1e-12.
+
+## Speed
+
+One time step with the Matlab winch takes about 1.5 ms (median) and 2.2 ms (99th percentile),
+measured with Matlab R2021a on Linux, compared to less than 0.1 ms with the Julia winch model.
+The simulation is real-time capable up to a time-lapse factor of about 18 (2.8 ms per time step of
+50 ms); 24 is too fast. The 3D viewer of autopilot.jl needs additional time.
+[`examples/matlab/bench_engine.jl`](https://github.com/OpenSourceAWE/KiteControllers.jl/blob/main/examples/matlab/bench_engine.jl)
+measures the cost of different ways of calling Matlab.
+
+## Your own Matlab winch model
+
+Replace the body of `async_machine_step.m` and keep its interface:
+
+```matlab
+async_machine_step(p)    % store the parameter struct p
+y = async_machine_step(x)
+% x = [speed, force, set_speed, use_brake, brake, last_set_speed]
+% y = [acc0, dacc_dv, dacc_dF, brake, last_set_speed]
+```
+
+`set_speed` is `NaN` if there is no set speed. The model is called once per time step, so it may
+keep its own state in persistent variables, but the brake and the rate limited set speed must be
+returned, because KiteModels saves and restores them. The partial derivatives must describe the
+acceleration within the time step well; the accuracy of the simulation depends on them.
+
+## Limitations
+
+- Only the asynchronous machine is available in Matlab, not the `TorqueControlledMachine`.
+- Matlab holds the parameters of one `MatlabMachine` at a time; they are sent again when another
+  `MatlabMachine` is used. Change the parameters of a `MatlabMachine` before its first time step.
+- The Matlab engine adds about 1.5 ms per time step, see [Speed](@ref).

examples/autopilot.jl

@@ -72,6 +72,8 @@ app::KiteApp = KiteApp(deepcopy(load_settings(PROJECT)), 0, 0, true, nothing, no
 let use_turbulence = KiteControllers.get_use_turbulence(PROJECT)
     isnothing(use_turbulence) || (app.set.use_turbulence = use_turbulence)
 end
+# autopilot_matlab.jl defines WINCH_MODEL to replace the winch model of the project
+@isdefined(WINCH_MODEL) && (app.set.winch_model = WINCH_MODEL)
 app.max_time      = app.set.sim_time
 app.next_max_time = app.max_time
 
@@ -574,6 +576,7 @@ on(viewer_menu(:menu_project).i_selected) do _
                 let use_turbulence = KiteControllers.get_use_turbulence(PROJECT)
                     isnothing(use_turbulence) || (app.set.use_turbulence = use_turbulence)
                 end
+                @isdefined(WINCH_MODEL) && (app.set.winch_model = WINCH_MODEL)
                 app.max_time      = app.set.sim_time
                 app.next_max_time = app.max_time
                 app.initialized = false

examples/autopilot_matlab.jl

@@ -0,0 +1,22 @@
+# Run autopilot.jl with the winch model implemented in Matlab (MatlabMachine, see examples/matlab)
+# instead of the winch model of the selected project. Needs a Matlab installation.
+#
+# MATLAB.jl cannot be part of the examples environment, because it cannot be installed without
+# Matlab. This script therefore runs in the environment examples/matlab. If another environment is
+# active, it starts a new Julia process for it, without system image: the packages in the system
+# image cannot be replaced by the versions the MatlabMachine needs.
+const MATLAB_ENV = joinpath(@__DIR__, "matlab")
+
+if dirname(Base.active_project()) != MATLAB_ENV
+    julia_exe = joinpath(Sys.BINDIR, "julia")
+    println("Starting a new Julia process for the environment examples/matlab...")
+    run(`$julia_exe --project=$MATLAB_ENV -t auto --gcthreads=4,1 $(@__FILE__)`)
+else
+    # a macro like mat"" cannot be used here: the whole if block is expanded before `using MATLAB`
+    # runs; invokelatest, because MATLAB is loaded while this block runs
+    @eval using MATLAB
+    Base.invokelatest(MATLAB.eval_string, "addpath('$(replace(MATLAB_ENV, "'" => "''"))');")
+    WINCH_MODEL = "MatlabMachine"
+    include("autopilot.jl")
+end
+nothing

examples/matlab/Project.toml

@@ -0,0 +1,29 @@
+[deps]
+FFMPEG_jll = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
+KiteControllers = "b19fd437-5b3a-4cca-b199-c3ded0f20ded"
+KiteModels = "b94af626-7959-4878-9336-2adc27959007"
+KitePlots = "0b16aec0-02e8-48df-811e-5ff0d82eaaf3"
+KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533"
+KiteViewers = "2e593061-95e7-45e4-95f4-df0491f2e601"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MATLAB = "10e44e05-a98a-55b3-a45b-ba969058deb6"
+MakieControlPlots = "6d616b69-6563-4f6e-8472-6f6c706c6f74"
+NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Timers = "21f18d07-b854-4dab-86f0-c15a3821819a"
+WinchModels = "7dcfa46b-7979-4771-bbf4-0aee0da42e1f"
+
+[sources]
+KiteControllers = {path = "../.."}
+KitePlots = {path = "../../lib/KitePlots"}
+
+[compat]
+KiteModels = "0.13"
+KiteViewers = "0.7"
+LaTeXStrings = "1.4"
+MATLAB = "0.10"
+MakieControlPlots = "0.1.5"
+NativeFileDialog = "0.2.1"

examples/matlab/README.md

@@ -0,0 +1,64 @@
+# Winch model in Matlab
+
+`async_machine_step.m` is a Matlab implementation of the model of the winch with an asynchronous
+machine (`AsyncMachine` of WinchModels.jl). With `winch_model: "MatlabMachine"` in the settings,
+KiteModels uses it instead of the Julia model, called through the Matlab engine once per time
+step. See [PlanMatlab.md](../../PlanMatlab.md) for the design.
+
+## Requirements
+
+- Matlab R2020b or newer, and a license that allows starting the Matlab engine
+- Linux: `csh` (the Matlab engine is started with it)
+
+## Setup
+
+```julia
+using Pkg
+Pkg.activate("examples/matlab")
+ENV["MATLAB_ROOT"] = "/usr/local/MATLAB/R2021a"   # only if the wrong Matlab is found
+Pkg.instantiate()
+Pkg.build("MATLAB")
+```
+
+Until the new versions of WinchModels.jl and KiteModels.jl are released, use the development
+versions, e.g. `Pkg.develop(path="../WinchModels.jl")`.
+
+## Test
+
+```julia
+include("examples/matlab/test_async_machine_step.jl")
+```
+
+compares the Matlab implementation with the Julia `AsyncMachine`: acceleration, partial
+derivatives, brake and rate limit of the set speed, and measures the time per call.
+
+## Full power cycles
+
+```julia
+include("examples/matlab/compare_full_cycle.jl")
+```
+
+flies the power cycles of `examples/autopilot.jl` twice, without viewer and turbulence, with the
+Julia `AsyncMachine` and with the `MatlabMachine`, and compares the energy, the tether force, the
+reel-out speed and the time per step. It needs KiteModels and WinchModels with support for the
+`MatlabMachine`; until they are released, develop them in this environment.
+
+By default, it uses the project `system.yaml`. To compare other projects, define `PROJECTS`
+before including the script, e.g. `PROJECTS = ["hydra20_426.yml", "hydra20_600.yml"]`, or pass
+them as arguments:
+
+```bash
+julia --project=examples/matlab examples/matlab/compare_full_cycle.jl hydra20_426.yml hydra20_600.yml
+```
+
+## Benchmark of the Matlab engine
+
+```julia
+include("examples/matlab/bench_engine.jl")
+```
+
+measures the time of one call of Matlab from Julia for different ways of calling it (step 0 of
+PlanMatlab.md), using the dummy functions `bench_step.m` and `bench_step_vec.m`. Each
+`engEvalString` costs milliseconds, and `mxcall` and the exception check of `eval_string` need
+several of them; the `MatlabMachine` therefore uses one unchecked `eval_string` with the inputs as
+text and reads back one vector.