From ac53b147fa3eab2057a1013072a498cb9d99fd2b Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sat, 10 Oct 2026 17:04:08 +0200 Subject: [PATCH 1/3] Make find_steady_state! for KPS3 find a real equilibrium nlsolve usually stopped because its steps became tiny, with the accelerations of the tether particles still at about 25 m/s^2, and on some machines it printed 'solver did not converge'. The Cartesian offsets made the system badly conditioned (stiff stretching modes, soft rotation modes), the start at zero stretch sat on the kink of the spring force, and the kite can only be fully balanced at its natural elevation. Now the unknowns are the angles and relative stretches of the segments, the solver starts from a slightly stretched tether, and, as for KPS4, the elevation of the kite is prescribed instead of its vertical force balance. The horizontal residuals are taken in the plane of the tether, so the result does not depend on upwind_dir. If the pass with reduced stiffness fails, the full-stiffness pass starts again from the straight tether. The test of find_steady_state! now checks the equilibrium, the elevation and that no warning is printed, and no longer skips failures on CI. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 24 ++++++++++++---- src/KPS3.jl | 70 +++++++++++++++++++++++++++++++++++++---------- test/test-kps3.jl | 45 ++++++++++++++++-------------- 3 files changed, 97 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c916fbb3a..ead600f72 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,12 +5,24 @@ SPDX-License-Identifier: MIT ### KiteModels v0.11.19 - unreleased #### Fixed - The state of the winch (brake on/off, rate limited set speed) is now updated once per time step - in `next_step!` (and reset in `init!` to the initial set speed `set.v_reel_out` and the brake state - after construction, so that repeated calls of `init!` give the same result) instead of on every call of the residual function. Before, - it depended on the number of residual evaluations of the DAE solver, so tiny numerical - differences (e.g. between Julia versions) could switch the brake at a different moment. With a - winch `v_min` of 0.15 m/s this made the reel-out speed of the hydra20 simulations of - KiteControllers.jl run away on Julia 1.13. The simulation results change slightly. + in `next_step!` (and reset in `init!` to the initial set speed `set.v_reel_out` and the brake + state after construction, so that repeated calls of `init!` give the same result) instead of on + every call of the residual function. Before, it depended on the number of residual evaluations + of the DAE solver, so tiny numerical differences (e.g. between Julia versions) could switch the + brake at a different moment. With a winch `v_min` of 0.15 m/s this made the reel-out speed of the + hydra20 simulations of KiteControllers.jl run away on Julia 1.13. The simulation results change + slightly. +- `find_steady_state!` (KPS3) now finds a real equilibrium. Before, `nlsolve` usually stopped + because its steps became tiny, with the accelerations of the tether particles still at about + 25 m/s², and on some machines it printed `find_steady_state!: solver did not converge!`. The + unknowns are now the angles and relative stretches of the tether segments instead of the offsets + of the particle positions, the solver starts from a slightly stretched tether (the spring force + has a kink at zero stretch), and, as for KPS4, the elevation of the kite is prescribed instead of + its vertical force balance, so `calc_elevation(s)` returns `set.elevation`. The tether particles + and the horizontal force balance of the kite are now solved to `ftol = 1e-6`; a warning is only + printed if this fails. The horizontal residuals are taken in the plane of the tether, so the + result no longer depends on `upwind_dir`. The initial state changes: for example, with a 392 m + tether the kite starts at 70° instead of 64.3°. #### Changed - requires WinchModels 0.3.12 (for `update_winch_state!` and `calc_acceleration(...; update_state)`) - support Julia 1.12 and 1.13 only, as WinchModels 0.3.12 and KiteUtils do: `julia` compat diff --git a/src/KPS3.jl b/src/KPS3.jl index 96118e679..7e3918ef2 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -531,46 +531,86 @@ function spring_forces(s::KPS3) forces end -function find_steady_state_inner(s::KPS3, X, prn=false; delta=0.0, upwind_dir=nothing) +# Convert the unknowns of the steady state solver, the angle offsets of the tether segments from the +# elevation angle [rad] (first half) and their relative stretches [-] (second half), to the offsets of +# the particle positions from the straight, unstretched tether in x and z direction, as used by init. +function steady_state_offsets(s::KPS3, p) + segments = s.set.segments + l0 = s.set.l_tether / segments + elevation = deg2rad(s.set.elevation) + X = zeros(SimFloat, 2segments) + x, z = 0.0, 0.0 + for i in 1:segments + angle = elevation + p[i] + len = l0 * (1 + p[segments+i]) + x += len * cos(angle) + z += len * sin(angle) + X[i] = x - i * l0 * cos(elevation) + X[segments+i] = z - i * l0 * sin(elevation) + end + X +end + +function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=nothing, warn=true) res = zeros(MVector{6*s.set.segments+2, SimFloat}) + segments = s.set.segments + # init turns the tether into the wind direction; turn the residuals and positions back + turnangle = something(upwind_dir, -pi/2) + pi/2 + horizontal(x, y) = cos(turnangle) * x - sin(turnangle) * y # helper function for the steady state finder + # Equations: the horizontal and vertical force balance of the tether particles, the horizontal + # force balance of the kite and its elevation angle. The kite can only be fully balanced at its + # natural elevation, so the elevation from the settings is prescribed instead of its vertical + # force balance. function test_initial_condition!(F, x::Vector) - y0, yd0 = init(s, x; delta, upwind_dir) + y0, yd0 = init(s, steady_state_offsets(s, x); delta, upwind_dir) residual!(res, yd0, y0, s) - for i in 1:s.set.segments - F[i] = res[1 + 3*(i-1) + 3*s.set.segments] - F[i+s.set.segments] = res[3 + 3*(i-1) + 3*s.set.segments] + for i in 1:segments + j = 3*(i-1) + 3*segments + F[i] = horizontal(res[j+1], res[j+2]) + F[i+segments] = res[j+3] end - return nothing + # replace the vertical force balance of the kite by the elevation angle, scaled to meters + j = 3*(segments-1) + x_kite = horizontal(y0[j+1], y0[j+2]) + z_kite = y0[j+3] + F[2segments] = (atan(z_kite, x_kite) - deg2rad(s.set.elevation)) * hypot(x_kite, z_kite) + return nothing end if prn println("\nStarted function test_nlsolve...") end - jac! = make_jac(test_initial_condition!, length(X)) - results = nlsolve(test_initial_condition!, jac!, X, xtol=1e-6, ftol=1e-6, autoscale=true, iterations=1000) + jac! = make_jac(test_initial_condition!, length(p)) + results = nlsolve(test_initial_condition!, jac!, p, xtol=1e-10, ftol=1e-6, autoscale=true, iterations=1000) if prn println("\nresult: $results") end - if !converged(results) + if !results.f_converged && warn @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations))" # Check if the solution contains finite values if !all(isfinite, results.zero) error("find_steady_state!: solver returned non-finite values. Cannot compute steady state.") end end - results.zero + results.zero, results.f_converged end """ - find_steady_state!(s::KPS3; prn=false, delta = 0.0, stiffness_factor=0.035, upwind_dir=-pi/2) + find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2) Find an initial equilibrium, based on the initial parameters `l_tether`, elevation and `v_reel_out`. + +The tether particles are in equilibrium. The kite is placed at the elevation angle from the +settings, so only its horizontal force balance is solved for; the kite can only be fully balanced +at its natural elevation. """ function find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2) set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) - zero = zeros(SimFloat, 2*s.set.segments) + # start with a straight, slightly stretched tether; at zero stretch the spring force has a kink + p0 = [zeros(SimFloat, s.set.segments); fill(SimFloat(1e-3), s.set.segments)] s.stiffness_factor=stiffness_factor - zero = find_steady_state_inner(s, zero, prn; delta, upwind_dir) + p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false) s.stiffness_factor=1.0 - zero = find_steady_state_inner(s, zero, prn; delta, upwind_dir) - init(s, zero; delta=delta, upwind_dir) + # if the solution with the reduced stiffness failed, start again from the straight tether + p, _ = find_steady_state_inner(s, converged ? p : p0, prn; delta, upwind_dir) + init(s, steady_state_offsets(s, p); delta=delta, upwind_dir) end \ No newline at end of file diff --git a/test/test-kps3.jl b/test/test-kps3.jl index fd1a2232a..d91b6c177 100644 --- a/test/test-kps3.jl +++ b/test/test-kps3.jl @@ -390,28 +390,31 @@ const SEGMENTS = load_settings("system.yaml").segments local kps = KPS3(KCU(load_settings("system.yaml"))) init_392(kps) KiteModels.set_depower_steering!(kps, 0.25, 0.0) - try - res1, res2 = find_steady_state!(kps; delta = 1e-6, prn = false) - @test norm(res2) < 1e-5 # velocity and acceleration must be near zero - pre_tension = KiteModels.calc_pre_tension(kps) - @test pre_tension > 1.0001 - @test pre_tension < 1.01 - @test unstretched_length(kps) ≈ 392.0 # initial, unstretched tether length - @test tether_length(kps) ≈ 392.1861381318156 rtol = 1e-5 # real, stretched tether length - @test winch_force(kps) ≈ 276.25751212763817 rtol = 3e-2 # initial force at the winch [N] - lift, drag = lift_drag(kps) - @test lift ≈ 443.63277537186394 rtol=2e-2 # initial lift force of the kite [N] - @test drag ≈ 94.25218065939362 rtol=2e-2 # initial drag force of the kite [N] - @test lift_over_drag(kps) ≈ 4.706870146326417 rtol=2e-3 # initial lift-over-drag - @test norm(v_wind_kite(kps)) ≈ 9.107670173739065 rtol=1e-2 # initial wind speed at the height of the kite [m/s] - catch e - if e isa ErrorException && contains(e.msg, "find_steady_state!") && get(ENV, "CI", "false") == "true" - @warn "Steady state solver failed to converge. Skipping test." - @test_broken false # Mark as known issue (CI flake only) - else - rethrow(e) - end + # the solver must converge without warning; the kite is placed at the elevation from the settings + res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; delta = 1e-6, prn = false) + @test rad2deg(calc_elevation(kps)) ≈ 70.0 atol=1e-6 + # the tether particles are in equilibrium (x and z components of their accelerations) + res = zeros(length(res1)) + KiteModels.residual!(res, res2, res1, kps) + segments = kps.set.segments + for i in 1:segments-1 + j = 3*(i-1) + 3*segments + @test abs(res[j+1]) < 1e-5 + @test abs(res[j+3]) < 1e-5 end + @test abs(res[3*(segments-1) + 3*segments + 1]) < 1e-5 # horizontal balance of the kite + @test norm(res2) < 1e-5 # velocity and acceleration must be near zero + pre_tension = KiteModels.calc_pre_tension(kps) + @test pre_tension > 1.0001 + @test pre_tension < 1.01 + @test unstretched_length(kps) ≈ 392.0 # initial, unstretched tether length + @test tether_length(kps) ≈ 392.2025410439266 rtol = 1e-5 # real, stretched tether length + @test winch_force(kps) ≈ 301.23463553992514 rtol = 3e-2 # initial force at the winch [N] + lift, drag = lift_drag(kps) + @test lift ≈ 327.30289678275346 rtol=2e-2 # initial lift force of the kite [N] + @test drag ≈ 73.01286087237223 rtol=2e-2 # initial drag force of the kite [N] + @test lift_over_drag(kps) ≈ 4.482811560485004 rtol=2e-3 # initial lift-over-drag + @test norm(v_wind_kite(kps)) ≈ 9.107670173739065 rtol=1e-2 # initial wind speed at the height of the kite [m/s] end function run_benchmarks() From a6544f70e6f2483777be00328abfcec933b87631 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sat, 10 Oct 2026 17:30:13 +0200 Subject: [PATCH 2/3] Test the KPS3 steady state for several upwind directions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Checks the equilibrium in the plane of the tether for upwind_dir = -π/2, 0, π/4 and π, so that a sign error in the coordinate transform of find_steady_state_inner is detected. Co-Authored-By: Claude Opus 5.5 (1M context) --- test/test-kps3.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/test-kps3.jl b/test/test-kps3.jl index d91b6c177..a26bd683d 100644 --- a/test/test-kps3.jl +++ b/test/test-kps3.jl @@ -417,6 +417,31 @@ const SEGMENTS = load_settings("system.yaml").segments @test norm(v_wind_kite(kps)) ≈ 9.107670173739065 rtol=1e-2 # initial wind speed at the height of the kite [m/s] end + @testset "test_find_steady_state upwind_dir = $(round(rad2deg(upwind_dir)))°" for upwind_dir in (-π/2, 0.0, π/4, π) + local kps = KPS3(KCU(load_settings("system.yaml"))) + init_392(kps) + KiteModels.set_depower_steering!(kps, 0.25, 0.0) + res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; delta = 1e-6, upwind_dir) + @test rad2deg(calc_elevation(kps)) ≈ 70.0 atol=1e-6 + @test winch_force(kps) ≈ 301.23463553992514 rtol = 1e-4 # independent of the wind direction + # the tether lies in the vertical plane of the wind direction; take the components of the + # accelerations in this plane (horizontal and vertical) and normal to it + turnangle = upwind_dir + π/2 + in_plane(x, y) = cos(turnangle) * x - sin(turnangle) * y + normal(x, y) = sin(turnangle) * x + cos(turnangle) * y + segments = kps.set.segments + j = 3*(segments-1) + @test abs(normal(res1[j+1], res1[j+2])) < 1e-3 # kite position in the plane + @test in_plane(res1[j+1], res1[j+2]) > 0 # downwind of the ground station + res = zeros(length(res1)) + KiteModels.residual!(res, res2, res1, kps) + for i in 1:segments + j = 3*(i-1) + 3*segments + @test abs(in_plane(res[j+1], res[j+2])) < 1e-5 # horizontal balance in the plane + i < segments && @test abs(res[j+3]) < 1e-5 # vertical balance of the tether particles + end + end + function run_benchmarks() global height println("\ncalc_rho:") From 8e785d0048b246e903ba58e7a5e710bfc492e3b8 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sat, 10 Oct 2026 17:33:15 +0200 Subject: [PATCH 3/3] Give the link checker more time for slow sites https://doi.org/10.5281/zenodo.13310253 redirects to zenodo.org, which takes about 20 s to answer, the default timeout of lychee, so the link checker failed with a timeout on most runs. Co-Authored-By: Claude Opus 5.5 (1M context) --- .lychee.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.lychee.toml b/.lychee.toml index 2f2b2bebd..940ba0f71 100644 --- a/.lychee.toml +++ b/.lychee.toml @@ -1,3 +1,7 @@ +# doi.org redirects to zenodo.org, which often needs about 20 s to answer (default timeout: 20 s) +timeout = 60 +max_retries = 3 + exclude = [ "@ref", "@cite",