Migration guide: NLsolve.nlsolve with manual Jacobian and autoscale for kite power system

[ChrisRackauckas-Claude]

Context

KiteModels.jl uses NLsolve.jl to find the steady-state equilibrium of a 4-point kite power system (3D mass-spring network with aerodynamic forces). We'd like to migrate to NonlinearSolve.jl but the current code relies on several NLsolve-specific features.

Current NLsolve usage

The core call pattern is:

using NLsolve

# Manual central finite-difference Jacobian with subnormal flushing
function make_jac(f!, n_vars)
    _F1 = zeros(Float64, n_vars)
    _F2 = zeros(Float64, n_vars)
    _xp = zeros(Float64, n_vars)
    _xm = zeros(Float64, n_vars)
    threshold = sqrt(floatmin(Float64))
    function jac!(J, x)
        h_factor = cbrt(eps(Float64))
        for j in 1:n_vars
            copyto!(_xp, x); copyto!(_xm, x)
            h = max(abs(x[j]), one(Float64)) * h_factor
            _xp[j] += h; _xm[j] -= h
            f!(_F1, _xp); f!(_F2, _xm)
            @views J[:, j] .= (_F1 .- _F2) ./ (2h)
        end
        @. J = ifelse(abs(J) < threshold, zero(Float64), J)
    end
    return jac!
end

n_unknowns = 20  # 2*(segments + kite_particles - 1) + 2
X00 = zeros(Float64, n_unknowns)
jac! = make_jac(f!, n_unknowns)
results = nlsolve(f!, jac!, X00;
                  autoscale=true, xtol=4e-7, ftol=4e-7, iterations=200)
if NLsolve.converged(results)
    solution = results.zero
end

NLsolve features we depend on

1. Manual Jacobian (jac!): Needed because NLSolversBase >= 7.10 produces subnormal Jacobian entries via DifferentiationInterface, which cause NaN in NLsolve's autoscale (subnormal column norms squared underflow to zero). The manual Jacobian flushes near-zero entries.

2. autoscale=true: Critical for convergence — the system variables have very different scales (positions in meters vs. small perturbations).

3. Separate xtol and ftol: Both set to 4e-7.

4. iterations limit: Set to 200.

5. results.zero and converged(results): Solution extraction and convergence check.

Questions

1. What is the equivalent NonlinearSolve.jl API for this call pattern?
2. Does NonlinearSolve handle autoscaling? (I see #425 is open for trust region autoscaling.)
3. Does NonlinearSolve's automatic differentiation avoid the subnormal Jacobian issue, or would we still need a manual Jacobian?
4. What's the recommended solver algorithm for this kind of problem (20-variable nonlinear system from physics, starting from a reasonable initial guess)?

Self-contained MWE

Below is a complete, self-contained MWE (~960 lines) that reproduces the KiteModels.jl steady-state test without any kite-specific package dependencies. It models a 4-point kite on a segmented tether with aerodynamic forces, spring/damping, and a winch. The nonlinear solve finds particle positions where all accelerations vanish (steady-state equilibrium).

Dependencies: StaticArrays, LinearAlgebra, NLsolve, Dierckx, Test

Click to expand full MWE (mwe_steady_state_kps4.jl)

$(cat /home/crackauc/sandbox/tmp_20260219_005256_56738/mwe_steady_state_kps4.jl)

Running the MWE

] add StaticArrays NLsolve Dierckx
include("mwe_steady_state_kps4.jl")

Expected output (all 15 tests pass):

  l=100.0: tether_length=100.92196416459542, pre_tension=1.0003226874576083 ✓
  l=200.0: tether_length=202.59541320501268, pre_tension=1.0004541973108771 ✓
  l=392.0: tether_length=399.15270436705833, pre_tension=1.0006386343184874 ✓
Test Summary:                | Pass  Total   Time
test_find_steady_state (MWE) |   15     15  15.3s

Environment

• Julia 1.11
• NLsolve v4.5.1
• StaticArrays v1.9.16

[ChrisRackauckas-Claude]

# KiteModels.jl-free MWE for test-steady-state-kps4.jl
#
# Self-contained reproduction of the KPS4 steady-state test without depending on
# KiteModels, KiteUtils, KitePodModels, AtmosphericModels, or WinchModels.
#
# External deps: StaticArrays, LinearAlgebra, NLsolve, Dierckx, Test

using StaticArrays, LinearAlgebra, NLsolve, Dierckx, Test

set_zero_subnormals(true)

# ═══════════════════════════════════════════════════════════════════════════════
# Constants
# ═══════════════════════════════════════════════════════════════════════════════
const SimFloat    = Float64
const KVec3       = MVector{3, SimFloat}
const SVec3       = SVector{3, SimFloat}
const G_EARTH     = 9.81
const KITE_PARTICLES = 4
const KITE_SPRINGS   = 9
const DRAG_CORR      = 0.93
const KITE_ANGLE     = 3.83   # angle offset between kite and last tether segment [deg]
const PRE_STRESS     = 0.9998 # multiplier for initial bridle spring lengths
const SEGMENTS       = 6
const NP = SEGMENTS + KITE_PARTICLES + 1  # 11 total particles (incl. origin)
const NS = SEGMENTS + KITE_SPRINGS        # 15 total springs

# Custom norm for speed (matches KiteModels)
import LinearAlgebra: norm
@inline @inbounds function norm(vec::SVector{3, SimFloat})
    sqrt(vec[1]*vec[1] + vec[2]*vec[2] + vec[3]*vec[3])
end
@inline @inbounds function norm(vec::MVector{3, SimFloat})
    sqrt(vec[1]*vec[1] + vec[2]*vec[2] + vec[3]*vec[3])
end

# Safe acos (clamps to [-1, 1])
acos2(x) = acos(clamp(x, -1.0, 1.0))

# ═══════════════════════════════════════════════════════════════════════════════
# Settings (hardcoded from KiteModels.jl/data/settings.yaml)
# ═══════════════════════════════════════════════════════════════════════════════
mutable struct Settings
    segments::Int64; sample_freq::Int64
    l_tethers::Vector{Float64}; v_reel_outs::Vector{Float64}
    elevations::Vector{Float64}; depowers::Vector{Float64}
    max_iter::Int64
    c0::Float64; c_s::Float64; k_ds::Float64; max_steering::Float64; cs_4p::Float64
    alpha_d_max::Float64; depower_offset::Float64
    version::Int64; mass::Float64; area::Float64; rel_side_area::Float64; height_k::Float64
    alpha_cl::Vector{Float64}; cl_list::Vector{Float64}
    alpha_cd::Vector{Float64}; cd_list::Vector{Float64}
    width::Float64; alpha_zero::Float64; alpha_ztip::Float64; m_k::Float64
    rel_nose_mass::Float64; rel_top_mass::Float64
    smc::Float64; cmq::Float64; cord_length::Float64
    d_line::Float64; l_bridle::Float64; h_bridle::Float64
    rel_compr_stiffness::Float64; rel_damping::Float64
    kcu_mass::Float64; kcu_diameter::Float64; cd_kcu::Float64
    power2steer_dist::Float64; depower_drum_diameter::Float64; tape_thickness::Float64
    d_tether::Float64; cd_tether::Float64
    axial_damping::Float64; axial_stiffness::Float64
    rho_tether::Float64; e_tether::Float64
    max_force::Float64; drum_radius::Float64; gear_ratio::Float64
    inertia_total::Float64; f_coulomb::Float64; c_vf::Float64; max_acc::Float64
    v_wind::Float64; temp_ref::Float64; height_gnd::Float64; h_ref::Float64
    rho_0::Float64; alpha::Float64; z0::Float64; profile_law::Int64
end

# Property accessors (mimic KiteUtils.Settings virtual properties)
function Base.getproperty(s::Settings, sym::Symbol)
    sym === :l_tether   && return getfield(s, :l_tethers)[1]
    sym === :v_reel_out && return getfield(s, :v_reel_outs)[1]
    sym === :elevation  && return getfield(s, :elevations)[1]
    sym === :depower    && return getfield(s, :depowers)[1]
    return getfield(s, sym)
end
function Base.setproperty!(s::Settings, sym::Symbol, val)
    sym === :l_tether   && (getfield(s, :l_tethers)[1] = val; return val)
    sym === :v_reel_out && (getfield(s, :v_reel_outs)[1] = val; return val)
    sym === :elevation  && (getfield(s, :elevations)[1] = val; return val)
    sym === :depower    && (getfield(s, :depowers)[1] = val; return val)
    if val isa Int && getfield(s, sym) isa Float64
        setfield!(s, sym, Float64(val))
    else
        setfield!(s, sym, val)
    end
    return val
end

function default_settings()
    Settings(
        6, 20,                              # segments, sample_freq
        [150.0], [0.0], [70.8], [25.0],    # l_tethers, v_reel_outs, elevations, depowers
        200,                                # max_iter
        0.0, 2.59, 1.5, 16.834, 1.0,       # c0, c_s, k_ds, max_steering, cs_4p
        31.0, 23.6,                         # alpha_d_max, depower_offset
        1, 6.2, 10.18, 30.6, 2.23,         # version, mass, area, rel_side_area, height_k
        [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0],
        [   0.0,    0.5,   0.0,  0.08, 0.125,  0.15, 0.2,  1.0,  1.0,  0.0,  -0.5,   0.0],
        [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0],
        [   0.5,    0.5,    0.5,   1.0,   0.2, 0.1,  0.2,  1.0,   0.5,   0.5,   0.5],
        5.77, 4.0, 10.0, 0.2,              # width, alpha_zero, alpha_ztip, m_k
        0.47, 0.4,                          # rel_nose_mass, rel_top_mass
        0.0, 0.0, 2.0,                     # smc, cmq, cord_length
        2.5, 33.4, 4.9,                    # d_line, l_bridle, h_bridle
        0.25, 6.0,                          # rel_compr_stiffness, rel_damping
        8.4, 0.4, 0.0,                     # kcu_mass, kcu_diameter, cd_kcu
        1.3, 0.069, 0.0006,                # power2steer_dist, depower_drum_diameter, tape_thickness
        4.0, 0.958,                         # d_tether, cd_tether
        473.0, 614600.0,                    # axial_damping, axial_stiffness
        724.0, 55e9,                        # rho_tether, e_tether
        4000.0, 0.1615, 6.2,               # max_force, drum_radius, gear_ratio
        0.204, 122.0, 30.6, 4.0,           # inertia_total, f_coulomb, c_vf, max_acc
        9.51, 15.0, 0.0, 6.0,              # v_wind, temp_ref, height_gnd, h_ref
        1.225, 0.08163, 0.0002, 3,         # rho_0, alpha, z0, profile_law
    )
end

# ═══════════════════════════════════════════════════════════════════════════════
# Atmospheric Model (from AtmosphericModels.jl)
# ═══════════════════════════════════════════════════════════════════════════════
struct AtmosphericModel
    set::Settings
    rho_zero_temp::Float64
end

function AtmosphericModel(set::Settings)
    rzt = (15.0 + 273.15) / (set.temp_ref + 273.15) * set.rho_0
    AtmosphericModel(set, rzt)
end

calc_rho(am::AtmosphericModel, height) = am.rho_zero_temp * exp(-(height + am.set.height_gnd) / 8550.0)

function calc_wind_factor(am::AtmosphericModel, height)
    s = am.set
    if s.profile_law == 1       # EXP
        exp(s.alpha * log(height / s.h_ref))
    elseif s.profile_law == 2   # LOG
        log(height / s.z0) / log(s.h_ref / s.z0)
    elseif s.profile_law == 3   # EXPLOG
        K = 1.0
        log1 = log(height / s.z0) / log(s.h_ref / s.z0)
        exp1 = exp(s.alpha * log(height / s.h_ref))
        log1 + K * (log1 - exp1)
    else
        error("invalid profile_law: $(s.profile_law)")
    end
end

# ═══════════════════════════════════════════════════════════════════════════════
# KCU — Kite Control Unit (from KitePodModels.jl)
# ═══════════════════════════════════════════════════════════════════════════════
mutable struct KCU
    set::Settings
    set_depower::Float64
    set_steering::Float64
    depower::Float64
    steering::Float64
end

function KCU(set::Settings)
    dp = set.depower_offset * 0.01
    KCU(set, dp, 0.0, dp, 0.0)
end

get_depower(kcu::KCU) = kcu.depower
get_steering(kcu::KCU) = kcu.steering

function calc_delta_l(kcu::KCU, rel_depower)
    u = kcu.set.depower_drum_diameter
    l_ro = 0.0
    rotations = (rel_depower - 0.01 * kcu.set.depower_offset) *
                10.0 * 11.0 / 3.0 * (3918.8 - 230.8) / 4096.0
    while rotations > 0.0
        l_ro += u * π
        rotations -= 1.0
        u -= kcu.set.tape_thickness
    end
    if rotations < 0.0
        l_ro += (-(u + kcu.set.tape_thickness) * rotations +
                 u * (rotations + 1.0)) * π * rotations
    end
    return l_ro
end

function calc_alpha_depower(kcu::KCU, rel_depower)
    a   = kcu.set.power2steer_dist
    b_0 = kcu.set.h_bridle + 0.5 * kcu.set.height_k
    b   = b_0 + 0.5 * calc_delta_l(kcu, rel_depower)
    c   = sqrt(a * a + b_0 * b_0)
    @assert c < a + b "ERROR in calc_alpha_depower!"
    tmp = clamp(1 / (2 * a * b) * (a * a + b * b - c * c), -1.0, 1.0)
    return π / 2.0 - acos(tmp)
end

# ═══════════════════════════════════════════════════════════════════════════════
# AsyncMachine — Winch Model (from WinchModels.jl)
# ═══════════════════════════════════════════════════════════════════════════════
mutable struct AsyncMachine
    set::Settings
    u_nom::Float64; omega_sn::Float64; omega_mn::Float64
    tau_n::Float64; tau_b::Float64
    v_min::Float64; brake_acc::Float64
    brake::Bool; last_set_speed::Float64
end

function AsyncMachine(set::Settings)
    AsyncMachine(set, 400.0/sqrt(3), 1500/60*2π, 1460/60*2π,
                 121.0, 363.0, 0.2, -25.0, true, 0.0)
end

smooth_sign(x) = x / sqrt(x * x + 36.0)   # EPSILON=6, 6^2=36

function calc_acceleration(wm::AsyncMachine, set_speed, speed, force, use_brake=false)
    dt = 1.0 / wm.set.sample_freq
    if use_brake
        if abs(set_speed) < 0.9 * wm.v_min;  wm.brake = true
        elseif abs(set_speed) > 1.1 * wm.v_min; wm.brake = false; end
        wm.brake && return wm.brake_acc * speed
    end
    MAX_ACC = wm.set.max_acc
    ls = set_speed
    if ls > wm.last_set_speed + MAX_ACC * dt;  ls = wm.last_set_speed + MAX_ACC * dt
    elseif ls < wm.last_set_speed - MAX_ACC * dt; ls = wm.last_set_speed - MAX_ACC * dt; end
    wm.last_set_speed = ls
    r = wm.set.drum_radius; n = wm.set.gear_ratio
    omega      = n / r * speed
    omega_sync = n / r * ls
    delta = omega_sync - omega
    R2 = (wm.u_nom^2 * (wm.omega_sn - wm.omega_mn)) / (2 * wm.omega_sn^2) *
         (1/wm.tau_n + sqrt(1/wm.tau_n^2 - 1/wm.tau_b^2))
    L  = (wm.u_nom^2 / (2 * wm.omega_sn * wm.tau_b)) / wm.omega_sn
    oref = abs(omega_sync) <= wm.omega_sn ? wm.omega_sn : omega_sync
    omega_dot = (wm.u_nom^2 * R2 * delta) / (oref^2 * (R2^2 + L^2 * delta^2))
    τ_c = wm.set.f_coulomb * r / n * smooth_sign(omega)
    τ_v = wm.set.c_vf * omega * r^2 / n^2
    omega_dot += r / n * force * 4000.0 / wm.set.max_force - τ_c - τ_v
    omega_dot /= wm.set.inertia_total
    return r / n * omega_dot
end

# ═══════════════════════════════════════════════════════════════════════════════
# Spring
# ═══════════════════════════════════════════════════════════════════════════════
struct Spring
    p1::Int16; p2::Int16
    length::SimFloat; axial_stiffness::SimFloat; axial_damping::SimFloat
end
Spring() = Spring(Int16(0), Int16(0), 0.0, 0.0, 0.0)

# ═══════════════════════════════════════════════════════════════════════════════
# Geometry helpers (from KiteUtils.jl)
# ═══════════════════════════════════════════════════════════════════════════════
function initial_kite_ref_frame(vec_c, v_app)
    z = normalize(vec_c)
    y = normalize(cross(v_app, vec_c))
    x = normalize(cross(y, vec_c))
    return (x, y, z)
end

function get_particles(height_k, height_b, width, m_k,
                       pos_pod  = [75.0, 0.0, 129.90381057],
                       vec_c    = [-15.0, 0.0, -25.98076211],
                       v_app    = [10.4855, 0.0, -3.08324])
    beta = π / 2.0
    x, y, z = initial_kite_ref_frame(vec_c, v_app)
    h_kx = height_k * cos(beta); h_kz = height_k * sin(beta)
    h_bx = height_b * cos(beta); h_bz = height_b * sin(beta)
    pos_kite = pos_pod - (h_kz + h_bz) * z + (h_kx + h_bx) * x
    pos_C    = pos_kite + h_kz * z + 0.5 * width * y + h_kx * x
    pos_A    = pos_kite + h_kz * z + (h_kx + width * m_k) * x
    pos_D    = pos_kite + h_kz * z - 0.5 * width * y + h_kx * x
    pos0     = pos_kite + (h_kz + h_bz) * z + (h_kx + h_bx) * x
    [zeros(3), pos0, pos_A, pos_kite, pos_C, pos_D]
end

function bridle_length(se::Settings)
    particles = get_particles(se.height_k, se.h_bridle, se.width, se.m_k)
    b = particles[2:end]  # [pos0, pos_A, pos_kite, pos_C, pos_D]
    norm(b[1]-b[2]) + norm(b[1]-b[4]) + norm(b[1]-b[5]) +
    norm(b[3]-b[2]) + norm(b[3]-b[4]) + norm(b[3]-b[5])
end

function rotate_around_y(vec, angle)
    T = typeof(angle)
    result = zeros(T, 3)
    result[1] =  cos(angle) * vec[1] + sin(angle) * vec[3]
    result[2] =  vec[2]
    result[3] = -sin(angle) * vec[1] + cos(angle) * vec[3]
    result
end

# ═══════════════════════════════════════════════════════════════════════════════
# KPS4 model struct
# ═══════════════════════════════════════════════════════════════════════════════
mutable struct KPS4
    set::Settings; kcu::KCU; am::AtmosphericModel; wm::AsyncMachine
    calc_cl::Spline1D; calc_cd::Spline1D
    iter::Int64
    v_wind::KVec3; v_wind_gnd::KVec3; v_wind_tether::KVec3; v_apparent::KVec3
    bridle_factor::SimFloat; side_cl::SimFloat
    drag_force::KVec3; side_force::KVec3; f_d::KVec3; f_s::KVec3
    ks::SimFloat; lift_force::KVec3; spring_force::KVec3; last_force::KVec3
    res1::SVector{NP, KVec3}; res2::SVector{NP, KVec3}
    pos::SVector{NP, KVec3}; vel::SVector{NP, KVec3}
    vel_kite::KVec3; segment_length::SimFloat
    param_cl::SimFloat; param_cd::SimFloat; psi::SimFloat
    alpha_depower::SimFloat; pitch::SimFloat; pitch_rate::SimFloat
    alpha_2::SimFloat; alpha_3::SimFloat; alpha_4::SimFloat; side_slip::SimFloat
    t_0::SimFloat; v_reel_out::SimFloat; last_v_reel_out::SimFloat
    l_tether::SimFloat; rho::SimFloat
    depower::SimFloat; steering::SimFloat; kcu_steering::SimFloat
    stiffness_factor::SimFloat
    initial_masses::MVector{NP, SimFloat}; masses::MVector{NP, SimFloat}
    springs::MVector{NS, Spring}; forces::SVector{NP, KVec3}
    sync_speed::Union{SimFloat, Nothing}; set_torque::Union{SimFloat, Nothing}
    x::KVec3; y::KVec3; z::KVec3
end

function new_svec_kvec3()
    SVector{NP}(ntuple(_ -> KVec3(0.0, 0.0, 0.0), NP))
end

function KPS4(kcu::KCU)
    set = kcu.set
    am = AtmosphericModel(set)
    wm = AsyncMachine(set)
    cl = Spline1D(set.alpha_cl, set.cl_list)
    cd = Spline1D(set.alpha_cd, set.cd_list)
    s = KPS4(
        set, kcu, am, wm, cl, cd, 0,
        KVec3(0,0,0), KVec3(0,0,0), KVec3(0,0,0), KVec3(0,0,0),  # winds
        1.0, 0.0,                                                   # bridle_factor, side_cl
        KVec3(0,0,0), KVec3(0,0,0), KVec3(0,0,0), KVec3(0,0,0),  # drag, side, f_d, f_s
        0.0, KVec3(0,0,0), KVec3(0,0,0), KVec3(0,0,0),           # ks, lift, spring, last
        new_svec_kvec3(), new_svec_kvec3(),                         # res1, res2
        new_svec_kvec3(), new_svec_kvec3(),                         # pos, vel
        KVec3(0,0,0), 0.0,                                         # vel_kite, segment_length
        0.2, 1.0, 0.0,                                             # param_cl, param_cd, psi
        0.0, 0.0, 0.0,                                             # alpha_depower, pitch, pitch_rate
        0.0, 0.0, 0.0, 0.0,                                       # alpha_2,3,4, side_slip
        0.0, 0.0, 0.0,                                             # t_0, v_reel_out, last_v_reel_out
        0.0, 0.0,                                                  # l_tether, rho
        0.0, 0.0, 0.0,                                             # depower, steering, kcu_steering
        1.0,                                                       # stiffness_factor
        ones(MVector{NP, SimFloat}), zeros(MVector{NP, SimFloat}), # initial_masses, masses
        MVector{NS}(ntuple(_ -> Spring(), NS)),                    # springs
        new_svec_kvec3(),                                          # forces
        0.0, nothing,                                              # sync_speed, set_torque
        KVec3(0,0,0), KVec3(0,0,0), KVec3(0,0,0),                # x, y, z
    )
    clear!(s)
    return s
end
KPS4(set::Settings) = KPS4(KCU(set))

# ═══════════════════════════════════════════════════════════════════════════════
# Initialization functions
# ═══════════════════════════════════════════════════════════════════════════════

# Bridle connection topology: [from_particle, to_particle, unstressed_length]
# -1 means computed from geometry
const SPRINGS_INPUT = [0.  1.  150.
                       1.  2.   -1.
                       4.  2.   -1.
                       4.  5.   -1.
                       3.  4.   -1.
                       5.  1.   -1.
                       4.  1.   -1.
                       3.  5.   -1.
                       5.  2.   -1.
                       2.  3.   -1.]

function init_springs!(s::KPS4)
    l_0 = s.set.l_tether / s.set.segments
    particles = get_particles(s.set.height_k, s.set.h_bridle, s.set.width, s.set.m_k)
    for j in 1:size(SPRINGS_INPUT, 1)
        if j == 1
            for i in 1:s.set.segments
                k = s.set.e_tether * (s.set.d_tether / 2000.0)^2 * π / l_0
                c = s.set.axial_damping / l_0
                s.springs[i] = Spring(Int16(i), Int16(i + 1), l_0, k, c)
            end
        else
            p0 = Int(SPRINGS_INPUT[j, 1]) + 1
            p1 = Int(SPRINGS_INPUT[j, 2]) + 1
            if SPRINGS_INPUT[j, 3] == -1
                l_b = norm(particles[p1] - particles[p0]) * PRE_STRESS
                k = s.set.e_tether * (s.set.d_line / 2000.0)^2 * π / l_b
                p0 += s.set.segments - 1
                p1 += s.set.segments - 1
                c = s.set.axial_damping / l_b
                s.springs[j + s.set.segments - 1] = Spring(Int16(p0), Int16(p1), l_b, k, c)
            end
        end
    end
    s.springs
end

function init_masses!(s::KPS4)
    s.masses .= 0.0
    l_0 = s.set.l_tether / s.set.segments
    for i in 1:s.set.segments
        m_half = 0.5 * l_0 * s.set.rho_tether * (s.set.d_tether / 2000.0)^2 * π
        s.masses[i]   += m_half
        s.masses[i+1] += m_half
    end
    s.masses[s.set.segments + 1] += s.set.kcu_mass
    k2 = s.set.rel_top_mass * (1.0 - s.set.rel_nose_mass)
    k3 = 0.5 * (1.0 - s.set.rel_top_mass) * (1.0 - s.set.rel_nose_mass)
    k4 = k3
    s.masses[s.set.segments + 2] += s.set.rel_nose_mass * s.set.mass
    s.masses[s.set.segments + 3] += k2 * s.set.mass
    s.masses[s.set.segments + 4] += k3 * s.set.mass
    s.masses[s.set.segments + 5] += k4 * s.set.mass
    s.masses
end

function init_pos_vel_acc(s::KPS4,
        X = zeros(2 * (s.set.segments + KITE_PARTICLES) + 1);
        old = false, delta = 0.0)
    seg = s.set.segments
    np  = seg + 1 + KITE_PARTICLES
    pos = SVector{np}(ntuple(_ -> KVec3(0.0, 0.0, 0.0), np))
    vel = MVector{np}(ntuple(_ -> KVec3(0.0, 0.0, 0.0), np))
    acc = SVector{np}(ntuple(_ -> KVec3(0.0, 0.0, 0.0), np))
    pos[1] .= [0.0, delta, 0.0]
    vel[1] .= [delta, delta, delta]
    acc[1] .= [delta, delta, delta]
    sin_el = sin(s.set.elevation / 180.0 * π)
    cos_el = cos(s.set.elevation / 180.0 * π)
    for i in 1:seg
        radius = -i * (s.set.l_tether / seg)
        pos[i+1] .= [-cos_el * radius + X[i], delta,
                      -sin_el * radius + X[seg + KITE_PARTICLES - 1 + i]]
        vel[i+1] .= [delta, delta, 0.0]
        acc[i+1] .= [delta, delta, -9.81]
    end
    vec_c = pos[seg] - pos[seg + 1]
    if old
        particles = get_particles(s.set.height_k, s.set.h_bridle, s.set.width, s.set.m_k)
    else
        particles = get_particles(s.set.height_k, s.set.h_bridle, s.set.width, s.set.m_k,
                                  pos[seg + 1],
                                  rotate_around_y(vec_c, deg2rad(KITE_ANGLE)),
                                  s.v_apparent)
    end
    j = 1
    for i in [1, 2, 3]   # p8=A, p9=B, p10=C
        pos[seg + 1 + i] .= particles[i + 2] .+
            [X[seg + j], 0.0, X[2 * seg + KITE_PARTICLES - 1 + j]]
        vel[seg + 1 + i] .= [delta, delta, delta]
        acc[seg + 1 + i] .= [delta, delta, -9.81]
        j += 1
    end
    # p11 = D (mirrors C in x,z; opposite y)
    acc[seg + 1 + 4] .= [delta, delta, -9.81]
    vel[seg + 1 + 4] .= [delta, delta, delta]
    pos[seg + 1 + 4][1] = pos[seg + 1 + 3][1]   # D.x = C.x
    pos[seg + 1 + 4][3] = pos[seg + 1 + 3][3]   # D.z = C.z
    pos[seg + 1 + 3][2] += X[end]                # Y of C
    pos[seg + 1 + 4][2]  = -pos[seg + 1 + 3][2]  # Y of D = -Y of C

    for i in eachindex(pos)
        s.pos[i] .= pos[i]
    end
    # Add reel-out velocity contribution
    for i in 2:seg + 1
        vel[i] .+= (pos[i + 1] - pos[i]) * (s.set.v_reel_out * (i - 1) / seg)
    end
    for i in seg + 2:seg + KITE_PARTICLES + 1
        vel[i] .+= vel[seg + 1]
    end
    pos, vel, acc
end

function set_initial_velocity!(s::KPS4)
    seg = s.set.segments
    for i in 2:seg + 1
        s.vel[i] .= (s.pos[i + 1] - s.pos[i]) * (s.set.v_reel_out * (i - 1) / seg)
    end
    for i in seg + 2:seg + KITE_PARTICLES + 1
        s.vel[i] .= s.vel[seg + 1]
    end
end

function init_inner(s::KPS4,
        X = zeros(2 * (s.set.segments + KITE_PARTICLES - 1) + 1);
        old = false, delta = 0.0)
    pos, vel, acc = init_pos_vel_acc(s, X; old, delta)
    vcat(pos[2:end], vel[2:end]), vcat(vel[2:end], acc[2:end])
end

function turn(res, upwind_dir)
    turnangle = upwind_dir + π / 2
    res2 = zeros(SimFloat, length(res))
    for i in 1:div(length(res), 3)
        x = res[3*(i-1)+1]; y = res[3*(i-1)+2]
        res2[3*(i-1)+1] = cos(turnangle) * x + sin(turnangle) * y
        res2[3*(i-1)+2] = cos(turnangle) * y - sin(turnangle) * x
        res2[3*(i-1)+3] = res[3*(i-1)+3]
    end
    res2
end

function init(s::KPS4,
        X = zeros(2 * (s.set.segments + KITE_PARTICLES - 1) + 1);
        old = false, delta = 0.0, upwind_dir = nothing)
    seg = s.set.segments; np = seg + KITE_PARTICLES
    res1_, res2_ = init_inner(s, X; old, delta)
    res1__ = reduce(vcat, res1_)
    res2__ = reduce(vcat, res2_)
    if !isnothing(upwind_dir)
        res1__ = turn(res1__, upwind_dir)
        res2__ = turn(res2__, upwind_dir)
    end
    set_initial_velocity!(s)
    for i in 1:div(length(res1__), 6)
        j = i + np
        res1__[3*(j-1)+1:3*j] .= s.vel[i+1]
        res2__[3*(i-1)+1:3*i] .= s.vel[i+1]
    end
    n_state = 6 * np + 2
    res1 = vcat(res1__, [s.l_tether, s.set.v_reel_out])
    res2 = vcat(res2__, [0.0, 0.0])
    MVector{n_state, SimFloat}(res1), MVector{n_state, SimFloat}(res2)
end

# ═══════════════════════════════════════════════════════════════════════════════
# clear! — reset model to initial state
# ═══════════════════════════════════════════════════════════════════════════════
function clear!(s::KPS4)
    s.t_0 = 0.0
    s.v_reel_out = s.set.v_reel_out
    s.last_v_reel_out = s.set.v_reel_out
    s.sync_speed = s.set.v_reel_out
    s.v_wind_gnd    .= [s.set.v_wind, 0, 0]
    s.v_wind_tether .= [s.set.v_wind, 0, 0]
    s.v_apparent    .= [s.set.v_wind, 0, 0]
    height = sin(deg2rad(s.set.elevation)) * s.set.l_tether
    s.v_wind .= s.v_wind_gnd * calc_wind_factor(s.am, height)
    s.l_tether = s.set.l_tether
    s.segment_length = s.l_tether / s.set.segments
    init_masses!(s)
    init_springs!(s)
    for i in 1:s.set.segments + KITE_PARTICLES + 1
        s.forces[i] .= zeros(3)
    end
    s.drag_force .= [0.0, 0, 0]
    s.lift_force .= [0.0, 0, 0]
    s.side_force .= [0.0, 0, 0]
    s.rho = s.set.rho_0
    s.bridle_factor = s.set.l_bridle / bridle_length(s.set)
    s.ks = deg2rad(s.set.max_steering)
    s.kcu.depower = s.set.depower / 100.0
    s.kcu.set_depower = s.kcu.depower
    # Simplified: skip orient_euler (pitch_rate=0 so cmq term vanishes)
    s.pitch = 0.0
    s.pitch_rate = 0.0
    set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu))
end

# ═══════════════════════════════════════════════════════════════════════════════
# set_depower_steering! / set_v_wind_ground!
# ═══════════════════════════════════════════════════════════════════════════════
function set_depower_steering!(s::KPS4, depower, steering)
    s.depower = depower
    s.kcu_steering = steering
    s.alpha_depower = calc_alpha_depower(s.kcu, depower)
    s.steering = (steering - s.set.c0) /
                 (1.0 + s.set.k_ds * (s.alpha_depower / deg2rad(s.set.alpha_d_max)))
    nothing
end

function set_v_wind_ground!(s::KPS4, height, v_wind_gnd = s.set.v_wind; upwind_dir = -π/2)
    height < 6.0 && (height = 6.0)
    wind_dir = -upwind_dir - π / 2
    s.v_wind .= v_wind_gnd * calc_wind_factor(s.am, height) .* [cos(wind_dir), sin(wind_dir), 0]
    s.v_wind_gnd .= [v_wind_gnd * cos(wind_dir), v_wind_gnd * sin(wind_dir), 0.0]
    s.v_wind_tether .= s.v_wind_gnd * calc_wind_factor(s.am, height / 2.0)
    s.rho = calc_rho(s.am, height)
    nothing
end

# ═══════════════════════════════════════════════════════════════════════════════
# Physics: force calculations
# ═══════════════════════════════════════════════════════════════════════════════

@inline function calc_particle_forces!(s::KPS4, pos1, pos2, vel1, vel2,
                                       spring, segments, d_tether, rho, i)
    l_0 = spring.length
    k = spring.axial_stiffness * s.stiffness_factor
    c = spring.axial_damping
    segment = pos1 - pos2
    rel_vel = vel1 - vel2
    av_vel = 0.5 * (vel1 + vel2)
    norm1 = norm(segment)
    unit_vector = segment / norm1

    k1 = s.set.rel_compr_stiffness * k
    k2_comp = 0.1 * k
    c1 = s.set.rel_damping * c
    spring_vel = unit_vector ⋅ rel_vel
    if (norm1 - l_0) > 0.0
        if i > segments
            s.spring_force .= (k * (norm1 - l_0) + c1 * spring_vel) * unit_vector
        else
            s.spring_force .= (k * (norm1 - l_0) + c * spring_vel) * unit_vector
        end
    elseif i > segments
        s.spring_force .= (k1 * (norm1 - l_0) + c * spring_vel) * unit_vector
    else
        s.spring_force .= (k2_comp * (norm1 - l_0) + c * spring_vel) * unit_vector
    end

    s.v_apparent .= s.v_wind_tether - av_vel
    v_app_kcu = s.v_wind_tether - vel2
    if s.set.version == 1
        area = norm1 * d_tether
    else
        if i > segments
            area = norm1 * s.set.d_line * 0.001 * s.bridle_factor
        else
            area = norm1 * d_tether
        end
    end

    v_app_perp = s.v_apparent - (s.v_apparent ⋅ unit_vector) * unit_vector
    half_drag = (-0.25 * rho * s.set.cd_tether * norm(v_app_perp) * area) * v_app_perp
    if i == segments
        v_app_perp_kcu = v_app_kcu - (v_app_kcu ⋅ unit_vector) * unit_vector
        kcu_area = π * (s.set.kcu_diameter / 2)^2
        kcu_drag = (-0.25 * rho * s.set.cd_kcu * norm(v_app_perp_kcu) * kcu_area) * v_app_perp_kcu
        s.forces[spring.p2] .+= kcu_drag
    end

    s.forces[spring.p1] .+= half_drag + s.spring_force
    s.forces[spring.p2] .+= half_drag - s.spring_force
    if i == 1; s.last_force .= s.forces[spring.p1]; end
    nothing
end

@inline function calc_aero_forces!(s::KPS4, pos, vel, rho, alpha_depower, rel_steering)
    seg = s.set.segments
    rel_side_area = s.set.rel_side_area / 100.0
    K = 1 - rel_side_area

    pos_B, pos_C, pos_D = pos[seg+3], pos[seg+4], pos[seg+5]
    v_A, v_B, v_C, v_D  = vel[seg+2], vel[seg+3], vel[seg+4], vel[seg+5]
    va_1 = s.v_wind - v_A; va_2 = s.v_wind - v_B
    va_3 = s.v_wind - v_C; va_4 = s.v_wind - v_D

    pos_centre = 0.5 * (pos_C + pos_D)
    delta_vec = pos_B - pos_centre
    z = -normalize(delta_vec)
    y = normalize(pos_C - pos_D)
    x = y × z
    s.x .= x; s.y .= y; s.z .= z

    va_xz2 = va_2 - (va_2 ⋅ y) * y
    va_xy3 = va_3 - (va_3 ⋅ z) * z
    va_xy4 = va_4 - (va_4 ⋅ z) * z
    va_xz1 = va_1 - (va_1 ⋅ y) * y

    R_k_w = SMatrix{3,3}([x y z])
    va_k = R_k_w' * SVector{3}(va_2)
    s.side_slip = atan(va_k[2], -va_k[1])

    alpha_2 = rad2deg(π - acos2(normalize(va_xz2) ⋅ x) - alpha_depower) + s.set.alpha_zero
    alpha_3 = rad2deg(π - acos2(normalize(va_xy3) ⋅ x) + rel_steering * s.ks) + s.set.alpha_ztip
    alpha_4 = rad2deg(π - acos2(normalize(va_xy4) ⋅ x) - rel_steering * s.ks) + s.set.alpha_ztip
    s.alpha_2 = alpha_2; s.alpha_3 = alpha_3; s.alpha_4 = alpha_4

    dc = s.set.version == 3 ? 1.0 : DRAG_CORR
    CL2, CD2 = s.calc_cl(alpha_2), dc * s.calc_cd(alpha_2)
    CL3, CD3 = s.calc_cl(alpha_3), dc * s.calc_cd(alpha_3)
    CL4, CD4 = s.calc_cl(alpha_4), dc * s.calc_cd(alpha_4)
    s.side_cl = CL4 - CL3
    L2 = (-0.5 * rho * norm(va_xz2)^2 * s.set.area * CL2) * normalize(va_2 × y)
    L3 = (-0.5 * rho * norm(va_xy3)^2 * s.set.area * rel_side_area * CL3) * normalize(va_3 × z)
    L4 = (-0.5 * rho * norm(va_xy4)^2 * s.set.area * rel_side_area * CL4) * normalize(z × va_4)
    D2 = (-0.5 * K * rho * norm(va_2) * s.set.area * CD2) * va_2
    D3 = (-0.5 * K * rho * norm(va_3) * s.set.area * rel_side_area * CD3) * va_3
    D4 = (-0.5 * K * rho * norm(va_4) * s.set.area * rel_side_area * CD4) * va_4

    s.lift_force .= L2
    if s.set.version == 3
        s.drag_force .= D2
        s.forces[seg + 3] .+= (L2 + D2 - D3 - D4)
    else
        s.drag_force .= D2 + D3 + D4
        s.forces[seg + 3] .+= (L2 + D2)
    end
    s.f_d .= (0.5 * rho * s.set.area * norm(va_xz1)^2 * s.set.cmq *
              s.pitch_rate * s.set.cord_length) * s.z
    s.f_s .= (0.5 * rho * s.set.area *
              (0.5 * (norm(va_xy3) + norm(va_xy4)))^2 * s.set.smc *
              rel_steering * s.ks) * s.x
    s.forces[seg + 2] .+= s.f_d
    s.forces[seg + 4] .+= (L3 + D3 - 0.5 * s.f_d + 0.5 * s.f_s)
    s.forces[seg + 5] .+= (L4 + D4 - 0.5 * s.f_d - 0.5 * s.f_s)
    s.side_force .= (L3 + L4)
end

@inline function inner_loop!(s::KPS4, pos, vel, v_wind_gnd, segments, d_tether)
    for i in eachindex(s.springs)
        p1 = s.springs[i].p1; p2 = s.springs[i].p2
        height = 0.5 * (pos[p1][3] + pos[p2][3])
        rho = calc_rho(s.am, height)
        @assert height > -1000
        height < 6 && (height = 6)
        s.v_wind_tether .= calc_wind_factor(s.am, height) * v_wind_gnd
        calc_particle_forces!(s, pos[p1], pos[p2], vel[p1], vel[p2],
                              s.springs[i], segments, d_tether, rho, i)
    end
    nothing
end

function loop!(s::KPS4, pos, vel, posd, veld)
    seg = s.set.segments
    L_0 = s.l_tether / seg
    mass_per_meter = (s.set.version == 1) ?
        0.011 : s.set.rho_tether * π * (s.set.d_tether / 2000.0)^2
    s.res1[1] .= pos[1]
    s.res2[1] .= vel[1]
    particles = seg + KITE_PARTICLES + 1
    for i in 2:particles
        s.res1[i] .= vel[i] - posd[i]
    end
    m_tether_particle = mass_per_meter * s.segment_length
    s.masses[seg + 1] = s.set.kcu_mass + 0.5 * m_tether_particle
    ad = s.set.axial_damping / L_0
    as = s.set.axial_stiffness / L_0
    for i in 1:seg
        s.masses[i] = m_tether_particle
        s.springs[i] = Spring(s.springs[i].p1, s.springs[i].p2, s.segment_length, as, ad)
    end
    inner_loop!(s, pos, vel, s.v_wind_gnd, seg, s.set.d_tether / 1000.0)
    for i in 2:particles
        s.res2[i] .= veld[i] - (SVector(0, 0, -G_EARTH) - s.forces[i] / s.masses[i])
    end
end

# ═══════════════════════════════════════════════════════════════════════════════
# DAE residual
# ═══════════════════════════════════════════════════════════════════════════════
function residual!(res, yd, y::Vector{SimFloat}, s::KPS4, time)
    S = length(y)
    residual!(res, MVector{S, SimFloat}(yd), MVector{S, SimFloat}(y), s, time)
end

function residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S
    T = S - 2
    segments = div(T, 6) - KITE_PARTICLES

    for i in 1:segments + KITE_PARTICLES + 1
        s.forces[i] .= SVector(0.0, 0, 0)
    end
    len, v_ro   = y[end-1],  y[end]
    lend, v_rod = yd[end-1], yd[end]
    y_  = @view y[1:end-2]
    yd_ = @view yd[1:end-2]
    part  = reshape(SVector{T}(y_),  Size(3, div(T,6), 2))
    partd = reshape(SVector{T}(yd_), Size(3, div(T,6), 2))
    pos1, vel1 = part[:,:,1], part[:,:,2]
    pos = SVector{div(T,6)+1}(
        i == 1 ? SVector(0.0, 0, 0) : SVector(pos1[:, i-1]) for i in 1:div(T,6)+1)
    vel = SVector{div(T,6)+1}(
        i == 1 ? SVector(0.0, 0, 0) : SVector(vel1[:, i-1]) for i in 1:div(T,6)+1)
    posd1, veld1 = partd[:,:,1], partd[:,:,2]
    posd = SVector{div(T,6)+1}(
        i == 1 ? SVector(0.0, 0, 0) : SVector(posd1[:, i-1]) for i in 1:div(T,6)+1)
    veld = SVector{div(T,6)+1}(
        i == 1 ? SVector(0.0, 0, 0) : SVector(veld1[:, i-1]) for i in 1:div(T,6)+1)
    @assert isfinite(pos[2][3])

    s.l_tether = len
    s.segment_length = len / segments
    calc_aero_forces!(s, pos, vel, s.rho, s.alpha_depower, s.steering)
    loop!(s, pos, vel, posd, veld)

    res[end-1] = lend - v_ro
    use_brake = true   # AsyncMachine + sync_speed == 0
    if !isnothing(s.sync_speed) && s.sync_speed == 0.0
        use_brake = true
    end
    res[end] = v_rod - calc_acceleration(s.wm, s.sync_speed, v_ro, norm(s.forces[1]), use_brake)

    for i in 2:div(T,6)+1
        for j in 1:3
            res[3*(i-2)+j]               = s.res1[i][j]
            res[3*(div(T,6))+3*(i-2)+j]  = s.res2[i][j]
        end
    end
    for i in 1:div(T,6)+1
        s.pos[i] .= pos[i]
        s.vel[i] .= vel[i]
    end
    s.vel_kite .= vel[end-2]
    s.v_reel_out = v_ro
    @assert !isnan(norm(res))
    s.iter += 1
    nothing
end

# ═══════════════════════════════════════════════════════════════════════════════
# Steady state solver
# ═══════════════════════════════════════════════════════════════════════════════
function make_jac(f!, n_vars)
    _F1 = zeros(SimFloat, n_vars)
    _F2 = zeros(SimFloat, n_vars)
    _xp = zeros(SimFloat, n_vars)
    _xm = zeros(SimFloat, n_vars)
    threshold = sqrt(floatmin(SimFloat))
    function jac!(J, x)
        h_factor = cbrt(eps(SimFloat))
        for j in 1:n_vars
            copyto!(_xp, x); copyto!(_xm, x)
            h = max(abs(x[j]), one(SimFloat)) * h_factor
            _xp[j] += h; _xm[j] -= h
            f!(_F1, _xp); f!(_F2, _xm)
            @views J[:, j] .= (_F1 .- _F2) ./ (2h)
        end
        @. J = ifelse(abs(J) < threshold, zero(SimFloat), J)
    end
    return jac!
end

function find_steady_state!(s::KPS4; prn=false, delta=0.001,
                            stiffness_factor=0.035, upwind_dir=-π/2)
    set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir=-π/2)
    s.stiffness_factor = stiffness_factor
    seg = s.set.segments
    n_state = 6 * (seg + KITE_PARTICLES) + 2
    res = zeros(MVector{n_state, SimFloat})
    iter = 0

    function test_initial_condition!(F, x::Vector)
        x1 = copy(x)
        y0, yd0 = init(s, x1; delta)
        try
            residual!(res, yd0, y0, s, 0.0)
        catch _
            @warn "Warning in test_initial_condition!"
            F .= 1.0e6
            iter += 1
            return nothing
        end
        np = seg + KITE_PARTICLES
        for i in 1:np - 1
            j = (i != np - 1) ? i : i + 1
            F[i]      = res[1 + 3*(j-1) + 3*np]  # x-accel
            F[i + np] = res[3 + 3*(j-1) + 3*np]  # z-accel
        end
        i_kcu = seg + 1
        F[end-1] = res[1 + 3*(i_kcu-1) + 3*np]   # KCU x-accel
        i_c = seg + 3
        F[end]   = res[2 + 3*(i_c-1) + 3*np]     # C y-accel
        iter += 1
        return nothing
    end

    n_unknowns = 2 * (seg + KITE_PARTICLES - 1) + 2
    if prn; println("\nStarted function test_nlsolve..."); end
    X00 = zeros(SimFloat, n_unknowns)
    jac! = make_jac(test_initial_condition!, n_unknowns)
    results = nlsolve(test_initial_condition!, jac!, X00;
                      autoscale=true, xtol=4e-7, ftol=4e-7, iterations=s.set.max_iter)
    if prn; println("\nresult: $results"); end
    if !NLsolve.converged(results)
        @warn "find_steady_state!: solver did not converge!"
    end
    y0, yd0 = init(s, results.zero; upwind_dir)
    set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir)
    residual!(res, yd0, y0, s, 0.0)
    y0, yd0
end

# ═══════════════════════════════════════════════════════════════════════════════
# Getter functions
# ═══════════════════════════════════════════════════════════════════════════════
calc_height(s::KPS4) = s.pos[end-2][3]
pos_kite(s::KPS4)    = s.pos[end-2]
unstretched_length(s::KPS4) = s.l_tether

function tether_length(s::KPS4)
    len = 0.0
    for i in 1:s.set.segments
        len += norm(s.pos[i+1] - s.pos[i])
    end
    len
end

function spring_forces(s::KPS4)
    seg = s.set.segments
    forces = zeros(SimFloat, seg + KITE_SPRINGS)
    for i in 1:seg
        forces[i] = s.springs[i].axial_stiffness *
                    (norm(s.pos[i+1] - s.pos[i]) - s.segment_length) * s.stiffness_factor
    end
    for i in 1:KITE_SPRINGS
        p1 = s.springs[i + seg].p1; p2 = s.springs[i + seg].p2
        sp = s.springs[i + seg]
        segment = s.pos[p1] - s.pos[p2]
        n1 = norm(segment)
        k = sp.axial_stiffness * s.stiffness_factor
        k1 = 0.25 * k
        forces[i + seg] = (n1 - sp.length) > 0 ? k * (n1 - sp.length) : k1 * (n1 - sp.length)
    end
    forces
end

function calc_pre_tension(s::KPS4)
    forces = spring_forces(s)
    av = sum(forces[1:s.set.segments]) / s.set.segments
    r = av / s.set.axial_stiffness
    r < 0 && (r = 0.0)
    isnan(r) && (r = 0.0)
    return r + 1.0
end

# ═══════════════════════════════════════════════════════════════════════════════
# Run the test
# ═══════════════════════════════════════════════════════════════════════════════
set = deepcopy(default_settings())

function init_392!(s, l, set)
    clear!(s)
    s.set.l_tethers[1] = l
    s.set.elevation = 70.0
    s.set.area = 10.0
    s.set.rel_side_area = 50.0
    s.set.v_wind = 9.1
    s.set.mass = 6.2
    s.set.c_s = 0.6
end

@testset "test_find_steady_state (MWE)" begin
    l_tethers = [100.0, 200.0, 392.0]
    for l in l_tethers
        kps4 = KPS4(KCU(deepcopy(set)))
        init_392!(kps4, l, set)
        clear!(kps4)
        set_depower_steering!(kps4, kps4.set.depower_offset / 100.0, 0.0)
        height = sin(deg2rad(kps4.set.elevation)) * kps4.set.l_tether
        kps4.v_wind .= kps4.v_wind_gnd * calc_wind_factor(kps4.am, height)

        res1, res2 = find_steady_state!(kps4; delta=0.001, stiffness_factor=0.035, prn=false)

        @test sum(res2) ≈ -9.81 * (set.segments + KITE_PARTICLES)
        pre_tension = calc_pre_tension(kps4)
        @test pre_tension > 1.0001
        @test pre_tension < 1.01
        @test unstretched_length(kps4) ≈ l
        @test isapprox(tether_length(kps4), 1.008954l, rtol=1e-2)
        println("  l=$l: tether_length=$(tether_length(kps4)), pre_tension=$pre_tension ✓")
    end
end

[ChrisRackauckas-Claude]

Resolution: NonlinearSolve TrustRegion works — and outperforms NLsolve

After testing both the MWE and the original KiteModels.jl code (bestie branch), TrustRegion() with the custom Jacobian passes all 15 tests — compared to NLsolve's 14/15.

Test Results: Original KiteModels.jl

l_tether	Retcode	Residual Norm	pre_tension	tether_len	tether_err	Result
100	Success	6.4e-8	1.000382	101.091	0.194%	PASS
200	Success	2.2e-9	1.000470	202.686	0.444%	PASS
392	Success	3.5e-10	1.000638	399.150	0.920%	PASS

15/15 pass vs NLsolve's 14/15 (NLsolve fails the l=392 tether_length test with 1.35% error, exceeding the 1% rtol threshold).

Head-to-Head: NLsolve vs NonlinearSolve TrustRegion (l=392)

Metric	NLsolve	NonlinearSolve TrustRegion
Tether error	1.35% (FAIL)	0.92% (PASS)
Residual norm	(not reported)	3.5e-10
Runtime (post-JIT)	0.154s	0.031s

Code Change Required

The swap in KPS4.jl is straightforward:

Before (NLsolve):

using NLsolve

jac! = make_jac(test_initial_condition!, length(X00))
results = nlsolve(test_initial_condition!, jac!, X00,
                  autoscale=true, xtol=4e-7, ftol=4e-7, iterations=s.set.max_iter)
if !converged(results)
    @warn "find_steady_state! did not converge"
end
y0, yd0 = init(s, results.zero; upwind_dir)

After (NonlinearSolve):

using NonlinearSolve

jac_2arg! = make_jac(test_initial_condition!, length(X00))

# Wrap for NonlinearSolve's 3-arg (out, u, p) convention
f_nl!(F, x, p) = test_initial_condition!(F, x)
jac_nl!(J, x, p) = jac_2arg!(J, x)

nf = NonlinearFunction(f_nl!; jac=jac_nl!)
prob = NonlinearProblem(nf, X00)
sol = solve(prob, TrustRegion(); abstol=4e-7, reltol=4e-7, maxiters=s.set.max_iter)
if !SciMLBase.successful_retcode(sol.retcode)
    @warn "find_steady_state! did not converge: $(sol.retcode)"
end
y0, yd0 = init(s, sol.u; upwind_dir)

Key Findings

1. autoscale=true is NOT needed. TrustRegion's built-in trust region management provides equivalent or better robustness without NLsolve's Moré/MINPACK diagonal scaling.

2. The custom make_jac Jacobian IS still required. Three reasons:

• NonlinearSolve's AutoForwardDiff fails because the model's MVector{3, Float64} internals can't hold ForwardDiff.Dual numbers (MethodError: no method matching Float64(::ForwardDiff.Dual{...}))
• NonlinearSolve's AutoFiniteDiff uses different step sizing and doesn't flush subnormals, causing convergence failures on l=200 and l=392 (hits MaxIters with large residual norms)
• The custom Jacobian's step size (h = max(abs(x[j]), 1.0) * cbrt(eps(Float64))) and subnormal flushing (@. J = ifelse(abs(J) < sqrt(floatmin(Float64)), zero(Float64), J)) are tuned for this physics model

3. NewtonRaphson() fails — converges to a wrong solution. This is expected; without trust region globalization, Newton's method can't handle the problem's basins of attraction.

4. KiteModels.jl already depends on NonlinearSolve (via OrdinaryDiffEqNonlinearSolve), so this swap doesn't add a new dependency.

API Migration Cheat Sheet

NLsolve	NonlinearSolve	Notes
f!(F, x)	f!(F, u, p)	Add unused p parameter
jac!(J, x)	jac!(J, u, p)	Add unused p parameter
ftol=4e-7	abstol=4e-7	Direct equivalent
xtol=4e-7	reltol=4e-7	Approximate equivalent
iterations=200	maxiters=200	Direct equivalent
autoscale=true	(not needed)	TrustRegion works without it here
NLsolve.converged(r)	SciMLBase.successful_retcode(sol)
results.zero	sol.u

Solver Comparison Summary

Solver	l=100	l=200	l=392	Notes
NLsolve (trust_region + autoscale)	PASS	PASS	FAIL (tether_length)	14/15
NonlinearSolve TrustRegion + custom jac	PASS	PASS	PASS	15/15
NonlinearSolve TrustRegion + AutoFiniteDiff	PASS	FAIL (MaxIters)	FAIL (MaxIters)	5/15
NonlinearSolve NewtonRaphson + custom jac	FAIL (wrong solution)	FAIL	FAIL	0/15

Bottom Line

The migration from NLsolve to NonlinearSolve TrustRegion is straightforward and improves results. The only requirement is retaining the custom finite-difference Jacobian (adapted to 3-arg signatures). No autoscaling support is needed for this use case.

[ChrisRackauckas-Claude]

Closing — NonlinearSolve.jl's TrustRegion() works out of the box for this use case (with the custom Jacobian retained). See the detailed findings above and the migration PR: OpenSourceAWE/KiteModels.jl#261