From a514aa6a45cd1dda9f231f0a844694fee548a670 Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 19:42:07 +0100
Subject: [PATCH 1/7] Add the Matlab implementation of the winch model

examples/matlab/async_machine_step.m implements the AsyncMachine of WinchModels.jl in
Matlab, used by the MatlabMachine (winch_model: "MatlabMachine"). It updates the brake
and the rate limited set speed and returns the acceleration and its partial derivatives.
test_async_machine_step.jl compares it with the Julia AsyncMachine and measures the time
per call; the folder has its own environment with MATLAB.jl.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .gitignore                                 |   1 +
 examples/matlab/Project.toml               |  12 +++
 examples/matlab/README.md                  |  33 +++++++
 examples/matlab/async_machine_step.m       | 105 +++++++++++++++++++++
 examples/matlab/test_async_machine_step.jl |  83 ++++++++++++++++
 5 files changed, 234 insertions(+)
 create mode 100644 examples/matlab/Project.toml
 create mode 100644 examples/matlab/README.md
 create mode 100644 examples/matlab/async_machine_step.m
 create mode 100644 examples/matlab/test_async_machine_step.jl

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/examples/matlab/Project.toml b/examples/matlab/Project.toml
new file mode 100644
index 00000000..bc860c00
--- /dev/null
+++ b/examples/matlab/Project.toml
@@ -0,0 +1,12 @@
+[deps]
+KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533"
+MATLAB = "10e44e05-a98a-55b3-a45b-ba969058deb6"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WinchModels = "7dcfa46b-7979-4771-bbf4-0aee0da42e1f"
+
+[sources]
+WinchModels = {path = "/home/ufechner/repos/WinchModels.jl"}
+
+[compat]
+MATLAB = "0.10"
diff --git a/examples/matlab/README.md b/examples/matlab/README.md
new file mode 100644
index 00000000..da6b94c9
--- /dev/null
+++ b/examples/matlab/README.md
@@ -0,0 +1,33 @@
+# 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.
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/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

From 709ebda154d66b7317cb13fb5e9c9839e89908f0 Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 19:45:11 +0100
Subject: [PATCH 2/7] examples/matlab: add the benchmark of the Matlab engine

bench_engine.jl measures the time of one call of Matlab from Julia for different ways of
calling it (step 0 of PlanMatlab.md), including the unchecked eval_string used by the
MatlabMachine.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 examples/matlab/README.md        | 12 ++++++
 examples/matlab/bench_engine.jl  | 65 ++++++++++++++++++++++++++++++++
 examples/matlab/bench_step.m     | 11 ++++++
 examples/matlab/bench_step_vec.m | 14 +++++++
 4 files changed, 102 insertions(+)
 create mode 100644 examples/matlab/bench_engine.jl
 create mode 100644 examples/matlab/bench_step.m
 create mode 100644 examples/matlab/bench_step_vec.m

diff --git a/examples/matlab/README.md b/examples/matlab/README.md
index da6b94c9..406429e0 100644
--- a/examples/matlab/README.md
+++ b/examples/matlab/README.md
@@ -31,3 +31,15 @@ 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.
+
+## 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/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

From d5fef87c955ff6dcbdb3db5395670138ebed65e9 Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 19:57:42 +0100
Subject: [PATCH 3/7] examples/matlab: compare full power cycles of the Julia
 and the Matlab winch model

compare_full_cycle.jl flies the power cycles of autopilot.jl, without viewer and
turbulence, with the AsyncMachine and with the MatlabMachine, and compares energy,
tether force, reel-out speed and time per step. The environment gets KiteControllers
(from this repository), KiteModels and Printf; the machine-specific path of WinchModels
is removed from its [sources].

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 examples/matlab/Project.toml          |   5 +-
 examples/matlab/README.md             |  11 +++
 examples/matlab/compare_full_cycle.jl | 106 ++++++++++++++++++++++++++
 3 files changed, 121 insertions(+), 1 deletion(-)
 create mode 100644 examples/matlab/compare_full_cycle.jl

diff --git a/examples/matlab/Project.toml b/examples/matlab/Project.toml
index bc860c00..b1befee4 100644
--- a/examples/matlab/Project.toml
+++ b/examples/matlab/Project.toml
@@ -1,12 +1,15 @@
 [deps]
+KiteControllers = "b19fd437-5b3a-4cca-b199-c3ded0f20ded"
+KiteModels = "b94af626-7959-4878-9336-2adc27959007"
 KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533"
 MATLAB = "10e44e05-a98a-55b3-a45b-ba969058deb6"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 WinchModels = "7dcfa46b-7979-4771-bbf4-0aee0da42e1f"
 
 [sources]
-WinchModels = {path = "/home/ufechner/repos/WinchModels.jl"}
+KiteControllers = {path = "../.."}
 
 [compat]
 MATLAB = "0.10"
diff --git a/examples/matlab/README.md b/examples/matlab/README.md
index 406429e0..2d5c3d18 100644
--- a/examples/matlab/README.md
+++ b/examples/matlab/README.md
@@ -32,6 +32,17 @@ 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.
+
 ## Benchmark of the Matlab engine
 
 ```julia
diff --git a/examples/matlab/compare_full_cycle.jl b/examples/matlab/compare_full_cycle.jl
new file mode 100644
index 00000000..1e952d1b
--- /dev/null
+++ b/examples/matlab/compare_full_cycle.jl
@@ -0,0 +1,106 @@
+# 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.
+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]
+# end of user parameter section #
+
+matlab_dir = @__DIR__
+mat"addpath($matlab_dir)"
+set_data_path(joinpath(@__DIR__, "..", "..", "data"))
+
+"""
+    fly(winch_model; 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; max_time=MAX_TIME)
+    set = deepcopy(load_settings("system.yaml"))
+    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
+
+println("Flying with the AsyncMachine (Julia)...")
+ref = fly("AsyncMachine")
+println("Flying with the MatlabMachine (Matlab)...")
+mat = fly("MatlabMachine")
+
+skip = Int(round(10 / (ref.time[2] - ref.time[1])))  # ignore the first 10 s (compilation)
+rel(a, b) = (a - b) / abs(b) * 100
+@printf("\n%-28s %14s %14s %10s\n", "", "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) : first_diff - 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])))
+nothing

From a4328fca52dcfbf625f220652d148e00f3185a8d Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 20:09:01 +0100
Subject: [PATCH 4/7] Add examples/autopilot_matlab.jl to fly autopilot.jl with
 the Matlab winch model

autopilot_matlab.jl runs autopilot.jl in the environment examples/matlab, with the
MatlabMachine instead of the winch model of the selected project; from another
environment, e.g. the menu, it starts a new Julia process without system image.
autopilot.jl uses WINCH_MODEL, if defined, as winch model. The environment
examples/matlab gets the packages autopilot.jl needs.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 examples/autopilot.jl        |  3 +++
 examples/autopilot_matlab.jl | 22 ++++++++++++++++++++++
 examples/matlab/Project.toml | 14 ++++++++++++++
 examples/menu.jl             |  1 +
 4 files changed, 40 insertions(+)
 create mode 100644 examples/autopilot_matlab.jl

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
index b1befee4..b42f69a9 100644
--- a/examples/matlab/Project.toml
+++ b/examples/matlab/Project.toml
@@ -1,15 +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/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\")",

From eb32b8bf1d50fb345948e536fd5f11b7be6ff2fe Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 20:16:21 +0100
Subject: [PATCH 5/7] compare_full_cycle.jl: compare the winch models for any
 project

The projects are passed as arguments, or as the global PROJECTS when the script is
included; system.yaml is the default.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 examples/matlab/README.md             |  8 +++
 examples/matlab/compare_full_cycle.jl | 83 ++++++++++++++++-----------
 2 files changed, 59 insertions(+), 32 deletions(-)

diff --git a/examples/matlab/README.md b/examples/matlab/README.md
index 2d5c3d18..b6a485be 100644
--- a/examples/matlab/README.md
+++ b/examples/matlab/README.md
@@ -43,6 +43,14 @@ Julia `AsyncMachine` and with the `MatlabMachine`, and compares the energy, 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
diff --git a/examples/matlab/compare_full_cycle.jl b/examples/matlab/compare_full_cycle.jl
index 1e952d1b..d9992ec7 100644
--- a/examples/matlab/compare_full_cycle.jl
+++ b/examples/matlab/compare_full_cycle.jl
@@ -1,6 +1,10 @@
 # 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__)
@@ -9,6 +13,9 @@ 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__
@@ -16,13 +23,13 @@ mat"addpath($matlab_dir)"
 set_data_path(joinpath(@__DIR__, "..", "..", "data"))
 
 """
-    fly(winch_model; max_time=MAX_TIME)
+    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; max_time=MAX_TIME)
-    set = deepcopy(load_settings("system.yaml"))
+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
@@ -72,35 +79,47 @@ function fly(winch_model; max_time=MAX_TIME)
     (; time, force, v_ro, length=length_, state, t_step, e_mech, cycles=ssc.fpp.fpca.cycle)
 end
 
-println("Flying with the AsyncMachine (Julia)...")
-ref = fly("AsyncMachine")
-println("Flying with the MatlabMachine (Matlab)...")
-mat = fly("MatlabMachine")
-
-skip = Int(round(10 / (ref.time[2] - ref.time[1])))  # ignore the first 10 s (compilation)
 rel(a, b) = (a - b) / abs(b) * 100
-@printf("\n%-28s %14s %14s %10s\n", "", "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))
+
+"""
+    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
-# 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) : first_diff - 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])))
 nothing

From e8c5eeb8473bb1bc2b790d3aa557e906b2d971a7 Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 20:19:20 +0100
Subject: [PATCH 6/7] Docs: add the page "Winch Model in Matlab"

Requirements, setup, usage, how the Matlab winch model is called, accuracy, speed, how
to write an own Matlab winch model and limitations. The examples page links to it
from a new section for autopilot_matlab.jl.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/make.jl         |   1 +
 docs/src/examples.md |   5 ++
 docs/src/matlab.md   | 143 +++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 149 insertions(+)
 create mode 100644 docs/src/matlab.md

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..e083203f
--- /dev/null
+++ b/docs/src/matlab.md
@@ -0,0 +1,143 @@
+# 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.
+
+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).

From 2220f95abc780404355998216a8fd7c9d283bf92 Mon Sep 17 00:00:00 2001
From: Uwe Fechner <fechner@aenarete.eu>
Date: Thu, 29 Oct 2026 20:25:47 +0100
Subject: [PATCH 7/7] Docs: explain the linearization of the Matlab winch model
 within the time step

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/matlab.md | 45 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/docs/src/matlab.md b/docs/src/matlab.md
index e083203f..61c423d9 100644
--- a/docs/src/matlab.md
+++ b/docs/src/matlab.md
@@ -86,6 +86,51 @@ per time step, in `update_winch_state!` of WinchModels.jl:
 - 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