Add trim option, deprecate delta and stiffness_factor, fix steady-state wind and reel-out velocities

.JETLSConfig.toml.default

@@ -22,4 +22,17 @@ severity = "off"
 pattern = "`integrator` is not defined"
 match_by = "message"
 match_type = "literal"
-severity = "off"
+severity = "off"
+
+[[diagnostic.patterns]]
+pattern = "`KiteModels.KPS4` has abstract field `wm::WinchModels.AbstractWinchModel`"
+match_by = "message"
+match_type = "literal"
+severity = "off"
+
+[[diagnostic.patterns]]
+pattern = "`solve` is captured and boxed"
+match_by = "message"
+match_type = "literal"
+severity = "off"
+path = "src/KPS4.jl"

.markdownlint.json

@@ -4,7 +4,7 @@
     "start_indented": false
   },
   "MD013": {
-    "line_length": 1000,
+    "line_length": 2000,
     "tables": false
   },
   "MD024": {

CHANGELOG.md

@@ -2,6 +2,21 @@
 SPDX-FileCopyrightText: 2025 Uwe Fechner, Bart van de Lint
 SPDX-License-Identifier: MIT
 -->
+### KiteModels v0.12.1 - unreleased
+#### Added
+- `find_steady_state!(s; trim=true)` (KPS3 and KPS4): instead of prescribing `set.elevation`, the vertical force balance of the kite is solved too, so the kite is placed at the elevation where its lift matches gravity (its natural elevation); `set.elevation` is only the start value. For KPS4 all forces including gravity are balanced in this case (with `trim=false`, the accelerations in `yd0` are `-g`, so gravity is not balanced), and the returned accelerations are zero. A second, fully balanced trim exists at low elevations with the kite in deep stall. Solutions with an angle of attack beyond the angle of maximum lift (from `alpha_cl`/`cl_list`) are therefore rejected, and if no valid trim is found from `set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 10°, with a cheaper solve (no perturbation, no last-resort continuation). If no trim with attached flow is found, an error is thrown instead of returning a deep stall or non-converged state; with `system.yaml` and the Maasvlakte wind profile at 4 m/s this takes 7 s (KPS3) and 11 s (KPS4). The default `trim=false` keeps the previous behaviour.
+- `init!(s; trim=true)` passes `trim` to `find_steady_state!`. The initial state is then an equilibrium of the simulation, so the elevation stays constant. If no trim with attached flow is found, `init!` returns `nothing`.
+- `find_steady_state!` with `trim=false` (KPS3 and KPS4) logs a warning if the kite is in deep stall (angle of attack beyond the angle of maximum lift). This happens if the prescribed elevation is far below the natural elevation (with `system.yaml`: below about 65° for KPS3 and 45° for KPS4), because then the kite can only be balanced in deep stall; before, such a state was returned silently.
+#### Changed
+- The keyword argument `delta` of `find_steady_state!` and `init!` is deprecated and ignored. `delta` perturbed the initial state (all particles shifted sideways by `delta` [m], with the velocity `delta` [m/s]), and the returned state was perturbed too, so it was not an equilibrium: the kite moved by up to 0.4° after `init!` even with `trim=true`. Now the steady state is solved without perturbation, with a perturbation of 0.001 only as an internal fallback, and it is always returned without perturbation. The initial states change slightly, also for `steady_state=false`. The examples, tests and `precompile.jl` no longer pass `delta`.
+- The keyword argument `stiffness_factor` of `find_steady_state!` and `init!` is deprecated and ignored. Before, the KPS4 steady state was solved for the reduced tether stiffness `stiffness_factor`, which `next_step!` then raised by 0.01 per time step to 1.0, so the initial state was not an equilibrium of the simulation (with a stretched tether and a transient during the first seconds), and the result depended on the chosen value. Now the steady state is always solved for the full tether stiffness, and `s.stiffness_factor` is 1.0 afterwards. A reduced stiffness is only used internally as start value: if the direct solve fails, the stiffness schedules 0.1 → 1.0, 0.1 → 0.3 → 1.0 and 0.3 → 1.0 are tried, and then the previous continuation. In a benchmark with `system.yaml` (wind speeds of 6 to 25 m/s, elevations of 40° to 85°, tether lengths of 50 to 1000 m, with and without trim), all 160 cases were solved, in 0.14 s on average. For KPS3, the first solve uses a fixed stiffness factor of 0.035; the result was already solved for the full stiffness. The initial states change. The examples, tests and `precompile.jl` no longer pass `stiffness_factor`.
+- With `trim=true`, `find_steady_state!` (KPS4) starts the search at 80° (or at `set.elevation`, if it is higher) instead of at `set.elevation`, because the solver converges from above the natural elevation, but often not from far below it. The restarts at other elevations are unchanged.
+- The cases of `test/test-steady-state-robustness.jl` at an elevation of 50° are moved to 60° and 65°, because at 50° the steady state is in deep stall.
+#### Fixed
+- `find_steady_state!` (KPS4) now solves for the wind at the height of the kite. Before, the wind was set once from the height of the kite before the solver started, which is 0 m after `clear!`, so the kite was balanced for the ground wind (9.51 m/s instead of 12.6 m/s with `system.yaml`); afterwards the wind at the real height was set, which left force residuals of about 130 N on the kite points. The initial state changes: for example, with `system.yaml` the winch force increases from 308 N to 456 N.
+- The initial velocities of the particles (KPS4, in `find_steady_state!` and `init!`) caused by reeling out were multiplied by the length of the tether segments (about 25 m with `system.yaml`), and the direction of the KCU was taken from the bridle. With `v_reel_out = 1 m/s` the KCU moved at 20.9 m/s and the kite at 5.0 m/s instead of 1.0 m/s, so the steady state was found for a much lower apparent wind speed. With `trim=true` the kite was therefore placed far below its natural elevation (54.8° instead of 68.2° in `examples/test_init_4p.jl`), and after `init!` it climbed for about 10 s, with a peak of the tether force of more than twice the steady value. Now the velocity of each tether particle points along its tether segment and grows linearly from zero at the winch to `v_reel_out` at the KCU; the kite particles move with the KCU. The initial states change for all simulations with `v_reel_out ≠ 0`.
+- The steady state of KPS3 (`find_steady_state!` and `init!`) ignored `v_reel_out`: all particles were at rest and the reel-out speed in the initial state vector was zero, so the initial elevation did not depend on `v_reel_out`, and in the simulation the winch first had to accelerate to the set speed. Now, as for KPS4, the velocity of each particle points along its tether segment and grows linearly from zero at the winch to `v_reel_out` at the kite, and the initial reel-out speed is `v_reel_out`. With `trim=true`, `system.yaml` and a wind speed of 14.5 m/s, the natural elevation decreases from 73.5° to 70.4° for `v_reel_out = 1 m/s` and to 60.9° for 4 m/s, and it stays constant during the simulation. The initial states change for all simulations with `v_reel_out ≠ 0`.
+
 ### KiteModels v0.12.0 2026-10-11
 #### Changed
 - BREAKING: the quaternion written to `SysState` is `KA` (aft-right-up, against ENU),

Manifest-v1.12.toml.default

@@ -2,7 +2,7 @@
 
 julia_version = "1.12.7"
 manifest_format = "2.0"
-project_hash = "02a519ef26de8d09c0837d7cff2d48fa2d64c8e7"
+project_hash = "9e386311c7e10678d81dbab1d2f61e3246ac3958"
 
 [[deps.ADTypes]]
 deps = ["PrecompileTools"]
@@ -1366,23 +1366,21 @@ version = "0.4.2"
 
 [[deps.KiteUtils]]
 deps = ["Arrow", "CSV", "DocStringExtensions", "LinearAlgebra", "OrderedCollections", "Parameters", "Parsers", "Pkg", "PrecompileTools", "RecursiveArrayTools", "ReferenceFrameRotations", "Rotations", "StaticArrays", "StructArrays", "StructTypes", "YAML"]
-git-tree-sha1 = "85b4f730bf098b5afc796daa11dd8fde66a1db2c"
+git-tree-sha1 = "feec003e7b99dac12af25693057647d38f48c7f0"
 uuid = "90980105-b163-44e5-ba9f-8b1c83bb0533"
-version = "0.13.1"
+version = "0.13.2"
 
 [[deps.KiteViewers]]
 deps = ["FileIO", "GLMakie", "GeometryBasics", "Joysticks", "KiteUtils", "LinearAlgebra", "Parameters", "Pkg", "PrecompileTools", "Printf", "Reexport", "Rotations", "StaticArrays", "StructArrays", "Timers"]
-git-tree-sha1 = "503f9eaaaae53811443a59a88f5bfd472bd6a66c"
-repo-rev = "agent/54-update-kiteutils"
-repo-url = "https://github.com/OpenSourceAWE/KiteViewers.jl"
+git-tree-sha1 = "017c164c65b5f6ca1efc97e7dfa3c5d9ef4c78e2"
 uuid = "2e593061-95e7-45e4-95f4-df0491f2e601"
-version = "0.6.0"
+version = "0.7.0"
 
 [[deps.Krylov]]
 deps = ["LinearAlgebra", "Printf", "SparseArrays"]
-git-tree-sha1 = "909890a0d1486c392128cf4bd6e038ee0c6a7af4"
+git-tree-sha1 = "e0574e7e2cd5c40d9a6ae444511fbadd51afbf30"
 uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
-version = "0.10.10"
+version = "0.10.11"
 
 [[deps.LAME_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1677,9 +1675,9 @@ version = "0.24.15"
 
 [[deps.MakieControlPlots]]
 deps = ["CairoMakie", "GLMakie", "GridLayoutBase", "JLD2", "LaTeXStrings", "Makie", "Pkg", "PrecompileTools", "Printf", "StaticArraysCore"]
-git-tree-sha1 = "aeb72de0f9bad76c4080d79b954864e6b090be2e"
+git-tree-sha1 = "1d45e9c07bb53aad7c9f666c221a07b3a9199e9b"
 uuid = "6d616b69-6563-4f6e-8472-6f6c706c6f74"
-version = "0.1.20"
+version = "0.1.21"
 weakdeps = ["ControlSystemsBase"]
 
     [deps.MakieControlPlots.extensions]
@@ -2085,9 +2083,9 @@ version = "1.58.2+0"
 
 [[deps.Parameters]]
 deps = ["OrderedCollections", "UnPack"]
-git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
+git-tree-sha1 = "0ed04c372da78ff5b98e35e3bd4f4ca9931299cc"
 uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-version = "0.12.3"
+version = "0.13.1"
 
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
@@ -3136,15 +3134,15 @@ version = "1.2.0+0"
 
 [[deps.Xorg_xkbcomp_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxkbfile_jll"]
-git-tree-sha1 = "801a858fc9fb90c11ffddee1801bb06a738bda9b"
+git-tree-sha1 = "9a5e6e8f1a7c40941465c04840145e527f0d6b4b"
 uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
-version = "1.4.7+0"
+version = "1.5.0+0"
 
 [[deps.Xorg_xkeyboard_config_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xkbcomp_jll"]
-git-tree-sha1 = "2e59214e017a55cb87474a00fa76035c82ac0e17"
+git-tree-sha1 = "df77b22eb77051b3d0ef30aab979684f7482c66f"
 uuid = "33bec58e-1273-512f-9401-5d533626f822"
-version = "2.47.0+2"
+version = "2.48.0+0"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -3253,9 +3251,9 @@ version = "1.64.0+1"
 
 [[deps.oneTBB_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
-git-tree-sha1 = "da8c1f6eee04831f14edcfa5dae611d309807e57"
+git-tree-sha1 = "ed2d8d4ad96c03d3abe3382cdbbdea9489b63636"
 uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e"
-version = "2022.3.0+0"
+version = "2023.1.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]

docs/src/examples_4p.md

@@ -38,8 +38,9 @@ Then we call the function `find_steady_state` which uses a non-linear solver to
 find_steady_state!(kps, prn=true)
 ```
 
-Finding the steady state of the 4-point model is difficult and it only works when we artificially reduce the stiffness by a factor
-of 0.035. In the function [`init!`](@ref) this factor is slowly increased to 1.0.
+The steady state is solved for the full tether stiffness. If the solver does not converge directly, a solution for a
+reduced stiffness is found first and used as start value. With `trim=true`, the kite is placed at its natural elevation
+instead of at the elevation from the settings.
 
 To plot the result in 2D we extract the vectors of the x and z coordinates of the tether particles with a for loop:
 
@@ -97,9 +98,6 @@ julia> calc_height(kps)
 152.18828499538841
 ```
 
-Because of the the stiffness_factor of 0.035 we have a longer tether-length then when using
-the 1 point kite model.
-
 Print the force at the winch (groundstation, in Newton) and at each tether segment:
 
 ```julia

examples/Project.toml

@@ -1,7 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
-MakieControlPlots = "6d616b69-6563-4f6e-8472-6f6c706c6f74"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DiscretePIDs = "c1363496-6848-4723-8758-079b737f6baf"
@@ -11,6 +10,7 @@ KiteModels = "b94af626-7959-4878-9336-2adc27959007"
 KitePodModels = "9de5dc81-f971-414a-927b-652b2f41c539"
 KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+MakieControlPlots = "6d616b69-6563-4f6e-8472-6f6c706c6f74"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -23,12 +23,12 @@ KiteModels = {path = ".."}
 [compat]
 BenchmarkTools = "1.6"
 Colors = "0"
-MakieControlPlots = "0.1.4"
 ControlSystemsBase = "1"
 DSP = "0.7.10"
 Documenter = "1.10.1"
 JLD2 = "0.5"
 LaTeXStrings = "1.4.0"
+MakieControlPlots = "0.1.4"
 SciMLBase = "2"
 StatsBase = "0.34"
 julia = "1.12, 1.13"

examples/bench.jl

@@ -39,7 +39,7 @@ function simulate(integrator, steps)
     iter / steps
 end
 
-integrator = KiteModels.init!(kps3; delta = 0.001, stiffness_factor = 0.1, prn = STATISTIC)
+integrator = KiteModels.init!(kps3; prn = STATISTIC)
 
 println("\nStarting simulation...")
 simulate(integrator, 100)

examples/bench_4p.jl

@@ -58,7 +58,7 @@ function simulate(integrator, steps, offset = 0)
     iter / steps
 end
 
-integrator = init!(kps4, delta = 0.001, stiffness_factor = 0.1, prn = STATISTIC)
+integrator = init!(kps4, prn = STATISTIC)
 
 println("\nStarting simulation...")
 simulate(integrator, 100, 100)

examples/calc_spectrum.jl

@@ -103,7 +103,7 @@ function sim_and_plot(set; depower=DEPOWER, f_ex)
     set.elevation = 67.0
     kcu = KCU(set)
     kps4::KPS4 = KPS4(kcu)
-    integrator = KiteModels.init!(kps4; delta=0.001, stiffness_factor=0.04, prn=STATISTIC)
+    integrator = KiteModels.init!(kps4; prn=STATISTIC)
     set_depower_steering(kps4.kcu, depower, 0.0)
     # scale excitation amplitude to the steady-state winch torque for any kite size
     sin_amplitude = 6.8 * norm(kps4.forces[1]) * set.drum_radius / set.gear_ratio

examples/compare_kps3_kps4.jl

@@ -47,12 +47,12 @@ function simulate(s, integrator, steps, plot = false; fig = "")
             s.pos,
             0.0;
             zoom = ZOOM,
-            xlim = (35, 55),
+            xlim = (40, 55),
             xy = (52.0, z_text),
             front = FRONT_VIEW,
             fig,
         )
-        sleep(0.025)
+        sleep(0.05)
     end
     for i in 1:steps
         if PRINT
@@ -72,19 +72,19 @@ function simulate(s, integrator, steps, plot = false; fig = "")
                     s.pos,
                     reltime;
                     zoom = ZOOM,
-                    xlim = (35, 75),
+                    xlim = (40, 55),
                     xy = (52.0, z_text),
                     front = FRONT_VIEW,
                     fig,
                 )
-                sleep(0.025)
+                sleep(0.05)
             end
         end
     end
     iter / steps
 end
 
-integrator = init!(kps3, delta = 0.001, stiffness_factor = 0.1, prn = STATISTIC)
+integrator = init!(kps3, prn = STATISTIC)
 av_steps = simulate(kps3, integrator, STEPS, true; fig = "kps3")
 reactivate_host_app()
 
@@ -95,7 +95,7 @@ println("winch_force [N]: $(round(winch_force(kps3), digits=2))")
 println("Average number of callbacks per time step: $(round(av_steps, digits=2))")
 
 kps4.set.alpha_zero = ALPHA_ZERO
-integrator = init!(kps4; delta = 0.003, stiffness_factor = 0.1, prn = STATISTIC)
+integrator = init!(kps4; prn = STATISTIC)
 av_steps = simulate(kps4, integrator, STEPS, true; fig = "kps4")
 reactivate_host_app()
 

examples/plot_alpha2.jl

@@ -59,7 +59,7 @@ function simulate(integrator, steps, offset=0)
     iter / steps
 end
 
-integrator = KiteModels.init!(kps4, stiffness_factor=0.5, prn=STATISTIC)
+integrator = KiteModels.init!(kps4, prn=STATISTIC)
 
 println("\nStarting simulation...")
 simulate(integrator, 100, 100)

examples/plot_cl_cd.jl

@@ -92,7 +92,7 @@ for depower in DEPOWER
 
     kcu = KCU(set)
     kps4 = KPS4(kcu)
-    integrator = KiteModels.init!(kps4; delta=0.03, stiffness_factor=0.05, prn=STATISTIC)
+    integrator = KiteModels.init!(kps4; prn=STATISTIC)
     if ! isnothing(integrator)
         try
             simulate(kps4, integrator, logger, STEPS)

examples/plot_parking_test.jl

@@ -120,7 +120,7 @@ for depower in DEPOWER
     kps4::KPS4 = KPS4(kcu)
     set.elevation += 5
     set.v_wind = V_WIND_200[i] / 1.348881340489221 * calc_wind_factor(kps4.am, HEIGHT[i])/1.348881340489221
-    integrator = KiteModels.init!(kps4; delta=0.001, stiffness_factor=0.04, prn=STATISTIC)
+    integrator = KiteModels.init!(kps4; prn=STATISTIC)
     if ! isnothing(integrator)
         try
             simulate(kps4, integrator, logger, STEPS)

examples/plot_pitch_stability.jl

@@ -109,7 +109,7 @@ kcu::KCU = KCU(set)
 kps4::KPS4 = KPS4(kcu)
 set.elevations[1] += 5
 set.v_wind = V_WIND_200
-integrator = KiteModels.init!(kps4; delta=0.001*0, stiffness_factor=1, prn=STATISTIC)
+integrator = KiteModels.init!(kps4; prn=STATISTIC)
 if ! isnothing(integrator)
     simulate(kps4, integrator, logger, STEPS)
     set_depower_steering(kps4.kcu, depower, 0.0)

examples/plot_side_cl.jl

@@ -85,7 +85,7 @@ STEERING::Vector{Float64} = zeros(length(STEERING_RANGE))
 SIDE_CL::Vector{Float64}  = zeros(length(STEERING_RANGE))
 for (i, set_steering) in pairs(STEERING_RANGE)
     local side_cl, integrator
-    integrator = KiteModels.init!(kps4;  delta=0.001, stiffness_factor=0.1, prn=STATISTIC)
+    integrator = KiteModels.init!(kps4; prn=STATISTIC)
     side_cl, steering = simulate(integrator, STEPS, set_steering, plot=false)
     SIDE_CL[i] = side_cl
     STEERING[i] = steering

examples/reel_out_1p.jl

@@ -72,7 +72,7 @@ function simulate(integrator, steps, plot=false)
     iter / steps
 end
 
-integrator = KiteModels.init!(kps3, delta=0.002, stiffness_factor=0.1, prn=STATISTIC)
+integrator = KiteModels.init!(kps3, prn=STATISTIC)
 kps3.sync_speed = 0.0
 
 av_steps = if PLOT

examples/reel_out_4p.jl

@@ -77,7 +77,7 @@ function simulate(integrator, steps, plot=false)
     iter / steps
 end
 
-integrator = KiteModels.init!(kps4; delta=0.002, stiffness_factor=0.1, prn=STATISTIC)
+integrator = KiteModels.init!(kps4; prn=STATISTIC)
 kps4.sync_speed = 0.0
 
 av_steps = if PLOT

examples/reel_out_4p_torque_control.jl

@@ -93,7 +93,7 @@ function simulate(integrator, steps, plot=false)
     iter / steps
 end
 
-integrator = KiteModels.init!(kps4; delta=0.003, stiffness_factor=0.1, prn=STATISTIC)
+integrator = KiteModels.init!(kps4; prn=STATISTIC)
 kps4.sync_speed = 0.0
 
 av_steps = if PLOT

examples/simulate_simple.jl

@@ -65,7 +65,7 @@ function simulate(integrator, steps, plot=false)
     iter / steps
 end
 kps4.set.upwind_dir = rad2deg(UPWIND_DIR)
-integrator = KiteModels.init!(kps4;  delta=0.001, stiffness_factor=0.1, prn=STATISTIC,
+integrator = KiteModels.init!(kps4; prn=STATISTIC,
                               steady_state=false)
 
 av_steps = if PLOT

examples/simulate_steering.jl

@@ -65,7 +65,7 @@ function simulate(integrator, steps, plot=false)
     iter / steps
 end
 
-integrator = KiteModels.init!(kps3; delta=0, stiffness_factor=0.04, prn=STATISTIC)
+integrator = KiteModels.init!(kps3; prn=STATISTIC)
 
 av_steps = if PLOT
     simulate(integrator, STEPS, true)

examples/steering_test_1p.jl

@@ -101,7 +101,7 @@ function simulate(integrator, steps; plot=false)
     end
 end
 
-integrator = KiteModels.init!(kps3;  delta=0.0, stiffness_factor=1, prn=STATISTIC)
+integrator = KiteModels.init!(kps3; prn=STATISTIC)
 simulate(integrator, STEPS; plot=true)
 
 function delay(x, y, t_max = 10)

examples/steering_test_4p.jl

@@ -105,7 +105,7 @@ function simulate(integrator, steps; plot=false)
 end
 
 
-integrator = KiteModels.init!(kps4;  delta=0.001, stiffness_factor=0.1, prn=STATISTIC)
+integrator = KiteModels.init!(kps4; prn=STATISTIC)
 simulate(integrator, STEPS; plot=true)
 
 function delay(x, y, t_max = 10)