From c86a7dd08a0efd768460c7512eafb4c2f47a3b8c Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 02:24:38 +0200 Subject: [PATCH 01/21] Add trim option to find_steady_state! and fix KPS4 steady-state wind find_steady_state!(s; trim=true) (KPS3 and KPS4) solves the vertical force balance of the kite instead of prescribing set.elevation, placing the kite at the elevation where its lift matches gravity. Deep stall trims (angle of attack beyond the angle of maximum lift) are rejected, and the solver is restarted at other elevations if no valid trim is found. KPS4: the wind is now set from the kite height of each solver iterate. Before, the kite was balanced for the ground wind (height 0 after clear!), leaving force residuals of about 130 N on the kite points. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 18 ++++++++ src/KPS3.jl | 68 +++++++++++++++++++++++------ src/KPS4.jl | 107 +++++++++++++++++++++++++++++++++++----------- test/test-trim.jl | 49 +++++++++++++++++++++ 4 files changed, 203 insertions(+), 39 deletions(-) create mode 100644 test/test-trim.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index c27f8b463..c70fcceff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,24 @@ SPDX-FileCopyrightText: 2025 Uwe Fechner, Bart van de Lint SPDX-License-Identifier: MIT --> +### KiteModels v0.11.20 - 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. A second, + fully balanced trim exists at low elevations with the kite in deep stall (e.g. at 25° with an + angle of attack of 75° for KPS4 with `system.yaml`). 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 5°. The default `trim=false` keeps the previous behaviour. +#### 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. + ### KiteModels v0.11.19 2026-10-10 #### Fixed - The state of the winch (brake on/off, rate limited set speed) is now updated once per time step diff --git a/src/KPS3.jl b/src/KPS3.jl index 7e3918ef2..7c093be28 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -551,7 +551,7 @@ function steady_state_offsets(s::KPS3, p) X end -function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=nothing, warn=true) +function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=nothing, warn=true, trim=false) 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 @@ -562,7 +562,7 @@ function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=no # 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. + # force balance, unless trim is true. function test_initial_condition!(F, x::Vector) y0, yd0 = init(s, steady_state_offsets(s, x); delta, upwind_dir) residual!(res, yd0, y0, s) @@ -571,6 +571,7 @@ function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=no F[i] = horizontal(res[j+1], res[j+2]) F[i+segments] = res[j+3] end + trim && 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]) @@ -594,7 +595,8 @@ function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=no end """ - find_steady_state!(s::KPS3; prn=false, delta = 0.002, 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, + trim=false) Find an initial equilibrium, based on the initial parameters `l_tether`, elevation and `v_reel_out`. @@ -602,15 +604,53 @@ Find an initial equilibrium, based on the initial parameters 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. + +If `trim` is `true`, the vertical force balance of the kite is solved for too and the elevation +is free: the kite is placed at the elevation where its lift matches gravity (its natural +elevation), starting the search at `s.set.elevation`. Trimmed steady states with an angle of +attack beyond the angle of maximum lift (deep stall) are rejected. If no valid trim is found from +`s.set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 5°. """ -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) - # 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 - p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false) - s.stiffness_factor=1.0 - # 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 +function find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2, + trim=false) + el_set = s.set.elevation + alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] + # Find the steady state, starting at s.set.elevation. Returns the state, whether the solver + # converged and whether the kite is in deep stall. + function solve(warn) + set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) + # 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 + p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false, trim) + s.stiffness_factor=1.0 + # if the solution with the reduced stiffness failed, start again from the straight tether + p, converged = find_steady_state_inner(s, converged ? p : p0, prn; delta, upwind_dir, warn, trim) + all(isfinite, p) || return nothing, nothing, false, false + y0, yd0 = init(s, steady_state_offsets(s, p); delta=delta, upwind_dir) + trim && residual!(zeros(length(y0)), yd0, y0, s) # updates the angle of attack s.alpha_2 + y0, yd0, converged, trim && s.alpha_2 > alpha_stall + end + y0, yd0, converged, stalled = solve(!trim) + (trim && (stalled || !converged)) || return y0, yd0 + try + # Restart from other elevations: from far below the natural elevation the solver can + # find a deep stall trim. + starts = [el_set .+ (5.0:5.0:85.0); el_set .- (5.0:5.0:85.0)] + for el in filter(el -> 5.0 <= el <= 89.0, starts) + s.set.elevation = el + y, yd, conv, st = solve(false) + if conv && !st + prn && println("find_steady_state!: trimmed steady state found, starting at $el°") + return y, yd + end + end + finally + s.set.elevation = el_set + end + # no valid trim found, return the result from el_set again + y0, yd0, _, stalled = solve(true) + stalled && @warn "find_steady_state!: only a trimmed steady state in deep stall was found " * + "(angle of attack $(round(s.alpha_2, digits=1))°)" + y0, yd0 +end diff --git a/src/KPS4.jl b/src/KPS4.jl index 99a431f38..0ae9bb9df 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -705,7 +705,8 @@ function turn(res, upwind_dir) end """ - find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor=0.035, upwind_dir=-pi/2)) + find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor=0.035, upwind_dir=-pi/2, + trim=false) Find an initial equilibrium, based on the initial parameters `l_tether`, elevation and `v_reel_out`. @@ -716,14 +717,82 @@ of the kite (`pos_kite`, point B) is kept at `s.set.elevation`, because the kite be fully balanced at its natural elevation. Therefore `calc_elevation(s)` returns `s.set.elevation` after the steady state was found. +If `trim` is `true`, the z components of the forces on C and D are balanced too and the +elevation is solved for instead: the kite is placed at the elevation where its lift +matches gravity (its natural elevation), starting the search at `s.set.elevation`. +Trimmed steady states with an angle of attack beyond the angle of maximum lift (deep +stall) are rejected. If no valid trim is found from `s.set.elevation`, the solver is +restarted at higher and then at lower elevations, in steps of 5°. + The nonlinear solver is sensitive to the starting point, so small floating-point differences (e.g. between Julia versions) can decide whether it converges. If the first attempt fails, it is repeated without autoscaling and then by continuation: a solution is first found for an easier `stiffness_factor` (and `delta`) and then used as start value while stepping towards the requested values. """ -function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor=0.035, upwind_dir=-pi/2) - set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir=-pi/2) +function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor=0.035, upwind_dir=-pi/2, + trim=false) + el_set = s.set.elevation + alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] + res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) + # the state for the offsets X, which are relative to the tether at s.set.elevation + function steady_state(X) + y_, yd_ = init(s, X; upwind_dir) + set_v_wind_ground!(s, y_[3*s.set.segments+6], s.set.v_wind; upwind_dir) + residual!(res, yd_, y_, s, 0.0) + y_, yd_ + end + try + results, converged = solve_steady_state!(s; prn, delta, stiffness_factor, trim) + converged && ((y00, yd00) = steady_state(results.zero)) + stalled = trim && converged && s.alpha_2 > alpha_stall + if trim && (stalled || !converged) + # Restart from other elevations: the solver converges from above the natural + # elevation, but diverges or finds a deep stall trim from far below it. + starts = [el_set .+ (5.0:5.0:85.0); el_set .- (5.0:5.0:85.0)] + for el in filter(el -> 5.0 <= el <= 89.0, starts) + s.set.elevation = el + r, conv = solve_steady_state!(s; prn, delta, stiffness_factor, trim) + conv || continue + y, yd = steady_state(r.zero) + if s.alpha_2 <= alpha_stall + prn && println("find_steady_state!: trimmed steady state found, starting at $el°") + results, converged, stalled = r, true, false + y00, yd00 = y, yd + break + end + end + if stalled + s.set.elevation = el_set + y00, yd00 = steady_state(results.zero) + @warn "find_steady_state!: only a trimmed steady state in deep stall was found " * + "(angle of attack $(round(s.alpha_2, digits=1))°)" + end + end + if !converged + s.set.elevation = el_set + @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations), residual_norm=$(results.residual_norm))" + # 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 + y00, yd00 = steady_state(results.zero) + end + return y00, yd00 + catch e + if e isa AssertionError + error("find_steady_state!: computation resulted in non-finite values. The solver did not find a valid steady state.") + else + rethrow(e) + end + finally + s.set.elevation = el_set + end +end + +# Solve for the offsets of the steady state, see find_steady_state!. Returns the result of nlsolve +# and whether it is usable. +function solve_steady_state!(s::KPS4; prn=false, delta=0.001, stiffness_factor=0.035, trim=false) s.stiffness_factor = stiffness_factor res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) delta_ref = Ref(delta) @@ -732,6 +801,8 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= function test_initial_condition!(F, x::Vector) x1 = copy(x) y0, yd0 = init(s, x1; delta=delta_ref[]) + # the wind at the height of the kite (z of point B) of this candidate state + set_v_wind_ground!(s, y0[3*s.set.segments+6], s.set.v_wind; upwind_dir=-pi/2) try residual!(res, yd0, y0, s, 0.0) catch e @@ -746,7 +817,7 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= # 2..segments+1, p8 and p9, the symmetric part of the x residual of C and D, the elevation # of the kite, the unused entry and the antisymmetric part of the y residual of C and D. # The kite can only be fully balanced at its natural elevation, therefore the elevation - # of the kite is prescribed instead of the z residual of C and D. + # of the kite is prescribed instead of the z residual of C and D, unless trim is true. segments = s.set.segments n_free = segments + KITE_PARTICLES - 1 # number of x (and of z) unknowns acc(i, k) = res[k + 3*(i-2) + 3*(segments+KITE_PARTICLES)] # residual of particle i, component k @@ -757,8 +828,12 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= C, D = segments+4, segments+5 F[n_free] = 0.5 * (acc(C, 1) + acc(D, 1)) # position of the kite as returned by pos_kite (point B = p9 = particle segments+3) - kite = SVector(y0[3*segments+4], y0[3*segments+5], y0[3*segments+6]) - F[2*n_free] = calc_elevation(kite) - deg2rad(s.set.elevation) + if trim + F[2*n_free] = 0.5 * (acc(C, 3) + acc(D, 3)) + else + kite = SVector(y0[3*segments+4], y0[3*segments+5], y0[3*segments+6]) + F[2*n_free] = calc_elevation(kite) - deg2rad(s.set.elevation) + end F[end-1] = x[end-1] F[end] = 0.5 * (acc(C, 2) - acc(D, 2)) return nothing @@ -810,25 +885,7 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= s.stiffness_factor = stiffness_factor delta_ref[] = delta if prn println("\nresult: $results") end - if !usable(results) - @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations), residual_norm=$(results.residual_norm))" - # 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 - try - y00, yd00 = init(s, results.zero; upwind_dir) - set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) - residual!(res, yd00, y00, s, 0.0) - return y00, yd00 - catch e - if e isa AssertionError - error("find_steady_state!: computation resulted in non-finite values. The solver did not find a valid steady state.") - else - rethrow(e) - end - end + results, usable(results) end const KITE_ANGLE = 3.83 # angle between the kite and the last tether segment due to the mass of the control pod diff --git a/test/test-trim.jl b/test/test-trim.jl new file mode 100644 index 000000000..c2c8d7adf --- /dev/null +++ b/test/test-trim.jl @@ -0,0 +1,49 @@ +# SPDX-FileCopyrightText: 2026 Uwe Fechner +# SPDX-License-Identifier: MIT + +using Pkg +if dirname(Pkg.project().path) != @__DIR__ + Pkg.activate(@__DIR__) +end +using Test +using KiteModels, KitePodModels + +set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) + +# Find the steady state, starting at the elevation el0 [°]. Returns the elevation [°] and the +# largest force residual (the residual of the accelerations of all particles). +function steady_state(Model, el0; trim, delta) + set = deepcopy(load_settings("system.yaml")) + set.elevation = el0 + s = Model(KCU(set)) + KiteModels.clear!(s) + y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(s; trim, delta) + res = zeros(length(y0)) + residual!(res, yd0, y0, s, 0.0) + if trim + @test s.alpha_2 < s.set.alpha_cl[argmax(s.set.cl_list)] # no deep stall + end + particles = Model == KPS4 ? s.set.segments + KiteModels.KITE_PARTICLES : s.set.segments + n = 3 * particles + rad2deg(calc_elevation(s)), maximum(abs, res[n+1:2n]) +end + +@testset "find_steady_state! with trim=true, $Model" for Model in (KPS3, KPS4) + el_set = load_settings("system.yaml").elevation + # without trim the elevation is prescribed, so the kite is not in equilibrium + el, max_res = steady_state(Model, el_set; trim=false, delta=0.0) + @test el ≈ el_set atol=1e-6 + @test max_res > 1.0 + # with trim all forces are balanced at the natural elevation of the kite + el_trim, max_res = steady_state(Model, el_set; trim=true, delta=0.0) + @test max_res < 1e-4 + @test abs(el_trim - el_set) > 0.5 + @test 60.0 < el_trim < 80.0 + # the trimmed elevation does not depend on the start value or on delta; from 30° the solver + # finds a deep stall trim first, and for KPS4 it diverges from 55°, so both need restarts + for el0 in (30.0, 55.0, 85.0) + @test steady_state(Model, el0; trim=true, delta=0.0)[1] ≈ el_trim atol=1e-3 + end + @test steady_state(Model, el_set; trim=true, delta=0.001)[1] ≈ el_trim atol=0.05 +end +nothing From 5bab4a73a7d732c9225651d389caedec899f413b Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 02:38:13 +0200 Subject: [PATCH 02/21] Add trim option to init! and balance gravity in KPS4 trim MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit init!(s; trim=true) passes trim to find_steady_state!. KPS4: init sets the vertical accelerations in yd0 to -g, so the steady state solver balanced the forces without gravity. With trim=true the accelerations are now zero, so the kite is trimmed with its weight and the returned state is an equilibrium of the simulation: with init!(s; trim=true, stiffness_factor=1.0) the elevation stays within 0.1° over 10 s. trim=false is unchanged. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 17 +++++++++++------ src/KPS4.jl | 18 +++++++++++++++--- src/KiteModels.jl | 10 +++++++--- test/test-trim.jl | 18 +++++++++++++++++- 4 files changed, 50 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c70fcceff..152fce07d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,12 +6,17 @@ SPDX-License-Identifier: MIT #### 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. A second, - fully balanced trim exists at low elevations with the kite in deep stall (e.g. at 25° with an - angle of attack of 75° for KPS4 with `system.yaml`). 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 5°. The default `trim=false` keeps the previous behaviour. + 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 5°. The default `trim=false` keeps the previous behaviour. +- `init!(s; trim=true)` passes `trim` to `find_steady_state!`. The trim depends on the stiffness + factor, which `next_step!` raises to 1.0, so use `init!(s; trim=true, stiffness_factor=1.0)` for an + initial state that is an equilibrium of the simulation; then the elevation stays constant + (within 0.1° over 10 s for KPS4 with `system.yaml`). #### 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!`, diff --git a/src/KPS4.jl b/src/KPS4.jl index 0ae9bb9df..fb5f887b8 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -717,9 +717,10 @@ of the kite (`pos_kite`, point B) is kept at `s.set.elevation`, because the kite be fully balanced at its natural elevation. Therefore `calc_elevation(s)` returns `s.set.elevation` after the steady state was found. -If `trim` is `true`, the z components of the forces on C and D are balanced too and the -elevation is solved for instead: the kite is placed at the elevation where its lift -matches gravity (its natural elevation), starting the search at `s.set.elevation`. +If `trim` is `true`, all forces including gravity are balanced, also the z components of +the forces on C and D, and the elevation is solved for instead: the kite is placed at the +elevation where its lift matches gravity (its natural elevation), starting the search at +`s.set.elevation`. Trimmed steady states with an angle of attack beyond the angle of maximum lift (deep stall) are rejected. If no valid trim is found from `s.set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 5°. @@ -738,6 +739,11 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= # the state for the offsets X, which are relative to the tether at s.set.elevation function steady_state(X) y_, yd_ = init(s, X; upwind_dir) + # the accelerations of a steady state that is balanced with gravity are zero + if trim + n = 3*(s.set.segments+KITE_PARTICLES) + yd_[n+1:2n] .= 0.0 + end set_v_wind_ground!(s, y_[3*s.set.segments+6], s.set.v_wind; upwind_dir) residual!(res, yd_, y_, s, 0.0) y_, yd_ @@ -801,6 +807,12 @@ function solve_steady_state!(s::KPS4; prn=false, delta=0.001, stiffness_factor=0 function test_initial_condition!(F, x::Vector) x1 = copy(x) y0, yd0 = init(s, x1; delta=delta_ref[]) + # init sets the vertical accelerations to -g, so that gravity is not balanced; with trim + # the accelerations are zero, so that the kite and the tether are balanced with gravity + if trim + n = 3*(s.set.segments+KITE_PARTICLES) + yd0[n+1:2n] .= 0.0 + end # the wind at the height of the kite (z of point B) of this candidate state set_v_wind_ground!(s, y0[3*s.set.segments+6], s.set.v_wind; upwind_dir=-pi/2) try diff --git a/src/KiteModels.jl b/src/KiteModels.jl index fdbc7fd79..24c54ff04 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -640,7 +640,7 @@ end """ init!(s::AKM; stiffness_factor=0.5, delta=0.005, - prn=false, steady_state=true) -> Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} + prn=false, steady_state=true, trim=false) -> Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} Initializes the integrator of the model (KPS3 and KPS4 only). @@ -655,11 +655,15 @@ Parameters: produces a non-equilibrium starting point, which is useful when the caller has already set up a consistent state or when the steady-state solver is known to struggle (e.g. unusual wind conditions). Returns `nothing` if steady-state solving fails. +- trim: if `true`, the steady state is found at the natural elevation of the kite, where its + lift matches gravity, instead of at `s.set.elevation`; see `find_steady_state!`. Ignored if + `steady_state` is `false`. Returns: An instance of an `ODEIntegrator` or `IDAIntegrator`, or `nothing` if initialization fails. """ -function init!(s::AKM; stiffness_factor=0.5, delta=0.005, prn=false, steady_state=true)::Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} +function init!(s::AKM; stiffness_factor=0.5, delta=0.005, prn=false, steady_state=true, + trim=false)::Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} clear!(s) init_winch_state!(s) upwind_dir = deg2rad(Float64(s.set.upwind_dir)) @@ -667,7 +671,7 @@ function init!(s::AKM; stiffness_factor=0.5, delta=0.005, prn=false, steady_stat y0, yd0 = if steady_state try - KiteModels.find_steady_state!(s; stiffness_factor, delta, upwind_dir, prn) + KiteModels.find_steady_state!(s; stiffness_factor, delta, upwind_dir, prn, trim) catch e if e isa AssertionError || (e isa ErrorException && startswith(e.msg, "find_steady_state!:")) println("ERROR: Failure to find initial steady state in find_steady_state! function!\n"* diff --git a/test/test-trim.jl b/test/test-trim.jl index c2c8d7adf..b501328f9 100644 --- a/test/test-trim.jl +++ b/test/test-trim.jl @@ -38,7 +38,7 @@ end el_trim, max_res = steady_state(Model, el_set; trim=true, delta=0.0) @test max_res < 1e-4 @test abs(el_trim - el_set) > 0.5 - @test 60.0 < el_trim < 80.0 + @test 50.0 < el_trim < 80.0 # the trimmed elevation does not depend on the start value or on delta; from 30° the solver # finds a deep stall trim first, and for KPS4 it diverges from 55°, so both need restarts for el0 in (30.0, 55.0, 85.0) @@ -46,4 +46,20 @@ end end @test steady_state(Model, el_set; trim=true, delta=0.001)[1] ≈ el_trim atol=0.05 end + +@testset "init! with trim=true, $Model" for Model in (KPS3, KPS4) + # at full stiffness the trimmed initial state is an equilibrium of the simulation, + # so the elevation stays constant + s = Model(KCU(deepcopy(load_settings("system.yaml")))) + integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true, stiffness_factor=1.0) + el0 = calc_elevation(s) + @test abs(rad2deg(el0) - s.set.elevation) > 0.5 + dt = 1/s.set.sample_freq + max_dev = 0.0 + for _ in 1:round(Int, 10/dt) + next_step!(s, integrator; set_speed=0.0, dt) + max_dev = max(max_dev, abs(rad2deg(calc_elevation(s) - el0))) + end + @test max_dev < 0.2 +end nothing From 46a1fec48c639f208670b19f12fad80eb169073a Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 03:01:23 +0200 Subject: [PATCH 03/21] Test trim for the Maasvlakte and Cabauw wind profiles MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit For KPS3 and KPS4 at 4, 7 and 10 m/s: after init!(s; trim=true, stiffness_factor=1.0, delta=0.0) the wind at the kite matches the profile, the kite is not stalled, and the elevation stays within 0.01° over 10 s. Maasvlakte at 4 m/s is skipped, because there the kite cannot be carried with attached flow. Co-Authored-By: Claude Opus 5.5 --- test/test-trim.jl | 41 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/test/test-trim.jl b/test/test-trim.jl index b501328f9..704d0c059 100644 --- a/test/test-trim.jl +++ b/test/test-trim.jl @@ -5,7 +5,7 @@ using Pkg if dirname(Pkg.project().path) != @__DIR__ Pkg.activate(@__DIR__) end -using Test +using LinearAlgebra, Test using KiteModels, KitePodModels set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) @@ -62,4 +62,43 @@ end end @test max_dev < 0.2 end + +# wind profiles of AtmosphericModels.jl, fitted to the near-shore location Maasvlakte, NL and +# to the inland location Cabauw, NL +const PROFILES = Dict("Maasvlakte" => (profile_law=3, alpha=0.08163, z0=0.0002), + "Cabauw" => (profile_law=1, alpha=0.234, z0=0.0002)) + +function profile_settings(profile, v_wind) + set = deepcopy(load_settings("system.yaml")) + p = PROFILES[profile] + set.profile_law, set.alpha, set.z0, set.v_wind = p.profile_law, p.alpha, p.z0, v_wind + set +end + +@testset "trim for wind profile $profile, $Model" for Model in (KPS3, KPS4), profile in ("Maasvlakte", "Cabauw") + elevations = Float64[] + for v_wind in (4.0, 7.0, 10.0) + # at 4 m/s, Maasvlakte has only 5.2 m/s at the height of the kite: the lift is too small + # to carry the kite, so no trim with attached flow exists + profile == "Maasvlakte" && v_wind == 4.0 && continue + s = Model(KCU(profile_settings(profile, v_wind))) + integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true, stiffness_factor=1.0, delta=0.0) + el0 = calc_elevation(s) + push!(elevations, el0) + # the kite is balanced for the wind of the profile at its height, with attached flow + v_profile = v_wind * KiteModels.calc_wind_factor(s.am, calc_height(s)) + @test norm(s.v_wind) ≈ v_profile rtol=1e-6 + @test s.alpha_2 < s.set.alpha_cl[argmax(s.set.cl_list)] + # the trimmed state is an equilibrium of the simulation, so the elevation stays constant + dt = 1/s.set.sample_freq + max_dev = 0.0 + for _ in 1:round(Int, 10/dt) + next_step!(s, integrator; set_speed=0.0, dt) + max_dev = max(max_dev, abs(rad2deg(calc_elevation(s) - el0))) + end + @test max_dev < 0.01 + end + # the kite flies higher at higher wind speeds + @test issorted(elevations) +end nothing From 9587cabb4b47f29bcc5ce33668de5a22a74cb076 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 03:23:50 +0200 Subject: [PATCH 04/21] Deprecate delta and warn about deep stall at a prescribed elevation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The keyword argument delta of find_steady_state! and init! is deprecated and ignored. 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, so it is an equilibrium of the simulation. KPS3 now converges in all 72 cases of a sweep over wind speed, elevation, tether length and stiffness factor. The examples, tests and precompile.jl no longer pass delta. With trim=false, find_steady_state! warns if the kite is in deep stall, which happens if the prescribed elevation is far below the natural elevation (with system.yaml below about 60°). The robustness test cases at 50° are moved to 60° and 65° for this reason. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 17 +++++++-- examples/bench.jl | 2 +- examples/bench_4p.jl | 2 +- examples/calc_spectrum.jl | 2 +- examples/compare_kps3_kps4.jl | 4 +-- examples/plot_cl_cd.jl | 2 +- examples/plot_parking_test.jl | 2 +- examples/plot_pitch_stability.jl | 2 +- examples/plot_side_cl.jl | 2 +- examples/reel_out_1p.jl | 2 +- examples/reel_out_4p.jl | 2 +- examples/reel_out_4p_torque_control.jl | 2 +- examples/simulate_simple.jl | 2 +- examples/simulate_steering.jl | 2 +- examples/steering_test_1p.jl | 2 +- examples/steering_test_4p.jl | 2 +- examples/test_init_1p.jl | 2 +- examples/test_init_4p.jl | 2 +- examples/test_steady_state.jl | 2 +- examples_3d/auto_parking_4p.jl | 2 +- examples_3d/parking_4p.jl | 2 +- examples_3d/parking_wind_dir.jl | 2 +- src/KPS3.jl | 50 +++++++++++++++++--------- src/KPS4.jl | 35 +++++++++++++----- src/KiteModels.jl | 28 +++++++++++---- src/precompile.jl | 4 +-- test/plot_kps4.jl | 2 +- test/test-interface.jl | 8 ++--- test/test-kps3.jl | 4 +-- test/test-kps4.jl | 2 +- test/test-simulate_4p.jl | 2 +- test/test-steady-state-kps4.jl | 2 +- test/test-steady-state-robustness.jl | 29 +++++++-------- test/test-trim.jl | 41 ++++++++++++++++----- test/test-update-sys-state.jl | 4 +-- test/test_for_precompile.jl | 2 +- 36 files changed, 182 insertions(+), 92 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 152fce07d..a629cd5db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,8 +15,21 @@ SPDX-License-Identifier: MIT lower elevations, in steps of 5°. The default `trim=false` keeps the previous behaviour. - `init!(s; trim=true)` passes `trim` to `find_steady_state!`. The trim depends on the stiffness factor, which `next_step!` raises to 1.0, so use `init!(s; trim=true, stiffness_factor=1.0)` for an - initial state that is an equilibrium of the simulation; then the elevation stays constant - (within 0.1° over 10 s for KPS4 with `system.yaml`). + initial state that is an equilibrium of the simulation; then the elevation stays constant. +- `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 60°), 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 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!`, diff --git a/examples/bench.jl b/examples/bench.jl index 46b2e2a54..971971792 100644 --- a/examples/bench.jl +++ b/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; stiffness_factor = 0.1, prn = STATISTIC) println("\nStarting simulation...") simulate(integrator, 100) diff --git a/examples/bench_4p.jl b/examples/bench_4p.jl index 69b10e47a..7a913b77f 100644 --- a/examples/bench_4p.jl +++ b/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, stiffness_factor = 0.1, prn = STATISTIC) println("\nStarting simulation...") simulate(integrator, 100, 100) diff --git a/examples/calc_spectrum.jl b/examples/calc_spectrum.jl index 266b0b247..6f0082db0 100644 --- a/examples/calc_spectrum.jl +++ b/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; stiffness_factor=0.04, 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 diff --git a/examples/compare_kps3_kps4.jl b/examples/compare_kps3_kps4.jl index a3bd06437..470a77e9e 100644 --- a/examples/compare_kps3_kps4.jl +++ b/examples/compare_kps3_kps4.jl @@ -84,7 +84,7 @@ function simulate(s, integrator, steps, plot = false; fig = "") iter / steps end -integrator = init!(kps3, delta = 0.001, stiffness_factor = 0.1, prn = STATISTIC) +integrator = init!(kps3, stiffness_factor = 0.1, 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; stiffness_factor = 0.1, prn = STATISTIC) av_steps = simulate(kps4, integrator, STEPS, true; fig = "kps4") reactivate_host_app() diff --git a/examples/plot_cl_cd.jl b/examples/plot_cl_cd.jl index 4e2731f0c..61d9b39a8 100644 --- a/examples/plot_cl_cd.jl +++ b/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; stiffness_factor=0.05, prn=STATISTIC) if ! isnothing(integrator) try simulate(kps4, integrator, logger, STEPS) diff --git a/examples/plot_parking_test.jl b/examples/plot_parking_test.jl index 4b0f70675..55dbeea8e 100644 --- a/examples/plot_parking_test.jl +++ b/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; stiffness_factor=0.04, prn=STATISTIC) if ! isnothing(integrator) try simulate(kps4, integrator, logger, STEPS) diff --git a/examples/plot_pitch_stability.jl b/examples/plot_pitch_stability.jl index ec258939f..59ca0421f 100644 --- a/examples/plot_pitch_stability.jl +++ b/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; stiffness_factor=1, prn=STATISTIC) if ! isnothing(integrator) simulate(kps4, integrator, logger, STEPS) set_depower_steering(kps4.kcu, depower, 0.0) diff --git a/examples/plot_side_cl.jl b/examples/plot_side_cl.jl index 9823fa3d8..34a718fb8 100644 --- a/examples/plot_side_cl.jl +++ b/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; stiffness_factor=0.1, prn=STATISTIC) side_cl, steering = simulate(integrator, STEPS, set_steering, plot=false) SIDE_CL[i] = side_cl STEERING[i] = steering diff --git a/examples/reel_out_1p.jl b/examples/reel_out_1p.jl index d49c8ed3f..f19e71441 100644 --- a/examples/reel_out_1p.jl +++ b/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, stiffness_factor=0.1, prn=STATISTIC) kps3.sync_speed = 0.0 av_steps = if PLOT diff --git a/examples/reel_out_4p.jl b/examples/reel_out_4p.jl index 4df8e0207..ad301029b 100644 --- a/examples/reel_out_4p.jl +++ b/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; stiffness_factor=0.1, prn=STATISTIC) kps4.sync_speed = 0.0 av_steps = if PLOT diff --git a/examples/reel_out_4p_torque_control.jl b/examples/reel_out_4p_torque_control.jl index e75d45ca6..e0996fd6b 100644 --- a/examples/reel_out_4p_torque_control.jl +++ b/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; stiffness_factor=0.1, prn=STATISTIC) kps4.sync_speed = 0.0 av_steps = if PLOT diff --git a/examples/simulate_simple.jl b/examples/simulate_simple.jl index 4b68ac33e..8fb079fd0 100644 --- a/examples/simulate_simple.jl +++ b/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; stiffness_factor=0.1, prn=STATISTIC, steady_state=false) av_steps = if PLOT diff --git a/examples/simulate_steering.jl b/examples/simulate_steering.jl index 116dda484..8744a3bc1 100644 --- a/examples/simulate_steering.jl +++ b/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; stiffness_factor=0.04, prn=STATISTIC) av_steps = if PLOT simulate(integrator, STEPS, true) diff --git a/examples/steering_test_1p.jl b/examples/steering_test_1p.jl index 5e9545a37..8ca941000 100644 --- a/examples/steering_test_1p.jl +++ b/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; stiffness_factor=1, prn=STATISTIC) simulate(integrator, STEPS; plot=true) function delay(x, y, t_max = 10) diff --git a/examples/steering_test_4p.jl b/examples/steering_test_4p.jl index 8e185439b..12c3d5b9d 100644 --- a/examples/steering_test_4p.jl +++ b/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; stiffness_factor=0.1, prn=STATISTIC) simulate(integrator, STEPS; plot=true) function delay(x, y, t_max = 10) diff --git a/examples/test_init_1p.jl b/examples/test_init_1p.jl index 785c48374..613d448b1 100644 --- a/examples/test_init_1p.jl +++ b/examples/test_init_1p.jl @@ -33,7 +33,7 @@ logger = Logger(set.segments + 1, STEPS) kcu::KCU = KCU(set) kps3::KPS3 = KPS3(kcu) -_ = KiteModels.init!(kps3; delta=0.03, stiffness_factor=0.1, prn=STATISTIC) +_ = KiteModels.init!(kps3; stiffness_factor=0.1, prn=STATISTIC) lift, drag = lift_drag(kps3) sys_state::SysState = KiteModels.SysState(kps3) log!(logger, sys_state) diff --git a/examples/test_init_4p.jl b/examples/test_init_4p.jl index d24c33727..844454b72 100644 --- a/examples/test_init_4p.jl +++ b/examples/test_init_4p.jl @@ -35,7 +35,7 @@ logger = Logger(set.segments + 5, STEPS) kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) -integrator = KiteModels.init!(kps4; delta=0.001, stiffness_factor=0.1, prn=STATISTIC) +integrator = KiteModels.init!(kps4; stiffness_factor=0.1, prn=STATISTIC) for _ in 1:80 next_step!(kps4, integrator; set_speed=kps4.set.v_reel_out, upwind_dir=deg2rad(UPWIND_DIR), dt) end diff --git a/examples/test_steady_state.jl b/examples/test_steady_state.jl index 2d86e3e25..5932fbf8a 100644 --- a/examples/test_steady_state.jl +++ b/examples/test_steady_state.jl @@ -13,7 +13,7 @@ kps4::KPS4 = KPS4(kcu) clear!(kps4) KiteModels.set_depower_steering!(kps4, kps4.set.depower_offset/100.0, 0.0) -@time KiteModels.init!(kps4; delta=0.001, stiffness_factor=0.5, prn=false) +@time KiteModels.init!(kps4; stiffness_factor=0.5, prn=false) println("\nlift, drag [N] : $(KiteModels.lift_drag(kps4))") println("\nSpring forces:") diff --git a/examples_3d/auto_parking_4p.jl b/examples_3d/auto_parking_4p.jl index 0d71e74b0..126f0ce4c 100644 --- a/examples_3d/auto_parking_4p.jl +++ b/examples_3d/auto_parking_4p.jl @@ -170,7 +170,7 @@ function simulate(integrator, sys_state) end function play() - integrator = KiteModels.init!(kps4; delta = 0.001, stiffness_factor = 0.1) + integrator = KiteModels.init!(kps4; stiffness_factor = 0.1) sys_state = SysState(kps4) toc() simulate(integrator, sys_state) diff --git a/examples_3d/parking_4p.jl b/examples_3d/parking_4p.jl index c6691be86..58726add8 100644 --- a/examples_3d/parking_4p.jl +++ b/examples_3d/parking_4p.jl @@ -98,7 +98,7 @@ function simulate(integrator, steps, plot = PLOT) end toc() -integrator = KiteModels.init!(kps4, delta = 0.001, stiffness_factor = 0.1, prn = STATISTIC) +integrator = KiteModels.init!(kps4, stiffness_factor = 0.1, prn = STATISTIC) toc() println("\nStarting simulation...") diff --git a/examples_3d/parking_wind_dir.jl b/examples_3d/parking_wind_dir.jl index 795f31a91..8359355ef 100644 --- a/examples_3d/parking_wind_dir.jl +++ b/examples_3d/parking_wind_dir.jl @@ -190,7 +190,7 @@ end function play_parking() saved_turbulence = kps4.set.use_turbulence kps4.set.use_turbulence = 0.0 - integrator = KiteModels.init!(kps4; delta=0.001, stiffness_factor=0.1) + integrator = KiteModels.init!(kps4; stiffness_factor=0.1) kps4.set.use_turbulence = saved_turbulence sys_state = SysState(kps4) if KiteUtils.PROJECT == "system.yaml" diff --git a/src/KPS3.jl b/src/KPS3.jl index 7c093be28..cc8d736f9 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -595,44 +595,62 @@ function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=no end """ - find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2, - trim=false) + find_steady_state!(s::KPS3; prn=false, stiffness_factor=0.035, upwind_dir=-pi/2, trim=false) 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. +at its natural elevation. If the prescribed elevation is far below the natural elevation, the kite +can only be balanced in deep stall, with an angle of attack beyond the angle of maximum lift; then +a warning is logged. If `trim` is `true`, the vertical force balance of the kite is solved for too and the elevation is free: the kite is placed at the elevation where its lift matches gravity (its natural elevation), starting the search at `s.set.elevation`. Trimmed steady states with an angle of attack beyond the angle of maximum lift (deep stall) are rejected. If no valid trim is found from `s.set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 5°. + +The keyword argument `delta` is deprecated and ignored. """ -function find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2, +function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor=0.035, upwind_dir=-pi/2, trim=false) + warn_delta_deprecated(delta, :find_steady_state!) el_set = s.set.elevation alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] # Find the steady state, starting at s.set.elevation. Returns the state, whether the solver # converged and whether the kite is in deep stall. function solve(warn) - set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) - # 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 - p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false, trim) - s.stiffness_factor=1.0 - # if the solution with the reduced stiffness failed, start again from the straight tether - p, converged = find_steady_state_inner(s, converged ? p : p0, prn; delta, upwind_dir, warn, trim) + local p, converged + for delta in STEADY_STATE_DELTAS + set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) + # 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 + p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false, trim) + s.stiffness_factor=1.0 + # if the solution with the reduced stiffness failed, start again from the straight tether + p, converged = find_steady_state_inner(s, converged ? p : p0, prn; delta, upwind_dir, warn=false, trim) + converged && break + end + if !converged && warn + @warn "find_steady_state!: solver did not converge!" + all(isfinite, p) || error("find_steady_state!: solver returned non-finite values. Cannot compute steady state.") + end all(isfinite, p) || return nothing, nothing, false, false - y0, yd0 = init(s, steady_state_offsets(s, p); delta=delta, upwind_dir) - trim && residual!(zeros(length(y0)), yd0, y0, s) # updates the angle of attack s.alpha_2 - y0, yd0, converged, trim && s.alpha_2 > alpha_stall + y0, yd0 = init(s, steady_state_offsets(s, p); upwind_dir) + residual!(zeros(length(y0)), yd0, y0, s) # updates the angle of attack s.alpha_2 + y0, yd0, converged, s.alpha_2 > alpha_stall end y0, yd0, converged, stalled = solve(!trim) - (trim && (stalled || !converged)) || return y0, yd0 + if !trim + converged && stalled && @warn "find_steady_state!: the kite is in deep stall at the prescribed " * + "elevation of $(el_set)° (angle of attack $(round(s.alpha_2, digits=1))°); use trim=true or a " * + "higher elevation" + return y0, yd0 + end + (stalled || !converged) || return y0, yd0 try # Restart from other elevations: from far below the natural elevation the solver can # find a deep stall trim. diff --git a/src/KPS4.jl b/src/KPS4.jl index fb5f887b8..e24be9eb6 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -705,8 +705,7 @@ function turn(res, upwind_dir) end """ - find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor=0.035, upwind_dir=-pi/2, - trim=false) + find_steady_state!(s::KPS4; prn=false, stiffness_factor=0.035, upwind_dir=-pi/2, trim=false) Find an initial equilibrium, based on the initial parameters `l_tether`, elevation and `v_reel_out`. @@ -715,7 +714,9 @@ The forces on the tether particles, the KCU and the kite particles are balanced gravity), except for the z component of the kite points C and D: instead, the elevation of the kite (`pos_kite`, point B) is kept at `s.set.elevation`, because the kite can only be fully balanced at its natural elevation. Therefore `calc_elevation(s)` returns -`s.set.elevation` after the steady state was found. +`s.set.elevation` after the steady state was found. If the prescribed elevation is far +below the natural elevation, the kite can only be balanced in deep stall, with an angle of +attack beyond the angle of maximum lift; then a warning is logged. If `trim` is `true`, all forces including gravity are balanced, also the z components of the forces on C and D, and the elevation is solved for instead: the kite is placed at the @@ -728,11 +729,24 @@ restarted at higher and then at lower elevations, in steps of 5°. The nonlinear solver is sensitive to the starting point, so small floating-point differences (e.g. between Julia versions) can decide whether it converges. If the first attempt fails, it is repeated without autoscaling and then by continuation: -a solution is first found for an easier `stiffness_factor` (and `delta`) and then -used as start value while stepping towards the requested values. +a solution is first found for an easier `stiffness_factor` and then used as start value +while stepping towards the requested value. The initial state is not perturbed, unless +this fails; the returned steady state is never perturbed. + +The keyword argument `delta` is deprecated and ignored. """ -function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor=0.035, upwind_dir=-pi/2, +function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor=0.035, upwind_dir=-pi/2, trim=false) + warn_delta_deprecated(delta, :find_steady_state!) + # solve without perturbation of the initial state, with a perturbation only as fallback + function solve() + local r, conv + for delta in STEADY_STATE_DELTAS + r, conv = solve_steady_state!(s; prn, delta, stiffness_factor, trim) + conv && break + end + r, conv + end el_set = s.set.elevation alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) @@ -749,7 +763,7 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= y_, yd_ end try - results, converged = solve_steady_state!(s; prn, delta, stiffness_factor, trim) + results, converged = solve() converged && ((y00, yd00) = steady_state(results.zero)) stalled = trim && converged && s.alpha_2 > alpha_stall if trim && (stalled || !converged) @@ -758,7 +772,7 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= starts = [el_set .+ (5.0:5.0:85.0); el_set .- (5.0:5.0:85.0)] for el in filter(el -> 5.0 <= el <= 89.0, starts) s.set.elevation = el - r, conv = solve_steady_state!(s; prn, delta, stiffness_factor, trim) + r, conv = solve() conv || continue y, yd = steady_state(r.zero) if s.alpha_2 <= alpha_stall @@ -784,6 +798,11 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= end y00, yd00 = steady_state(results.zero) end + if !trim && converged && s.alpha_2 > alpha_stall + @warn "find_steady_state!: the kite is in deep stall at the prescribed elevation of " * + "$(el_set)° (angle of attack $(round(s.alpha_2, digits=1))°); use trim=true or a " * + "higher elevation" + end return y00, yd00 catch e if e isa AssertionError diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 24c54ff04..6c7dfa892 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -77,6 +77,18 @@ Short alias for the AbstractKiteModel. """ const AKM = AbstractKiteModel +# Perturbations of the initial state, tried in this order by the steady state solvers: delta +# shifts all particles sideways by delta [m] and gives them the velocity delta [m/s]. A perturbed +# start can help the nonlinear solver, but the returned steady state is never perturbed. +const STEADY_STATE_DELTAS = (0.0, 0.001) + +function warn_delta_deprecated(delta, f) + isnothing(delta) || Base.depwarn("The keyword argument `delta` of `$f` is deprecated and " * + "ignored: the steady state is solved without perturbation (with a small perturbation " * + "only as fallback) and returned without perturbation.", f) + nothing +end + function __init__() if isdir(joinpath(pwd(), "data")) && isfile(joinpath(pwd(), "data", "system.yaml")) set_data_path(joinpath(pwd(), "data")) @@ -639,15 +651,14 @@ function init_winch_state!(s::AKM) end """ - init!(s::AKM; stiffness_factor=0.5, delta=0.005, - prn=false, steady_state=true, trim=false) -> Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} + init!(s::AKM; stiffness_factor=0.5, prn=false, steady_state=true, + trim=false) -> Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} Initializes the integrator of the model (KPS3 and KPS4 only). Parameters: - s: an instance of an abstract kite model - stiffness_factor: factor applied to the tether stiffness during initialization -- delta: initial stretch of the tether during the steady state calculation - prn: if set to true, print the detailed solver results - steady_state: if `true` (default), call `find_steady_state!` to solve for an equilibrium before constructing the integrator. If `false`, skip the steady-state solver and initialise @@ -659,11 +670,14 @@ Parameters: lift matches gravity, instead of at `s.set.elevation`; see `find_steady_state!`. Ignored if `steady_state` is `false`. +The keyword argument `delta` is deprecated and ignored. + Returns: An instance of an `ODEIntegrator` or `IDAIntegrator`, or `nothing` if initialization fails. """ -function init!(s::AKM; stiffness_factor=0.5, delta=0.005, prn=false, steady_state=true, +function init!(s::AKM; stiffness_factor=0.5, delta=nothing, prn=false, steady_state=true, trim=false)::Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} + warn_delta_deprecated(delta, :init!) clear!(s) init_winch_state!(s) upwind_dir = deg2rad(Float64(s.set.upwind_dir)) @@ -671,18 +685,18 @@ function init!(s::AKM; stiffness_factor=0.5, delta=0.005, prn=false, steady_stat y0, yd0 = if steady_state try - KiteModels.find_steady_state!(s; stiffness_factor, delta, upwind_dir, prn, trim) + KiteModels.find_steady_state!(s; stiffness_factor, upwind_dir, prn, trim) catch e if e isa AssertionError || (e isa ErrorException && startswith(e.msg, "find_steady_state!:")) println("ERROR: Failure to find initial steady state in find_steady_state! function!\n"* - "Try to increase the delta parameter or to decrease the stiffness_factor of the init! function.") + "Try to decrease the stiffness_factor of the init! function.") return nothing else rethrow(e) end end else - y0, yd0 = KiteModels.init(s; delta, upwind_dir) + y0, yd0 = KiteModels.init(s; upwind_dir) set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) y0, yd0 diff --git a/src/precompile.jl b/src/precompile.jl index 62009fbf5..7a1667a2f 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -18,8 +18,8 @@ @compile_workload begin # all calls in this block will be precompiled, regardless of whether( # they belong to your package or not (on Julia 1.8 and higher) - integrator = KiteModels.init!(kps3_; delta=0.001, stiffness_factor=0.035, prn=false) - integrator = KiteModels.init!(kps4_; delta=0.001, stiffness_factor=0.035, prn=false) + integrator = KiteModels.init!(kps3_; stiffness_factor=0.035, prn=false) + integrator = KiteModels.init!(kps4_; stiffness_factor=0.035, prn=false) nothing end end \ No newline at end of file diff --git a/test/plot_kps4.jl b/test/plot_kps4.jl index a71c702f4..0227b3d3e 100644 --- a/test/plot_kps4.jl +++ b/test/plot_kps4.jl @@ -23,7 +23,7 @@ kps4.stiffness_factor = 0.04 y0, yd0 = KiteModels.init(kps4) -find_steady_state!(kps4; delta=0.001, stiffness_factor=1, prn=false) +find_steady_state!(kps4; stiffness_factor=1, prn=false) println("kite distance: $(norm(kps4.pos[end-2]))") println(KiteModels.spring_forces(kps4)) diff --git a/test/test-interface.jl b/test/test-interface.jl index 4dbe9de62..0876fb1ee 100644 --- a/test/test-interface.jl +++ b/test/test-interface.jl @@ -32,28 +32,28 @@ end @testset "KPS3 init! interface" begin kps4.set.upwind_dir = rad2deg(-pi/2) - integ = init!(kps3; stiffness_factor=0.5, delta=0.0001, prn=false) + integ = init!(kps3; stiffness_factor=0.5, prn=false) @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator @test kps3.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end @testset "KPS4 init! interface" begin kps4.set.upwind_dir = rad2deg(-pi/2) - integ = init!(kps4; stiffness_factor=0.5, delta=0.005, prn=false) + integ = init!(kps4; stiffness_factor=0.5, prn=false) @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator @test kps4.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end @testset "KPS3 init! steady_state=false" begin kps3.set.upwind_dir = rad2deg(-pi/2) - integ = init!(kps3; delta=0.0001, steady_state=false) + integ = init!(kps3; steady_state=false) @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator @test kps3.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end @testset "KPS4 init! steady_state=false" begin kps4.set.upwind_dir = rad2deg(-pi/2) - integ = init!(kps4; delta=0.005, steady_state=false) + integ = init!(kps4; steady_state=false) @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator @test kps4.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end diff --git a/test/test-kps3.jl b/test/test-kps3.jl index a26bd683d..2fd757d38 100644 --- a/test/test-kps3.jl +++ b/test/test-kps3.jl @@ -391,7 +391,7 @@ const SEGMENTS = load_settings("system.yaml").segments init_392(kps) KiteModels.set_depower_steering!(kps, 0.25, 0.0) # 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) + res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; 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)) @@ -421,7 +421,7 @@ 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) - res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; delta = 1e-6, upwind_dir) + res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; 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 diff --git a/test/test-kps4.jl b/test/test-kps4.jl index dbdc5b583..f02095e39 100644 --- a/test/test-kps4.jl +++ b/test/test-kps4.jl @@ -687,7 +687,7 @@ set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) # and the updraft shows up in v_wind_kite (and therefore in the logged system state) local kps4_ = KPS4(KCU(deepcopy(load_settings("system.yaml")))) kps4_.set.kcu_diameter = 0.0 - integrator = KiteModels.init!(kps4_; stiffness_factor=0.5, delta=0.005, prn=false) + integrator = KiteModels.init!(kps4_; stiffness_factor=0.5, prn=false) dt = 1/kps4_.set.sample_freq for _ in 1:5 next_step!(kps4_, integrator; set_speed=0, v_wind_vert=1.5, dt) diff --git a/test/test-simulate_4p.jl b/test/test-simulate_4p.jl index f96871b81..9cc3c9060 100644 --- a/test/test-simulate_4p.jl +++ b/test/test-simulate_4p.jl @@ -44,7 +44,7 @@ end kps4.set.solver = "IDA" kps4.set.profile_law = Int(EXPLOG) kps4.set.alpha = 0.08163 - integrator = KiteModels.init!(kps4; stiffness_factor=0.4, delta=0.0015, prn=false) + integrator = KiteModels.init!(kps4; stiffness_factor=0.4, prn=false) simulate(kps4, integrator, 100) av_steps = simulate(kps4, integrator, STEPS-100) if Sys.isapple() diff --git a/test/test-steady-state-kps4.jl b/test/test-steady-state-kps4.jl index 7bff5a298..e7159fa49 100644 --- a/test/test-steady-state-kps4.jl +++ b/test/test-steady-state-kps4.jl @@ -35,7 +35,7 @@ end init_392(kps4_local, l) KiteModels.set_depower_steering!(kps4_local, kps4_local.set.depower_offset/100.0, 0.0) try - _, res2 = find_steady_state!(kps4_local; delta=0.001, stiffness_factor=0.07, prn=false) + _, res2 = find_steady_state!(kps4_local; stiffness_factor=0.07, prn=false) @test sum(res2) ≈ -9.81*(set.segments+ KiteModels.KITE_PARTICLES) # velocity and acceleration must be near zero pre_tension = KiteModels.calc_pre_tension(kps4_local) @test pre_tension > 1.0001 diff --git a/test/test-steady-state-robustness.jl b/test/test-steady-state-robustness.jl index d85b6cf22..83ae1c53e 100644 --- a/test/test-steady-state-robustness.jl +++ b/test/test-steady-state-robustness.jl @@ -8,20 +8,21 @@ end using Test using KiteModels, KitePodModels -# Cases for which a single nlsolve attempt in find_steady_state! does not converge -# (non-finite result or a residual norm far above the tolerance), so that the retries -# without autoscaling and the continuation in stiffness_factor and delta are needed. -# Fields: v_wind [m/s], elevation [°], tether length [m], stiffness_factor, delta +# Cases for which a single nlsolve attempt in find_steady_state! did not converge when the +# initial state was perturbed by delta, so that the retries without autoscaling and the +# continuation in stiffness_factor were needed. The elevations are high enough for a steady +# state with attached flow: below about 60°, the kite can only be balanced in deep stall. +# Fields: v_wind [m/s], elevation [°], tether length [m], stiffness_factor set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) const ROBUSTNESS_CASES = [ - (12.0, 50.0, 392.0, 0.5, 0.001), - (9.51, 50.0, 150.0, 0.1, 0.005), - (9.51, 70.0, 150.0, 0.5, 0.005), - (6.0, 50.0, 150.0, 0.04, 0.005), - (12.0, 70.0, 392.0, 0.1, 0.005), - (6.0, 70.0, 392.0, 0.1, 0.001), + (12.0, 60.0, 392.0, 0.5), + (9.51, 65.0, 150.0, 0.1), + (9.51, 70.0, 150.0, 0.5), + (6.0, 65.0, 150.0, 0.04), + (12.0, 70.0, 392.0, 0.1), + (6.0, 70.0, 392.0, 0.1), ] function robustness_kps4(v_wind, elevation, l_tether) @@ -34,11 +35,11 @@ function robustness_kps4(v_wind, elevation, l_tether) end @testset "find_steady_state! robustness" begin - for (v_wind, elevation, l_tether, stiffness_factor, delta) in ROBUSTNESS_CASES - @testset "v_wind=$v_wind, elevation=$elevation, l_tether=$l_tether, stiffness_factor=$stiffness_factor, delta=$delta" begin + for (v_wind, elevation, l_tether, stiffness_factor) in ROBUSTNESS_CASES + @testset "v_wind=$v_wind, elevation=$elevation, l_tether=$l_tether, stiffness_factor=$stiffness_factor" begin kps4 = robustness_kps4(v_wind, elevation, l_tether) # a warning is logged if no usable steady state was found - y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps4; delta, stiffness_factor) + y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps4; stiffness_factor) @test all(isfinite, y0) @test all(isfinite, yd0) @test kps4.stiffness_factor == stiffness_factor @@ -51,7 +52,7 @@ end @test rad2deg(calc_elevation(kps4)) ≈ elevation atol=1e-6 kps4 = robustness_kps4(v_wind, elevation, l_tether) - integrator = @test_logs min_level=Base.CoreLogging.Warn KiteModels.init!(kps4; delta, stiffness_factor) + integrator = @test_logs min_level=Base.CoreLogging.Warn KiteModels.init!(kps4; stiffness_factor) @test !isnothing(integrator) end end diff --git a/test/test-trim.jl b/test/test-trim.jl index 704d0c059..fc4386915 100644 --- a/test/test-trim.jl +++ b/test/test-trim.jl @@ -12,12 +12,12 @@ set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) # Find the steady state, starting at the elevation el0 [°]. Returns the elevation [°] and the # largest force residual (the residual of the accelerations of all particles). -function steady_state(Model, el0; trim, delta) +function steady_state(Model, el0; trim) set = deepcopy(load_settings("system.yaml")) set.elevation = el0 s = Model(KCU(set)) KiteModels.clear!(s) - y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(s; trim, delta) + y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(s; trim) res = zeros(length(y0)) residual!(res, yd0, y0, s, 0.0) if trim @@ -31,20 +31,19 @@ end @testset "find_steady_state! with trim=true, $Model" for Model in (KPS3, KPS4) el_set = load_settings("system.yaml").elevation # without trim the elevation is prescribed, so the kite is not in equilibrium - el, max_res = steady_state(Model, el_set; trim=false, delta=0.0) + el, max_res = steady_state(Model, el_set; trim=false) @test el ≈ el_set atol=1e-6 @test max_res > 1.0 # with trim all forces are balanced at the natural elevation of the kite - el_trim, max_res = steady_state(Model, el_set; trim=true, delta=0.0) + el_trim, max_res = steady_state(Model, el_set; trim=true) @test max_res < 1e-4 @test abs(el_trim - el_set) > 0.5 @test 50.0 < el_trim < 80.0 - # the trimmed elevation does not depend on the start value or on delta; from 30° the solver + # the trimmed elevation does not depend on the start value; from 30° the solver # finds a deep stall trim first, and for KPS4 it diverges from 55°, so both need restarts for el0 in (30.0, 55.0, 85.0) - @test steady_state(Model, el0; trim=true, delta=0.0)[1] ≈ el_trim atol=1e-3 + @test steady_state(Model, el0; trim=true)[1] ≈ el_trim atol=1e-3 end - @test steady_state(Model, el_set; trim=true, delta=0.001)[1] ≈ el_trim atol=0.05 end @testset "init! with trim=true, $Model" for Model in (KPS3, KPS4) @@ -82,7 +81,7 @@ end # to carry the kite, so no trim with attached flow exists profile == "Maasvlakte" && v_wind == 4.0 && continue s = Model(KCU(profile_settings(profile, v_wind))) - integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true, stiffness_factor=1.0, delta=0.0) + integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true, stiffness_factor=1.0) el0 = calc_elevation(s) push!(elevations, el0) # the kite is balanced for the wind of the profile at its height, with attached flow @@ -101,4 +100,30 @@ end # the kite flies higher at higher wind speeds @test issorted(elevations) end + +@testset "deprecated keyword delta, $Model" for Model in (KPS3, KPS4) + # delta is still accepted, but ignored: the steady state is the same as without it + s = Model(KCU(deepcopy(load_settings("system.yaml")))) + KiteModels.clear!(s) + y_ref, _ = find_steady_state!(s) + s = Model(KCU(deepcopy(load_settings("system.yaml")))) + KiteModels.clear!(s) + y, _ = find_steady_state!(s; delta=0.005) + @test y ≈ y_ref + @test !isnothing(init!(Model(KCU(deepcopy(load_settings("system.yaml")))); delta=0.005)) +end + +@testset "deep stall at a low prescribed elevation, $Model" for Model in (KPS3, KPS4) + # 50° is far below the natural elevation: without trim, the kite can only be balanced in deep + # stall, which is reported; with trim, the natural elevation is found + set = deepcopy(load_settings("system.yaml")) + set.elevation = 50.0 + s = Model(KCU(set)) + KiteModels.clear!(s) + @test_logs (:warn, r"deep stall at the prescribed elevation") find_steady_state!(s) + s = Model(KCU(deepcopy(set))) + KiteModels.clear!(s) + @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(s; trim=true) + @test s.alpha_2 < s.set.alpha_cl[argmax(s.set.cl_list)] +end nothing diff --git a/test/test-update-sys-state.jl b/test/test-update-sys-state.jl index cb609534c..71c746084 100644 --- a/test/test-update-sys-state.jl +++ b/test/test-update-sys-state.jl @@ -17,7 +17,7 @@ wrap(delta) = atan(sin(delta), cos(delta)) set = deepcopy(load_settings("system.yaml")) set.upwind_dir = rad2deg(-pi/2) kps4::KPS4 = KPS4(deepcopy(set)) -integrator = init!(kps4; stiffness_factor=0.5, delta=0.005, prn=false) +integrator = init!(kps4; stiffness_factor=0.5, prn=false) DT = 1/set.sample_freq # `next_step!` stores the time of the *previous* step in `t_0`, so two steps are needed # before `s.t_0 > 0` and the rate calculation sees a non-zero time difference @@ -124,7 +124,7 @@ end @testset "update_sys_state! KPS3" begin kps3 = KPS3(deepcopy(set)) - integrator3 = init!(kps3; stiffness_factor=0.5, delta=0.0001, prn=false) + integrator3 = init!(kps3; stiffness_factor=0.5, prn=false) next_step!(kps3, integrator3; set_speed=0, dt=DT) ss = SysState(kps3) @test length(ss.X) == length(kps3.pos) diff --git a/test/test_for_precompile.jl b/test/test_for_precompile.jl index f1590ea1a..e808b9f76 100644 --- a/test/test_for_precompile.jl +++ b/test/test_for_precompile.jl @@ -99,7 +99,7 @@ if "KiteViewers" ∈ keys(Pkg.project().dependencies) iter / steps end - integrator = KiteModels.init!(kps4, delta=0.001, stiffness_factor=0.1, prn=STATISTIC) + integrator = KiteModels.init!(kps4, stiffness_factor=0.1, prn=STATISTIC) println("\nStarting simulation...") simulate(integrator, STEPS) From 46f696eebdd32467c0669dc21e41c3983502d293 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 09:41:50 +0200 Subject: [PATCH 05/21] Fix the initial reel-out velocities of the KPS4 particles The velocities were multiplied by the length of the tether segments, so with trim=true the kite was placed far below its natural elevation. Use trim=true in test_init_4p.jl and print v_reel_out. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 10 ++++++++++ examples/test_init_4p.jl | 3 ++- src/KPS4.jl | 4 ++-- 3 files changed, 14 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a629cd5db..dbbda68d5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,16 @@ SPDX-License-Identifier: MIT 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`. ### KiteModels v0.11.19 2026-10-10 #### Fixed diff --git a/examples/test_init_4p.jl b/examples/test_init_4p.jl index 844454b72..32eeb772b 100644 --- a/examples/test_init_4p.jl +++ b/examples/test_init_4p.jl @@ -35,7 +35,7 @@ logger = Logger(set.segments + 5, STEPS) kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) -integrator = KiteModels.init!(kps4; stiffness_factor=0.1, prn=STATISTIC) +integrator = KiteModels.init!(kps4; stiffness_factor=1, prn=STATISTIC, trim=true) for _ in 1:80 next_step!(kps4, integrator; set_speed=kps4.set.v_reel_out, upwind_dir=deg2rad(UPWIND_DIR), dt) end @@ -57,6 +57,7 @@ println("--> orient_quat (viewer): roll: ", roll, " pitch: ", pitch, " yaw: ", println("x:", kps4.x) # from trailing edge to leading edge in ENU reference frame println("y:", kps4.y) # to the right looking in flight direction println("z:", kps4.z) # down +println("v_reel_out: ", round(kps4.v_reel_out, digits = 3), " m/s") azimuth = calc_azimuth(kps4) println("azimuth: ", round(rad2deg(azimuth), digits = 2), "°") println("azimuth_north: ", round(rad2deg(KiteUtils.azimuth_north(pos_kite(kps4))), digits = 2), "°") diff --git a/src/KPS4.jl b/src/KPS4.jl index e24be9eb6..7e4fb545c 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -1010,7 +1010,7 @@ function init_pos_vel_acc(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES)+1 s.pos[i] .= pos[i] end for i in 2:s.set.segments+1 - vel[i] .+= (pos[i+1] - pos[i]) * (s.set.v_reel_outs[1]*(i-1)/s.set.segments) + vel[i] .+= normalize(pos[i] - pos[i-1]) * (s.set.v_reel_outs[1]*(i-1)/s.set.segments) end # the velocity vector of the kite particles is the same as the velocity of the last tether point for i in s.set.segments+2:s.set.segments+KITE_PARTICLES+1 @@ -1021,7 +1021,7 @@ end function set_initial_velocity!(s::KPS4) for i in 2:s.set.segments+1 - s.vel[i] .= (s.pos[i+1] - s.pos[i]) * (s.set.v_reel_outs[1]*(i-1)/s.set.segments) + s.vel[i] .= normalize(s.pos[i] - s.pos[i-1]) * (s.set.v_reel_outs[1]*(i-1)/s.set.segments) end # the velocity vector of the kite particles is the same as the velocity of the last tether point for i in s.set.segments+2:s.set.segments+KITE_PARTICLES+1 From fc8fa927775ced73546a6cf583c6ebd3843cafaf Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 10:30:42 +0200 Subject: [PATCH 06/21] Deprecate the keyword stiffness_factor of init! and find_steady_state! The steady state is now always solved for the full tether stiffness, so the initial state is an equilibrium of the simulation. A reduced stiffness is only used internally as start value: if the direct solve fails, the schedules 0.1 -> 1, 0.1 -> 0.3 -> 1 and 0.3 -> 1 are tried, then the previous continuation. With trim=true, the KPS4 search starts at 80 degrees, from where the solver converges. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 26 ++++++++--- docs/src/examples_4p.md | 3 -- examples/bench.jl | 2 +- examples/bench_4p.jl | 2 +- examples/calc_spectrum.jl | 2 +- examples/compare_kps3_kps4.jl | 4 +- examples/plot_alpha2.jl | 2 +- examples/plot_cl_cd.jl | 2 +- examples/plot_parking_test.jl | 2 +- examples/plot_pitch_stability.jl | 2 +- examples/plot_side_cl.jl | 2 +- examples/reel_out_1p.jl | 2 +- examples/reel_out_4p.jl | 2 +- examples/reel_out_4p_torque_control.jl | 2 +- examples/simulate_simple.jl | 2 +- examples/simulate_steering.jl | 2 +- examples/steering_test_1p.jl | 2 +- examples/steering_test_4p.jl | 2 +- examples/test_init_1p.jl | 2 +- examples/test_init_4p.jl | 2 +- examples/test_steady_state.jl | 2 +- src/KPS3.jl | 15 +++++-- src/KPS4.jl | 62 ++++++++++++++++++-------- src/KiteModels.jl | 28 ++++++++---- src/precompile.jl | 4 +- test/plot_kps4.jl | 2 +- test/test-interface.jl | 4 +- test/test-kps4.jl | 9 ++-- test/test-simulate_4p.jl | 2 +- test/test-steady-state-kps4.jl | 2 +- test/test-steady-state-robustness.jl | 29 ++++++------ test/test-trim.jl | 24 ++++++++-- test/test-update-sys-state.jl | 4 +- test/test-winch-state.jl | 8 ++-- test/test_for_precompile.jl | 2 +- 35 files changed, 170 insertions(+), 94 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index dbbda68d5..5421dac5d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,13 +13,13 @@ SPDX-License-Identifier: MIT 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 5°. The default `trim=false` keeps the previous behaviour. -- `init!(s; trim=true)` passes `trim` to `find_steady_state!`. The trim depends on the stiffness - factor, which `next_step!` raises to 1.0, so use `init!(s; trim=true, stiffness_factor=1.0)` for an - initial state that is an equilibrium of the simulation; then the elevation stays constant. +- `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. - `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 60°), because then the kite - can only be balanced in deep stall; before, such a state was returned silently. + 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 @@ -28,6 +28,22 @@ SPDX-License-Identifier: MIT 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 diff --git a/docs/src/examples_4p.md b/docs/src/examples_4p.md index cded24615..dcd683937 100644 --- a/docs/src/examples_4p.md +++ b/docs/src/examples_4p.md @@ -97,9 +97,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 diff --git a/examples/bench.jl b/examples/bench.jl index 971971792..6edf18f9e 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -39,7 +39,7 @@ function simulate(integrator, steps) iter / steps end -integrator = KiteModels.init!(kps3; stiffness_factor = 0.1, prn = STATISTIC) +integrator = KiteModels.init!(kps3; prn = STATISTIC) println("\nStarting simulation...") simulate(integrator, 100) diff --git a/examples/bench_4p.jl b/examples/bench_4p.jl index 7a913b77f..01ea9a426 100644 --- a/examples/bench_4p.jl +++ b/examples/bench_4p.jl @@ -58,7 +58,7 @@ function simulate(integrator, steps, offset = 0) iter / steps end -integrator = init!(kps4, stiffness_factor = 0.1, prn = STATISTIC) +integrator = init!(kps4, prn = STATISTIC) println("\nStarting simulation...") simulate(integrator, 100, 100) diff --git a/examples/calc_spectrum.jl b/examples/calc_spectrum.jl index 6f0082db0..60d5f5b3f 100644 --- a/examples/calc_spectrum.jl +++ b/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; 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 diff --git a/examples/compare_kps3_kps4.jl b/examples/compare_kps3_kps4.jl index 470a77e9e..694d638b7 100644 --- a/examples/compare_kps3_kps4.jl +++ b/examples/compare_kps3_kps4.jl @@ -84,7 +84,7 @@ function simulate(s, integrator, steps, plot = false; fig = "") iter / steps end -integrator = init!(kps3, 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; stiffness_factor = 0.1, prn = STATISTIC) +integrator = init!(kps4; prn = STATISTIC) av_steps = simulate(kps4, integrator, STEPS, true; fig = "kps4") reactivate_host_app() diff --git a/examples/plot_alpha2.jl b/examples/plot_alpha2.jl index 095ca4dfd..bb729ebf8 100644 --- a/examples/plot_alpha2.jl +++ b/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) diff --git a/examples/plot_cl_cd.jl b/examples/plot_cl_cd.jl index 61d9b39a8..27710cc3f 100644 --- a/examples/plot_cl_cd.jl +++ b/examples/plot_cl_cd.jl @@ -92,7 +92,7 @@ for depower in DEPOWER kcu = KCU(set) kps4 = KPS4(kcu) - integrator = KiteModels.init!(kps4; stiffness_factor=0.05, prn=STATISTIC) + integrator = KiteModels.init!(kps4; prn=STATISTIC) if ! isnothing(integrator) try simulate(kps4, integrator, logger, STEPS) diff --git a/examples/plot_parking_test.jl b/examples/plot_parking_test.jl index 55dbeea8e..232010433 100644 --- a/examples/plot_parking_test.jl +++ b/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; stiffness_factor=0.04, prn=STATISTIC) + integrator = KiteModels.init!(kps4; prn=STATISTIC) if ! isnothing(integrator) try simulate(kps4, integrator, logger, STEPS) diff --git a/examples/plot_pitch_stability.jl b/examples/plot_pitch_stability.jl index 59ca0421f..f9ddac19f 100644 --- a/examples/plot_pitch_stability.jl +++ b/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; 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) diff --git a/examples/plot_side_cl.jl b/examples/plot_side_cl.jl index 34a718fb8..e96c513c8 100644 --- a/examples/plot_side_cl.jl +++ b/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; 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 diff --git a/examples/reel_out_1p.jl b/examples/reel_out_1p.jl index f19e71441..28f3df832 100644 --- a/examples/reel_out_1p.jl +++ b/examples/reel_out_1p.jl @@ -72,7 +72,7 @@ function simulate(integrator, steps, plot=false) iter / steps end -integrator = KiteModels.init!(kps3, stiffness_factor=0.1, prn=STATISTIC) +integrator = KiteModels.init!(kps3, prn=STATISTIC) kps3.sync_speed = 0.0 av_steps = if PLOT diff --git a/examples/reel_out_4p.jl b/examples/reel_out_4p.jl index ad301029b..696b63aa5 100644 --- a/examples/reel_out_4p.jl +++ b/examples/reel_out_4p.jl @@ -77,7 +77,7 @@ function simulate(integrator, steps, plot=false) iter / steps end -integrator = KiteModels.init!(kps4; stiffness_factor=0.1, prn=STATISTIC) +integrator = KiteModels.init!(kps4; prn=STATISTIC) kps4.sync_speed = 0.0 av_steps = if PLOT diff --git a/examples/reel_out_4p_torque_control.jl b/examples/reel_out_4p_torque_control.jl index e0996fd6b..dd7dfb8c6 100644 --- a/examples/reel_out_4p_torque_control.jl +++ b/examples/reel_out_4p_torque_control.jl @@ -93,7 +93,7 @@ function simulate(integrator, steps, plot=false) iter / steps end -integrator = KiteModels.init!(kps4; stiffness_factor=0.1, prn=STATISTIC) +integrator = KiteModels.init!(kps4; prn=STATISTIC) kps4.sync_speed = 0.0 av_steps = if PLOT diff --git a/examples/simulate_simple.jl b/examples/simulate_simple.jl index 8fb079fd0..f72ef475f 100644 --- a/examples/simulate_simple.jl +++ b/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; stiffness_factor=0.1, prn=STATISTIC, +integrator = KiteModels.init!(kps4; prn=STATISTIC, steady_state=false) av_steps = if PLOT diff --git a/examples/simulate_steering.jl b/examples/simulate_steering.jl index 8744a3bc1..94c9ea851 100644 --- a/examples/simulate_steering.jl +++ b/examples/simulate_steering.jl @@ -65,7 +65,7 @@ function simulate(integrator, steps, plot=false) iter / steps end -integrator = KiteModels.init!(kps3; stiffness_factor=0.04, prn=STATISTIC) +integrator = KiteModels.init!(kps3; prn=STATISTIC) av_steps = if PLOT simulate(integrator, STEPS, true) diff --git a/examples/steering_test_1p.jl b/examples/steering_test_1p.jl index 8ca941000..bc93ba040 100644 --- a/examples/steering_test_1p.jl +++ b/examples/steering_test_1p.jl @@ -101,7 +101,7 @@ function simulate(integrator, steps; plot=false) end end -integrator = KiteModels.init!(kps3; stiffness_factor=1, prn=STATISTIC) +integrator = KiteModels.init!(kps3; prn=STATISTIC) simulate(integrator, STEPS; plot=true) function delay(x, y, t_max = 10) diff --git a/examples/steering_test_4p.jl b/examples/steering_test_4p.jl index 12c3d5b9d..84da30061 100644 --- a/examples/steering_test_4p.jl +++ b/examples/steering_test_4p.jl @@ -105,7 +105,7 @@ function simulate(integrator, steps; plot=false) end -integrator = KiteModels.init!(kps4; stiffness_factor=0.1, prn=STATISTIC) +integrator = KiteModels.init!(kps4; prn=STATISTIC) simulate(integrator, STEPS; plot=true) function delay(x, y, t_max = 10) diff --git a/examples/test_init_1p.jl b/examples/test_init_1p.jl index 613d448b1..8bb7f38a5 100644 --- a/examples/test_init_1p.jl +++ b/examples/test_init_1p.jl @@ -33,7 +33,7 @@ logger = Logger(set.segments + 1, STEPS) kcu::KCU = KCU(set) kps3::KPS3 = KPS3(kcu) -_ = KiteModels.init!(kps3; stiffness_factor=0.1, prn=STATISTIC) +_ = KiteModels.init!(kps3; prn=STATISTIC) lift, drag = lift_drag(kps3) sys_state::SysState = KiteModels.SysState(kps3) log!(logger, sys_state) diff --git a/examples/test_init_4p.jl b/examples/test_init_4p.jl index 32eeb772b..260880907 100644 --- a/examples/test_init_4p.jl +++ b/examples/test_init_4p.jl @@ -35,7 +35,7 @@ logger = Logger(set.segments + 5, STEPS) kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) -integrator = KiteModels.init!(kps4; stiffness_factor=1, prn=STATISTIC, trim=true) +integrator = KiteModels.init!(kps4; prn=STATISTIC, trim=true) for _ in 1:80 next_step!(kps4, integrator; set_speed=kps4.set.v_reel_out, upwind_dir=deg2rad(UPWIND_DIR), dt) end diff --git a/examples/test_steady_state.jl b/examples/test_steady_state.jl index 5932fbf8a..ab7c44c7c 100644 --- a/examples/test_steady_state.jl +++ b/examples/test_steady_state.jl @@ -13,7 +13,7 @@ kps4::KPS4 = KPS4(kcu) clear!(kps4) KiteModels.set_depower_steering!(kps4, kps4.set.depower_offset/100.0, 0.0) -@time KiteModels.init!(kps4; stiffness_factor=0.5, prn=false) +@time KiteModels.init!(kps4; prn=false) println("\nlift, drag [N] : $(KiteModels.lift_drag(kps4))") println("\nSpring forces:") diff --git a/src/KPS3.jl b/src/KPS3.jl index cc8d736f9..75218a743 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -594,8 +594,12 @@ function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=no results.zero, results.f_converged end +# stiffness factor of the first solve of find_steady_state! (KPS3); its result is the start value +# for the solve with the full tether stiffness +const KPS3_START_STIFFNESS_FACTOR = 0.035 + """ - find_steady_state!(s::KPS3; prn=false, stiffness_factor=0.035, upwind_dir=-pi/2, trim=false) + find_steady_state!(s::KPS3; prn=false, upwind_dir=-pi/2, trim=false) Find an initial equilibrium, based on the initial parameters `l_tether`, elevation and `v_reel_out`. @@ -612,11 +616,14 @@ elevation), starting the search at `s.set.elevation`. Trimmed steady states with attack beyond the angle of maximum lift (deep stall) are rejected. If no valid trim is found from `s.set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 5°. -The keyword argument `delta` is deprecated and ignored. +The steady state is always solved for the full tether stiffness; a solution for a reduced +stiffness is only used as start value. The keyword arguments `delta` and `stiffness_factor` +are deprecated and ignored. """ -function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor=0.035, upwind_dir=-pi/2, +function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor=nothing, upwind_dir=-pi/2, trim=false) warn_delta_deprecated(delta, :find_steady_state!) + warn_stiffness_factor_deprecated(stiffness_factor, :find_steady_state!) el_set = s.set.elevation alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] # Find the steady state, starting at s.set.elevation. Returns the state, whether the solver @@ -627,7 +634,7 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) # 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 + s.stiffness_factor=KPS3_START_STIFFNESS_FACTOR p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false, trim) s.stiffness_factor=1.0 # if the solution with the reduced stiffness failed, start again from the straight tether diff --git a/src/KPS4.jl b/src/KPS4.jl index 7e4fb545c..492f51655 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -705,7 +705,7 @@ function turn(res, upwind_dir) end """ - find_steady_state!(s::KPS4; prn=false, stiffness_factor=0.035, upwind_dir=-pi/2, trim=false) + find_steady_state!(s::KPS4; prn=false, upwind_dir=-pi/2, trim=false) Find an initial equilibrium, based on the initial parameters `l_tether`, elevation and `v_reel_out`. @@ -721,33 +721,38 @@ attack beyond the angle of maximum lift; then a warning is logged. If `trim` is `true`, all forces including gravity are balanced, also the z components of the forces on C and D, and the elevation is solved for instead: the kite is placed at the elevation where its lift matches gravity (its natural elevation), starting the search at -`s.set.elevation`. +80° or at `s.set.elevation`, if it is higher: the solver converges from above the natural +elevation, but often not from far below it. Trimmed steady states with an angle of attack beyond the angle of maximum lift (deep -stall) are rejected. If no valid trim is found from `s.set.elevation`, the solver is -restarted at higher and then at lower elevations, in steps of 5°. +stall) are rejected. If no valid trim is found from the start elevation, the solver is +restarted at higher and then at lower elevations than `s.set.elevation`, in steps of 5°. -The nonlinear solver is sensitive to the starting point, so small floating-point +The steady state is always solved for the full tether stiffness (`s.stiffness_factor` is 1.0 +afterwards). The nonlinear solver is sensitive to the starting point, so small floating-point differences (e.g. between Julia versions) can decide whether it converges. If the first attempt fails, it is repeated without autoscaling and then by continuation: -a solution is first found for an easier `stiffness_factor` and then used as start value -while stepping towards the requested value. The initial state is not perturbed, unless -this fails; the returned steady state is never perturbed. +a solution is first found for a reduced tether stiffness and then used as start value +for the full stiffness. The initial state is not perturbed, unless this fails; the returned +steady state is never perturbed. -The keyword argument `delta` is deprecated and ignored. +The keyword arguments `delta` and `stiffness_factor` are deprecated and ignored. """ -function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor=0.035, upwind_dir=-pi/2, +function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor=nothing, upwind_dir=-pi/2, trim=false) warn_delta_deprecated(delta, :find_steady_state!) + warn_stiffness_factor_deprecated(stiffness_factor, :find_steady_state!) # solve without perturbation of the initial state, with a perturbation only as fallback function solve() local r, conv for delta in STEADY_STATE_DELTAS - r, conv = solve_steady_state!(s; prn, delta, stiffness_factor, trim) + r, conv = solve_steady_state!(s; prn, delta, trim) conv && break end r, conv end el_set = s.set.elevation + # elevation of the first attempt; the offsets of its result are relative to it + el_start = trim ? max(el_set, TRIM_START_ELEVATION) : el_set alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) # the state for the offsets X, which are relative to the tether at s.set.elevation @@ -763,6 +768,7 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= y_, yd_ end try + s.set.elevation = el_start results, converged = solve() converged && ((y00, yd00) = steady_state(results.zero)) stalled = trim && converged && s.alpha_2 > alpha_stall @@ -770,7 +776,7 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= # Restart from other elevations: the solver converges from above the natural # elevation, but diverges or finds a deep stall trim from far below it. starts = [el_set .+ (5.0:5.0:85.0); el_set .- (5.0:5.0:85.0)] - for el in filter(el -> 5.0 <= el <= 89.0, starts) + for el in filter(el -> 5.0 <= el <= 89.0 && el != el_start, starts) s.set.elevation = el r, conv = solve() conv || continue @@ -783,14 +789,14 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= end end if stalled - s.set.elevation = el_set + s.set.elevation = el_start y00, yd00 = steady_state(results.zero) @warn "find_steady_state!: only a trimmed steady state in deep stall was found " * "(angle of attack $(round(s.alpha_2, digits=1))°)" end end if !converged - s.set.elevation = el_set + s.set.elevation = el_start @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations), residual_norm=$(results.residual_norm))" # Check if the solution contains finite values if !all(isfinite, results.zero) @@ -815,9 +821,19 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= end end -# Solve for the offsets of the steady state, see find_steady_state!. Returns the result of nlsolve -# and whether it is usable. -function solve_steady_state!(s::KPS4; prn=false, delta=0.001, stiffness_factor=0.035, trim=false) +# Stiffness factors for the continuation in solve_steady_state!: if the solve for the full tether +# stiffness fails, the solution is first found for the reduced stiffness factors of a schedule, +# each starting from the previous solution. 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) showed +# that the solve for the full stiffness converges directly in most cases and is the fastest +# option, starting at a reduced stiffness is needed mainly with trim, and the schedules fail in +# different cases, so they are tried in turn. +const STIFFNESS_SCHEDULES = ((0.1, 1.0), (0.1, 0.3, 1.0), (0.3, 1.0)) + +# Solve for the offsets of the steady state for the full tether stiffness, see +# find_steady_state!. Returns the result of nlsolve and whether it is usable. +function solve_steady_state!(s::KPS4; prn=false, delta=0.001, trim=false) + stiffness_factor = 1.0 s.stiffness_factor = stiffness_factor res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) delta_ref = Ref(delta) @@ -892,8 +908,18 @@ function solve_steady_state!(s::KPS4; prn=false, delta=0.001, stiffness_factor=0 if !usable(results) r = solve(stiffness_factor, delta, X00; autoscale=false) usable(r) || (r = nothing) + # continuation with a fixed schedule of stiffness factors + for schedule in STIFFNESS_SCHEDULES + isnothing(r) || break + r = try_solve(schedule[1], delta, X00) + for sf in schedule[2:end] + isnothing(r) && break + r = try_solve(sf, delta, r.zero) + end + !isnothing(r) && prn && println("find_steady_state!: converged with the stiffness schedule $schedule") + end if isnothing(r) - # continuation: solve an easier problem, then step towards the requested parameters + # continuation: solve an easier problem, then step towards the full stiffness starts = [(sf0, delta) for sf0 in (0.1, 0.2, 0.3, 0.05) if sf0 != stiffness_factor] if delta != 0.001 push!(starts, (stiffness_factor, 0.001)) diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 6c7dfa892..53ddac4d2 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -82,6 +82,11 @@ const AKM = AbstractKiteModel # start can help the nonlinear solver, but the returned steady state is never perturbed. const STEADY_STATE_DELTAS = (0.0, 0.001) +# With trim=true, the KPS4 steady state solver starts at this elevation [°], or at +# set.elevation if it is higher: the solver converges from above the natural elevation, +# but often diverges or finds a deep stall trim from far below it. +const TRIM_START_ELEVATION = 80.0 + function warn_delta_deprecated(delta, f) isnothing(delta) || Base.depwarn("The keyword argument `delta` of `$f` is deprecated and " * "ignored: the steady state is solved without perturbation (with a small perturbation " * @@ -89,6 +94,13 @@ function warn_delta_deprecated(delta, f) nothing end +function warn_stiffness_factor_deprecated(stiffness_factor, f) + isnothing(stiffness_factor) || Base.depwarn("The keyword argument `stiffness_factor` of `$f` is " * + "deprecated and ignored: the steady state is always solved for the full tether stiffness " * + "(with a reduced stiffness only as an internal start value).", f) + nothing +end + function __init__() if isdir(joinpath(pwd(), "data")) && isfile(joinpath(pwd(), "data", "system.yaml")) set_data_path(joinpath(pwd(), "data")) @@ -651,14 +663,13 @@ function init_winch_state!(s::AKM) end """ - init!(s::AKM; stiffness_factor=0.5, prn=false, steady_state=true, + init!(s::AKM; prn=false, steady_state=true, trim=false) -> Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} Initializes the integrator of the model (KPS3 and KPS4 only). Parameters: - s: an instance of an abstract kite model -- stiffness_factor: factor applied to the tether stiffness during initialization - prn: if set to true, print the detailed solver results - steady_state: if `true` (default), call `find_steady_state!` to solve for an equilibrium before constructing the integrator. If `false`, skip the steady-state solver and initialise @@ -670,26 +681,27 @@ Parameters: lift matches gravity, instead of at `s.set.elevation`; see `find_steady_state!`. Ignored if `steady_state` is `false`. -The keyword argument `delta` is deprecated and ignored. +The keyword arguments `delta` and `stiffness_factor` are deprecated and ignored: the initial +state is always an equilibrium for the full tether stiffness. Returns: An instance of an `ODEIntegrator` or `IDAIntegrator`, or `nothing` if initialization fails. """ -function init!(s::AKM; stiffness_factor=0.5, delta=nothing, prn=false, steady_state=true, +function init!(s::AKM; stiffness_factor=nothing, delta=nothing, prn=false, steady_state=true, trim=false)::Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.IDAIntegrator, Nothing} warn_delta_deprecated(delta, :init!) + warn_stiffness_factor_deprecated(stiffness_factor, :init!) clear!(s) init_winch_state!(s) upwind_dir = deg2rad(Float64(s.set.upwind_dir)) - s.stiffness_factor = stiffness_factor + s.stiffness_factor = 1.0 y0, yd0 = if steady_state try - KiteModels.find_steady_state!(s; stiffness_factor, upwind_dir, prn, trim) + KiteModels.find_steady_state!(s; upwind_dir, prn, trim) catch e if e isa AssertionError || (e isa ErrorException && startswith(e.msg, "find_steady_state!:")) - println("ERROR: Failure to find initial steady state in find_steady_state! function!\n"* - "Try to decrease the stiffness_factor of the init! function.") + println("ERROR: Failure to find initial steady state in find_steady_state! function!") return nothing else rethrow(e) diff --git a/src/precompile.jl b/src/precompile.jl index 7a1667a2f..2a4b2f32d 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -18,8 +18,8 @@ @compile_workload begin # all calls in this block will be precompiled, regardless of whether( # they belong to your package or not (on Julia 1.8 and higher) - integrator = KiteModels.init!(kps3_; stiffness_factor=0.035, prn=false) - integrator = KiteModels.init!(kps4_; stiffness_factor=0.035, prn=false) + integrator = KiteModels.init!(kps3_; prn=false) + integrator = KiteModels.init!(kps4_; prn=false) nothing end end \ No newline at end of file diff --git a/test/plot_kps4.jl b/test/plot_kps4.jl index 0227b3d3e..017ee3422 100644 --- a/test/plot_kps4.jl +++ b/test/plot_kps4.jl @@ -23,7 +23,7 @@ kps4.stiffness_factor = 0.04 y0, yd0 = KiteModels.init(kps4) -find_steady_state!(kps4; stiffness_factor=1, prn=false) +find_steady_state!(kps4; prn=false) println("kite distance: $(norm(kps4.pos[end-2]))") println(KiteModels.spring_forces(kps4)) diff --git a/test/test-interface.jl b/test/test-interface.jl index 0876fb1ee..5b12f4109 100644 --- a/test/test-interface.jl +++ b/test/test-interface.jl @@ -32,14 +32,14 @@ end @testset "KPS3 init! interface" begin kps4.set.upwind_dir = rad2deg(-pi/2) - integ = init!(kps3; stiffness_factor=0.5, prn=false) + integ = init!(kps3; prn=false) @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator @test kps3.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end @testset "KPS4 init! interface" begin kps4.set.upwind_dir = rad2deg(-pi/2) - integ = init!(kps4; stiffness_factor=0.5, prn=false) + integ = init!(kps4; prn=false) @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator @test kps4.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end diff --git a/test/test-kps4.jl b/test/test-kps4.jl index f02095e39..7a91eba6e 100644 --- a/test/test-kps4.jl +++ b/test/test-kps4.jl @@ -583,7 +583,7 @@ set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) height = sin(deg2rad(kps4.set.elevations[1])) * kps4.set.l_tethers[1] kps4.v_wind .= kps4.v_wind_gnd * calc_wind_factor(kps4.am, height) try - res1, res2 = find_steady_state!(kps4; stiffness_factor = 0.07, prn = false) + res1, res2 = find_steady_state!(kps4; prn = false) # TODO check why -9.81 appears in the residual @test sum(res2) ≈ -9.81*(kps4.set.segments + KiteModels.KITE_PARTICLES) # velocity and acceleration must be near zero pre_tension = KiteModels.calc_pre_tension(kps4) @@ -687,7 +687,7 @@ set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) # and the updraft shows up in v_wind_kite (and therefore in the logged system state) local kps4_ = KPS4(KCU(deepcopy(load_settings("system.yaml")))) kps4_.set.kcu_diameter = 0.0 - integrator = KiteModels.init!(kps4_; stiffness_factor=0.5, prn=false) + integrator = KiteModels.init!(kps4_; prn=false) dt = 1/kps4_.set.sample_freq for _ in 1:5 next_step!(kps4_, integrator; set_speed=0, v_wind_vert=1.5, dt) @@ -704,9 +704,10 @@ set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) kps4_ = KPS4(KCU(deepcopy(load_settings("system.yaml")))) kps4_.set.kcu_diameter = 0.0 kps4_.set.max_iter = 500 - integrator = KiteModels.init!(kps4_; stiffness_factor = 0.07, prn = false) + integrator = KiteModels.init!(kps4_; prn = false) kps4_.set.version = 2 - kps4_.stiffness_factor = 3 + # the steady state is solved for the full stiffness; increase it to break the bridle + kps4_.stiffness_factor = 20 @test maximum(spring_forces(kps4_; prn = false)) > 9000 end diff --git a/test/test-simulate_4p.jl b/test/test-simulate_4p.jl index 9cc3c9060..28e5256fa 100644 --- a/test/test-simulate_4p.jl +++ b/test/test-simulate_4p.jl @@ -44,7 +44,7 @@ end kps4.set.solver = "IDA" kps4.set.profile_law = Int(EXPLOG) kps4.set.alpha = 0.08163 - integrator = KiteModels.init!(kps4; stiffness_factor=0.4, prn=false) + integrator = KiteModels.init!(kps4; prn=false) simulate(kps4, integrator, 100) av_steps = simulate(kps4, integrator, STEPS-100) if Sys.isapple() diff --git a/test/test-steady-state-kps4.jl b/test/test-steady-state-kps4.jl index e7159fa49..23b6ba7e8 100644 --- a/test/test-steady-state-kps4.jl +++ b/test/test-steady-state-kps4.jl @@ -35,7 +35,7 @@ end init_392(kps4_local, l) KiteModels.set_depower_steering!(kps4_local, kps4_local.set.depower_offset/100.0, 0.0) try - _, res2 = find_steady_state!(kps4_local; stiffness_factor=0.07, prn=false) + _, res2 = find_steady_state!(kps4_local; prn=false) @test sum(res2) ≈ -9.81*(set.segments+ KiteModels.KITE_PARTICLES) # velocity and acceleration must be near zero pre_tension = KiteModels.calc_pre_tension(kps4_local) @test pre_tension > 1.0001 diff --git a/test/test-steady-state-robustness.jl b/test/test-steady-state-robustness.jl index 83ae1c53e..4e76cb1f2 100644 --- a/test/test-steady-state-robustness.jl +++ b/test/test-steady-state-robustness.jl @@ -10,19 +10,19 @@ using KiteModels, KitePodModels # Cases for which a single nlsolve attempt in find_steady_state! did not converge when the # initial state was perturbed by delta, so that the retries without autoscaling and the -# continuation in stiffness_factor were needed. The elevations are high enough for a steady -# state with attached flow: below about 60°, the kite can only be balanced in deep stall. -# Fields: v_wind [m/s], elevation [°], tether length [m], stiffness_factor +# continuation in the stiffness factor were needed. The elevations are high enough for a steady +# state with attached flow: below about 45°, the kite can only be balanced in deep stall. +# Fields: v_wind [m/s], elevation [°], tether length [m] set_data_path(joinpath(dirname(dirname(pathof(KiteModels)::String)), "data")) const ROBUSTNESS_CASES = [ - (12.0, 60.0, 392.0, 0.5), - (9.51, 65.0, 150.0, 0.1), - (9.51, 70.0, 150.0, 0.5), - (6.0, 65.0, 150.0, 0.04), - (12.0, 70.0, 392.0, 0.1), - (6.0, 70.0, 392.0, 0.1), + (12.0, 60.0, 392.0), + (9.51, 65.0, 150.0), + (9.51, 70.0, 150.0), + (6.0, 65.0, 150.0), + (12.0, 70.0, 392.0), + (6.0, 70.0, 392.0), ] function robustness_kps4(v_wind, elevation, l_tether) @@ -35,14 +35,15 @@ function robustness_kps4(v_wind, elevation, l_tether) end @testset "find_steady_state! robustness" begin - for (v_wind, elevation, l_tether, stiffness_factor) in ROBUSTNESS_CASES - @testset "v_wind=$v_wind, elevation=$elevation, l_tether=$l_tether, stiffness_factor=$stiffness_factor" begin + for (v_wind, elevation, l_tether) in ROBUSTNESS_CASES + @testset "v_wind=$v_wind, elevation=$elevation, l_tether=$l_tether" begin kps4 = robustness_kps4(v_wind, elevation, l_tether) # a warning is logged if no usable steady state was found - y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps4; stiffness_factor) + y0, yd0 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps4) @test all(isfinite, y0) @test all(isfinite, yd0) - @test kps4.stiffness_factor == stiffness_factor + # the steady state is solved for the full tether stiffness + @test kps4.stiffness_factor == 1.0 @test unstretched_length(kps4) ≈ l_tether # in a steady state the tether must be under tension; a result where the solver # stalled far away from the solution has a (nearly) unstretched tether @@ -52,7 +53,7 @@ end @test rad2deg(calc_elevation(kps4)) ≈ elevation atol=1e-6 kps4 = robustness_kps4(v_wind, elevation, l_tether) - integrator = @test_logs min_level=Base.CoreLogging.Warn KiteModels.init!(kps4; stiffness_factor) + integrator = @test_logs min_level=Base.CoreLogging.Warn KiteModels.init!(kps4) @test !isnothing(integrator) end end diff --git a/test/test-trim.jl b/test/test-trim.jl index fc4386915..25d578225 100644 --- a/test/test-trim.jl +++ b/test/test-trim.jl @@ -50,7 +50,7 @@ end # at full stiffness the trimmed initial state is an equilibrium of the simulation, # so the elevation stays constant s = Model(KCU(deepcopy(load_settings("system.yaml")))) - integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true, stiffness_factor=1.0) + integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true) el0 = calc_elevation(s) @test abs(rad2deg(el0) - s.set.elevation) > 0.5 dt = 1/s.set.sample_freq @@ -81,7 +81,7 @@ end # to carry the kite, so no trim with attached flow exists profile == "Maasvlakte" && v_wind == 4.0 && continue s = Model(KCU(profile_settings(profile, v_wind))) - integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true, stiffness_factor=1.0) + integrator = @test_logs min_level=Base.CoreLogging.Warn init!(s; trim=true) el0 = calc_elevation(s) push!(elevations, el0) # the kite is balanced for the wind of the profile at its height, with attached flow @@ -113,11 +113,27 @@ end @test !isnothing(init!(Model(KCU(deepcopy(load_settings("system.yaml")))); delta=0.005)) end +@testset "deprecated keyword stiffness_factor, $Model" for Model in (KPS3, KPS4) + # stiffness_factor is still accepted, but ignored: the steady state is solved for the full + # tether stiffness + s = Model(KCU(deepcopy(load_settings("system.yaml")))) + KiteModels.clear!(s) + y_ref, _ = find_steady_state!(s) + s = Model(KCU(deepcopy(load_settings("system.yaml")))) + KiteModels.clear!(s) + y, _ = find_steady_state!(s; stiffness_factor=0.1) + @test y ≈ y_ref + @test s.stiffness_factor == 1.0 + s = Model(KCU(deepcopy(load_settings("system.yaml")))) + @test !isnothing(init!(s; stiffness_factor=0.1)) + @test s.stiffness_factor == 1.0 +end + @testset "deep stall at a low prescribed elevation, $Model" for Model in (KPS3, KPS4) - # 50° is far below the natural elevation: without trim, the kite can only be balanced in deep + # 40° is far below the natural elevation: without trim, the kite can only be balanced in deep # stall, which is reported; with trim, the natural elevation is found set = deepcopy(load_settings("system.yaml")) - set.elevation = 50.0 + set.elevation = 40.0 s = Model(KCU(set)) KiteModels.clear!(s) @test_logs (:warn, r"deep stall at the prescribed elevation") find_steady_state!(s) diff --git a/test/test-update-sys-state.jl b/test/test-update-sys-state.jl index 71c746084..77b0e2f54 100644 --- a/test/test-update-sys-state.jl +++ b/test/test-update-sys-state.jl @@ -17,7 +17,7 @@ wrap(delta) = atan(sin(delta), cos(delta)) set = deepcopy(load_settings("system.yaml")) set.upwind_dir = rad2deg(-pi/2) kps4::KPS4 = KPS4(deepcopy(set)) -integrator = init!(kps4; stiffness_factor=0.5, prn=false) +integrator = init!(kps4; prn=false) DT = 1/set.sample_freq # `next_step!` stores the time of the *previous* step in `t_0`, so two steps are needed # before `s.t_0 > 0` and the rate calculation sees a non-zero time difference @@ -124,7 +124,7 @@ end @testset "update_sys_state! KPS3" begin kps3 = KPS3(deepcopy(set)) - integrator3 = init!(kps3; stiffness_factor=0.5, prn=false) + integrator3 = init!(kps3; prn=false) next_step!(kps3, integrator3; set_speed=0, dt=DT) ss = SysState(kps3) @test length(ss.X) == length(kps3.pos) diff --git a/test/test-winch-state.jl b/test/test-winch-state.jl index 504b53e49..71c25f555 100644 --- a/test/test-winch-state.jl +++ b/test/test-winch-state.jl @@ -15,7 +15,7 @@ winch_state(s) = (s.sync_speed, s.wm.brake, s.wm.last_set_speed) @testset "init! resets the winch state, $(nameof(Model))" for Model in (KPS3, KPS4) s = Model(KCU(load_settings("system.yaml"))) - integrator = KiteModels.init!(s; stiffness_factor=0.04) + integrator = KiteModels.init!(s) @test !isnothing(integrator) initial = winch_state(s) @test initial[1] == s.set.v_reel_out @@ -28,14 +28,14 @@ winch_state(s) = (s.sync_speed, s.wm.brake, s.wm.last_set_speed) # set the brake to a state inside the hysteresis band that differs from the initial one s.sync_speed = s.wm.v_min s.wm.brake = !initial[2] - integrator = KiteModels.init!(s; stiffness_factor=0.04) + integrator = KiteModels.init!(s) @test !isnothing(integrator) @test winch_state(s) == initial end @testset "residual! does not change the winch state, $(nameof(Model))" for Model in (KPS3, KPS4) s = Model(KCU(load_settings("system.yaml"))) - integrator = KiteModels.init!(s; stiffness_factor=0.04) + integrator = KiteModels.init!(s) @test !isnothing(integrator) for set_speed in (0.0, 0.5s.wm.v_min, 10.0, -10.0) # a set speed that would switch the brake or move the rate limited set speed @@ -51,7 +51,7 @@ end @testset "next_step! advances the winch state once, $(nameof(Model))" for Model in (KPS3, KPS4) s = Model(KCU(load_settings("system.yaml"))) - integrator = KiteModels.init!(s; stiffness_factor=0.04) + integrator = KiteModels.init!(s) @test !isnothing(integrator) @test s.wm.brake # the initial set speed is zero step = s.set.max_acc / s.set.sample_freq diff --git a/test/test_for_precompile.jl b/test/test_for_precompile.jl index e808b9f76..58c5d9dcf 100644 --- a/test/test_for_precompile.jl +++ b/test/test_for_precompile.jl @@ -99,7 +99,7 @@ if "KiteViewers" ∈ keys(Pkg.project().dependencies) iter / steps end - integrator = KiteModels.init!(kps4, stiffness_factor=0.1, prn=STATISTIC) + integrator = KiteModels.init!(kps4, prn=STATISTIC) println("\nStarting simulation...") simulate(integrator, STEPS) From 1ed2c7c34ca7bb68fb1c7688232a5cec7fedc93b Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 10:42:51 +0200 Subject: [PATCH 07/21] Use v_reel_out in the KPS3 steady state The particles were at rest and the initial reel-out speed was zero, so the initial elevation did not depend on v_reel_out. Now the particle velocities grow linearly along the tether to v_reel_out at the kite, and the initial state contains the reel-out speed, as for KPS4. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 9 +++++++++ examples/test_init_1p.jl | 4 +++- src/KPS3.jl | 9 +++++++-- test/test-trim.jl | 15 +++++++++++++++ 4 files changed, 34 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5421dac5d..b8e48f029 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -63,6 +63,15 @@ SPDX-License-Identifier: MIT 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.11.19 2026-10-10 #### Fixed diff --git a/examples/test_init_1p.jl b/examples/test_init_1p.jl index 8bb7f38a5..16b579fea 100644 --- a/examples/test_init_1p.jl +++ b/examples/test_init_1p.jl @@ -23,6 +23,7 @@ V_WIND = 14.5 # the following values can be changed to match your interest dt::Float64 = 0.05 set.solver="DFBDF" # IDA or DFBDF +set.v_reel_out = 1.0 # initial reel-out speed [m/s] STEPS = 1 const PLOT = true STATISTIC = false @@ -33,7 +34,7 @@ logger = Logger(set.segments + 1, STEPS) kcu::KCU = KCU(set) kps3::KPS3 = KPS3(kcu) -_ = KiteModels.init!(kps3; prn=STATISTIC) +_ = KiteModels.init!(kps3; prn=STATISTIC, trim=true) lift, drag = lift_drag(kps3) sys_state::SysState = KiteModels.SysState(kps3) log!(logger, sys_state) @@ -53,6 +54,7 @@ x, y, z = kite_ref_frame(kps3) println("x:", x) # from trailing edge to leading edge in ENU reference frame println("y:", y) # to the right looking in flight direction println("z:", z) # down +println("v_reel_out: ", round(kps4.v_reel_out, digits = 3), " m/s") azimuth = calc_azimuth(kps3) println("azimuth: ", round(rad2deg(azimuth), digits = 2), "°") diff --git a/src/KPS3.jl b/src/KPS3.jl index 75218a743..df3482342 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -479,6 +479,11 @@ function init_inner(s::KPS3, X=zeros(2 * s.set.segments); delta=0.0, upwind_dir= vel[i+1] .= SVec3(DELTA, DELTA, DELTA) acc[i+1] .= SVec3(DELTA, DELTA, DELTA) end + # reeling out: 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 + for i in 1:s.set.segments + vel[i+1] .+= normalize(pos[i+1] - pos[i]) * (s.set.v_reel_out * i / s.set.segments) + end for i in eachindex(pos) @inbounds s.pos[i] .= pos[i] end @@ -502,10 +507,10 @@ end # same as above, but returns a tuple of two one dimensional arrays function init(s::KPS3, X=zeros(2 * (s.set.segments)); delta = 0.0, upwind_dir=nothing) res1_, res2_ = init_inner(s, X; delta, upwind_dir=something(upwind_dir, -pi/2)) - res1 = vcat(reduce(vcat, res1_), [s.l_tether, 0]) + res1 = vcat(reduce(vcat, res1_), [s.l_tether, s.set.v_reel_out]) res2 = vcat(reduce(vcat, res2_), [0, 0]) if !isnothing(upwind_dir) - n = 6 * s.set.segments # only rotate the 3D vector part; the last 2 scalars (l_tether, 0) must not be touched + n = 6 * s.set.segments # only rotate the 3D vector part; the last 2 scalars (l_tether, v_reel_out) must not be touched res1 = vcat(turn(res1[1:n], upwind_dir), res1[n+1:end]) res2 = vcat(turn(res2[1:n], upwind_dir), res2[n+1:end]) # Keep s.pos consistent with the rotated state vectors so callers that diff --git a/test/test-trim.jl b/test/test-trim.jl index 25d578225..dc8dd55c8 100644 --- a/test/test-trim.jl +++ b/test/test-trim.jl @@ -142,4 +142,19 @@ end @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(s; trim=true) @test s.alpha_2 < s.set.alpha_cl[argmax(s.set.cl_list)] end + +@testset "reel-out speed in the steady state, $Model" for Model in (KPS3, KPS4) + # reeling out reduces the apparent wind speed, so the natural elevation decreases, and the + # initial state contains the reel-out speed + elevations = map((0.0, 2.0)) do v_reel_out + set = deepcopy(load_settings("system.yaml")) + set.v_reel_out = v_reel_out + set.v_reel_outs[1] = v_reel_out + s = Model(KCU(set)) + integrator = init!(s; trim=true) + @test integrator.u[end] ≈ v_reel_out + rad2deg(calc_elevation(s)) + end + @test elevations[2] < elevations[1] - 1.0 +end nothing From 20a8588dd63ee942ecc5aa7f8e1801b9c04c09ec Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 10:43:53 +0200 Subject: [PATCH 08/21] Bugfix --- examples/test_init_1p.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/test_init_1p.jl b/examples/test_init_1p.jl index 16b579fea..93af6d836 100644 --- a/examples/test_init_1p.jl +++ b/examples/test_init_1p.jl @@ -54,7 +54,7 @@ x, y, z = kite_ref_frame(kps3) println("x:", x) # from trailing edge to leading edge in ENU reference frame println("y:", y) # to the right looking in flight direction println("z:", z) # down -println("v_reel_out: ", round(kps4.v_reel_out, digits = 3), " m/s") +println("v_reel_out: ", round(kps3.v_reel_out, digits = 3), " m/s") azimuth = calc_azimuth(kps3) println("azimuth: ", round(rad2deg(azimuth), digits = 2), "°") From 59fdb75eb6422a56ab4ada78e0b1198faab8af95 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 11:19:03 +0200 Subject: [PATCH 09/21] Use the same x limits for all frames in compare_kps3_kps4.jl The first frame used xlim = (35, 55) and the later frames (35, 75), so the kite seemed to jump after the first frame. Both now use (40, 55). Also use MakieControlPlots from its main branch in the examples. Co-Authored-By: Claude Opus 5.5 (1M context) --- examples/Project.toml | 5 +++-- examples/compare_kps3_kps4.jl | 8 ++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/Project.toml b/examples/Project.toml index 7a42444bf..a641c1bf4 100644 --- a/examples/Project.toml +++ b/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" @@ -19,16 +19,17 @@ Timers = "21f18d07-b854-4dab-86f0-c15a3821819a" [sources] KiteModels = {path = ".."} +MakieControlPlots = {rev = "main", url = "https://github.com/OpenSourceAWE/MakieControlPlots.jl.git"} [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" diff --git a/examples/compare_kps3_kps4.jl b/examples/compare_kps3_kps4.jl index 694d638b7..38c13320d 100644 --- a/examples/compare_kps3_kps4.jl +++ b/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,12 +72,12 @@ 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 From 5b3f8e51e5a409a1e106c7344acf47c3713a88f1 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 12:13:29 +0200 Subject: [PATCH 10/21] Throw an error if no trim with attached flow is found With trim=true, find_steady_state! (KPS3 and KPS4) returned a deep stall or non-converged state if all restarts failed. Now it throws an error, so init! returns nothing instead of initializing from an invalid trim. The restarts use 10 degree steps and a cheaper solve, so that failing takes 7 s (KPS3) and 11 s (KPS4) instead of 25 s and 78 s. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 8 ++++++-- src/KPS3.jl | 23 ++++++++++++----------- src/KPS4.jl | 30 ++++++++++++++++-------------- src/KiteModels.jl | 9 +++++++-- test/test-trim.jl | 11 +++++++++++ 5 files changed, 52 insertions(+), 29 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b8e48f029..6038a22ef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,9 +12,13 @@ SPDX-License-Identifier: MIT 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 5°. The default `trim=false` keeps the previous behaviour. + 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. + 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 diff --git a/src/KPS3.jl b/src/KPS3.jl index df3482342..88ab34ee2 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -619,7 +619,8 @@ If `trim` is `true`, the vertical force balance of the kite is solved for too an is free: the kite is placed at the elevation where its lift matches gravity (its natural elevation), starting the search at `s.set.elevation`. Trimmed steady states with an angle of attack beyond the angle of maximum lift (deep stall) are rejected. If no valid trim is found from -`s.set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 5°. +`s.set.elevation`, the solver is restarted at higher and then at lower elevations, in steps of 10°. +If this fails too, an error is thrown, and `init!` returns `nothing`. The steady state is always solved for the full tether stiffness; a solution for a reduced stiffness is only used as start value. The keyword arguments `delta` and `stiffness_factor` @@ -632,10 +633,11 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= el_set = s.set.elevation alpha_stall = s.set.alpha_cl[argmax(s.set.cl_list)] # angle of attack of maximum lift [°] # Find the steady state, starting at s.set.elevation. Returns the state, whether the solver - # converged and whether the kite is in deep stall. - function solve(warn) + # converged and whether the kite is in deep stall. The restarts at other elevations use only + # the unperturbed start value, so that failing to find a trim is fast. + function solve(warn; restart=false) local p, converged - for delta in STEADY_STATE_DELTAS + for delta in (restart ? (0.0,) : STEADY_STATE_DELTAS) set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) # 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)] @@ -666,10 +668,11 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= try # Restart from other elevations: from far below the natural elevation the solver can # find a deep stall trim. - starts = [el_set .+ (5.0:5.0:85.0); el_set .- (5.0:5.0:85.0)] + starts = [el_set .+ (TRIM_RESTART_STEP:TRIM_RESTART_STEP:85.0); + el_set .- (TRIM_RESTART_STEP:TRIM_RESTART_STEP:85.0)] for el in filter(el -> 5.0 <= el <= 89.0, starts) s.set.elevation = el - y, yd, conv, st = solve(false) + y, yd, conv, st = solve(false; restart=true) if conv && !st prn && println("find_steady_state!: trimmed steady state found, starting at $el°") return y, yd @@ -678,9 +681,7 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= finally s.set.elevation = el_set end - # no valid trim found, return the result from el_set again - y0, yd0, _, stalled = solve(true) - stalled && @warn "find_steady_state!: only a trimmed steady state in deep stall was found " * - "(angle of attack $(round(s.alpha_2, digits=1))°)" - y0, yd0 + # no valid trim found: a deep stall or non-converged state is not a valid initial state + error("find_steady_state!: no trimmed steady state with attached flow was found; the " * + "solver did not converge, or the kite can only be balanced in deep stall") end diff --git a/src/KPS4.jl b/src/KPS4.jl index 492f51655..b614b65eb 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -725,7 +725,8 @@ elevation where its lift matches gravity (its natural elevation), starting the s elevation, but often not from far below it. Trimmed steady states with an angle of attack beyond the angle of maximum lift (deep stall) are rejected. If no valid trim is found from the start elevation, the solver is -restarted at higher and then at lower elevations than `s.set.elevation`, in steps of 5°. +restarted at higher and then at lower elevations than `s.set.elevation`, in steps of 10°. +If this fails too, an error is thrown, and `init!` returns `nothing`. The steady state is always solved for the full tether stiffness (`s.stiffness_factor` is 1.0 afterwards). The nonlinear solver is sensitive to the starting point, so small floating-point @@ -741,11 +742,12 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= trim=false) warn_delta_deprecated(delta, :find_steady_state!) warn_stiffness_factor_deprecated(stiffness_factor, :find_steady_state!) - # solve without perturbation of the initial state, with a perturbation only as fallback - function solve() + # solve without perturbation of the initial state, with a perturbation only as fallback; the + # restarts at other elevations use a cheaper solve, so that failing to find a trim is fast + function solve(; restart=false) local r, conv - for delta in STEADY_STATE_DELTAS - r, conv = solve_steady_state!(s; prn, delta, trim) + for delta in (restart ? (0.0,) : STEADY_STATE_DELTAS) + r, conv = solve_steady_state!(s; prn, delta, trim, continuation=!restart) conv && break end r, conv @@ -775,10 +777,11 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= if trim && (stalled || !converged) # Restart from other elevations: the solver converges from above the natural # elevation, but diverges or finds a deep stall trim from far below it. - starts = [el_set .+ (5.0:5.0:85.0); el_set .- (5.0:5.0:85.0)] + starts = [el_set .+ (TRIM_RESTART_STEP:TRIM_RESTART_STEP:85.0); + el_set .- (TRIM_RESTART_STEP:TRIM_RESTART_STEP:85.0)] for el in filter(el -> 5.0 <= el <= 89.0 && el != el_start, starts) s.set.elevation = el - r, conv = solve() + r, conv = solve(restart=true) conv || continue y, yd = steady_state(r.zero) if s.alpha_2 <= alpha_stall @@ -788,11 +791,10 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= break end end - if stalled - s.set.elevation = el_start - y00, yd00 = steady_state(results.zero) - @warn "find_steady_state!: only a trimmed steady state in deep stall was found " * - "(angle of attack $(round(s.alpha_2, digits=1))°)" + # no valid trim found: a deep stall or non-converged state is not a valid initial state + if stalled || !converged + error("find_steady_state!: no trimmed steady state with attached flow was found; " * + "the solver did not converge, or the kite can only be balanced in deep stall") end end if !converged @@ -832,7 +834,7 @@ const STIFFNESS_SCHEDULES = ((0.1, 1.0), (0.1, 0.3, 1.0), (0.3, 1.0)) # Solve for the offsets of the steady state for the full tether stiffness, see # find_steady_state!. Returns the result of nlsolve and whether it is usable. -function solve_steady_state!(s::KPS4; prn=false, delta=0.001, trim=false) +function solve_steady_state!(s::KPS4; prn=false, delta=0.001, trim=false, continuation=true) stiffness_factor = 1.0 s.stiffness_factor = stiffness_factor res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) @@ -918,7 +920,7 @@ function solve_steady_state!(s::KPS4; prn=false, delta=0.001, trim=false) end !isnothing(r) && prn && println("find_steady_state!: converged with the stiffness schedule $schedule") end - if isnothing(r) + if isnothing(r) && continuation # continuation: solve an easier problem, then step towards the full stiffness starts = [(sf0, delta) for sf0 in (0.1, 0.2, 0.3, 0.05) if sf0 != stiffness_factor] if delta != 0.001 diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 53ddac4d2..3d5c08809 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -87,6 +87,10 @@ const STEADY_STATE_DELTAS = (0.0, 0.001) # but often diverges or finds a deep stall trim from far below it. const TRIM_START_ELEVATION = 80.0 +# If no trim with attached flow is found from the start elevation, the steady state solvers +# restart at higher and lower elevations, in steps of TRIM_RESTART_STEP [°] +const TRIM_RESTART_STEP = 10.0 + function warn_delta_deprecated(delta, f) isnothing(delta) || Base.depwarn("The keyword argument `delta` of `$f` is deprecated and " * "ignored: the steady state is solved without perturbation (with a small perturbation " * @@ -678,8 +682,8 @@ Parameters: set up a consistent state or when the steady-state solver is known to struggle (e.g. unusual wind conditions). Returns `nothing` if steady-state solving fails. - trim: if `true`, the steady state is found at the natural elevation of the kite, where its - lift matches gravity, instead of at `s.set.elevation`; see `find_steady_state!`. Ignored if - `steady_state` is `false`. + lift matches gravity, instead of at `s.set.elevation`; see `find_steady_state!`. If no trim + with attached flow is found, `nothing` is returned. Ignored if `steady_state` is `false`. The keyword arguments `delta` and `stiffness_factor` are deprecated and ignored: the initial state is always an equilibrium for the full tether stiffness. @@ -702,6 +706,7 @@ function init!(s::AKM; stiffness_factor=nothing, delta=nothing, prn=false, stead catch e if e isa AssertionError || (e isa ErrorException && startswith(e.msg, "find_steady_state!:")) println("ERROR: Failure to find initial steady state in find_steady_state! function!") + e isa ErrorException && println(e.msg) return nothing else rethrow(e) diff --git a/test/test-trim.jl b/test/test-trim.jl index dc8dd55c8..156ebecdf 100644 --- a/test/test-trim.jl +++ b/test/test-trim.jl @@ -143,6 +143,17 @@ end @test s.alpha_2 < s.set.alpha_cl[argmax(s.set.cl_list)] end +@testset "no trim with attached flow, $Model" for Model in (KPS3, KPS4) + # at 4 m/s, Maasvlakte has only 5.2 m/s at the height of the kite, so no trim with attached + # flow exists: find_steady_state! throws an error and init! returns nothing instead of + # initializing from a deep stall or non-converged state + s = Model(KCU(profile_settings("Maasvlakte", 4.0))) + KiteModels.clear!(s) + @test_throws r"no trimmed steady state with attached flow" find_steady_state!(s; trim=true) + s = Model(KCU(profile_settings("Maasvlakte", 4.0))) + @test isnothing(init!(s; trim=true)) +end + @testset "reel-out speed in the steady state, $Model" for Model in (KPS3, KPS4) # reeling out reduces the apparent wind speed, so the natural elevation decreases, and the # initial state contains the reel-out speed From 40ba14c6b7aa456ec1cf360e0665f03a3726a8cc Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 12:17:50 +0200 Subject: [PATCH 11/21] Remove stiffness_factor from examples_3d, mwes and the docs The keyword is deprecated and ignored. Also update the paragraph in examples_4p.md, which still described the reduced stiffness. Co-Authored-By: Claude Opus 5.5 (1M context) --- docs/src/examples_4p.md | 5 +++-- examples_3d/auto_parking_4p.jl | 2 +- examples_3d/parking_4p.jl | 2 +- examples_3d/parking_wind_dir.jl | 2 +- mwes/mwe_01.jl | 2 +- mwes/mwe_02.jl | 2 +- mwes/mwe_11.jl | 2 +- mwes/mwe_11b.jl | 2 +- mwes/mwe_12.jl | 2 +- mwes/mwe_12b.jl | 2 +- mwes/mwe_26.jl | 2 +- 11 files changed, 13 insertions(+), 12 deletions(-) diff --git a/docs/src/examples_4p.md b/docs/src/examples_4p.md index dcd683937..5abb8bc70 100644 --- a/docs/src/examples_4p.md +++ b/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: diff --git a/examples_3d/auto_parking_4p.jl b/examples_3d/auto_parking_4p.jl index 126f0ce4c..f52051b1a 100644 --- a/examples_3d/auto_parking_4p.jl +++ b/examples_3d/auto_parking_4p.jl @@ -170,7 +170,7 @@ function simulate(integrator, sys_state) end function play() - integrator = KiteModels.init!(kps4; stiffness_factor = 0.1) + integrator = KiteModels.init!(kps4) sys_state = SysState(kps4) toc() simulate(integrator, sys_state) diff --git a/examples_3d/parking_4p.jl b/examples_3d/parking_4p.jl index 58726add8..436c7a15a 100644 --- a/examples_3d/parking_4p.jl +++ b/examples_3d/parking_4p.jl @@ -98,7 +98,7 @@ function simulate(integrator, steps, plot = PLOT) end toc() -integrator = KiteModels.init!(kps4, stiffness_factor = 0.1, prn = STATISTIC) +integrator = KiteModels.init!(kps4, prn = STATISTIC) toc() println("\nStarting simulation...") diff --git a/examples_3d/parking_wind_dir.jl b/examples_3d/parking_wind_dir.jl index 8359355ef..04a3ad112 100644 --- a/examples_3d/parking_wind_dir.jl +++ b/examples_3d/parking_wind_dir.jl @@ -190,7 +190,7 @@ end function play_parking() saved_turbulence = kps4.set.use_turbulence kps4.set.use_turbulence = 0.0 - integrator = KiteModels.init!(kps4; stiffness_factor=0.1) + integrator = KiteModels.init!(kps4) kps4.set.use_turbulence = saved_turbulence sys_state = SysState(kps4) if KiteUtils.PROJECT == "system.yaml" diff --git a/mwes/mwe_01.jl b/mwes/mwe_01.jl index 8e38aa9e3..b42385756 100644 --- a/mwes/mwe_01.jl +++ b/mwes/mwe_01.jl @@ -6,5 +6,5 @@ using KiteModels kps4::KPS4 = KPS4(KCU(se())) -integrator = KiteModels.init!(kps4; stiffness_factor=0.035, prn=false) +integrator = KiteModels.init!(kps4; prn=false) nothing \ No newline at end of file diff --git a/mwes/mwe_02.jl b/mwes/mwe_02.jl index 5a9882abc..50e9e01ec 100644 --- a/mwes/mwe_02.jl +++ b/mwes/mwe_02.jl @@ -11,7 +11,7 @@ set.solver="DFBDF" kps4::KPS4 = KPS4(KCU(set)) STEPS = 200 -integrator = KiteModels.init!(kps4; stiffness_factor=0.035, prn=false) +integrator = KiteModels.init!(kps4; prn=false) next_step!(kps4, integrator, set_speed = 0, dt=0.05) @timev next_step!(kps4, integrator; set_speed = 0, dt=0.05) function nsteps(kps4, integrator) diff --git a/mwes/mwe_11.jl b/mwes/mwe_11.jl index 6cdf8692f..18859cf66 100644 --- a/mwes/mwe_11.jl +++ b/mwes/mwe_11.jl @@ -13,5 +13,5 @@ kps3::KPS3 = KPS3(kcu) reltime = 0.0 -integrator = KiteModels.init!(kps3, stiffness_factor=0.04) +integrator = KiteModels.init!(kps3) plot2d(kps3.pos, reltime; zoom=false, front=false, segments=set.segments) \ No newline at end of file diff --git a/mwes/mwe_11b.jl b/mwes/mwe_11b.jl index 6979aa490..1f9f4cff9 100644 --- a/mwes/mwe_11b.jl +++ b/mwes/mwe_11b.jl @@ -13,5 +13,5 @@ kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) reltime = 0.0 -integrator = KiteModels.init!(kps4, stiffness_factor=0.5) +integrator = KiteModels.init!(kps4) plot2d(kps4.pos, reltime; zoom=true, front=false, segments=set.segments) \ No newline at end of file diff --git a/mwes/mwe_12.jl b/mwes/mwe_12.jl index a62a568d0..2d80c4f3c 100644 --- a/mwes/mwe_12.jl +++ b/mwes/mwe_12.jl @@ -13,5 +13,5 @@ kcu::KCU = KCU(set) kps3::KPS3 = KPS3(kcu) reltime = 0.0 -integrator = KiteModels.init!(kps3, stiffness_factor=0.04) +integrator = KiteModels.init!(kps3) plot2d(kps3.pos, reltime; zoom=false, front=false, segments=set.segments) \ No newline at end of file diff --git a/mwes/mwe_12b.jl b/mwes/mwe_12b.jl index b72d4a593..e65ace7ff 100644 --- a/mwes/mwe_12b.jl +++ b/mwes/mwe_12b.jl @@ -13,7 +13,7 @@ kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) reltime = 0.0 -integrator = KiteModels.init!(kps4, stiffness_factor=0.5) +integrator = KiteModels.init!(kps4) plot2d(kps4.pos, reltime; zoom=true, front=false, segments=set.segments) nothing \ No newline at end of file diff --git a/mwes/mwe_26.jl b/mwes/mwe_26.jl index 1c5c6414e..59ac0d62b 100644 --- a/mwes/mwe_26.jl +++ b/mwes/mwe_26.jl @@ -13,7 +13,7 @@ set.v_wind = 50 kps4::KPS4 = KPS4(KCU(set)) STEPS = 200 -integrator = KiteModels.init!(kps4; stiffness_factor=0.5, prn=false) +integrator = KiteModels.init!(kps4; prn=false) next_step!(kps4, integrator, set_speed = 0, dt=0.05) function nsteps(kps4, integrator) From 1359b7d9161662d4a4fddc5a08fd2ae090538c2f Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 12:19:28 +0200 Subject: [PATCH 12/21] Qualify the equilibrium guarantee in the init! docstring With steady_state=false, the initial state is not an equilibrium. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/KiteModels.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 3d5c08809..54f0936db 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -685,8 +685,8 @@ Parameters: lift matches gravity, instead of at `s.set.elevation`; see `find_steady_state!`. If no trim with attached flow is found, `nothing` is returned. Ignored if `steady_state` is `false`. -The keyword arguments `delta` and `stiffness_factor` are deprecated and ignored: the initial -state is always an equilibrium for the full tether stiffness. +The keyword arguments `delta` and `stiffness_factor` are deprecated and ignored: with +`steady_state=true`, the initial state is an equilibrium for the full tether stiffness. Returns: An instance of an `ODEIntegrator` or `IDAIntegrator`, or `nothing` if initialization fails. From 2a8395485eee9dd1ce1267862baf26aa9d920075 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 12:25:18 +0200 Subject: [PATCH 13/21] Update default manifest --- Manifest-v1.13.toml.default | 8 ++++---- examples/Project.toml | 1 - 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/Manifest-v1.13.toml.default b/Manifest-v1.13.toml.default index 63968033d..5dd6bfd86 100644 --- a/Manifest-v1.13.toml.default +++ b/Manifest-v1.13.toml.default @@ -1861,10 +1861,10 @@ 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" registries = "General" uuid = "6d616b69-6563-4f6e-8472-6f6c706c6f74" -version = "0.1.20" +version = "0.1.21" weakdeps = ["ControlSystemsBase"] [deps.MakieControlPlots.extensions] @@ -2715,10 +2715,10 @@ version = "1.0.0" [[deps.SIMD]] deps = ["PrecompileTools"] -git-tree-sha1 = "1ed93b3c4ffd33cc0c5d17eeb46f1b504ce42731" +git-tree-sha1 = "e24dc23107d426a096d3eae6c165b921e74c18e4" registries = "General" uuid = "fdea26ae-647d-5447-a871-4b548cad5224" -version = "3.8.0" +version = "3.7.2" [[deps.SIMDTypes]] git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" diff --git a/examples/Project.toml b/examples/Project.toml index a641c1bf4..efd07c889 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -19,7 +19,6 @@ Timers = "21f18d07-b854-4dab-86f0-c15a3821819a" [sources] KiteModels = {path = ".."} -MakieControlPlots = {rev = "main", url = "https://github.com/OpenSourceAWE/MakieControlPlots.jl.git"} [compat] BenchmarkTools = "1.6" From fe035d1cb63666f2c8a2730967845217e1361a16 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 21:32:16 +0200 Subject: [PATCH 14/21] Write each v0.12.1 changelog entry on one line Unwrap the bullets of the unreleased section so that diffs of single entries stay readable, and raise the markdownlint line-length limit to 2000 so the longest entries pass. Co-Authored-By: Claude Opus 5.5 --- .markdownlint.json | 2 +- CHANGELOG.md | 82 +++++++--------------------------------------- 2 files changed, 12 insertions(+), 72 deletions(-) diff --git a/.markdownlint.json b/.markdownlint.json index a965f00a9..a06efd782 100644 --- a/.markdownlint.json +++ b/.markdownlint.json @@ -4,7 +4,7 @@ "start_indented": false }, "MD013": { - "line_length": 1000, + "line_length": 2000, "tables": false }, "MD024": { diff --git a/CHANGELOG.md b/CHANGELOG.md index f5f9b085a..1a1332b27 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,78 +4,18 @@ 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. +- `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`. +- `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 From 26167d3bfce2ec887313786ac2229ff77230abfe Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 23:04:07 +0200 Subject: [PATCH 15/21] Fix string interpolation of el in find_steady_state! message MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit "$el°" interpolated the undefined variable el° because ° is a valid identifier character; use "$(el)°" instead. Co-Authored-By: Claude Opus 5.5 --- src/KPS3.jl | 2 +- src/KPS4.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KPS3.jl b/src/KPS3.jl index 88ab34ee2..1c289927d 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -674,7 +674,7 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= s.set.elevation = el y, yd, conv, st = solve(false; restart=true) if conv && !st - prn && println("find_steady_state!: trimmed steady state found, starting at $el°") + prn && println("find_steady_state!: trimmed steady state found, starting at $(el)°") return y, yd end end diff --git a/src/KPS4.jl b/src/KPS4.jl index b614b65eb..c8600f5ea 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -785,7 +785,7 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= conv || continue y, yd = steady_state(r.zero) if s.alpha_2 <= alpha_stall - prn && println("find_steady_state!: trimmed steady state found, starting at $el°") + prn && println("find_steady_state!: trimmed steady state found, starting at $(el)°") results, converged, stalled = r, true, false y00, yd00 = y, yd break From f82e4c6857faa19a4eb583823fd76e62a59e7403 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 23:14:29 +0200 Subject: [PATCH 16/21] Initialise the solve results before the loop in find_steady_state! The loops always run at least once, but JETLS cannot know this and reported the results as possibly undefined (lowering/undef-local-var). Co-Authored-By: Claude Opus 5.5 --- src/KPS3.jl | 2 +- src/KPS4.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KPS3.jl b/src/KPS3.jl index 1c289927d..d8b4e2dc5 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -636,7 +636,7 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= # converged and whether the kite is in deep stall. The restarts at other elevations use only # the unperturbed start value, so that failing to find a trim is fast. function solve(warn; restart=false) - local p, converged + p, converged = nothing, false for delta in (restart ? (0.0,) : STEADY_STATE_DELTAS) set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) # start with a straight, slightly stretched tether; at zero stretch the spring force has a kink diff --git a/src/KPS4.jl b/src/KPS4.jl index c8600f5ea..d01d01eb3 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -745,7 +745,7 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= # solve without perturbation of the initial state, with a perturbation only as fallback; the # restarts at other elevations use a cheaper solve, so that failing to find a trim is fast function solve(; restart=false) - local r, conv + r, conv = nothing, false for delta in (restart ? (0.0,) : STEADY_STATE_DELTAS) r, conv = solve_steady_state!(s; prn, delta, trim, continuation=!restart) conv && break From 3ff049c8855c6317938bd867b7d458e6ede036e8 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 23:19:55 +0200 Subject: [PATCH 17/21] Silence two JETLS false positives in KPS4 Initialise y00, yd00 in find_steady_state!: every path sets them before the return, but JETLS reported them as possibly undefined. Suppress the abstract-field diagnostic for KPS4.wm, which is abstract on purpose. Co-Authored-By: Claude Opus 5.5 --- .JETLSConfig.toml.default | 8 +++++++- src/KPS4.jl | 1 + 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/.JETLSConfig.toml.default b/.JETLSConfig.toml.default index d8b081065..fb9df7758 100644 --- a/.JETLSConfig.toml.default +++ b/.JETLSConfig.toml.default @@ -22,4 +22,10 @@ severity = "off" pattern = "`integrator` is not defined" match_by = "message" match_type = "literal" -severity = "off" \ No newline at end of file +severity = "off" + +[[diagnostic.patterns]] +pattern = "`KiteModels.KPS4` has abstract field `wm::WinchModels.AbstractWinchModel`" +match_by = "message" +match_type = "literal" +severity = "off" diff --git a/src/KPS4.jl b/src/KPS4.jl index d01d01eb3..663716f19 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -769,6 +769,7 @@ function find_steady_state!(s::KPS4; prn=false, delta=nothing, stiffness_factor= residual!(res, yd_, y_, s, 0.0) y_, yd_ end + y00, yd00 = nothing, nothing try s.set.elevation = el_start results, converged = solve() From 25902f8c35136fe0eaa64652f49751682e2a0eba Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 23:23:27 +0200 Subject: [PATCH 18/21] Avoid boxed captured variables in solve_steady_state! Make r local in try_solve, which wrote to the r of the enclosing function. Suppress the JETLS warning that solve is captured and boxed: the boxing comes from the keyword argument of the local function and costs nothing compared with nlsolve. Co-Authored-By: Claude Opus 5.5 --- .JETLSConfig.toml.default | 7 +++++++ src/KPS4.jl | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/.JETLSConfig.toml.default b/.JETLSConfig.toml.default index fb9df7758..a1a9b51bd 100644 --- a/.JETLSConfig.toml.default +++ b/.JETLSConfig.toml.default @@ -29,3 +29,10 @@ pattern = "`KiteModels.KPS4` has abstract field `wm::WinchModels.AbstractWinchMo 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" diff --git a/src/KPS4.jl b/src/KPS4.jl index 663716f19..ea27793ba 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -902,7 +902,7 @@ function solve_steady_state!(s::KPS4; prn=false, delta=0.001, trim=false, contin # returns the result for (sf, d), trying both scaling options, or nothing function try_solve(sf, d, x0) for autoscale in (true, false) - r = solve(sf, d, x0; autoscale) + local r = solve(sf, d, x0; autoscale) usable(r) && return r end nothing From 3cc5af9496b2617f0f150856ca9b8768f20d768e Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 23:25:43 +0200 Subject: [PATCH 19/21] Add type assertions to Settings fields in KPS3 Settings has a custom getproperty, so JETLS cannot always infer the field types and reported method errors for String arguments. Co-Authored-By: Claude Opus 5.5 --- src/KPS3.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/KPS3.jl b/src/KPS3.jl index d8b4e2dc5..985d08fcb 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -482,7 +482,7 @@ function init_inner(s::KPS3, X=zeros(2 * s.set.segments); delta=0.0, upwind_dir= # reeling out: 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 for i in 1:s.set.segments - vel[i+1] .+= normalize(pos[i+1] - pos[i]) * (s.set.v_reel_out * i / s.set.segments) + vel[i+1] .+= normalize(pos[i+1] - pos[i]) * (s.set.v_reel_out::Float64 * i / s.set.segments) end for i in eachindex(pos) @inbounds s.pos[i] .= pos[i] @@ -540,9 +540,9 @@ end # 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) + segments = s.set.segments::Int64 + l0 = s.set.l_tether::Float64 / segments + elevation = deg2rad(s.set.elevation::Float64) X = zeros(SimFloat, 2segments) x, z = 0.0, 0.0 for i in 1:segments From f83f332389cd5292c6e2a0a8eedcc4a59403c7c1 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Sun, 11 Oct 2026 23:34:14 +0200 Subject: [PATCH 20/21] Make the results of solve local in find_steady_state! of KPS3 The local function solve assigned y0, yd0 and converged of the enclosing function, so these variables were boxed. Co-Authored-By: Claude Opus 5.5 --- src/KPS3.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KPS3.jl b/src/KPS3.jl index 985d08fcb..4796f9d5a 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -636,7 +636,7 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= # converged and whether the kite is in deep stall. The restarts at other elevations use only # the unperturbed start value, so that failing to find a trim is fast. function solve(warn; restart=false) - p, converged = nothing, false + local p, converged = nothing, false for delta in (restart ? (0.0,) : STEADY_STATE_DELTAS) set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) # start with a straight, slightly stretched tether; at zero stretch the spring force has a kink @@ -653,7 +653,7 @@ function find_steady_state!(s::KPS3; prn=false, delta=nothing, stiffness_factor= all(isfinite, p) || error("find_steady_state!: solver returned non-finite values. Cannot compute steady state.") end all(isfinite, p) || return nothing, nothing, false, false - y0, yd0 = init(s, steady_state_offsets(s, p); upwind_dir) + local y0, yd0 = init(s, steady_state_offsets(s, p); upwind_dir) residual!(zeros(length(y0)), yd0, y0, s) # updates the angle of attack s.alpha_2 y0, yd0, converged, s.alpha_2 > alpha_stall end From 914e8aee63afc8b901dfefa69790e9010f3792ec Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Mon, 12 Oct 2026 00:47:58 +0200 Subject: [PATCH 21/21] Update default manifest --- Manifest-v1.12.toml.default | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/Manifest-v1.12.toml.default b/Manifest-v1.12.toml.default index a4166ac89..d1f098d43 100644 --- a/Manifest-v1.12.toml.default +++ b/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"]