diff --git a/.gitignore b/.gitignore index 07f59919..a167ec49 100644 --- a/.gitignore +++ b/.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 diff --git a/docs/make.jl b/docs/make.jl index f403d740..84556b4f 100644 --- a/docs/make.jl +++ b/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, diff --git a/docs/src/examples.md b/docs/src/examples.md index eec48ee0..c4711f09 100644 --- a/docs/src/examples.md +++ b/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 diff --git a/docs/src/matlab.md b/docs/src/matlab.md new file mode 100644 index 00000000..61c423d9 --- /dev/null +++ b/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). diff --git a/examples/autopilot.jl b/examples/autopilot.jl index edf83385..344c897e 100644 --- a/examples/autopilot.jl +++ b/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 diff --git a/examples/autopilot_matlab.jl b/examples/autopilot_matlab.jl new file mode 100644 index 00000000..08fecf3c --- /dev/null +++ b/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 diff --git a/examples/matlab/Project.toml b/examples/matlab/Project.toml new file mode 100644 index 00000000..b42f69a9 --- /dev/null +++ b/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" diff --git a/examples/matlab/README.md b/examples/matlab/README.md new file mode 100644 index 00000000..b6a485be --- /dev/null +++ b/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. diff --git a/examples/matlab/async_machine_step.m b/examples/matlab/async_machine_step.m new file mode 100644 index 00000000..320d0af6 --- /dev/null +++ b/examples/matlab/async_machine_step.m @@ -0,0 +1,105 @@ +function y = async_machine_step(x) +%ASYNC_MACHINE_STEP One time step of the model of a winch with an asynchronous machine. +% +% Matlab implementation of the AsyncMachine of WinchModels.jl, called from Julia (MatlabMachine) +% once per time step. +% +% async_machine_step(p) stores the parameter struct p. Required fields: +% gear_ratio, drum_radius [m], inertia_total [kgm²], max_force [N], f_coulomb [N], +% c_vf [Ns/m], max_acc [m/s²], sample_freq [Hz], u_nom [V], omega_sn [rad/s], +% omega_mn [rad/s], tau_n [Nm], tau_b [Nm], v_min [m/s], brake_acc [1/s²] +% +% y = async_machine_step(x) with +% x = [speed, force, set_speed, use_brake, brake, last_set_speed] +% y = [acc0, dacc_dv, dacc_dF, brake, last_set_speed] +% speed: reel-out speed at the start of the time step [m/s] +% force: tether force at the start of the time step [N] +% set_speed: set value of the reel-out speed [m/s], NaN if there is none +% use_brake, brake: 0 or 1 +% last_set_speed: rate limited set speed of the last time step [m/s] +% acc0: acceleration [m/s²] at speed and force; dacc_dv [1/s] and dacc_dF [m/(s²N)] are its +% partial derivatives. brake and last_set_speed are the new state of the winch. +% +% Like update_winch_state! of WinchModels.jl, the brake and the rate limited set speed are +% updated once, before the acceleration is calculated. + + persistent p + if isstruct(x) + fields = {'gear_ratio', 'drum_radius', 'inertia_total', 'max_force', 'f_coulomb', ... + 'c_vf', 'max_acc', 'sample_freq', 'u_nom', 'omega_sn', 'omega_mn', 'tau_n', ... + 'tau_b', 'v_min', 'brake_acc'}; + missing_fields = fields(~isfield(x, fields)); + if ~isempty(missing_fields) + error('WinchModels:params', 'missing parameters: %s', strjoin(missing_fields, ', ')); + end + p = x; + y = []; + return + end + if isempty(p) + error('WinchModels:params', 'call async_machine_step(p) with the parameters first'); + end + + speed = x(1); + force = x(2); + set_speed = x(3); + use_brake = x(4) ~= 0; + brake = x(5) ~= 0; + last_set_speed = x(6); + + % update_winch_state!: brake with hysteresis, then rate limit of the set speed + if ~isnan(set_speed) + if use_brake + if abs(set_speed) < 0.9 * p.v_min + brake = true; + elseif abs(set_speed) > 1.1 * p.v_min + brake = false; + end + end + if ~(use_brake && brake) + max_step = p.max_acc / p.sample_freq; + last_set_speed = min(max(set_speed, last_set_speed - max_step), last_set_speed + max_step); + end + end + + if use_brake && brake + % if the brake is active the acceleration is proportional to the speed + acc0 = p.brake_acc * speed; + dacc_dv = p.brake_acc; + dacc_dF = 0.0; + else + r = p.drum_radius; + n = p.gear_ratio; + % motor resistance R2 [Ohm] and inductance L [H] from the nameplate data + R2 = (p.u_nom^2 * (p.omega_sn - p.omega_mn)) / (2 * p.omega_sn^2) * ... + (1 / p.tau_n + sqrt(1 / p.tau_n^2 - 1 / p.tau_b^2)); + L = p.u_nom^2 / (2 * p.omega_sn * p.tau_b) / p.omega_sn; + omega = n / r * speed; + omega_sync = n / r * last_set_speed; + delta = omega_sync - omega; + if abs(omega_sync) <= p.omega_sn + c = p.omega_sn^2; + else + c = omega_sync^2; + end + den = R2^2 + L^2 * delta^2; + tau_motor = (p.u_nom^2 * R2 * delta) / (c * den); + dtau_motor_ddelta = p.u_nom^2 * R2 * (R2^2 - L^2 * delta^2) / (c * den^2); + % friction: Coulomb (smooth sign) and viscous + EPSILON = 6; + tau_c = p.f_coulomb * r / n; + c_v = p.c_vf * r^2 / n^2; + smooth_sign = omega / sqrt(omega^2 + EPSILON^2); + d_smooth_sign = EPSILON^2 / (omega^2 + EPSILON^2)^1.5; + tau_friction = tau_c * smooth_sign + c_v * omega; + dtau_friction_domega = tau_c * d_smooth_sign + c_v; + + k_force = r / n * 4000.0 / p.max_force; + omega_dot_m = (tau_motor + k_force * force - tau_friction) / p.inertia_total; + acc0 = r / n * omega_dot_m; + % d(delta)/d(speed) = d(omega)/d(speed) * -1 = -n/r + dacc_dv = r / n / p.inertia_total * (-dtau_motor_ddelta - dtau_friction_domega) * n / r; + dacc_dF = r / n / p.inertia_total * k_force; + end + y = [acc0, dacc_dv, dacc_dF, double(brake), last_set_speed]; +end diff --git a/examples/matlab/bench_engine.jl b/examples/matlab/bench_engine.jl new file mode 100644 index 00000000..a7cf58d5 --- /dev/null +++ b/examples/matlab/bench_engine.jl @@ -0,0 +1,65 @@ +# Measure the cost of calling Matlab from Julia through the Matlab engine (MATLAB.jl), as done in +# step 0 of PlanMatlab.md. Every engEvalString costs milliseconds, so the number of engine calls +# per Matlab call decides the speed. eval_string checks for Matlab exceptions by default, which +# needs another eval and two reads. +using Pkg +if dirname(Pkg.project().path) != @__DIR__ + Pkg.activate(@__DIR__) +end +using MATLAB, Statistics, Printf + +const N = 500 + +matlab_dir = @__DIR__ +mat"addpath($matlab_dir)" +println(mat"version") + +p = Dict("gear_ratio" => 6.2, "drum_radius" => 0.1615, "inertia_total" => 0.024, + "max_force" => 4000.0, "brake_acc" => -25.0) +x = [2.0, 1000.0, 2.1, 1.0, 0.0, 2.0] +mxcall(:bench_step_vec, 0, p) # store the parameters once + +function measure(f, n) + f() # warm-up + t = zeros(n) + for i in 1:n + t[i] = @elapsed f() + end + t .* 1e3 +end + +# one eval without the exception check of MATLAB.jl +function eval_unchecked(cmd) + session = get_default_msession() + MATLAB.disable_exception_check!(session) + try + eval_string(session, cmd) + finally + MATLAB.enable_exception_check!(session) + end +end + +literal = "[" * join(repr.(x), " ") * "]" +results = [ + "mxcall plus" => measure(() -> mxcall(:plus, 1, 1.0, 2.0), N), + "mxcall 7 in, 5 out" => measure(() -> mxcall(:bench_step, 5, 2.0, 1000.0, 2.1, true, false, 2.0, p), N), + "mxcall packed" => measure(() -> mxcall(:bench_step_vec, 1, x), N), + "eval_string only" => measure(() -> eval_string("1;"), N), + "eval_string literal + get" => measure(() -> begin + eval_string("y = bench_step_vec($literal);") + jarray(get_mvariable(:y)) + end, N), + # what the MatlabMachine does + "unchecked literal + get" => measure(() -> begin + eval_unchecked("y = bench_step_vec($literal);") + jarray(get_mvariable(:y)) + end, N), +] + +for (name, t) in results + @printf("%-26s median %7.3f ms p99 %7.3f ms max %7.3f ms\n", + name, median(t), quantile(t, 0.99), maximum(t)) +end +@printf("budget per time step at 4x real time (dt = 0.05 s): %.1f ms\n", 50 / 4) +eval_string("clear y;") +nothing diff --git a/examples/matlab/bench_step.m b/examples/matlab/bench_step.m new file mode 100644 index 00000000..dea62d78 --- /dev/null +++ b/examples/matlab/bench_step.m @@ -0,0 +1,11 @@ +function [acc0, dacc_dv, dacc_dF, brake, last_set_speed] = bench_step(speed, force, set_speed, use_brake, brake, last_set_speed, p) +% Dummy with seven inputs and five outputs, for timing the Matlab engine call only (bench_engine.jl). + if use_brake && abs(set_speed) < 0.1 + brake = true; + end + last_set_speed = set_speed; + omega = p.gear_ratio / p.drum_radius * speed; + acc0 = p.drum_radius / p.gear_ratio * (force * 4000.0 / p.max_force - 0.01 * omega) / p.inertia_total; + dacc_dv = -0.01 / p.inertia_total; + dacc_dF = p.drum_radius / p.gear_ratio * 4000.0 / p.max_force / p.inertia_total; +end diff --git a/examples/matlab/bench_step_vec.m b/examples/matlab/bench_step_vec.m new file mode 100644 index 00000000..7566c124 --- /dev/null +++ b/examples/matlab/bench_step_vec.m @@ -0,0 +1,14 @@ +function y = bench_step_vec(x) +% Packed variant of bench_step, with the interface of async_machine_step: +% x = [speed, force, set_speed, use_brake, brake, last_set_speed], +% y = [acc0, dacc_dv, dacc_dF, brake, last_set_speed]. bench_step_vec(p) stores the parameter +% struct p in a persistent variable. For timing the Matlab engine call only (bench_engine.jl). + persistent p + if isstruct(x) + p = x; + y = []; + return + end + [acc0, dacc_dv, dacc_dF, brake, last_set_speed] = bench_step(x(1), x(2), x(3), x(4) ~= 0, x(5) ~= 0, x(6), p); + y = [acc0, dacc_dv, dacc_dF, double(brake), last_set_speed]; +end diff --git a/examples/matlab/compare_full_cycle.jl b/examples/matlab/compare_full_cycle.jl new file mode 100644 index 00000000..d9992ec7 --- /dev/null +++ b/examples/matlab/compare_full_cycle.jl @@ -0,0 +1,125 @@ +# Fly the same power cycles twice, with the Julia AsyncMachine and with the MatlabMachine (the same +# model, implemented in Matlab and linearized within each time step), and compare the results. +# Without turbulence, so that both runs are deterministic. Needs a Matlab installation. +# +# Usage: julia --project=examples/matlab examples/matlab/compare_full_cycle.jl [project.yml ...] +# without projects, system.yaml is used. When included, the projects are taken from the global +# PROJECTS, if defined. +using Pkg +if dirname(Pkg.project().path) != @__DIR__ + Pkg.activate(@__DIR__) +end +using MATLAB, KiteControllers, KiteModels, KiteUtils, Statistics, Printf + +# the following values can be changed to match your interest +MAX_TIME = 460.0 # simulation time [s] +if !@isdefined(PROJECTS) + PROJECTS = isempty(ARGS) ? ["system.yaml"] : ARGS +end +# end of user parameter section # + +matlab_dir = @__DIR__ +mat"addpath($matlab_dir)" +set_data_path(joinpath(@__DIR__, "..", "..", "data")) + +""" + fly(winch_model, project; max_time=MAX_TIME) + +Simulate with the autopilot, as examples/autopilot.jl does, but without viewer and turbulence. +Returns a named tuple with the time series and the time per step. +""" +function fly(winch_model, project; max_time=MAX_TIME) + set = deepcopy(load_settings(project)) + set.winch_model = winch_model + set.use_turbulence = 0.0 + dt = 1/set.sample_freq + kcu = KCU(set) + kps4 = KPS4(kcu) + kps4.wm.v_min = 0.1 + wcs = WCSettings(true; dt) + fcs = FPCSettings(true; dt) + fpps = FPPSettings(true) + u_d0 = 0.01 * set.depower_offset + u_d = 0.01 * set.depowers[1] + ssc = SystemStateControl(wcs, fcs, fpps; u_d0, u_d, v_wind=set.v_wind) + on_parking(ssc) + integrator = KiteModels.init!(kps4; trim=true) + isnothing(integrator) && error("Could not find the initial steady state.") + sys_state = SysState(kps4) + on_new_systate(ssc, sys_state) + + steps = Int(round(max_time/dt)) + time = zeros(steps) + force = zeros(steps) + v_ro = zeros(steps) + length_ = zeros(steps) + state = zeros(Int, steps) + t_step = zeros(steps) + e_mech = 0.0 + for i in 1:steps + if i > 100 + dp = max(KiteControllers.get_depower(ssc), 0.22) + heading = calc_heading(kps4; neg_azimuth=true, one_point=false) + azimuth = -calc_azimuth(kps4) + steering = -calc_steering(ssc; heading, azimuth) + set_depower_steering(kps4.kcu, dp, steering) + end + i == 200 && on_autopilot(ssc) + v_set = calc_v_set(ssc) + t_step[i] = @elapsed KiteModels.next_step!(kps4, integrator; set_speed=v_set, dt) + KiteModels.update_sys_state!(sys_state, kps4) + on_new_systate(ssc, sys_state) + e_mech += sys_state.winch_force[1] * sys_state.v_reelout[1] / 3600 * dt + time[i] = i * dt + force[i] = sys_state.winch_force[1] + v_ro[i] = sys_state.v_reelout[1] + length_[i] = sys_state.l_tether[1] + state[i] = Int(ssc.fpp._state) + end + (; time, force, v_ro, length=length_, state, t_step, e_mech, cycles=ssc.fpp.fpca.cycle) +end + +rel(a, b) = (a - b) / abs(b) * 100 + +""" + compare(project) + +Fly with both winch models and print the comparison. +""" +function compare(project) + println("\n$project: flying with the AsyncMachine (Julia)...") + ref = fly("AsyncMachine", project) + println("$project: flying with the MatlabMachine (Matlab)...") + mat = fly("MatlabMachine", project) + + skip = Int(round(10 / (ref.time[2] - ref.time[1]))) # ignore the first 10 s (compilation) + @printf("\n%-28s %14s %14s %10s\n", project, "AsyncMachine", "MatlabMachine", "diff") + @printf("%-28s %14.1f %14.1f %9.3f %%\n", "energy e_mech [Wh]", ref.e_mech, mat.e_mech, + rel(mat.e_mech, ref.e_mech)) + @printf("%-28s %14d %14d\n", "cycles", ref.cycles, mat.cycles) + @printf("%-28s %14.1f %14.1f %9.3f %%\n", "mean force [N]", mean(ref.force), mean(mat.force), + rel(mean(mat.force), mean(ref.force))) + @printf("%-28s %14.1f %14.1f %9.3f %%\n", "max force [N]", maximum(ref.force), maximum(mat.force), + rel(maximum(mat.force), maximum(ref.force))) + @printf("%-28s %14.3f %14.3f\n", "min reel-out speed [m/s]", minimum(ref.v_ro), minimum(mat.v_ro)) + @printf("%-28s %14.3f %14.3f\n", "max reel-out speed [m/s]", maximum(ref.v_ro), maximum(mat.v_ro)) + for (name, t) in (("AsyncMachine", ref.t_step), ("MatlabMachine", mat.t_step)) + t = t[skip:end] .* 1e3 + @printf("time per step, %-14s median %.3f ms p99 %.3f ms max %.3f ms\n", + name, median(t), quantile(t, 0.99), maximum(t)) + end + # the difference of the time series, until the first change of the state of the flight path + # planner that happens at a different time step + first_diff = findfirst(ref.state .!= mat.state) + n = isnothing(first_diff) ? length(ref.time) : max(first_diff - 1, 1) + @printf("identical sequence of the flight path planner states until t = %.2f s\n", ref.time[n]) + @printf("until then: max |Δforce| = %.1f N, max |Δv_ro| = %.4f m/s, max |Δlength| = %.4f m\n", + maximum(abs.(mat.force[1:n] .- ref.force[1:n])), + maximum(abs.(mat.v_ro[1:n] .- ref.v_ro[1:n])), + maximum(abs.(mat.length[1:n] .- ref.length[1:n]))) +end + +for project in PROJECTS + compare(project) +end +nothing diff --git a/examples/matlab/test_async_machine_step.jl b/examples/matlab/test_async_machine_step.jl new file mode 100644 index 00000000..1bd0833a --- /dev/null +++ b/examples/matlab/test_async_machine_step.jl @@ -0,0 +1,83 @@ +# Compare the Matlab implementation of the winch model (async_machine_step.m, used by the +# MatlabMachine of WinchModels.jl) with the Julia AsyncMachine. Needs a Matlab installation. +using Pkg +if dirname(Pkg.project().path) != @__DIR__ + Pkg.activate(@__DIR__) +end +using MATLAB, WinchModels, KiteUtils, Statistics, Test + +matlab_dir = @__DIR__ +mat"addpath($matlab_dir)" +set_data_path(joinpath(@__DIR__, "..", "..", "data")) +set = load_settings("system.yaml") + +ext = Base.get_extension(WinchModels, :WinchModelsMATLABExt) + +# acceleration of the Julia model for the state of am, without changing it +acc_julia(am, speed, force, use_brake) = + calc_acceleration(am, speed, force; use_brake, update_state=false) + +@testset "MatlabMachine" begin + @testset "numbers arrive in Matlab unchanged" begin + x = [0.1, 1/3, -2.5e-17, 1.0e300, 5.0e-324, prevfloat(1.0), NaN, Inf, -Inf] + eval_string("wm_t = [" * join(ext.mat_num.(x), " ") * "];") + y = vec(jarray(get_mvariable(:wm_t))) + @test isequal(y, x) + eval_string("clear wm_t;") + end + + @testset "Matlab model equals the Julia AsyncMachine" begin + mm = MatlabMachine(set) + n_cases = 0 + for use_brake in (false, true), brake in (false, true), + last_set_speed in (-4.0, -0.1, 0.0, 0.15, 2.0, 7.5), + set_speed in (nothing, -6.0, -0.15, 0.0, 0.1, 0.3, 2.2, 9.0), + speed in (-5.0, -0.5, 0.0, 0.05, 1.9, 6.0), + force in (0.0, 300.0, 4000.0) + + am = AsyncMachine(set) + am.brake, am.last_set_speed = brake, last_set_speed + mm.brake, mm.last_set_speed = brake, last_set_speed + update_winch_state!(am, set_speed; use_brake) + update_winch_state!(mm, set_speed; use_brake, speed, force) + @test mm.brake == am.brake + @test mm.last_set_speed == am.last_set_speed + acc = acc_julia(am, speed, force, use_brake) + @test mm.acc0 ≈ acc rtol=1e-12 atol=1e-12 + # partial derivatives against central differences of the Julia model + h = 1e-6 + dacc_dv = (acc_julia(am, speed + h, force, use_brake) - + acc_julia(am, speed - h, force, use_brake)) / 2h + dacc_dF = (acc_julia(am, speed, force + 1e-3, use_brake) - + acc_julia(am, speed, force - 1e-3, use_brake)) / 2e-3 + @test mm.dacc_dv ≈ dacc_dv rtol=1e-6 atol=1e-6 + @test mm.dacc_dF ≈ dacc_dF rtol=1e-6 atol=1e-9 + # the first order expansion is exact at the point of the linearization + @test calc_acceleration(mm, speed, force; use_brake, update_state=false) ≈ mm.acc0 + n_cases += 1 + end + println("compared $n_cases cases") + end + + @testset "errors of Matlab are reported" begin + mm = MatlabMachine(set) + update_winch_state!(mm, 2.0; speed=2.0, force=1000.0) + eval_string("clear async_machine_step;") # forget the parameters + @test_throws ErrorException update_winch_state!(mm, 2.0; speed=2.0, force=1000.0) + mm.initialized = false # sends the parameters again + update_winch_state!(mm, 2.0; speed=2.0, force=1000.0) + @test mm.initialized + end + + @testset "time per call" begin + mm = MatlabMachine(set) + mm.brake = false + update_winch_state!(mm, 2.0; use_brake=true, speed=2.0, force=1000.0) + t = [(@elapsed update_winch_state!(mm, 2.0; use_brake=true, speed=2.0, force=1000.0)) + for _ in 1:500] .* 1e3 + println("update_winch_state!(::MatlabMachine): median $(round(median(t), digits=3)) ms, " * + "max $(round(maximum(t), digits=3)) ms") + @test median(t) < 50.0 # less than one time step of 50 ms + end +end +nothing diff --git a/examples/menu.jl b/examples/menu.jl index 7073ad8d..de2bb54f 100644 --- a/examples/menu.jl +++ b/examples/menu.jl @@ -11,6 +11,7 @@ include("select_project.jl") options = ["select_project()", "set_turbulence = set_default_turbulence()", "autopilot_4p = include(\"autopilot.jl\")", + "autopilot_matlab = include(\"autopilot_matlab.jl\")", "batch_pilot = include(\"batch_pilot.jl\")", "batch_plot = include(\"batch_plot.jl\")", "create_wind_fields = include(\"create_wind_fields.jl\")",