Make find_steady_state! for KPS3 find a real equilibrium

.lychee.toml

@@ -1,3 +1,7 @@
+# doi.org redirects to zenodo.org, which often needs about 20 s to answer (default timeout: 20 s)
+timeout = 60
+max_retries = 3
+
 exclude = [
   "@ref",
   "@cite",

CHANGELOG.md

@@ -5,12 +5,24 @@ SPDX-License-Identifier: MIT
 ### KiteModels v0.11.19 - unreleased
 #### Fixed
 - The state of the winch (brake on/off, rate limited set speed) is now updated once per time step
-  in `next_step!` (and reset in `init!` to the initial set speed `set.v_reel_out` and the brake state
-  after construction, so that repeated calls of `init!` give the same result) instead of on every call of the residual function. Before,
-  it depended on the number of residual evaluations of the DAE solver, so tiny numerical
-  differences (e.g. between Julia versions) could switch the brake at a different moment. With a
-  winch `v_min` of 0.15 m/s this made the reel-out speed of the hydra20 simulations of
-  KiteControllers.jl run away on Julia 1.13. The simulation results change slightly.
+  in `next_step!` (and reset in `init!` to the initial set speed `set.v_reel_out` and the brake
+  state after construction, so that repeated calls of `init!` give the same result) instead of on
+  every call of the residual function. Before, it depended on the number of residual evaluations
+  of the DAE solver, so tiny numerical differences (e.g. between Julia versions) could switch the
+  brake at a different moment. With a winch `v_min` of 0.15 m/s this made the reel-out speed of the
+  hydra20 simulations of KiteControllers.jl run away on Julia 1.13. The simulation results change
+  slightly.
+- `find_steady_state!` (KPS3) now finds a real equilibrium. Before, `nlsolve` usually stopped
+  because its steps became tiny, with the accelerations of the tether particles still at about
+  25 m/s², and on some machines it printed `find_steady_state!: solver did not converge!`. The
+  unknowns are now the angles and relative stretches of the tether segments instead of the offsets
+  of the particle positions, the solver starts from a slightly stretched tether (the spring force
+  has a kink at zero stretch), and, as for KPS4, the elevation of the kite is prescribed instead of
+  its vertical force balance, so `calc_elevation(s)` returns `set.elevation`. The tether particles
+  and the horizontal force balance of the kite are now solved to `ftol = 1e-6`; a warning is only
+  printed if this fails. The horizontal residuals are taken in the plane of the tether, so the
+  result no longer depends on `upwind_dir`. The initial state changes: for example, with a 392 m
+  tether the kite starts at 70° instead of 64.3°.
 #### Changed
 - requires WinchModels 0.3.12 (for `update_winch_state!` and `calc_acceleration(...; update_state)`)
 - support Julia 1.12 and 1.13 only, as WinchModels 0.3.12 and KiteUtils do: `julia` compat

src/KPS3.jl

@@ -531,46 +531,86 @@ function spring_forces(s::KPS3)
     forces
 end
 
-function find_steady_state_inner(s::KPS3, X, prn=false; delta=0.0, upwind_dir=nothing)
+# Convert the unknowns of the steady state solver, the angle offsets of the tether segments from the
+# elevation angle [rad] (first half) and their relative stretches [-] (second half), to the offsets of
+# the particle positions from the straight, unstretched tether in x and z direction, as used by init.
+function steady_state_offsets(s::KPS3, p)
+    segments = s.set.segments
+    l0 = s.set.l_tether / segments
+    elevation = deg2rad(s.set.elevation)
+    X = zeros(SimFloat, 2segments)
+    x, z = 0.0, 0.0
+    for i in 1:segments
+        angle = elevation + p[i]
+        len = l0 * (1 + p[segments+i])
+        x += len * cos(angle)
+        z += len * sin(angle)
+        X[i]          = x - i * l0 * cos(elevation)
+        X[segments+i] = z - i * l0 * sin(elevation)
+    end
+    X
+end
+
+function find_steady_state_inner(s::KPS3, p, prn=false; delta=0.0, upwind_dir=nothing, warn=true)
     res = zeros(MVector{6*s.set.segments+2, SimFloat})
+    segments = s.set.segments
+    # init turns the tether into the wind direction; turn the residuals and positions back
+    turnangle = something(upwind_dir, -pi/2) + pi/2
+    horizontal(x, y) = cos(turnangle) * x - sin(turnangle) * y
 
     # helper function for the steady state finder
+    # Equations: the horizontal and vertical force balance of the tether particles, the horizontal
+    # force balance of the kite and its elevation angle. The kite can only be fully balanced at its
+    # natural elevation, so the elevation from the settings is prescribed instead of its vertical
+    # force balance.
     function test_initial_condition!(F, x::Vector)
-        y0, yd0 = init(s, x; delta, upwind_dir)
+        y0, yd0 = init(s, steady_state_offsets(s, x); delta, upwind_dir)
         residual!(res, yd0, y0, s)
-        for i in 1:s.set.segments
-            F[i]                = res[1 + 3*(i-1) + 3*s.set.segments]
-            F[i+s.set.segments] = res[3 + 3*(i-1) + 3*s.set.segments]
+        for i in 1:segments
+            j = 3*(i-1) + 3*segments
+            F[i]          = horizontal(res[j+1], res[j+2])
+            F[i+segments] = res[j+3]
         end
-        return nothing 
+        # replace the vertical force balance of the kite by the elevation angle, scaled to meters
+        j = 3*(segments-1)
+        x_kite = horizontal(y0[j+1], y0[j+2])
+        z_kite = y0[j+3]
+        F[2segments] = (atan(z_kite, x_kite) - deg2rad(s.set.elevation)) * hypot(x_kite, z_kite)
+        return nothing
     end
 
     if prn println("\nStarted function test_nlsolve...") end
-    jac! = make_jac(test_initial_condition!, length(X))
-    results = nlsolve(test_initial_condition!, jac!, X, xtol=1e-6, ftol=1e-6, autoscale=true, iterations=1000)
+    jac! = make_jac(test_initial_condition!, length(p))
+    results = nlsolve(test_initial_condition!, jac!, p, xtol=1e-10, ftol=1e-6, autoscale=true, iterations=1000)
     if prn println("\nresult: $results") end
-    if !converged(results)
+    if !results.f_converged && warn
         @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations))"
         # Check if the solution contains finite values
         if !all(isfinite, results.zero)
             error("find_steady_state!: solver returned non-finite values. Cannot compute steady state.")
         end
     end
-    results.zero
+    results.zero, results.f_converged
  end
 
 """
-    find_steady_state!(s::KPS3; prn=false, delta = 0.0, stiffness_factor=0.035, upwind_dir=-pi/2)
+    find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2)
 
 Find an initial equilibrium, based on the initial parameters
 `l_tether`, elevation and `v_reel_out`.
+
+The tether particles are in equilibrium. The kite is placed at the elevation angle from the
+settings, so only its horizontal force balance is solved for; the kite can only be fully balanced
+at its natural elevation.
 """
 function find_steady_state!(s::KPS3; prn=false, delta = 0.002, stiffness_factor=0.035, upwind_dir=-pi/2)
     set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir)
-    zero = zeros(SimFloat, 2*s.set.segments)
+    # start with a straight, slightly stretched tether; at zero stretch the spring force has a kink
+    p0 = [zeros(SimFloat, s.set.segments); fill(SimFloat(1e-3), s.set.segments)]
     s.stiffness_factor=stiffness_factor
-    zero = find_steady_state_inner(s, zero, prn; delta, upwind_dir)
+    p, converged = find_steady_state_inner(s, p0, prn; delta, upwind_dir, warn=false)
     s.stiffness_factor=1.0
-    zero = find_steady_state_inner(s, zero, prn; delta, upwind_dir)
-    init(s, zero; delta=delta, upwind_dir)
+    # if the solution with the reduced stiffness failed, start again from the straight tether
+    p, _ = find_steady_state_inner(s, converged ? p : p0, prn; delta, upwind_dir)
+    init(s, steady_state_offsets(s, p); delta=delta, upwind_dir)
 end

test/test-kps3.jl

@@ -390,27 +390,55 @@ const SEGMENTS = load_settings("system.yaml").segments
         local kps = KPS3(KCU(load_settings("system.yaml")))
         init_392(kps)
         KiteModels.set_depower_steering!(kps, 0.25, 0.0)
-        try
-            res1, res2 = find_steady_state!(kps; delta = 1e-6, prn = false)
-            @test norm(res2) < 1e-5                            # velocity and acceleration must be near zero
-            pre_tension = KiteModels.calc_pre_tension(kps)
-            @test pre_tension > 1.0001
-            @test pre_tension < 1.01
-            @test unstretched_length(kps) ≈ 392.0                        # initial, unstretched tether length
-            @test tether_length(kps) ≈ 392.1861381318156 rtol = 1e-5     # real, stretched tether length
-            @test winch_force(kps) ≈ 276.25751212763817 rtol = 3e-2      # initial force at the winch [N]
-            lift, drag = lift_drag(kps)
-            @test lift ≈ 443.63277537186394 rtol=2e-2                # initial lift force of the kite [N]
-            @test drag ≈ 94.25218065939362 rtol=2e-2               # initial drag force of the kite [N]
-            @test lift_over_drag(kps) ≈ 4.706870146326417 rtol=2e-3     # initial lift-over-drag
-            @test norm(v_wind_kite(kps)) ≈ 9.107670173739065 rtol=1e-2   # initial wind speed at the height of the kite [m/s]
-        catch e
-            if e isa ErrorException && contains(e.msg, "find_steady_state!") && get(ENV, "CI", "false") == "true"
-                @warn "Steady state solver failed to converge. Skipping test."
-                @test_broken false  # Mark as known issue (CI flake only)
-            else
-                rethrow(e)
-            end
+        # the solver must converge without warning; the kite is placed at the elevation from the settings
+        res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; delta = 1e-6, prn = false)
+        @test rad2deg(calc_elevation(kps)) ≈ 70.0 atol=1e-6
+        # the tether particles are in equilibrium (x and z components of their accelerations)
+        res = zeros(length(res1))
+        KiteModels.residual!(res, res2, res1, kps)
+        segments = kps.set.segments
+        for i in 1:segments-1
+            j = 3*(i-1) + 3*segments
+            @test abs(res[j+1]) < 1e-5
+            @test abs(res[j+3]) < 1e-5
+        end
+        @test abs(res[3*(segments-1) + 3*segments + 1]) < 1e-5     # horizontal balance of the kite
+        @test norm(res2) < 1e-5                            # velocity and acceleration must be near zero
+        pre_tension = KiteModels.calc_pre_tension(kps)
+        @test pre_tension > 1.0001
+        @test pre_tension < 1.01
+        @test unstretched_length(kps) ≈ 392.0                        # initial, unstretched tether length
+        @test tether_length(kps) ≈ 392.2025410439266 rtol = 1e-5     # real, stretched tether length
+        @test winch_force(kps) ≈ 301.23463553992514 rtol = 3e-2      # initial force at the winch [N]
+        lift, drag = lift_drag(kps)
+        @test lift ≈ 327.30289678275346 rtol=2e-2                # initial lift force of the kite [N]
+        @test drag ≈ 73.01286087237223 rtol=2e-2                 # initial drag force of the kite [N]
+        @test lift_over_drag(kps) ≈ 4.482811560485004 rtol=2e-3      # initial lift-over-drag
+        @test norm(v_wind_kite(kps)) ≈ 9.107670173739065 rtol=1e-2   # initial wind speed at the height of the kite [m/s]
+    end
+
+    @testset "test_find_steady_state upwind_dir = $(round(rad2deg(upwind_dir)))°" for upwind_dir in (-π/2, 0.0, π/4, π)
+        local kps = KPS3(KCU(load_settings("system.yaml")))
+        init_392(kps)
+        KiteModels.set_depower_steering!(kps, 0.25, 0.0)
+        res1, res2 = @test_logs min_level=Base.CoreLogging.Warn find_steady_state!(kps; delta = 1e-6, upwind_dir)
+        @test rad2deg(calc_elevation(kps)) ≈ 70.0 atol=1e-6
+        @test winch_force(kps) ≈ 301.23463553992514 rtol = 1e-4   # independent of the wind direction
+        # the tether lies in the vertical plane of the wind direction; take the components of the
+        # accelerations in this plane (horizontal and vertical) and normal to it
+        turnangle = upwind_dir + π/2
+        in_plane(x, y) = cos(turnangle) * x - sin(turnangle) * y
+        normal(x, y)   = sin(turnangle) * x + cos(turnangle) * y
+        segments = kps.set.segments
+        j = 3*(segments-1)
+        @test abs(normal(res1[j+1], res1[j+2])) < 1e-3          # kite position in the plane
+        @test in_plane(res1[j+1], res1[j+2]) > 0                 # downwind of the ground station
+        res = zeros(length(res1))
+        KiteModels.residual!(res, res2, res1, kps)
+        for i in 1:segments
+            j = 3*(i-1) + 3*segments
+            @test abs(in_plane(res[j+1], res[j+2])) < 1e-5       # horizontal balance in the plane
+            i < segments && @test abs(res[j+3]) < 1e-5            # vertical balance of the tether particles
         end
     end