From 5f38f464f1385d269a930c961c1c613d920aa37b Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Fri, 11 Apr 2025 14:05:17 +0200 Subject: [PATCH 01/25] Add KPS5 --- src/KPS5.jl | 478 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 478 insertions(+) create mode 100644 src/KPS5.jl diff --git a/src/KPS5.jl b/src/KPS5.jl new file mode 100644 index 000000000..1c230aba5 --- /dev/null +++ b/src/KPS5.jl @@ -0,0 +1,478 @@ +#= MIT License + +Copyright (c) 2020, 2021, 2022, 2024 Uwe Fechner + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. =# + +#= Model of a kite-power system in implicit form: residual = f(y, yd) + +This model implements a 3D mass-spring system with reel-out. It uses six tether segments (the number can be +configured in the file data/settings.yaml). The kite is modelled using 4 point masses and 3 aerodynamic +surfaces. The spring constant and the damping decrease with the segment length. The aerodynamic kite forces +are acting on three of the four kite point masses. + +Four point kite model, included from KiteModels.jl. + +Scientific background: http://arxiv.org/abs/1406.6218 =# + +# Array of connections of bridlepoints. +# First point, second point, unstressed length. +const SPRINGS_INPUT = [0. 1. 150. + 1. 2. -1. # s1, p7, p8 + 4. 2. -1. # s2, p10, p8 + 4. 5. -1. # s3, p10, p11 + 3. 4. -1. # s4, p9, p10 + 5. 1. -1. # s5, p11, p7 + 4. 1. -1. # s6, p10, p7 + 3. 5. -1. # s7, p9, p11 + 5. 2. -1. # s8, p11, p8 + 2. 3. -1.] # s9, p8, p9 + +# KCU = p7, A = p8, B = p9, C = p10, D = p11 + +# struct, defining the phyical parameters of one spring +@with_kw struct Spring{I, S} + p1::I = 1 # number of the first point + p2::I = 2 # number of the second point + length::S = 1.0 # current unstressed spring length + c_spring::S = 1.0 # spring constant [N/m] + damping::S = 0.1 # damping coefficent [Ns/m] +end + +const SP = Spring{Int16, SimFloat} +const KITE_PARTICLES = 4 +const KITE_SPRINGS = 9 +const DRAG_CORR = 0.93 # correction of the drag for the 4-point model +function zero(::Type{SP}) + SP(0,0,0,0,0) +end + +""" + mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel + +State of the kite power system, using a 4 point kite model. Parameters: +- S: Scalar type, e.g. SimFloat + In the documentation mentioned as Any, but when used in this module it is always SimFloat and not Any. +- T: Vector type, e.g. MVector{3, SimFloat} +- P: number of points of the system, segments+1 +- Q: number of springs in the system, P-1 +- SP: struct type, describing a spring +Normally a user of this package will not have to access any of the members of this type directly, +use the input and output functions instead. + +$(TYPEDFIELDS) +""" +@with_kw mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel + "Reference to the settings struct" + set::Settings + "Reference to the KCU model (Kite Control Unit as implemented in the package KitePodModels" + kcu::KCU + "Reference to the atmospheric model as implemented in the package AtmosphericModels" + am::AtmosphericModel = AtmosphericModel() + "Reference to winch model as implemented in the package WinchModels" + wm::AbstractWinchModel + "Iterations, number of calls to the function residual!" + iter:: Int64 = 0 + "Function for calculation the lift coefficent, using a spline based on the provided value pairs." + calc_cl::Spline1D + "Function for calculation the drag coefficent, using a spline based on the provided value pairs." + calc_cd::Spline1D + "wind vector at the height of the kite" + v_wind::T = zeros(S, 3) + "wind vector at reference height" + v_wind_gnd::T = zeros(S, 3) + "wind vector used for the calculation of the tether drag" + v_wind_tether::T = zeros(S, 3) + "apparent wind vector at the kite" + v_apparent::T = zeros(S, 3) + "bridle_factor = set.l_bridle/bridle_length(set)" + bridle_factor::S = 1.0 + "side lift coefficient, the difference of the left and right lift coefficients" + side_cl::S = 0.0 + "drag force of kite and bridle; output of calc_aero_forces!" + drag_force::T = zeros(S, 3) + "side_force acting on the kite" + side_force::T = zeros(S, 3) + "max_steering angle in radian" + ks::S = 0.0 + "lift force of the kite; output of calc_aero_forces!" + lift_force::T = zeros(S, 3) + "spring force of the current tether segment, output of calc_particle_forces!" + spring_force::T = zeros(S, 3) + "last winch force" + last_force::T = zeros(S, 3) + "a copy of the residual one (pos,vel) for debugging and unit tests" + res1::SVector{P, KVec3} = zeros(SVector{P, KVec3}) + "a copy of the residual two (vel,acc) for debugging and unit tests" + res2::SVector{P, KVec3} = zeros(SVector{P, KVec3}) + "a copy of the actual positions as output for the user" + pos::SVector{P, KVec3} = zeros(SVector{P, KVec3}) + "a copy of the actual velocities as output for the user" + vel::SVector{P, KVec3} = zeros(SVector{P, KVec3}) + "velocity vector of the kite" + vel_kite::T = zeros(S, 3) + "unstressed segment length [m]" + segment_length::S = 0.0 + "lift coefficient of the kite, depending on the angle of attack" + param_cl::S = 0.2 + "drag coefficient of the kite, depending on the angle of attack" + param_cd::S = 1.0 + "azimuth angle in radian; initial value is zero" + psi::S = zero(S) + "depower angle [deg]" + alpha_depower::S = 0.0 + "pitch angle [rad]" + pitch::S = 0.0 + "pitch rate [rad/s]" + pitch_rate::S = 0.0 + "aoa at particle B" + alpha_2::S = 0.0 + "aoa at particle B, corrected formula" + alpha_2b::S = 0.0 + "aoa at particle C" + alpha_3::S = 0.0 + alpha_3b::S = 0.0 + "aoa at particle D" + alpha_4::S = 0.0 + alpha_4b::S = 0.0 + "relative start time of the current time interval" + t_0::S = 0.0 + "reel out speed of the winch" + v_reel_out::S = 0.0 + "reel out speed at the last time step" + last_v_reel_out::S = 0.0 + "unstretched tether length" + l_tether::S = 0.0 + "air density at the height of the kite" + rho::S = 0.0 + "actual relative depower setting, must be between 0 .. 1.0" + depower::S = 0.0 + "actual relative steering setting, must be between -1.0 .. 1.0" + steering::S = 0.0 + "steering after the kcu, before applying offset and depower sensitivity, -1.0 .. 1.0" + kcu_steering::S = 0.0 + "multiplier for the stiffniss of tether and bridle" + stiffness_factor::S = 1.0 + "initial masses of the point masses" + initial_masses::MVector{P, S} = ones(P) + "current masses, depending on the total tether length" + masses::MVector{P, S} = zeros(P) + "vector of the springs, defined as struct" + springs::MVector{Q, SP} = zeros(SP, Q) + "vector of the forces, acting on the particles" + forces::SVector{P, KVec3} = zeros(SVector{P, KVec3}) + "synchronous speed of the motor/ generator" + sync_speed::Union{S, Nothing} = 0.0 + "set_torque of the motor/generator" + set_torque::Union{S, Nothing} = nothing + "set value of the force at the winch, for logging only" + set_force::Union{S, Nothing} = nothing + "set value of the bearing angle in radians, for logging only" + bearing::Union{S, Nothing} = nothing + "coordinates of the attractor point [azimuth, elevation] in radian, for logging only" + attractor::Union{SVector{2, S}, Nothing} = nothing + "x vector of kite reference frame" + x::T = zeros(S, 3) + "y vector of kite reference frame" + y::T = zeros(S, 3) + "z vector of kite reference frame" + z::T = zeros(S, 3) +end + +""" + clear!(s::KPS5) + +Initialize the kite power model. +""" +function clear!(s::KPS5) + s.t_0 = 0.0 # relative start time of the current time interval + 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] # wind vector at reference height + 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 + KiteModels.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 + roll, pitch, yaw = orient_euler(s) + s.pitch = pitch + s.pitch_rate = 0.0 + KiteModels.set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu)) +end + +function KPS5(kcu::KCU) + if kcu.set.winch_model == "AsyncMachine" + wm = AsyncMachine(kcu.set) + elseif kcu.set.winch_model == "TorqueControlledMachine" + wm = TorqueControlledMachine(kcu.set) + end + # wm.last_set_speed = kcu.set.v_reel_out + s = KPS5{SimFloat, KVec3, kcu.set.segments+KITE_PARTICLES+1, kcu.set.segments+KITE_SPRINGS, SP}(set=kcu.set, + kcu=kcu, wm=wm, calc_cl = Spline1D(kcu.set.alpha_cl, kcu.set.cl_list), + calc_cd=Spline1D(kcu.set.alpha_cd, kcu.set.cd_list) ) + clear!(s) + return s +end + +""" + calc_particle_forces!(s::KPS5, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i) + +Calculate the drag force of the tether segment, defined by the parameters pos1, pos2, vel1 and vel2 +and distribute it equally on the two particles, that are attached to the segment. +The result is stored in the array s.forces. +""" +@inline function calc_particle_forces!(s::KPS5, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i) + l_0 = spring.length # Unstressed length + k = spring.c_spring * s.stiffness_factor # Spring constant + c = spring.damping # Damping coefficient + 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 # compression stiffness kite springs + k2 = 0.1 * k # compression stiffness tether springs + c1 = s.set.rel_damping * c # damping kite springs + spring_vel = unit_vector ⋅ rel_vel + if (norm1 - l_0) > 0.0 + if i > segments # kite springs + 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 # kite spring + s.spring_force .= (k1 * (norm1 - l_0) + (c * spring_vel)) * unit_vector + else + s.spring_force .= (k2 * (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 # 6.0 = A_real/A_simulated + else + area = norm1 * d_tether + end + end + + v_app_perp = s.v_apparent - s.v_apparent ⋅ unit_vector * unit_vector + half_drag_force = (-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_force = (-0.25 * rho * s.set.cd_kcu * norm(v_app_perp_kcu) * kcu_area) * v_app_perp_kcu + @inbounds s.forces[spring.p2] .+= kcu_drag_force + end + + @inbounds s.forces[spring.p1] .+= half_drag_force + s.spring_force + @inbounds s.forces[spring.p2] .+= half_drag_force - s.spring_force + if i == 1 s.last_force .= s.forces[spring.p1] end + nothing +end + +""" + calc_aero_forces!(s::KPS5, pos, vel, rho, alpha_depower, rel_steering) + +Calculates the aerodynamic forces acting on the kite particles. + +Parameters: +- pos: vector of the particle positions +- vel: vector of the particle velocities +- rho: air density [kg/m^3] +- rel_depower: value between 0.0 and 1.0 +- alpha_depower: depower angle [degrees] +- rel_steering: value between -1.0 and +1.0 + +Updates the vector s.forces of the first parameter. +""" +function calc_aero_forces!(s::KPS5, pos, vel, rho, alpha_depower, rel_steering) +end + + +# =================== getter functions ==================================================== + +""" + calc_height(s::KPS5) + +Determine the height of the topmost kite particle above ground. +""" +function calc_height(s::KPS5) + pos_kite(s)[3] +end + +""" + pos_kite(s::KPS5) + +Return the position of the kite (top particle). +""" +function pos_kite(s::KPS5) + s.pos[end-2] +end + +""" + kite_ref_frame(s::KPS5; one_point=false) + +Returns a tuple of the x, y, and z vectors of the kite reference frame. +""" +function kite_ref_frame(s::KPS5; one_point=false) + if one_point + c = s.z + y = normalize(s.v_apparent × c) + x = normalize(y × c) + return x, y, c + else + return s.x, s.y, s.z + end +end + +""" + winch_force(s::KPS5) + +Return the absolute value of the force at the winch as calculated during the last timestep. +""" +function winch_force(s::KPS5) norm(s.last_force) end + +""" + cl_cd(s::KPS5) + +Calculate the lift and drag coefficients of the kite, based on the current angles of attack. +""" +function cl_cd(s::KPS5) + rel_side_area = s.set.rel_side_area/100.0 # defined in percent + K = 1 - rel_side_area # correction factor for the drag + if s.set.version == 3 + drag_corr = 1.0 + else + drag_corr = DRAG_CORR + end + CL2, CD2 = s.calc_cl(s.alpha_2), drag_corr * s.calc_cd(s.alpha_2) + CL3, CD3 = s.calc_cl(s.alpha_3), drag_corr * s.calc_cd(s.alpha_3) + CL4, CD4 = s.calc_cl(s.alpha_4), drag_corr * s.calc_cd(s.alpha_4) + if s.set.version == 3 + return CL2, CD2 + else + return CL2, K*(CD2+CD3+CD4) + end +end + +# ==================== end of getter functions ================================================ + +function spring_forces(s::KPS5; prn=true) + forces = zeros(SimFloat, s.set.segments+KITE_SPRINGS) + for i in 1:s.set.segments + forces[i] = s.springs[i].c_spring * (norm(s.pos[i+1] - s.pos[i]) - s.segment_length) * s.stiffness_factor + if forces[i] > s.set.max_force && prn + println("Tether raptures for segment $i !") + end + end + for i in 1:KITE_SPRINGS + p1 = s.springs[i+s.set.segments].p1 # First point nr. + p2 = s.springs[i+s.set.segments].p2 # Second point nr. + pos1, pos2 = s.pos[p1], s.pos[p2] + spring = s.springs[i+s.set.segments] + l_0 = spring.length # Unstressed length + k = spring.c_spring * s.stiffness_factor # Spring constant + segment = pos1 - pos2 + norm1 = norm(segment) + k1 = 0.25 * k # compression stiffness kite segments + if (norm1 - l_0) > 0.0 + spring_force = k * (norm1 - l_0) + else + spring_force = k1 * (norm1 - l_0) + end + forces[i+s.set.segments] = spring_force + if norm(s.spring_force) > 4000.0 + println("Bridle brakes for spring $i !") + end + end + forces +end + +""" + find_steady_state!(s::KPS5; prn=false, delta = 0.01, stiffness_factor=0.035, upwind_dir=-pi/2)) + +Find an initial equilibrium, based on the initial parameters +`l_tether`, elevation and `v_reel_out`. +""" +function find_steady_state!(s::KPS5; prn=false, delta = 0.01, stiffness_factor=0.035, upwind_dir=-pi/2) + set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir=-pi/2) + s.stiffness_factor = stiffness_factor + res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) + iter = 0 + + # helper function for the steady state finder + 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 e + println("Warning in test_initial_condition!") + end + for i in 1:s.set.segments+KITE_PARTICLES-1 + if i != s.set.segments+KITE_PARTICLES-1 + j = i + else + j = i + 1 + end + # copy the x-component of the residual res2 (acceleration) + F[i] = res[1 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)] + # copy the z-component of the residual res2 + F[i+s.set.segments+KITE_PARTICLES] = res[3 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)] + end + # copy the acceleration of point KCU in x direction + i = s.set.segments+1 + F[end-1] = res[1 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] + # copy the acceleration of point C in y direction + i = s.set.segments+3 + x = res[1 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] + y = res[2 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] + F[end] = res[2 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] + iter += 1 + return nothing + end + if prn println("\nStarted function test_nlsolve...") end + X00 = zeros(SimFloat, 2*(s.set.segments+KITE_PARTICLES-1)+2) + results = nlsolve(test_initial_condition!, X00, autoscale=true, xtol=4e-7, ftol=4e-7, iterations=s.set.max_iter) + if prn println("\nresult: $results") 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 From d57eec9ba8412dcfaaed6daf2c82aa1013b785e7 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Fri, 11 Apr 2025 14:09:35 +0200 Subject: [PATCH 02/25] Bugfixes --- src/KPS5.jl | 37 ++++++++++--------------------------- src/KiteModels.jl | 3 ++- 2 files changed, 12 insertions(+), 28 deletions(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index 1c230aba5..df7dff65a 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -33,36 +33,19 @@ Scientific background: http://arxiv.org/abs/1406.6218 =# # Array of connections of bridlepoints. # First point, second point, unstressed length. -const SPRINGS_INPUT = [0. 1. 150. - 1. 2. -1. # s1, p7, p8 - 4. 2. -1. # s2, p10, p8 - 4. 5. -1. # s3, p10, p11 - 3. 4. -1. # s4, p9, p10 - 5. 1. -1. # s5, p11, p7 - 4. 1. -1. # s6, p10, p7 - 3. 5. -1. # s7, p9, p11 - 5. 2. -1. # s8, p11, p8 - 2. 3. -1.] # s9, p8, p9 +const SPRINGS_INPUT_5P = [0. 1. 150. + 1. 2. -1. # s1, p7, p8 + 4. 2. -1. # s2, p10, p8 + 4. 5. -1. # s3, p10, p11 + 3. 4. -1. # s4, p9, p10 + 5. 1. -1. # s5, p11, p7 + 4. 1. -1. # s6, p10, p7 + 3. 5. -1. # s7, p9, p11 + 5. 2. -1. # s8, p11, p8 + 2. 3. -1.] # s9, p8, p9 # KCU = p7, A = p8, B = p9, C = p10, D = p11 -# struct, defining the phyical parameters of one spring -@with_kw struct Spring{I, S} - p1::I = 1 # number of the first point - p2::I = 2 # number of the second point - length::S = 1.0 # current unstressed spring length - c_spring::S = 1.0 # spring constant [N/m] - damping::S = 0.1 # damping coefficent [Ns/m] -end - -const SP = Spring{Int16, SimFloat} -const KITE_PARTICLES = 4 -const KITE_SPRINGS = 9 -const DRAG_CORR = 0.93 # correction of the drag for the 4-point model -function zero(::Type{SP}) - SP(0,0,0,0,0) -end - """ mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel diff --git a/src/KiteModels.jl b/src/KiteModels.jl index c3d020146..14ac39899 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -56,7 +56,7 @@ using ADTypes: AutoFiniteDiff using UnPack import ModelingToolkit.SciMLBase: successful_retcode -export KPS3, KPS4, RamAirKite, KVec3, SimFloat, Measurement, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG # constants and types +export KPS3, KPS4, KPS5, RamAirKite, KVec3, SimFloat, Measurement, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG # constants and types export calc_set_cl_cd!, copy_examples, copy_bin, update_sys_state! # helper functions export clear!, find_steady_state!, residual! # low level workers export init_sim!, init!, reinit!, next_step!, init_pos_vel, init_pos, model! # high level workers @@ -119,6 +119,7 @@ function __init__() end include("KPS4.jl") # include code, specific for the four point kite model +include("KPS5.jl") # include code, specific for the five point kite model include("ram_air_kite.jl") # include code, specific for the four point 3 line kite model include("mtk_model.jl") include("KPS3.jl") # include code, specific for the one point kite model From b2ff2c8c48170edaedb34fe619eae55ddba24410 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Fri, 11 Apr 2025 14:25:30 +0200 Subject: [PATCH 03/25] Add test_kps5.jl --- src/KPS5.jl | 4 +-- test/runtests.jl | 23 ++++++------ test/test_kps5.jl | 89 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 103 insertions(+), 13 deletions(-) create mode 100644 test/test_kps5.jl diff --git a/src/KPS5.jl b/src/KPS5.jl index df7dff65a..74ea647fe 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -196,8 +196,8 @@ function clear!(s::KPS5) s.l_tether = s.set.l_tether s.segment_length = s.l_tether / s.set.segments - init_masses!(s) - init_springs!(s) + # init_masses!(s) + # init_springs!(s) for i in 1:s.set.segments + KiteModels.KITE_PARTICLES + 1 s.forces[i] .= zeros(3) end diff --git a/test/runtests.jl b/test/runtests.jl index a64d776bf..e5be9f77b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,16 +17,17 @@ end::Bool cd("..") KiteUtils.set_data_path("") @testset verbose = true "Testing KiteModels..." begin - include("test_orientation.jl") - include("test_kps3.jl") + # include("test_orientation.jl") + # include("test_kps3.jl") include("test_kps4.jl") - if build_is_production_build - include("bench3.jl") - include("bench4.jl") - end - if ! haskey(ENV, "NO_MTK") - include("test_ram_air_kite.jl") - end - include("test_helpers.jl") - include("test_inertia_calculation.jl") + include("test_kps5.jl") + # if build_is_production_build + # include("bench3.jl") + # include("bench4.jl") + # end + # if ! haskey(ENV, "NO_MTK") + # include("test_ram_air_kite.jl") + # end + # include("test_helpers.jl") + # include("test_inertia_calculation.jl") end \ No newline at end of file diff --git a/test/test_kps5.jl b/test/test_kps5.jl new file mode 100644 index 000000000..c3aef27f9 --- /dev/null +++ b/test/test_kps5.jl @@ -0,0 +1,89 @@ +using Pkg +if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() + Pkg.update() +end +using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils +using KiteModels, KitePodModels + +set_data_path(joinpath(dirname(dirname(pathof(KiteModels))), "data")) +set = deepcopy(load_settings("system.yaml")) +kcu::KCU = KCU(set) +kps5::KPS5 = KPS5(kcu) + +pos, vel = nothing, nothing +poss, vels = nothing, nothing + +@testset verbose = true "KPS5 tests...." begin + +function set_defaults() + KiteModels.clear!(kps5) + kps5.set.l_tether = 150.0 + kps5.set.elevation = 60.0 + kps5.set.area = 20.0 + kps5.set.rel_side_area = 50.0 + kps5.set.v_wind = 8.0 + kps5.set.mass = 11.4 + kps5.set.damping = 2 * 473.0 + kps5.set.alpha = 1.0/7 + kps5.set.c_s = 0.6 + kps5.set.kcu_diameter = 0 +end + +function init_392() + KiteModels.clear!(kps5) + kps5.set.l_tether = 392.0 + kps5.set.elevation = 70.0 + kps5.set.area = 10.0 + kps5.set.rel_side_area = 50.0 + kps5.set.v_wind = 9.1 + kps5.set.mass = 6.2 + kps5.set.c_s = 0.6 +end + +function init_150() + KiteModels.clear!(kps5) + kps5.set.l_tether = 150.0 + kps5.set.elevation = 70.0 + kps5.set.area = 10.18 + kps5.set.rel_side_area = 30.6 + kps5.set.v_wind = 9.1 + kps5.set.mass = 6.21 + kps5.set.c_s = 0.6 + kps5.set.damping = 473.0 # unit damping + kps5.set.c_spring = 614600.0 # unit spring coefficent + kps5.set.width = 4.9622 +end + +function init3() + kps5.set.alpha = 0.08163 + KiteModels.clear!(kps5) + kps5.set.l_tether = 150.0 # - kps5.set.height_k - kps5.set.h_bridle + kps5.set.area = 10.18 + kps5.set.rel_side_area = 30.6 + kps5.set.mass = 6.21 + kps5.set.c_s = 0.6 + kps5.set.damping = 473.0 # unit damping coefficient + kps5.set.c_spring = 614600.0 # unit spring coefficent + kps5.set.width = 4.9622 + kps5.set.elevation = 70.7 + kps5.set.profile_law = Int(EXPLOG) + pos, vel = KiteModels.init_pos_vel(kps5) + posd = copy(vel) + veld = zero(vel) + height = 134.14733504839947 + kps5.v_wind .= kps5.v_wind_gnd * calc_wind_factor(kps5.am, height) + kps5.stiffness_factor = 1.0 + KiteModels.init_springs!(kps5) + return pos, vel, posd, veld +end + +set_defaults() + +@testset "calc_rho " begin + @test isapprox(calc_rho(kps5.am, 0.0), 1.225, atol=1e-5) + @test isapprox(calc_rho(kps5.am, 100.0), 1.210756, atol=1e-5) +end + +end +nothing \ No newline at end of file From b8e0625813cda193317b0133af6b9952d3508094 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Fri, 11 Apr 2025 14:54:55 +0200 Subject: [PATCH 04/25] Cleanup --- data/settings.yaml | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/data/settings.yaml b/data/settings.yaml index a4b00b7f7..5aee4f242 100644 --- a/data/settings.yaml +++ b/data/settings.yaml @@ -60,18 +60,6 @@ kps4: smc: 0.0 # steering moment coefficient [-] cmq: 0.0 # pitch rate dependant moment coefficient [-] cord_length: 2.0 # average aerodynamic cord length of the kite [m] - -kps4_3l: - radius: 2.0 # the radius of the circle shape on which the kite lines, viewed - # from the front [m] - bridle_center_distance: 4.0 # the distance from point the center bridle connection point of - # the middle line to the kite [m] - middle_length: 1.5 # the cord length of the kite in the middle [m] - tip_length: 0.62 # the cord length of the kite at the tips [m] - min_steering_line_distance: 1.0 # the distance between the left and right steering bridle [m] - # line connections on the kite that are closest to each other [m] - width_3l: 4.1 # width of the kite [m] - aero_surfaces: 3 # the number of aerodynamic surfaces to use per mass point [-] bridle: d_line: 2.5 # bridle line diameter [mm] From d6dbd011cb64775a17978d8fe5ca87001fbc39f6 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Mon, 14 Apr 2025 15:12:19 +0200 Subject: [PATCH 05/25] Importing packages works --- examples/kps5_example.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 examples/kps5_example.jl diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl new file mode 100644 index 000000000..29badd565 --- /dev/null +++ b/examples/kps5_example.jl @@ -0,0 +1,15 @@ +using KiteModels +using Timers +using Pkg +if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end +using KiteUtils +using Dierckx +using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkit: Symbolics, @register_symbolic +using OrdinaryDiffEqCore +using ModelingToolkit, LinearAlgebra, Timers, Parameters, ControlPlots +using Pkg + +set = load_settings("system_kps5.yaml") From 7934a7ccb02e644859565bf4cd9c69494281e57d Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Mon, 14 Apr 2025 17:27:23 +0200 Subject: [PATCH 06/25] Efforts in structuring code 14/04 --- examples/kps5_example.jl | 124 +++++++++- src/KPS5.jl | 510 ++++++++++++++++++++++++--------------- 2 files changed, 424 insertions(+), 210 deletions(-) diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index 29badd565..97790349e 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -4,12 +4,118 @@ using Pkg if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) using TestEnv; TestEnv.activate() end -using KiteUtils -using Dierckx -using ModelingToolkit: t_nounits as t, D_nounits as D -using ModelingToolkit: Symbolics, @register_symbolic -using OrdinaryDiffEqCore -using ModelingToolkit, LinearAlgebra, Timers, Parameters, ControlPlots -using Pkg - -set = load_settings("system_kps5.yaml") +# using ModelingToolkit: Symbolics, @register_symbolic +# using ModelingToolkit, ControlPlots, Parameters +s = KPS5(kcu=kcu) +# kps4::KPS4 = KPS4(kcu) +# dt = 1/set.sample_freq +# time_range = 0:dt:set.sim_time-dt +# steps = length(time_range) +# logger = Logger(KPS5.points(s), steps) +# init_sim!(s) +# generate_getters!(s) +# # simulate(s, logger) +# save_log(logger, "tmp") +# lg = load_log("tmp") +# play(s, lg) +# dt = 1/set.sample_freq + + +# @with_kw mutable struct Settings_not_there_yet @deftype Float64 +# g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] +# save::Bool = false # save animation frames +# # spring constants +# springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] +# springconstant_bridle::Float64 = 10000.0 # BRIDLE unit spring constant [N] +# springconstant_kite::Float64 = 60000.0 # KITE unit spring constant [N] +# # damping +# damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] +# rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] +# rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] +# end +# @with_kw mutable struct KPS5 +# "Reference to the settings2 struct" +# #set::Settings2 = Settings2() +# sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing +# t_0::Float64 = 0.0 +# iter::Int64 = 0 +# prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing +# integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing +# get_state::Function = () -> nothing +# end + +# # ----------------------------- +# # Simulation Function +# # ----------------------------- +# function simulate(s, logger) +# dt = 1/s.set.sample_freq +# tol = s.set.tol +# tspan = (0.0, dt) +# time_range = 0:dt:s.set.sim_time-dt +# steps = length(time_range) +# iter = 0 +# for i in 1:steps +# next_step!(s; dt=s.set.dt) +# u = s.get_state(s.integrator) +# x = u[1][1, :] +# y = u[1][2, :] +# z = u[1][3, :] +# iter += s.iter +# sys_state = SysState{points(s)}() +# sys_state.X .= x +# sys_state.Y .= y +# sys_state.Z .= z +# println("iter: $iter", " steps: $steps") +# log!(logger, sys_state) +# println(x[end], ", ", sys_state.X[end]) +# end +# println("iter: $iter", " steps: $steps") +# return nothing +# end + +# function play(s, lg) +# dt = 1/s.set.sample_freq +# conn = getconnections(s) +# sl = lg.syslog +# total_segmentsvector = Vector{Int64}[] +# for conn_pair in conn +# push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) +# end +# # Add tether segments +# for i in 0:(s.set.segments-2) +# push!(total_segmentsvector, [6+i, 6+i+1]) +# end +# # Add final connection from last tether point to bridle point +# push!(total_segmentsvector, [6+s.set.segments-1, 1]) +# for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt +# # Get positions at this time step +# x = sl.X[step] +# y = sl.Y[step] +# z = sl.Z[step] +# # Create points array for all points in the system +# pointsvector = Vector{Float64}[] +# for i in 1:points(s) #FIX THIS! +# push!(pointsvector, Float64[x[i], y[i], z[i]]) +# end +# # Calculate appropriate limits for the plot +# x_min, x_max = 0, 40 +# z_min, z_max = 0, 60 +# t = dt * (step-1) +# # Plot the kite system at this time step +# plot2d(pointsvector, total_segmentsvector, t; +# zoom = false, +# xlim = (x_min, x_max), +# ylim = (z_min, z_max) +# ) +# # Add a small delay to control animation speed +# sleep(0.05) +# end +# nothing +# end +# function plot_front_view3(lg) +# display(plotxy(lg.y, lg.z; +# xlabel="pos_y [m]", +# ylabel="height [m]", +# fig="front_view")) +# nothing +# end \ No newline at end of file diff --git a/src/KPS5.jl b/src/KPS5.jl index 74ea647fe..e7ff68037 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -30,22 +30,6 @@ are acting on three of the four kite point masses. Four point kite model, included from KiteModels.jl. Scientific background: http://arxiv.org/abs/1406.6218 =# - -# Array of connections of bridlepoints. -# First point, second point, unstressed length. -const SPRINGS_INPUT_5P = [0. 1. 150. - 1. 2. -1. # s1, p7, p8 - 4. 2. -1. # s2, p10, p8 - 4. 5. -1. # s3, p10, p11 - 3. 4. -1. # s4, p9, p10 - 5. 1. -1. # s5, p11, p7 - 4. 1. -1. # s6, p10, p7 - 3. 5. -1. # s7, p9, p11 - 5. 2. -1. # s8, p11, p8 - 2. 3. -1.] # s9, p8, p9 - -# KCU = p7, A = p8, B = p9, C = p10, D = p11 - """ mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel @@ -72,10 +56,6 @@ $(TYPEDFIELDS) wm::AbstractWinchModel "Iterations, number of calls to the function residual!" iter:: Int64 = 0 - "Function for calculation the lift coefficent, using a spline based on the provided value pairs." - calc_cl::Spline1D - "Function for calculation the drag coefficent, using a spline based on the provided value pairs." - calc_cd::Spline1D "wind vector at the height of the kite" v_wind::T = zeros(S, 3) "wind vector at reference height" @@ -91,9 +71,9 @@ $(TYPEDFIELDS) "drag force of kite and bridle; output of calc_aero_forces!" drag_force::T = zeros(S, 3) "side_force acting on the kite" - side_force::T = zeros(S, 3) + side_force::T = zeros(S, 3) # not yet "max_steering angle in radian" - ks::S = 0.0 + ks::S = 0.0 # not yet "lift force of the kite; output of calc_aero_forces!" lift_force::T = zeros(S, 3) "spring force of the current tether segment, output of calc_particle_forces!" @@ -117,23 +97,23 @@ $(TYPEDFIELDS) "drag coefficient of the kite, depending on the angle of attack" param_cd::S = 1.0 "azimuth angle in radian; initial value is zero" - psi::S = zero(S) + # psi::S = zero(S) "depower angle [deg]" - alpha_depower::S = 0.0 + alpha_depower::S = 0.0 # not yet "pitch angle [rad]" pitch::S = 0.0 "pitch rate [rad/s]" - pitch_rate::S = 0.0 + # pitch_rate::S = 0.0 "aoa at particle B" alpha_2::S = 0.0 - "aoa at particle B, corrected formula" - alpha_2b::S = 0.0 - "aoa at particle C" - alpha_3::S = 0.0 - alpha_3b::S = 0.0 - "aoa at particle D" - alpha_4::S = 0.0 - alpha_4b::S = 0.0 + #"aoa at particle B, corrected formula" + # alpha_2b::S = 0.0 + # "aoa at particle C" + # alpha_3::S = 0.0 + # alpha_3b::S = 0.0 + # "aoa at particle D" + # alpha_4::S = 0.0 + # alpha_4b::S = 0.0 "relative start time of the current time interval" t_0::S = 0.0 "reel out speed of the winch" @@ -144,14 +124,14 @@ $(TYPEDFIELDS) l_tether::S = 0.0 "air density at the height of the kite" rho::S = 0.0 - "actual relative depower setting, must be between 0 .. 1.0" - depower::S = 0.0 + #"actual relative depower setting, must be between 0 .. 1.0" + # depower::S = 0.0 #not yet "actual relative steering setting, must be between -1.0 .. 1.0" - steering::S = 0.0 + #steering::S = 0.0 "steering after the kcu, before applying offset and depower sensitivity, -1.0 .. 1.0" - kcu_steering::S = 0.0 + # kcu_steering::S = 0.0 "multiplier for the stiffniss of tether and bridle" - stiffness_factor::S = 1.0 + #stiffness_factor::S = 1.0 "initial masses of the point masses" initial_masses::MVector{P, S} = ones(P) "current masses, depending on the total tether length" @@ -177,44 +157,6 @@ $(TYPEDFIELDS) "z vector of kite reference frame" z::T = zeros(S, 3) end - -""" - clear!(s::KPS5) - -Initialize the kite power model. -""" -function clear!(s::KPS5) - s.t_0 = 0.0 # relative start time of the current time interval - 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] # wind vector at reference height - 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 + KiteModels.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 - roll, pitch, yaw = orient_euler(s) - s.pitch = pitch - s.pitch_rate = 0.0 - KiteModels.set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu)) -end - function KPS5(kcu::KCU) if kcu.set.winch_model == "AsyncMachine" wm = AsyncMachine(kcu.set) @@ -223,88 +165,306 @@ function KPS5(kcu::KCU) end # wm.last_set_speed = kcu.set.v_reel_out s = KPS5{SimFloat, KVec3, kcu.set.segments+KITE_PARTICLES+1, kcu.set.segments+KITE_SPRINGS, SP}(set=kcu.set, - kcu=kcu, wm=wm, calc_cl = Spline1D(kcu.set.alpha_cl, kcu.set.cl_list), - calc_cd=Spline1D(kcu.set.alpha_cd, kcu.set.cd_list) ) - clear!(s) + kcu=kcu, wm=wm) return s end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------ +# deriving a constants for points, total segments, and connections +# --------- +function points(s) + return 5 + s.set.segments +end +function total_segments(s) + return 9 + s.set.segments +end +function getconnections(s) + conn = [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] # connections between KITE points + conn = vcat(conn, [(6+i, 6+i+1) for i in 0:(s.set.segments-2)]...) # connection between tether points + conn = vcat(conn, [(6+s.set.segments-1, 1)]) # connection final tether point to bridle point + return conn +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +# Interpolating polars using Dierckx +# ------------------------------------ +alpha_cl = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] +cl_list = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] +alpha_cd = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] +cd_list = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] +function cl_interp(alpha) + cl_spline = Spline1D(alpha_cl, cl_list) + return cl_spline(alpha) +end +function cd_interp(alpha) + cd_spline = Spline1D(alpha_cd, cd_list) + return cd_spline(alpha) +end +@register_symbolic cl_interp(alpha) +@register_symbolic cd_interp(alpha) +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Initialize the simulation +# ----------------------------------------- +function init_sim!(s::KPS5) + pos, vel = calc_initial_state(s) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) + s.sys = simple_sys + tspan = (0.0, s.set.sim_time) + s.prob = ODEProblem(simple_sys, nothing, tspan) + s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) +end -""" - calc_particle_forces!(s::KPS5, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i) - -Calculate the drag force of the tether segment, defined by the parameters pos1, pos2, vel1 and vel2 -and distribute it equally on the two particles, that are attached to the segment. -The result is stored in the array s.forces. -""" -@inline function calc_particle_forces!(s::KPS5, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i) - l_0 = spring.length # Unstressed length - k = spring.c_spring * s.stiffness_factor # Spring constant - c = spring.damping # Damping coefficient - 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 # compression stiffness kite springs - k2 = 0.1 * k # compression stiffness tether springs - c1 = s.set.rel_damping * c # damping kite springs - spring_vel = unit_vector ⋅ rel_vel - if (norm1 - l_0) > 0.0 - if i > segments # kite springs - 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 +# ------------------------------ +# Calculate Initial State +# ------------------------------ +function calc_initial_state(s::KPS5) + p1location = [s.set.l_tether*cos(s.set.elevation) 0 s.set.l_tether*sin(s.set.elevation)] + kitepos0rot = get_kite_points(s) + POS0 = kitepos0rot .+ p1location' + POS0 = hcat(POS0, zeros(3, 1)) + if s.set.segments > 1 + extra_nodes = [POS0[:,6] + (POS0[:,1] - POS0[:,6]) * i / s.set.segments for i in 1:(s.set.segments-1)] + POS0 = hcat(POS0, extra_nodes...) + end + VEL0 = zeros(3, points(s)) + return POS0, VEL0 +end +# ----------------------------------------------------- +# initializing kite points +# ----------------------------------------------------- +function get_kite_points(s::KPS5) + # Original kite points in local reference frame + kitepos0 = # KITE points + # P1 Bridle P2 P3 P4 P5 + [0.000 s.set.chord_length/2 -s.set.chord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height+s.set.h_bridle s.set.height+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + + beta = s.set.elevation + Y_r = [sin(beta) 0 cos(beta); + 0 1 0; + -cos(beta) 0 sin(beta)] + # Apply rotation to all points + kitepos0rot = Y_r * kitepos0 + + return kitepos0rot +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Define a function to create reference frame update equations +# --------------------------------------------------------------- +function create_reference_frame_equations(pos, e_x, e_y, e_z) + # Calculate vectors for the reference frame + X = pos[:, 2] - pos[:, 3] # Vector from P3 to P2 + Y = pos[:, 5] - pos[:, 4] # Vector from P4 to P5 + Z = cross(X, Y) # Cross product for Z axis + # Normalize these vectors to get unit vectors + norm_X = sqrt(sum(X .^ 2)) + norm_Y = sqrt(sum(Y .^ 2)) + norm_Z = sqrt(sum(Z .^ 2)) + # Create equations to update the reference frame + ref_frame_eqs = [ + e_x[1] ~ X[1] / norm_X, + e_x[2] ~ X[2] / norm_X, + e_x[3] ~ X[3] / norm_X, + e_y[1] ~ Y[1] / norm_Y, + e_y[2] ~ Y[2] / norm_Y, + e_y[3] ~ Y[3] / norm_Y, + e_z[1] ~ Z[1] / norm_Z, + e_z[2] ~ Z[2] / norm_Z, + e_z[3] ~ Z[3] / norm_Z + ] + return ref_frame_eqs +end +# -------------------------------------------------------------------------- +# computing angle of attack +# -------------------------------------------------------------------------- +function compute_alpha1p(v_a, e_z, e_x) + # Calculate the angle Alpha1p + v_a_z = v_a ⋅ e_z # Use the symbolic dot product + v_a_x = v_a ⋅ e_x # Use the symbolic dot product + + # Use atan for the symbolic computation + alpha1p =rad2deg.(atan(v_a_z, v_a_x)) + return alpha1p +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Calculate Rest Lengths +# --------------------------- +function calc_rest_lengths(s::KPS5) + conn = getconnections(s) + POS0, VEL0 = calc_initial_state(s) + lengths = [norm(POS0[:,conn[i][2]] - POS0[:,conn[i][1]]) for i in 1:9] + l10 = norm(POS0[:,1] - POS0[:,6]) + lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) + return lengths, l10 +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ +#Define the Model +# ----------------------------------------------- +function model(s::KPS5, pos, vel) + POS0, VEL0 = pos, vel + rest_lengths, l_tether = calc_rest_lengths(s) + @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether + @parameters rho=s.set.rho cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.Area + @parameters kcu_cd=s.set.kcu_cd kcu_diameter=s.set.kcu_diameter + @variables pos(t)[1:3, 1:points(s)] = POS0 + @variables vel(t)[1:3, 1:points(s)] = VEL0 + @variables acc(t)[1:3, 1:points(s)] + @variables v_app_point(t)[1:3, 1:points(s)] + @variables segment(t)[1:3, 1:total_segments(s)] + @variables unit_vector(t)[1:3, 1:total_segments(s)] + @variables norm1(t)[1:total_segments(s)] + @variables rel_vel(t)[1:3, 1:total_segments(s)] + @variables spring_vel(t)[1:total_segments(s)] + @variables k_spring(t)[1:total_segments(s)] + @variables spring_force(t)[1:3, 1:total_segments(s)] + @variables v_apparent(t)[1:3, 1:total_segments(s)] + @variables v_app_perp(t)[1:3, 1:total_segments(s)] + @variables norm_v_app(t)[1:total_segments(s)] + @variables half_drag_force(t)[1:3, 1:total_segments(s)] + @variables drag_force(t)[1:3, 1:total_segments(s)] + @variables total_force(t)[1:3, 1:points(s)] + # local kite reference frame + @variables e_x(t)[1:3] + @variables e_y(t)[1:3] + @variables e_z(t)[1:3] + @variables alpha1p(t)[1:1] + + eqs1 = vcat(D.(pos) .~ vel, + D.(vel) .~ acc) + eqs2 = vcat(eqs1...) + eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded + # ----------------------------- + # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # ----------------------------- + # connections adding segment connections, from origin to bridle + conn = getconnections(s) + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) + k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] + k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) + # unit damping constants (C1 tether, C2 bridle, C3 kite) + c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] + c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) + # masses + mass_bridlelines = ((s.set.d_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint + mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether + mass_tetherpoints = mass_tether/(s.set.segments+1) + mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether + m_kitepoints = (m_kite/4) + mass_halfbridleline + PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] + PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # ----------------------------- + # Equations for Each Segment (Spring Forces, Drag, etc.) + # ----------------------------- + for i in 1:total_segments(s) + eqs = [ + segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], + norm1[i] ~ norm(segment[:, i]), + unit_vector[:, i] ~ -segment[:, i] / norm1[i], + rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], + spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), + spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], + v_apparent[:, i] ~ s.set.v_wind_tether .- (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, + v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], + norm_v_app[i] ~ norm(v_app_perp[:, i]) + ] + if i > 9 # tether segments + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*d_tether/1000) * v_app_perp[:, i]) + elseif i in [1, 3, 4, 7] # bridle lines, try to find Cd_bridlelines later + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*(s.set.d_line/1000)) * v_app_perp[:, i]) + else # kite + push!(eqs, half_drag_force[:, i] ~ zeros(3)) end - elseif i > segments # kite spring - s.spring_force .= (k1 * (norm1 - l_0) + (c * spring_vel)) * unit_vector - else - s.spring_force .= (k2 * (norm1 - l_0) + (c * spring_vel)) * unit_vector + eqs2 = vcat(eqs2, reduce(vcat, eqs)) 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 # 6.0 = A_real/A_simulated - else - area = norm1 * d_tether + # ----------------------------- + # Reference Frame and Aerodynamic Coefficients + # ----------------------------- + ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) + eqs2 = vcat(eqs2, ref_frame_eqs) + # only 1 AOA + v_a_kite = s.set.v_wind_tether - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step + eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) + # getting Cl and Cd + Cl = cl_interp(s, alpha1p) + Cd = cd_interp(alpha1p) + + # ----------------------------- + # Force Balance at Each Point + # ----------------------------- + for i in 1:points(s) + eqs = [] + force = sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - + sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) + v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i] + if i == 1 # KCU drag at bridle point + area_kcu = pi * ((kcu_diameter / 2) ^ 2) + Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) + Dy_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[2, i]*v_app_point[2, i]) + Dz_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[3, i]*v_app_point[3, i]) + D = [Dx_kcu, Dy_kcu, Dz_kcu] + push!(eqs, total_force[:, i] ~ force + D) + elseif i in 2:5 # the kite points that get Aero Forces + + v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 + # Lift calculation + L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) + # Cross product and normalization + cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) + normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) + L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y + # Final lift force vector + L = L_perpoint * L_direction + + # Drag calculation + D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared + # Create drag direction components + D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] + # Final drag force vector components + D = D_perpoint * D_direction + + # Total aerodynamic force + Fa = [L[1]+ D[1], L[2]+ D[2], L[3]+ D[3]] + push!(eqs, total_force[1, i] ~ force[1] + Fa[1]) + push!(eqs, total_force[2, i] ~ force[2] + Fa[2]) + push!(eqs, total_force[3, i] ~ force[3] + Fa[3]) + elseif i != 6 + push!(eqs, total_force[:, i] ~ force) end + push!(eqs, acc[:, i] ~ s.set.g_earth + total_force[:, i] / PointMasses[i]) + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i]) end - - v_app_perp = s.v_apparent - s.v_apparent ⋅ unit_vector * unit_vector - half_drag_force = (-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_force = (-0.25 * rho * s.set.cd_kcu * norm(v_app_perp_kcu) * kcu_area) * v_app_perp_kcu - @inbounds s.forces[spring.p2] .+= kcu_drag_force - end - - @inbounds s.forces[spring.p1] .+= half_drag_force + s.spring_force - @inbounds s.forces[spring.p2] .+= half_drag_force - s.spring_force - if i == 1 s.last_force .= s.forces[spring.p1] end - nothing + @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) + simple_sys = structural_simplify(sys) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p +end +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# next step function +# ----------------------------- +function next_step!(s::KPS5; dt=(1/s.set.sample_freq)) + s.t_0 = s.integrator.t + steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) + s.iter += 1 + s.integrator.t, steptime end -""" - calc_aero_forces!(s::KPS5, pos, vel, rho, alpha_depower, rel_steering) - -Calculates the aerodynamic forces acting on the kite particles. - -Parameters: -- pos: vector of the particle positions -- vel: vector of the particle velocities -- rho: air density [kg/m^3] -- rel_depower: value between 0.0 and 1.0 -- alpha_depower: depower angle [degrees] -- rel_steering: value between -1.0 and +1.0 - -Updates the vector s.forces of the first parameter. -""" -function calc_aero_forces!(s::KPS5, pos, vel, rho, alpha_depower, rel_steering) +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Generate Getters +# ----------------------------- +function generate_getters!(s::KPS5) + sys = s.sys + c = collect + get_state = ModelingToolkit.getu(sys, + [c(sys.pos)] + ) + s.get_state = (integ) -> get_state(integ) + return nothing end @@ -407,55 +567,3 @@ function spring_forces(s::KPS5; prn=true) forces end -""" - find_steady_state!(s::KPS5; prn=false, delta = 0.01, stiffness_factor=0.035, upwind_dir=-pi/2)) - -Find an initial equilibrium, based on the initial parameters -`l_tether`, elevation and `v_reel_out`. -""" -function find_steady_state!(s::KPS5; prn=false, delta = 0.01, stiffness_factor=0.035, upwind_dir=-pi/2) - set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir=-pi/2) - s.stiffness_factor = stiffness_factor - res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) - iter = 0 - - # helper function for the steady state finder - 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 e - println("Warning in test_initial_condition!") - end - for i in 1:s.set.segments+KITE_PARTICLES-1 - if i != s.set.segments+KITE_PARTICLES-1 - j = i - else - j = i + 1 - end - # copy the x-component of the residual res2 (acceleration) - F[i] = res[1 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)] - # copy the z-component of the residual res2 - F[i+s.set.segments+KITE_PARTICLES] = res[3 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)] - end - # copy the acceleration of point KCU in x direction - i = s.set.segments+1 - F[end-1] = res[1 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] - # copy the acceleration of point C in y direction - i = s.set.segments+3 - x = res[1 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] - y = res[2 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] - F[end] = res[2 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] - iter += 1 - return nothing - end - if prn println("\nStarted function test_nlsolve...") end - X00 = zeros(SimFloat, 2*(s.set.segments+KITE_PARTICLES-1)+2) - results = nlsolve(test_initial_condition!, X00, autoscale=true, xtol=4e-7, ftol=4e-7, iterations=s.set.max_iter) - if prn println("\nresult: $results") 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 From c08b81326946f7667461a319d46a4e2a73b98390 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Mon, 14 Apr 2025 17:28:34 +0200 Subject: [PATCH 07/25] Additional changes --- Project.toml | 1 + data/settings_kps5.yaml | 140 ++++++++++++++++++++++++++++++++++++++++ data/system_kps5.yaml | 3 + test/runtests.jl | 2 +- test/test_kps5.jl | 2 + 5 files changed, 147 insertions(+), 1 deletion(-) create mode 100644 data/settings_kps5.yaml create mode 100644 data/system_kps5.yaml diff --git a/Project.toml b/Project.toml index 44aa56e4b..dce237a42 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.6.17" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" AtmosphericModels = "c59cac55-771d-4f45-b14d-1c681463a295" +ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/data/settings_kps5.yaml b/data/settings_kps5.yaml new file mode 100644 index 000000000..28874b1ad --- /dev/null +++ b/data/settings_kps5.yaml @@ -0,0 +1,140 @@ +system: + log_file: "data/log_8700W_8ms" # filename without extension [replay only] + # use / as path delimiter, even on Windows + log_level: 2 # 0: no logging + time_lapse: 1.0 # relative replay speed + sim_time: 100.0 # simulation time [sim only] + segments: 6 # number of tether segments + sample_freq: 20 # sample frequency in Hz + zoom: 0.03 # zoom factor for the system view + kite_scale: 3.0 # relative zoom factor for the 4 point kite + fixed_font: "" # name or filepath+filename of alternative fixed pitch font, e.g. Liberation Mono + +initial: + l_tether: 150.0 # initial tether length [m] + elevation: 70.8 # initial elevation angle [deg] + v_reel_out: 0.0 # initial reel out speed [m/s] + depower: 25.0 # initial depower settings [%] + +solver: + abs_tol: 0.0006 # absolute tolerance of the DAE solver [m, m/s] + rel_tol: 0.001 # relative tolerance of the DAE solver [-] + solver: "DFBDF" # DAE solver, IDA or DFBDF or DImplicitEuler + linear_solver: "GMRES" # can be GMRES or LapackDense or Dense (only for IDA) + max_order: 4 # maximal order, usually between 3 and 5 (IDA and DFBDF) + max_iter: 200 # max number of iterations of the steady-state-solver + +steering: + c0: 0.0 # steering offset -0.0032 [-] + c_s: 2.59 # steering coefficient one point model; 2.59 was 0.6; TODO: check if it must be divided by kite_area + c2_cor: 0.93 # correction factor one point model + k_ds: 1.5 # influence of the depower angle on the steering sensitivity + delta_st: 0.02 # steering increment (when pressing RIGHT) + max_steering: 16.834 # max. steering angle of the side planes for four point model [degrees] + cs_4p: 1.0 # correction factor for the steering coefficient of the four point model + +depower: + alpha_d_max: 31.0 # max depower angle [deg] + depower_offset: 23.6 # at rel_depower=0.236 the kite is fully powered [%] + +kite: + model: "data/kite.obj" # 3D model of the kite + physical_model: "KPS4" # name of the kite model to use (KPS3 or KPS4) + version: 1 # version of the model to use + mass: 6.2 # kite mass incl. sensor unit [kg] + area: 10.18 # projected kite area [m²] + rel_side_area: 30.6 # relative side area [%] + height: 2.23 # height of the kite [m] + alpha_cl: [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] + cl_list: [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] + alpha_cd: [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] + cd_list: [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] + +kps4: + width: 5.77 # width of the kite [m] + alpha_zero: 4.0 # should be 4 .. 10 [degrees] + alpha_ztip: 10.0 # [degrees] + m_k: 0.2 # relative nose distance; increasing m_k increases C2 of the turn-rate law + rel_nose_mass: 0.47 # relative nose mass + rel_top_mass: 0.4 # mass of the top particle relative to the sum of top and side particles + smc: 0.0 # steering moment coefficient [-] + cmq: 0.0 # pitch rate dependant moment coefficient [-] + cord_length: 2.0 # average aerodynamic cord length of the kite [m] + +bridle: + d_line: 2.5 # bridle line diameter [mm] + l_bridle: 33.4 # sum of the lengths of the bridle lines [m] + h_bridle: 4.9 # height of bridle [m] + rel_compr_stiffness: 0.25 # relative compression stiffness of the kite springs [-] + rel_damping: 6.0 # relative damping of the kite spring (relative to main tether) [-] + +kcu: + kcu_model: "KCU1" # name of the kite control unit model, KCU1 or KCU2 + kcu_mass: 8.4 # mass of the kite control unit [kg] + kcu_diameter: 0.4 # diameter of the KCU for drag calculation [m] + cd_kcu: 0.47 #nonzero, changed form 0.0 tp 0.47 # drag coefficient of the KCU [-] + power2steer_dist: 1.3 # [m] + depower_drum_diameter: 0.069 # [m] + tape_thickness: 0.0006 # [m] + v_depower: 0.075 # max velocity of depowering in units per second (full range: 1 unit) + v_steering: 0.2 # max velocity of steering in units per second (full range: 2 units) + depower_gain: 3.0 # 3.0 means: more than 33% error -> full speed + steering_gain: 3.0 + +tether: + d_tether: 4 # tether diameter [mm] + cd_tether: 0.958 # drag coefficient of the tether + damping: 473.0 # unit damping coefficient [Ns] + c_spring: 614600.0 # unit spring constant coefficient [N] + rho_tether: 724.0 # density of Dyneema [kg/m³] + e_tether: 55000000000.0 # axial tensile modulus of Dyneema (M.B. Ruppert) [Pa] + # SK75: 109 to 132 GPa according to datasheet + +winch: + winch_model: "AsyncMachine" # or TorqueControlledMachine + max_force: 4000 # maximal (nominal) tether force; short overload allowed [N] + v_ro_max: 8.0 # maximal reel-out speed [m/s] + v_ro_min: -8.0 # minimal reel-out speed (=max reel-in speed) [m/s] + drum_radius: 0.1615 # radius of the drum [m] + max_acc: 4.0 # maximal acceleration of the winch [m/s²] + gear_ratio: 6.2 # gear ratio of the winch [-] + inertia_total: 0.204 # total inertia, as seen from the motor/generator [kgm²] + f_coulomb: 122.0 # coulomb friction [N] + c_vf: 30.6 # coefficient for the viscous friction [Ns/m] + p_speed: 10000.0 # proportional gain of the winch speed controller [-] + i_speed: 2500.0 # integral gain of the winch speed controller [-] + +environment: + v_wind: 9.51 + # [9.51, 0.0, 0.0] #changed to vec # wind speed at reference height [m/s] + upwind_dir: -90.0 # upwind direction [deg] + temp_ref: 15.0 # temperature at reference height [°C] + height_gnd: 0.0 # height of groundstation above see level [m] + h_ref: 6.0 # reference height for the wind speed [m] + + rho_0: 1.225 # air density at zero height and 15 °C [kg/m³] + alpha: 0.08163 # exponent of the wind profile law + z0: 0.0002 # surface roughness [m] + profile_law: 3 # 1=EXP, 2=LOG, 3=EXPLOG, 4=FAST_EXP, 5=FAST_LOG, 6=FAST_EXPLOG + # the following parameters are for calculating the turbulent wind field using the Mann model + use_turbulence: 0.0 # turbulence intensity relative to Cabauw, NL + v_wind_gnds: [3.483, 5.324, 8.163] # wind speeds at ref height for calculating the turbulent wind field [m/s] + avg_height: 200.0 # average height during reel out [m] + rel_turbs: [0.342, 0.465, 0.583] # relative turbulence at the v_wind_gnds + i_ref: 0.14 # is the expected value of the turbulence intensity at 15 m/s. + v_ref: 42.9 # five times the average wind speed in m/s at hub height over the full year [m/s] + # Cabauw: 8.5863 m/s * 5.0 = 42.9 m/s + height_step: 2.0 # use a grid with 2m resolution in z direction [m] + grid_step: 2.0 # grid resolution in x and y direction [m] + + + + + # needed: + # dt = 1/sample_freq # time step of the simulation [s] + # conn: [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] + # springcomstamt_bridle + # springcomstamt_tether + # rel_damping_kite: + # total_segments: 9 + segments # total segments [-] + # points: 5 + segments \ No newline at end of file diff --git a/data/system_kps5.yaml b/data/system_kps5.yaml new file mode 100644 index 000000000..171ec6e75 --- /dev/null +++ b/data/system_kps5.yaml @@ -0,0 +1,3 @@ +system: + sim_settings: "settings_kps5.yaml" # simulator settings + \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index e5be9f77b..25704197e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,7 +19,7 @@ KiteUtils.set_data_path("") @testset verbose = true "Testing KiteModels..." begin # include("test_orientation.jl") # include("test_kps3.jl") - include("test_kps4.jl") + # include("test_kps4.jl") include("test_kps5.jl") # if build_is_production_build # include("bench3.jl") diff --git a/test/test_kps5.jl b/test/test_kps5.jl index c3aef27f9..7525e0d35 100644 --- a/test/test_kps5.jl +++ b/test/test_kps5.jl @@ -83,6 +83,8 @@ set_defaults() @testset "calc_rho " begin @test isapprox(calc_rho(kps5.am, 0.0), 1.225, atol=1e-5) @test isapprox(calc_rho(kps5.am, 100.0), 1.210756, atol=1e-5) + kps5_= KPS5(kcu) + @test kps5_ isa KPS5 end end From cc34be98994ed54b3ef530dc6e04a0fdd51119db Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Mon, 14 Apr 2025 17:47:46 +0200 Subject: [PATCH 08/25] Cleaner debugging --- examples/kps5_example.jl | 6 ++- src/KPS5.jl | 98 +--------------------------------------- 2 files changed, 6 insertions(+), 98 deletions(-) diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index 97790349e..c30861771 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -6,8 +6,10 @@ if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) end # using ModelingToolkit: Symbolics, @register_symbolic # using ModelingToolkit, ControlPlots, Parameters -s = KPS5(kcu=kcu) -# kps4::KPS4 = KPS4(kcu) +set = deepcopy(load_settings("system_kps5.yaml")) +kcu::KCU = KCU(set) +s = KPS5(kcu) + # dt = 1/set.sample_freq # time_range = 0:dt:set.sim_time-dt # steps = length(time_range) diff --git a/src/KPS5.jl b/src/KPS5.jl index e7ff68037..cbcafa331 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -45,7 +45,7 @@ use the input and output functions instead. $(TYPEDFIELDS) """ -@with_kw mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel +@with_kw mutable struct KPS5{S, T} <: AbstractKiteModel "Reference to the settings struct" set::Settings "Reference to the KCU model (Kite Control Unit as implemented in the package KitePodModels" @@ -56,100 +56,6 @@ $(TYPEDFIELDS) wm::AbstractWinchModel "Iterations, number of calls to the function residual!" iter:: Int64 = 0 - "wind vector at the height of the kite" - v_wind::T = zeros(S, 3) - "wind vector at reference height" - v_wind_gnd::T = zeros(S, 3) - "wind vector used for the calculation of the tether drag" - v_wind_tether::T = zeros(S, 3) - "apparent wind vector at the kite" - v_apparent::T = zeros(S, 3) - "bridle_factor = set.l_bridle/bridle_length(set)" - bridle_factor::S = 1.0 - "side lift coefficient, the difference of the left and right lift coefficients" - side_cl::S = 0.0 - "drag force of kite and bridle; output of calc_aero_forces!" - drag_force::T = zeros(S, 3) - "side_force acting on the kite" - side_force::T = zeros(S, 3) # not yet - "max_steering angle in radian" - ks::S = 0.0 # not yet - "lift force of the kite; output of calc_aero_forces!" - lift_force::T = zeros(S, 3) - "spring force of the current tether segment, output of calc_particle_forces!" - spring_force::T = zeros(S, 3) - "last winch force" - last_force::T = zeros(S, 3) - "a copy of the residual one (pos,vel) for debugging and unit tests" - res1::SVector{P, KVec3} = zeros(SVector{P, KVec3}) - "a copy of the residual two (vel,acc) for debugging and unit tests" - res2::SVector{P, KVec3} = zeros(SVector{P, KVec3}) - "a copy of the actual positions as output for the user" - pos::SVector{P, KVec3} = zeros(SVector{P, KVec3}) - "a copy of the actual velocities as output for the user" - vel::SVector{P, KVec3} = zeros(SVector{P, KVec3}) - "velocity vector of the kite" - vel_kite::T = zeros(S, 3) - "unstressed segment length [m]" - segment_length::S = 0.0 - "lift coefficient of the kite, depending on the angle of attack" - param_cl::S = 0.2 - "drag coefficient of the kite, depending on the angle of attack" - param_cd::S = 1.0 - "azimuth angle in radian; initial value is zero" - # psi::S = zero(S) - "depower angle [deg]" - alpha_depower::S = 0.0 # not yet - "pitch angle [rad]" - pitch::S = 0.0 - "pitch rate [rad/s]" - # pitch_rate::S = 0.0 - "aoa at particle B" - alpha_2::S = 0.0 - #"aoa at particle B, corrected formula" - # alpha_2b::S = 0.0 - # "aoa at particle C" - # alpha_3::S = 0.0 - # alpha_3b::S = 0.0 - # "aoa at particle D" - # alpha_4::S = 0.0 - # alpha_4b::S = 0.0 - "relative start time of the current time interval" - t_0::S = 0.0 - "reel out speed of the winch" - v_reel_out::S = 0.0 - "reel out speed at the last time step" - last_v_reel_out::S = 0.0 - "unstretched tether length" - l_tether::S = 0.0 - "air density at the height of the kite" - rho::S = 0.0 - #"actual relative depower setting, must be between 0 .. 1.0" - # depower::S = 0.0 #not yet - "actual relative steering setting, must be between -1.0 .. 1.0" - #steering::S = 0.0 - "steering after the kcu, before applying offset and depower sensitivity, -1.0 .. 1.0" - # kcu_steering::S = 0.0 - "multiplier for the stiffniss of tether and bridle" - #stiffness_factor::S = 1.0 - "initial masses of the point masses" - initial_masses::MVector{P, S} = ones(P) - "current masses, depending on the total tether length" - masses::MVector{P, S} = zeros(P) - "vector of the springs, defined as struct" - springs::MVector{Q, SP} = zeros(SP, Q) - "vector of the forces, acting on the particles" - forces::SVector{P, KVec3} = zeros(SVector{P, KVec3}) - "synchronous speed of the motor/ generator" - sync_speed::Union{S, Nothing} = 0.0 - "set_torque of the motor/generator" - set_torque::Union{S, Nothing} = nothing - "set value of the force at the winch, for logging only" - set_force::Union{S, Nothing} = nothing - "set value of the bearing angle in radians, for logging only" - bearing::Union{S, Nothing} = nothing - "coordinates of the attractor point [azimuth, elevation] in radian, for logging only" - attractor::Union{SVector{2, S}, Nothing} = nothing "x vector of kite reference frame" x::T = zeros(S, 3) "y vector of kite reference frame" @@ -164,7 +70,7 @@ function KPS5(kcu::KCU) wm = TorqueControlledMachine(kcu.set) end # wm.last_set_speed = kcu.set.v_reel_out - s = KPS5{SimFloat, KVec3, kcu.set.segments+KITE_PARTICLES+1, kcu.set.segments+KITE_SPRINGS, SP}(set=kcu.set, + s = KPS5{SimFloat, KVec3}(set=kcu.set, kcu=kcu, wm=wm) return s end From 0b8c6ec37c8e5b258efded9d1dc8e74615a79ed9 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Tue, 15 Apr 2025 21:59:00 +0200 Subject: [PATCH 09/25] Works until integrator part --- data/settings_kps5.yaml | 12 ++++-------- examples/kps5_example.jl | 21 +++++++++++++++------ src/KPS5.jl | 38 +++++++++++++++++++++++++------------- 3 files changed, 44 insertions(+), 27 deletions(-) diff --git a/data/settings_kps5.yaml b/data/settings_kps5.yaml index 28874b1ad..279d9f3dd 100644 --- a/data/settings_kps5.yaml +++ b/data/settings_kps5.yaml @@ -3,7 +3,7 @@ system: # use / as path delimiter, even on Windows log_level: 2 # 0: no logging time_lapse: 1.0 # relative replay speed - sim_time: 100.0 # simulation time [sim only] + sim_time: 10.0 # simulation time [sim only] segments: 6 # number of tether segments sample_freq: 20 # sample frequency in Hz zoom: 0.03 # zoom factor for the system view @@ -44,7 +44,7 @@ kite: mass: 6.2 # kite mass incl. sensor unit [kg] area: 10.18 # projected kite area [m²] rel_side_area: 30.6 # relative side area [%] - height: 2.23 # height of the kite [m] + height_k: 2.23 # height of the kite [m] alpha_cl: [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] cl_list: [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] alpha_cd: [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] @@ -130,11 +130,7 @@ environment: - # needed: - # dt = 1/sample_freq # time step of the simulation [s] - # conn: [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] + # needed: # springcomstamt_bridle # springcomstamt_tether - # rel_damping_kite: - # total_segments: 9 + segments # total segments [-] - # points: 5 + segments \ No newline at end of file + # rel_damping_kite: \ No newline at end of file diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index c30861771..f837443fd 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -1,21 +1,30 @@ using KiteModels using Timers using Pkg +using DAEProblemLibrary if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) using TestEnv; TestEnv.activate() end +import OrdinaryDiffEqCore.init +import OrdinaryDiffEqCore.step! +using Dierckx, Interpolations, Serialization, StaticArrays, LinearAlgebra, Statistics, Parameters, NLsolve, + DocStringExtensions, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, NonlinearSolve, FiniteDiff, DifferentiationInterface # using ModelingToolkit: Symbolics, @register_symbolic # using ModelingToolkit, ControlPlots, Parameters set = deepcopy(load_settings("system_kps5.yaml")) kcu::KCU = KCU(set) s = KPS5(kcu) -# dt = 1/set.sample_freq -# time_range = 0:dt:set.sim_time-dt -# steps = length(time_range) -# logger = Logger(KPS5.points(s), steps) -# init_sim!(s) -# generate_getters!(s) +dt = 1/s.set.sample_freq +time_range = 0:dt:set.sim_time-dt +steps = length(time_range) +println("points: ", KiteModels.points(s)) +logger = Logger(KiteModels.points(s), steps) +solver = DImplicitEuler(autodiff=false)#@AutoFiniteDiff()) +KiteModels.get_kite_points(s) +KiteModels.calc_initial_state(s) +KiteModels.init_sim!(s) #(integrate gen getter in initsim) +# KiteModels.generate_getters!(s) # # simulate(s, logger) # save_log(logger, "tmp") # lg = load_log("tmp") diff --git a/src/KPS5.jl b/src/KPS5.jl index cbcafa331..d84b37a2a 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -30,6 +30,7 @@ are acting on three of the four kite point masses. Four point kite model, included from KiteModels.jl. Scientific background: http://arxiv.org/abs/1406.6218 =# +KiteModels.stiffnea=1000 """ mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel @@ -62,6 +63,12 @@ $(TYPEDFIELDS) y::T = zeros(S, 3) "z vector of kite reference frame" z::T = zeros(S, 3) + sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + #t_0::Float64 = 0.0 + #iter::Int64 = 0 + prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing + get_state::Function = () -> nothing end function KPS5(kcu::KCU) if kcu.set.winch_model == "AsyncMachine" @@ -115,9 +122,12 @@ function init_sim!(s::KPS5) s.sys = simple_sys tspan = (0.0, s.set.sim_time) s.prob = ODEProblem(simple_sys, nothing, tspan) + #differential_vars = ones(Bool, length(y0)) + #prob = DAEProblem{true}(residual!, yd0, y0, tspan, s; differential_vars) + #solver = DFBDF(autodiff=AutoFiniteDiff(), max_order=Val{s.set.max_order}()) s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) + #s.integrator = OrdinaryDiffEqCore.init(prob, solver; abstol=abstol, reltol=s.set.rel_tol, save_everystep=false, initializealg=OrdinaryDiffEqCore.NoInit()) end - # ------------------------------ # Calculate Initial State # ------------------------------ @@ -140,9 +150,9 @@ function get_kite_points(s::KPS5) # Original kite points in local reference frame kitepos0 = # KITE points # P1 Bridle P2 P3 P4 P5 - [0.000 s.set.chord_length/2 -s.set.chord_length/2 0 0; + [0.000 s.set.cord_length/2 -s.set.cord_length/2 0 0; 0.000 0 0 -s.set.width/2 s.set.width/2; - 0.000 s.set.height+s.set.h_bridle s.set.height+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridle s.set.h_bridle] beta = s.set.elevation Y_r = [sin(beta) 0 cos(beta); @@ -208,10 +218,12 @@ end function model(s::KPS5, pos, vel) POS0, VEL0 = pos, vel rest_lengths, l_tether = calc_rest_lengths(s) - @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + # @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + # same as Uwe here, only rel c and k , wrt tether, can/should be improved, to be for kite/bridle separately + @parameters K1=s.set.c_spring K2=s.set.rel_compr_stiffness*s.set.c_spring K3=s.set.rel_compr_stiffness*s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether - @parameters rho=s.set.rho cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.Area - @parameters kcu_cd=s.set.kcu_cd kcu_diameter=s.set.kcu_diameter + @parameters rho=s.set.rho_tether g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area + @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter @variables pos(t)[1:3, 1:points(s)] = POS0 @variables vel(t)[1:3, 1:points(s)] = VEL0 @variables acc(t)[1:3, 1:points(s)] @@ -263,7 +275,7 @@ function model(s::KPS5, pos, vel) # ----------------------------- # Equations for Each Segment (Spring Forces, Drag, etc.) # ----------------------------- - for i in 1:total_segments(s) + for i in 1:total_segments(s) eqs = [ segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], norm1[i] ~ norm(segment[:, i]), @@ -272,7 +284,7 @@ function model(s::KPS5, pos, vel) spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], - v_apparent[:, i] ~ s.set.v_wind_tether .- (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, + v_apparent[:, i] ~ [s.set.v_wind, 0.0, 0.0] - (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, # change wind for winddirection! integrate v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], norm_v_app[i] ~ norm(v_app_perp[:, i]) ] @@ -291,11 +303,11 @@ function model(s::KPS5, pos, vel) ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) eqs2 = vcat(eqs2, ref_frame_eqs) # only 1 AOA - v_a_kite = s.set.v_wind_tether - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + v_a_kite = [s.set.v_wind, 0.0, 0.0] - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) # getting Cl and Cd - Cl = cl_interp(s, alpha1p) + Cl = cl_interp(alpha1p) Cd = cd_interp(alpha1p) # ----------------------------- @@ -307,7 +319,7 @@ function model(s::KPS5, pos, vel) sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i] + v_app_point[:, i] ~ [s.set.v_wind, 0.0, 0.0] - vel[:, i] if i == 1 # KCU drag at bridle point area_kcu = pi * ((kcu_diameter / 2) ^ 2) Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) @@ -342,9 +354,9 @@ function model(s::KPS5, pos, vel) elseif i != 6 push!(eqs, total_force[:, i] ~ force) end - push!(eqs, acc[:, i] ~ s.set.g_earth + total_force[:, i] / PointMasses[i]) + push!(eqs, acc[:, i] ~ [0.0, 0.0, g_earth] + total_force[:, i] / PointMasses[i]) eqs2 = vcat(eqs2, reduce(vcat, eqs)) - eqs2 = vcat(eqs2, v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i]) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ [s.set.v_wind, 0.0, 0.0] - vel[:, i]) end @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) simple_sys = structural_simplify(sys) From 6f68256a1ec8fe7be1c79db656f6eefb2a22fe65 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Wed, 16 Apr 2025 12:22:52 +0200 Subject: [PATCH 10/25] Video runs now --- data/settings_kps5.yaml | 4 +- examples/kps5_example.jl | 116 ++++++++++++++++++++------------------- src/KPS5.jl | 37 +++++++++++-- 3 files changed, 93 insertions(+), 64 deletions(-) diff --git a/data/settings_kps5.yaml b/data/settings_kps5.yaml index 279d9f3dd..33c2f3539 100644 --- a/data/settings_kps5.yaml +++ b/data/settings_kps5.yaml @@ -3,7 +3,7 @@ system: # use / as path delimiter, even on Windows log_level: 2 # 0: no logging time_lapse: 1.0 # relative replay speed - sim_time: 10.0 # simulation time [sim only] + sim_time: 10 # simulation time [sim only] segments: 6 # number of tether segments sample_freq: 20 # sample frequency in Hz zoom: 0.03 # zoom factor for the system view @@ -11,7 +11,7 @@ system: fixed_font: "" # name or filepath+filename of alternative fixed pitch font, e.g. Liberation Mono initial: - l_tether: 150.0 # initial tether length [m] + l_tether: 50 # initial tether length [m] elevation: 70.8 # initial elevation angle [deg] v_reel_out: 0.0 # initial reel out speed [m/s] depower: 25.0 # initial depower settings [%] diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index f837443fd..27530663f 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -2,11 +2,14 @@ using KiteModels using Timers using Pkg using DAEProblemLibrary -if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end +# if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) +# using TestEnv; TestEnv.activate() +# end +using ControlPlots import OrdinaryDiffEqCore.init import OrdinaryDiffEqCore.step! +import OrdinaryDiffEqCore.solve +import OrdinaryDiffEqCore using Dierckx, Interpolations, Serialization, StaticArrays, LinearAlgebra, Statistics, Parameters, NLsolve, DocStringExtensions, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, NonlinearSolve, FiniteDiff, DifferentiationInterface # using ModelingToolkit: Symbolics, @register_symbolic @@ -18,18 +21,55 @@ s = KPS5(kcu) dt = 1/s.set.sample_freq time_range = 0:dt:set.sim_time-dt steps = length(time_range) -println("points: ", KiteModels.points(s)) logger = Logger(KiteModels.points(s), steps) -solver = DImplicitEuler(autodiff=false)#@AutoFiniteDiff()) +#solver = DImplicitEuler(autodiff=false)#@AutoFiniteDiff()) KiteModels.get_kite_points(s) KiteModels.calc_initial_state(s) KiteModels.init_sim!(s) #(integrate gen getter in initsim) -# KiteModels.generate_getters!(s) -# # simulate(s, logger) -# save_log(logger, "tmp") -# lg = load_log("tmp") -# play(s, lg) -# dt = 1/set.sample_freq +KiteModels.generate_getters!(s) +KiteModels.simulate(s, logger) +save_log(logger, "tmp") +lg = load_log("tmp") +function play(s, lg) + dt = 1/s.set.sample_freq + conn = KiteModels.getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + for i in 1:KiteModels.points(s) #FIX THIS! + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + # Calculate appropriate limits for the plot + x_min, x_max = 0, 40 + z_min, z_max = 0, 60 + t = (dt) * (step-1) + # Plot the kite system at this time step + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +play(s, lg) # @with_kw mutable struct Settings_not_there_yet @deftype Float64 @@ -84,49 +124,11 @@ KiteModels.init_sim!(s) #(integrate gen getter in initsim) # return nothing # end -# function play(s, lg) -# dt = 1/s.set.sample_freq -# conn = getconnections(s) -# sl = lg.syslog -# total_segmentsvector = Vector{Int64}[] -# for conn_pair in conn -# push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) -# end -# # Add tether segments -# for i in 0:(s.set.segments-2) -# push!(total_segmentsvector, [6+i, 6+i+1]) -# end -# # Add final connection from last tether point to bridle point -# push!(total_segmentsvector, [6+s.set.segments-1, 1]) -# for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt -# # Get positions at this time step -# x = sl.X[step] -# y = sl.Y[step] -# z = sl.Z[step] -# # Create points array for all points in the system -# pointsvector = Vector{Float64}[] -# for i in 1:points(s) #FIX THIS! -# push!(pointsvector, Float64[x[i], y[i], z[i]]) -# end -# # Calculate appropriate limits for the plot -# x_min, x_max = 0, 40 -# z_min, z_max = 0, 60 -# t = dt * (step-1) -# # Plot the kite system at this time step -# plot2d(pointsvector, total_segmentsvector, t; -# zoom = false, -# xlim = (x_min, x_max), -# ylim = (z_min, z_max) -# ) -# # Add a small delay to control animation speed -# sleep(0.05) -# end -# nothing -# end -# function plot_front_view3(lg) -# display(plotxy(lg.y, lg.z; -# xlabel="pos_y [m]", -# ylabel="height [m]", -# fig="front_view")) -# nothing -# end \ No newline at end of file + +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end \ No newline at end of file diff --git a/src/KPS5.jl b/src/KPS5.jl index d84b37a2a..d463e8003 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -64,7 +64,7 @@ $(TYPEDFIELDS) "z vector of kite reference frame" z::T = zeros(S, 3) sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing - #t_0::Float64 = 0.0 + t_0::Float64 = 0.0 #iter::Int64 = 0 prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing @@ -118,6 +118,7 @@ end # ----------------------------------------- function init_sim!(s::KPS5) pos, vel = calc_initial_state(s) + dt = 1/s.set.sample_freq simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) s.sys = simple_sys tspan = (0.0, s.set.sim_time) @@ -125,14 +126,15 @@ function init_sim!(s::KPS5) #differential_vars = ones(Bool, length(y0)) #prob = DAEProblem{true}(residual!, yd0, y0, tspan, s; differential_vars) #solver = DFBDF(autodiff=AutoFiniteDiff(), max_order=Val{s.set.max_order}()) - s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) + #s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) #s.integrator = OrdinaryDiffEqCore.init(prob, solver; abstol=abstol, reltol=s.set.rel_tol, save_everystep=false, initializealg=OrdinaryDiffEqCore.NoInit()) + s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); dt, abstol=s.set.abs_tol, save_on=false) end # ------------------------------ # Calculate Initial State # ------------------------------ function calc_initial_state(s::KPS5) - p1location = [s.set.l_tether*cos(s.set.elevation) 0 s.set.l_tether*sin(s.set.elevation)] + p1location = [s.set.l_tether*cos(deg2rad(s.set.elevation)) 0 s.set.l_tether*sin(deg2rad(s.set.elevation))] kitepos0rot = get_kite_points(s) POS0 = kitepos0rot .+ p1location' POS0 = hcat(POS0, zeros(3, 1)) @@ -154,7 +156,7 @@ function get_kite_points(s::KPS5) 0.000 0 0 -s.set.width/2 s.set.width/2; 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridle s.set.h_bridle] - beta = s.set.elevation + beta = deg2rad(s.set.elevation) Y_r = [sin(beta) 0 cos(beta); 0 1 0; -cos(beta) 0 sin(beta)] @@ -220,7 +222,7 @@ function model(s::KPS5, pos, vel) rest_lengths, l_tether = calc_rest_lengths(s) # @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether # same as Uwe here, only rel c and k , wrt tether, can/should be improved, to be for kite/bridle separately - @parameters K1=s.set.c_spring K2=s.set.rel_compr_stiffness*s.set.c_spring K3=s.set.rel_compr_stiffness*s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping + @parameters K1=s.set.c_spring K2=s.set.c_spring K3=s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether @parameters rho=s.set.rho_tether g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter @@ -371,6 +373,31 @@ function next_step!(s::KPS5; dt=(1/s.set.sample_freq)) s.iter += 1 s.integrator.t, steptime end +function simulate(s, logger) + dt = 1/s.set.sample_freq + tol = s.set.abs_tol + tspan = (0.0, dt) + time_range = 0:dt:s.set.sim_time-dt + steps = length(time_range) + iter = 0 + for i in 1:steps + next_step!(s; dt=dt) + u = s.get_state(s.integrator) + x = u[1][1, :] + y = u[1][2, :] + z = u[1][3, :] + iter += s.iter + sys_state = SysState{points(s)}() + sys_state.X .= x + sys_state.Y .= y + sys_state.Z .= z + println("iter: $iter", " steps: $steps") + log!(logger, sys_state) + println(x[end], ", ", sys_state.X[end]) + end + println("iter: $iter", " steps: $steps") + return nothing +end # ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Generate Getters From 678fde854bbdc601fd3bb12d77dcc98f21151ebb Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 10:58:26 +0200 Subject: [PATCH 11/25] Side view and front view added as a boolean --- examples/kps5_example.jl | 64 ++++++++++++++++++++++++++++++---------- 1 file changed, 49 insertions(+), 15 deletions(-) diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index 27530663f..a4a11b480 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -30,7 +30,7 @@ KiteModels.generate_getters!(s) KiteModels.simulate(s, logger) save_log(logger, "tmp") lg = load_log("tmp") -function play(s, lg) +function play(s, lg, front_view = false, side_view = true) dt = 1/s.set.sample_freq conn = KiteModels.getconnections(s) sl = lg.syslog @@ -44,31 +44,65 @@ function play(s, lg) end # Add final connection from last tether point to bridle point push!(total_segmentsvector, [6+s.set.segments-1, 1]) + for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt # Get positions at this time step x = sl.X[step] y = sl.Y[step] - z = sl.Z[step] + z = sl.Z[step] + # Create points array for all points in the system pointsvector = Vector{Float64}[] - for i in 1:KiteModels.points(s) #FIX THIS! - push!(pointsvector, Float64[x[i], y[i], z[i]]) - end - # Calculate appropriate limits for the plot - x_min, x_max = 0, 40 - z_min, z_max = 0, 60 - t = (dt) * (step-1) - # Plot the kite system at this time step - plot2d(pointsvector, total_segmentsvector, t; - zoom = false, - xlim = (x_min, x_max), - ylim = (z_min, z_max) - ) + + # Use comparison operator == instead of assignment = + if front_view == true + for i in 1:KiteModels.points(s) + # For front view: Y on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[y[i], z[i], x[i]]) + end + + # Calculate appropriate limits for front view + y_min, y_max = -10, 10 + z_min, z_max = 0, 50 + t = (dt) * (step-1) + + # Plot the kite system at this time step (front view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (y_min, y_max), + ylim = (z_min, z_max) + ) + elseif side_view == true + for i in 1:KiteModels.points(s) + # For side view: X on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + + # Calculate appropriate limits for side view + x_min, x_max = 0, 60 + z_min, z_max = 0, 60 + t = (dt) * (step-1) + + # Plot the kite system at this time step (side view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + end + # Add a small delay to control animation speed sleep(0.05) end nothing end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end play(s, lg) From e7f224ae6d9fd94da11701418b149bdb5dcaa561 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 10:58:58 +0200 Subject: [PATCH 12/25] error resolved concerning rho of air instead of tether in aerodynamic model --- src/KPS5.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index d463e8003..70a6b40c9 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -128,7 +128,7 @@ function init_sim!(s::KPS5) #solver = DFBDF(autodiff=AutoFiniteDiff(), max_order=Val{s.set.max_order}()) #s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) #s.integrator = OrdinaryDiffEqCore.init(prob, solver; abstol=abstol, reltol=s.set.rel_tol, save_everystep=false, initializealg=OrdinaryDiffEqCore.NoInit()) - s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); dt, abstol=s.set.abs_tol, save_on=false) + s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); dt, abstol=s.set.abs_tol, reltol = s.set.rel_tol, save_on=false) end # ------------------------------ # Calculate Initial State @@ -154,7 +154,7 @@ function get_kite_points(s::KPS5) # P1 Bridle P2 P3 P4 P5 [0.000 s.set.cord_length/2 -s.set.cord_length/2 0 0; 0.000 0 0 -s.set.width/2 s.set.width/2; - 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridl s.set.h_bridle] beta = deg2rad(s.set.elevation) Y_r = [sin(beta) 0 cos(beta); @@ -224,7 +224,7 @@ function model(s::KPS5, pos, vel) # same as Uwe here, only rel c and k , wrt tether, can/should be improved, to be for kite/bridle separately @parameters K1=s.set.c_spring K2=s.set.c_spring K3=s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether - @parameters rho=s.set.rho_tether g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area + @parameters rho=s.set.rho_0 g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter @variables pos(t)[1:3, 1:points(s)] = POS0 @variables vel(t)[1:3, 1:points(s)] = VEL0 @@ -332,20 +332,20 @@ function model(s::KPS5, pos, vel) elseif i in 2:5 # the kite points that get Aero Forces v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 - # Lift calculation + # # Lift calculation L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) # Cross product and normalization cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y - # Final lift force vector + # Final lift force vector L = L_perpoint * L_direction # Drag calculation D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared - # Create drag direction components + # # Create drag direction components D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] - # Final drag force vector components + # # Final drag force vector components D = D_perpoint * D_direction # Total aerodynamic force From a577bab3479ccd70066c679f5fb84bcaa79a27a1 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 11:04:55 +0200 Subject: [PATCH 13/25] Typo fixed bridle --- src/KPS5.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index 70a6b40c9..0269117b6 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -152,9 +152,9 @@ function get_kite_points(s::KPS5) # Original kite points in local reference frame kitepos0 = # KITE points # P1 Bridle P2 P3 P4 P5 - [0.000 s.set.cord_length/2 -s.set.cord_length/2 0 0; - 0.000 0 0 -s.set.width/2 s.set.width/2; - 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridl s.set.h_bridle] + [0.000 s.set.cord_length/2 -s.set.cord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridle s.set.h_bridle] beta = deg2rad(s.set.elevation) Y_r = [sin(beta) 0 cos(beta); From 2884899f09e43c7ac9c021272ed5f725b5cc0f4d Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 12:44:25 +0200 Subject: [PATCH 14/25] Relative compression stiffness added from settings, along with wind direction working --- src/KPS5.jl | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index 0269117b6..1bc2a89f2 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -214,6 +214,12 @@ function calc_rest_lengths(s::KPS5) lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) return lengths, l10 end +function get_wind_vector(s::KPS5) + v_wind_magnitude = s.set.v_wind + wind_angle = deg2rad(s.set.upwind_dir) # angle measuring form North (pos x), CW(+) + wind_vector = v_wind_magnitude*[-cos(wind_angle),-sin(wind_angle), 0.0] + return wind_vector +end # ------------------------------------------------------------------------------------------------------------------------------------------------------------------ #Define the Model # ----------------------------------------------- @@ -223,7 +229,7 @@ function model(s::KPS5, pos, vel) # @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether # same as Uwe here, only rel c and k , wrt tether, can/should be improved, to be for kite/bridle separately @parameters K1=s.set.c_spring K2=s.set.c_spring K3=s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping - @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether rel_compr_k=s.set.rel_compr_stiffness @parameters rho=s.set.rho_0 g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter @variables pos(t)[1:3, 1:points(s)] = POS0 @@ -274,6 +280,8 @@ function model(s::KPS5, pos, vel) m_kitepoints = (m_kite/4) + mass_halfbridleline PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # getting wind vector + wind_vector = get_wind_vector(s) # ----------------------------- # Equations for Each Segment (Spring Forces, Drag, etc.) # ----------------------------- @@ -284,9 +292,9 @@ function model(s::KPS5, pos, vel) unit_vector[:, i] ~ -segment[:, i] / norm1[i], rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], - k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (rel_compr_k + (1-rel_compr_k)*(norm1[i] > rest_lengths[i])), spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], - v_apparent[:, i] ~ [s.set.v_wind, 0.0, 0.0] - (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, # change wind for winddirection! integrate + v_apparent[:, i] ~ wind_vector - (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, # change wind for winddirection! integrate v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], norm_v_app[i] ~ norm(v_app_perp[:, i]) ] @@ -305,7 +313,7 @@ function model(s::KPS5, pos, vel) ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) eqs2 = vcat(eqs2, ref_frame_eqs) # only 1 AOA - v_a_kite = [s.set.v_wind, 0.0, 0.0] - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + v_a_kite = wind_vector - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) # getting Cl and Cd @@ -321,7 +329,7 @@ function model(s::KPS5, pos, vel) sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - v_app_point[:, i] ~ [s.set.v_wind, 0.0, 0.0] - vel[:, i] + v_app_point[:, i] ~ wind_vector - vel[:, i] if i == 1 # KCU drag at bridle point area_kcu = pi * ((kcu_diameter / 2) ^ 2) Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) @@ -358,7 +366,7 @@ function model(s::KPS5, pos, vel) end push!(eqs, acc[:, i] ~ [0.0, 0.0, g_earth] + total_force[:, i] / PointMasses[i]) eqs2 = vcat(eqs2, reduce(vcat, eqs)) - eqs2 = vcat(eqs2, v_app_point[:, i] ~ [s.set.v_wind, 0.0, 0.0] - vel[:, i]) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ wind_vector - vel[:, i]) end @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) simple_sys = structural_simplify(sys) From fec93c22a24dcedcf08531cbca8776dc3b33cb29 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 15:34:34 +0200 Subject: [PATCH 15/25] Atmospheric model added --- src/KPS5.jl | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index 1bc2a89f2..ebc4aabb3 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -119,7 +119,7 @@ end function init_sim!(s::KPS5) pos, vel = calc_initial_state(s) dt = 1/s.set.sample_freq - simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p, height, wind_vector = model(s, pos, vel) s.sys = simple_sys tspan = (0.0, s.set.sim_time) s.prob = ODEProblem(simple_sys, nothing, tspan) @@ -129,6 +129,7 @@ function init_sim!(s::KPS5) #s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) #s.integrator = OrdinaryDiffEqCore.init(prob, solver; abstol=abstol, reltol=s.set.rel_tol, save_everystep=false, initializealg=OrdinaryDiffEqCore.NoInit()) s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); dt, abstol=s.set.abs_tol, reltol = s.set.rel_tol, save_on=false) + generate_getters!(s) end # ------------------------------ # Calculate Initial State @@ -214,8 +215,8 @@ function calc_rest_lengths(s::KPS5) lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) return lengths, l10 end -function get_wind_vector(s::KPS5) - v_wind_magnitude = s.set.v_wind +function get_wind_vector(s::KPS5, v_wind_height) + v_wind_magnitude = v_wind_height # add function wind_angle = deg2rad(s.set.upwind_dir) # angle measuring form North (pos x), CW(+) wind_vector = v_wind_magnitude*[-cos(wind_angle),-sin(wind_angle), 0.0] return wind_vector @@ -231,7 +232,7 @@ function model(s::KPS5, pos, vel) @parameters K1=s.set.c_spring K2=s.set.c_spring K3=s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether rel_compr_k=s.set.rel_compr_stiffness @parameters rho=s.set.rho_0 g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area - @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter + @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter v_wind_ground=s.set.v_wind @variables pos(t)[1:3, 1:points(s)] = POS0 @variables vel(t)[1:3, 1:points(s)] = VEL0 @variables acc(t)[1:3, 1:points(s)] @@ -239,6 +240,9 @@ function model(s::KPS5, pos, vel) @variables segment(t)[1:3, 1:total_segments(s)] @variables unit_vector(t)[1:3, 1:total_segments(s)] @variables norm1(t)[1:total_segments(s)] + @variables height(t)[1:total_segments(s)] + @variables v_wind_height(t)[1:total_segments(s)] + @variables wind_vector(t)[1:3, 1:total_segments(s)] @variables rel_vel(t)[1:3, 1:total_segments(s)] @variables spring_vel(t)[1:total_segments(s)] @variables k_spring(t)[1:total_segments(s)] @@ -280,8 +284,8 @@ function model(s::KPS5, pos, vel) m_kitepoints = (m_kite/4) + mass_halfbridleline PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) - # getting wind vector - wind_vector = get_wind_vector(s) + # # getting wind vector + # wind_vector = get_wind_vector(s) # ----------------------------- # Equations for Each Segment (Spring Forces, Drag, etc.) # ----------------------------- @@ -289,12 +293,15 @@ function model(s::KPS5, pos, vel) eqs = [ segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], norm1[i] ~ norm(segment[:, i]), + height[i] ~ max((pos[3, conn[i][2]] + pos[3, conn[i][1]])/2, 6.0), # 6m is the minimum height + v_wind_height[i] ~ calc_wind_factor(s.am, height[i])*v_wind_ground, + wind_vector[:, i] ~ get_wind_vector(s, v_wind_height[i]), unit_vector[:, i] ~ -segment[:, i] / norm1[i], rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (rel_compr_k + (1-rel_compr_k)*(norm1[i] > rest_lengths[i])), spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], - v_apparent[:, i] ~ wind_vector - (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, # change wind for winddirection! integrate + v_apparent[:, i] ~ wind_vector[:, i] - (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, # change wind for winddirection! integrate v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], norm_v_app[i] ~ norm(v_app_perp[:, i]) ] @@ -313,7 +320,7 @@ function model(s::KPS5, pos, vel) ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) eqs2 = vcat(eqs2, ref_frame_eqs) # only 1 AOA - v_a_kite = wind_vector - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + v_a_kite = wind_vector[:, 2] - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) # getting Cl and Cd @@ -329,7 +336,7 @@ function model(s::KPS5, pos, vel) sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - v_app_point[:, i] ~ wind_vector - vel[:, i] + v_app_point[:, i] ~ wind_vector[:, i] - vel[:, i] if i == 1 # KCU drag at bridle point area_kcu = pi * ((kcu_diameter / 2) ^ 2) Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) @@ -366,11 +373,11 @@ function model(s::KPS5, pos, vel) end push!(eqs, acc[:, i] ~ [0.0, 0.0, g_earth] + total_force[:, i] / PointMasses[i]) eqs2 = vcat(eqs2, reduce(vcat, eqs)) - eqs2 = vcat(eqs2, v_app_point[:, i] ~ wind_vector - vel[:, i]) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ wind_vector[:, i] - vel[:, i]) end @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) simple_sys = structural_simplify(sys) - simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p, height, wind_vector end # ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # next step function @@ -379,6 +386,14 @@ function next_step!(s::KPS5; dt=(1/s.set.sample_freq)) s.t_0 = s.integrator.t steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) s.iter += 1 + # Get the current state including wind vector + state = s.get_state(s.integrator) + if length(state) >= 2 # Make sure wind_vector is included in state + height = state[2] # The height from the state + wind_vec = state[3] # The wind vector from the state + #println("Height at each segment : $(height[:, :])") + #println("Time: $(s.integrator.t), Wind Vector at each segment: $(wind_vec[:, :])") + end s.integrator.t, steptime end function simulate(s, logger) @@ -414,7 +429,7 @@ function generate_getters!(s::KPS5) sys = s.sys c = collect get_state = ModelingToolkit.getu(sys, - [c(sys.pos)] + [c(sys.pos), c(sys.height) , c(sys.wind_vector)] ) s.get_state = (integ) -> get_state(integ) return nothing From 90e27705102669eb7563c6869183efc0764362dc Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 17:34:30 +0200 Subject: [PATCH 16/25] Final code 17/04 not fully working yet --- src/KPS5.jl | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index ebc4aabb3..e18f2dcf6 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -229,9 +229,9 @@ function model(s::KPS5, pos, vel) rest_lengths, l_tether = calc_rest_lengths(s) # @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether # same as Uwe here, only rel c and k , wrt tether, can/should be improved, to be for kite/bridle separately - @parameters K1=s.set.c_spring K2=s.set.c_spring K3=s.set.c_spring C1=s.set.damping C2=s.set.damping*s.set.rel_damping C3=s.set.damping*s.set.rel_damping + @parameters E_modulus_tether=s.set.e_tether K3=s.set.c_spring_kite C1=s.set.damping C3=s.set.damping_kite_springs @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether rel_compr_k=s.set.rel_compr_stiffness - @parameters rho=s.set.rho_0 g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.area + @parameters rho=s.set.rho_0 g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether d_bridle_line=s.set.d_line S=s.set.area l_bridle=s.set.l_bridle @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter v_wind_ground=s.set.v_wind @variables pos(t)[1:3, 1:points(s)] = POS0 @variables vel(t)[1:3, 1:points(s)] = VEL0 @@ -269,23 +269,26 @@ function model(s::KPS5, pos, vel) # connections adding segment connections, from origin to bridle conn = getconnections(s) # final connection last tether point to bridle point - # unit spring constants (K1 tether, K2 bridle, K3 kite) + # unit spring constants (K1 tether, K2 bridle, K3 kite) K3 from settings + K1 = E_modulus_tether*pi*(d_tether/2000)^2 + K2 = E_modulus_tether*pi*(d_bridle_line/2000)^2 k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) - # unit damping constants (C1 tether, C2 bridle, C3 kite) + # unit damping constants (C1 tether, C2 bridle, C3 kite) C1 and C3 from settings + totalbridlelengths_model = rest_lengths[1]+rest_lengths[3]+rest_lengths[4]+rest_lengths[7] + lines_per_segment = l_bridle/totalbridlelengths_model + C2 = lines_per_segment*C1*(d_bridle_line/d_tether)^2 c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) # masses - mass_bridlelines = ((s.set.d_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_bridlelines = ((d_bridle_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether mass_tetherpoints = mass_tether/(s.set.segments+1) mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether - m_kitepoints = (m_kite/4) + mass_halfbridleline - PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] + m_kitepoints = [(m_kite*s.set.rel_mass_p2) + mass_halfbridleline, (m_kite*s.set.rel_mass_p3) + mass_halfbridleline, (m_kite*s.set.rel_mass_p4) + mass_halfbridleline, (m_kite*s.set.rel_mass_p4) + mass_halfbridleline] # mass P4=P5 + PointMasses = [mass_bridlepoint, m_kitepoints[1], m_kitepoints[2], m_kitepoints[3], m_kitepoints[4]] PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) - # # getting wind vector - # wind_vector = get_wind_vector(s) # ----------------------------- # Equations for Each Segment (Spring Forces, Drag, etc.) # ----------------------------- @@ -377,7 +380,7 @@ function model(s::KPS5, pos, vel) end @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) simple_sys = structural_simplify(sys) - simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p, height, wind_vector + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p, height, wind_vector, norm1 end # ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # next step function @@ -388,11 +391,17 @@ function next_step!(s::KPS5; dt=(1/s.set.sample_freq)) s.iter += 1 # Get the current state including wind vector state = s.get_state(s.integrator) + if !successful_retcode(s.integrator.sol) + println("Return code for solution: ", s.integrator.sol.retcode) + end if length(state) >= 2 # Make sure wind_vector is included in state height = state[2] # The height from the state - wind_vec = state[3] # The wind vector from the state + wind_vec = state[3] + norm1 = state[4] # The norm of the segments from the state + # The wind vector from the state #println("Height at each segment : $(height[:, :])") #println("Time: $(s.integrator.t), Wind Vector at each segment: $(wind_vec[:, :])") + println("Norm of the segments at each segment: $(norm1[:, :])") end s.integrator.t, steptime end @@ -429,7 +438,7 @@ function generate_getters!(s::KPS5) sys = s.sys c = collect get_state = ModelingToolkit.getu(sys, - [c(sys.pos), c(sys.height) , c(sys.wind_vector)] + [c(sys.pos), c(sys.height) , c(sys.wind_vector), c(sys.norm1)] ) s.get_state = (integ) -> get_state(integ) return nothing From 543dd11c784c2471e14e0edd46ad749a5c81e085 Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Thu, 17 Apr 2025 17:35:10 +0200 Subject: [PATCH 17/25] final commit 17/04 --- Project.toml | 3 +- data/settings_kps5.yaml | 22 +- examples/kps5_example.jl | 4 +- src/KPS5Friso.jl | 452 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 470 insertions(+), 11 deletions(-) create mode 100644 src/KPS5Friso.jl diff --git a/Project.toml b/Project.toml index dce237a42..4deff4ad7 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.6.17" ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" AtmosphericModels = "c59cac55-771d-4f45-b14d-1c681463a295" ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" +DAEProblemLibrary = "dfb8ca35-80a1-48ba-a605-84916a45b4f8" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" @@ -61,7 +62,7 @@ DocStringExtensions = "0.9.4" Documenter = "1.10.1" Interpolations = "0.15.1" KitePodModels = "0.3.8" -KiteUtils = "0.10" +KiteUtils = "0.10.2" LaTeXStrings = "1.4.0" LinearSolve = "~2.39.0" ModelingToolkit = "~9.71.0" diff --git a/data/settings_kps5.yaml b/data/settings_kps5.yaml index 33c2f3539..5d2e6c300 100644 --- a/data/settings_kps5.yaml +++ b/data/settings_kps5.yaml @@ -3,7 +3,7 @@ system: # use / as path delimiter, even on Windows log_level: 2 # 0: no logging time_lapse: 1.0 # relative replay speed - sim_time: 10 # simulation time [sim only] + sim_time: 5 # simulation time [sim only] segments: 6 # number of tether segments sample_freq: 20 # sample frequency in Hz zoom: 0.03 # zoom factor for the system view @@ -11,9 +11,9 @@ system: fixed_font: "" # name or filepath+filename of alternative fixed pitch font, e.g. Liberation Mono initial: - l_tether: 50 # initial tether length [m] - elevation: 70.8 # initial elevation angle [deg] - v_reel_out: 0.0 # initial reel out speed [m/s] + l_tether: 10 # initial tether length [m] + elevation: 45 # initial elevation angle [deg] + v_reel_out: 0.5 # initial reel out speed [m/s] depower: 25.0 # initial depower settings [%] solver: @@ -60,7 +60,14 @@ kps4: smc: 0.0 # steering moment coefficient [-] cmq: 0.0 # pitch rate dependant moment coefficient [-] cord_length: 2.0 # average aerodynamic cord length of the kite [m] - + +kps5: + c_spring_kite: 614600.0 # unit spring constant coefficient of the kite springs [N] + damping_kite_springs: 473.0 # unit damping coefficient [Ns] + rel_mass_p2: 0.25 # rel mass of p2 + rel_mass_p3: 0.25 # rel mass of p3 + rel_mass_p4: 0.25 # rel mass of p4 and p5 + bridle: d_line: 2.5 # bridle line diameter [mm] l_bridle: 33.4 # sum of the lengths of the bridle lines [m] @@ -105,9 +112,8 @@ winch: i_speed: 2500.0 # integral gain of the winch speed controller [-] environment: - v_wind: 9.51 - # [9.51, 0.0, 0.0] #changed to vec # wind speed at reference height [m/s] - upwind_dir: -90.0 # upwind direction [deg] + v_wind: 9.5 # wind speed at reference height [m/s] + upwind_dir: -180 # upwind direction (0 is form North, clockwise positive) [deg] temp_ref: 15.0 # temperature at reference height [°C] height_gnd: 0.0 # height of groundstation above see level [m] h_ref: 6.0 # reference height for the wind speed [m] diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index a4a11b480..960fd4ecc 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -79,8 +79,8 @@ function play(s, lg, front_view = false, side_view = true) end # Calculate appropriate limits for side view - x_min, x_max = 0, 60 - z_min, z_max = 0, 60 + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 t = (dt) * (step-1) # Plot the kite system at this time step (side view) diff --git a/src/KPS5Friso.jl b/src/KPS5Friso.jl new file mode 100644 index 000000000..001b9a566 --- /dev/null +++ b/src/KPS5Friso.jl @@ -0,0 +1,452 @@ +# average AOA is used, so now only 1 value, not 4, making editable +# clean up code +# New Plot +# Using logger +# segments: now all segments +using Timers +tic() +using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters, ControlPlots +using ModelingToolkit: Symbolics, @register_symbolic +using OrdinaryDiffEqCore +using Dierckx +using ModelingToolkit: t_nounits as t, D_nounits as D +using KiteUtils +include("plots.jl") +toc() +@with_kw mutable struct Settings2 @deftype Float64 + g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] + v_wind_tether::Vector{Float64} = [13.5, 0.0, 0.0] # wind velocity [m/s] + rho::Float64 = 1.225 # air density [kg/m³] + sim_time::Float64 = 5 # simulation duration [s] + dt = 0.05 # time step [s] + tol = 1e-6 # tolerance for the solver + save::Bool = false # save animation frames + # kite + mass::Float64 = 10.58 # mass of kite [kg] + Area::Float64 = 20.36 # surface area [m²] + width::Float64 = 8.16 # width of kite [m] + height::Float64 = 3.15 # height of kite [m] + chord_length::Float64 = 2.0 # chord length [m] + # KCU + bridle + h_bridle = 4.9 # height of bridle [m] + d_line::Float64 = 2.5 # bridle line diameter [mm] + l_bridle::Float64 = 33.4 # sum of the lengths of the bridle lines [m] + kcu_cd::Float64 = 0.47 # KCU drag coefficient + kcu_diameter::Float64 = 0.4 # KCU diameter [m] + kcu_mass::Float64 = 11 # mass of KCU (at bridle point)[kg] + # spring constants + springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] + springconstant_bridle::Float64 = 10000.0 # BRIDLE unit spring constant [N] + springconstant_kite::Float64 = 60000.0 # KITE unit spring constant [N] + # damping + damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] + rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] + rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] + # tether + l_tether::Float64 = 50.0 # tether length [m] + v_reel_out ::Float64 = 0.5 # reel-out speed [m/s] + rho_tether::Float64 = 724.0 # density of tether [kg/m³] + cd_tether::Float64 = 0.958 # drag coefficient of tether + d_tether::Float64 = 4 # tether diameter [mm] + elevation::Float64 = 1.22173048 # Elevation angle, angle XZ plane, between origin and bridle point [rad] + segments::Int64 = 14 # number of tether segments [-] +end +@with_kw mutable struct KPS5 + "Reference to the settings2 struct" + set::Settings2 = Settings2() + sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + t_0::Float64 = 0.0 + iter::Int64 = 0 + prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing + get_state::Function = () -> nothing +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------ +# deriving a constants for points, total segments, and connections +# --------- +function points(s) + return 5 + s.set.segments +end +function total_segments(s) + return 9 + s.set.segments +end +function getconnections(s) + conn = [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] # connections between KITE points + conn = vcat(conn, [(6+i, 6+i+1) for i in 0:(s.set.segments-2)]...) # connection between tether points + conn = vcat(conn, [(6+s.set.segments-1, 1)]) # connection final tether point to bridle point + return conn +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +# Interpolating polars using Dierckx +# ------------------------------------ +alpha_cl = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] +cl_list = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] +alpha_cd = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] +cd_list = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] +function cl_interp(alpha) + cl_spline = Spline1D(alpha_cl, cl_list) + return cl_spline(alpha) +end +function cd_interp(alpha) + cd_spline = Spline1D(alpha_cd, cd_list) + return cd_spline(alpha) +end +@register_symbolic cl_interp(alpha) +@register_symbolic cd_interp(alpha) +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Initialize the simulation +# ----------------------------------------- +function init_sim!(s::KPS5) + pos, vel = calc_initial_state(s) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) + s.sys = simple_sys + tspan = (0.0, s.set.sim_time) + s.prob = ODEProblem(simple_sys, nothing, tspan) + s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) +end +# ------------------------------ +# Calculate Initial State +# ------------------------------ +function calc_initial_state(s) + p1location = [s.set.l_tether*cos(s.set.elevation) 0 s.set.l_tether*sin(s.set.elevation)] + kitepos0rot = get_kite_points(s) + POS0 = kitepos0rot .+ p1location' + POS0 = hcat(POS0, zeros(3, 1)) + if s.set.segments > 1 + extra_nodes = [POS0[:,6] + (POS0[:,1] - POS0[:,6]) * i / s.set.segments for i in 1:(s.set.segments-1)] + POS0 = hcat(POS0, extra_nodes...) + end + VEL0 = zeros(3, points(s)) + return POS0, VEL0 +end +# ----------------------------------------------------- +# initializing kite points +# ----------------------------------------------------- +function get_kite_points(s) + # Original kite points in local reference frame + kitepos0 = # KITE points + # P1 Bridle P2 P3 P4 P5 + [0.000 s.set.chord_length/2 -s.set.chord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height+s.set.h_bridle s.set.height+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + + beta = s.set.elevation + Y_r = [sin(beta) 0 cos(beta); + 0 1 0; + -cos(beta) 0 sin(beta)] + # Apply rotation to all points + kitepos0rot = Y_r * kitepos0 + + return kitepos0rot +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Define a function to create reference frame update equations +# --------------------------------------------------------------- +function create_reference_frame_equations(pos, e_x, e_y, e_z) + # Calculate vectors for the reference frame + X = pos[:, 2] - pos[:, 3] # Vector from P3 to P2 + Y = pos[:, 5] - pos[:, 4] # Vector from P4 to P5 + Z = cross(X, Y) # Cross product for Z axis + # Normalize these vectors to get unit vectors + norm_X = sqrt(sum(X .^ 2)) + norm_Y = sqrt(sum(Y .^ 2)) + norm_Z = sqrt(sum(Z .^ 2)) + # Create equations to update the reference frame + ref_frame_eqs = [ + e_x[1] ~ X[1] / norm_X, + e_x[2] ~ X[2] / norm_X, + e_x[3] ~ X[3] / norm_X, + e_y[1] ~ Y[1] / norm_Y, + e_y[2] ~ Y[2] / norm_Y, + e_y[3] ~ Y[3] / norm_Y, + e_z[1] ~ Z[1] / norm_Z, + e_z[2] ~ Z[2] / norm_Z, + e_z[3] ~ Z[3] / norm_Z + ] + return ref_frame_eqs +end +# -------------------------------------------------------------------------- +# computing angle of attack +# -------------------------------------------------------------------------- +function compute_alpha1p(v_a, e_z, e_x) + # Calculate the angle Alpha1p + v_a_z = v_a ⋅ e_z # Use the symbolic dot product + v_a_x = v_a ⋅ e_x # Use the symbolic dot product + + # Use atan for the symbolic computation + alpha1p =rad2deg.(atan(v_a_z, v_a_x)) + return alpha1p +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Calculate Rest Lengths +# --------------------------- +function calc_rest_lengths(s) + conn = getconnections(s) + POS0, VEL0 = calc_initial_state(s) + lengths = [norm(POS0[:,conn[i][2]] - POS0[:,conn[i][1]]) for i in 1:9] + l10 = norm(POS0[:,1] - POS0[:,6]) + lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) + return lengths, l10 +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ +#Define the Model +# ----------------------------------------------- +function model(s, pos, vel) + POS0, VEL0 = pos, vel + rest_lengths, l_tether = calc_rest_lengths(s) + @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether + @parameters rho=s.set.rho cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.Area + @parameters kcu_cd=s.set.kcu_cd kcu_diameter=s.set.kcu_diameter + @variables pos(t)[1:3, 1:points(s)] = POS0 + @variables vel(t)[1:3, 1:points(s)] = VEL0 + @variables acc(t)[1:3, 1:points(s)] + @variables v_app_point(t)[1:3, 1:points(s)] + @variables segment(t)[1:3, 1:total_segments(s)] + @variables unit_vector(t)[1:3, 1:total_segments(s)] + @variables norm1(t)[1:total_segments(s)] + @variables rel_vel(t)[1:3, 1:total_segments(s)] + @variables spring_vel(t)[1:total_segments(s)] + @variables k_spring(t)[1:total_segments(s)] + @variables spring_force(t)[1:3, 1:total_segments(s)] + @variables v_apparent(t)[1:3, 1:total_segments(s)] + @variables v_app_perp(t)[1:3, 1:total_segments(s)] + @variables norm_v_app(t)[1:total_segments(s)] + @variables half_drag_force(t)[1:3, 1:total_segments(s)] + @variables drag_force(t)[1:3, 1:total_segments(s)] + @variables total_force(t)[1:3, 1:points(s)] + # local kite reference frame + @variables e_x(t)[1:3] + @variables e_y(t)[1:3] + @variables e_z(t)[1:3] + @variables alpha1p(t)[1:1] + + eqs1 = vcat(D.(pos) .~ vel, + D.(vel) .~ acc) + eqs2 = vcat(eqs1...) + eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded + # ----------------------------- + # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # ----------------------------- + # connections adding segment connections, from origin to bridle + conn = getconnections(s) + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) + k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] + k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) + # unit damping constants (C1 tether, C2 bridle, C3 kite) + c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] + c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) + # masses + mass_bridlelines = ((s.set.d_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint + mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether + mass_tetherpoints = mass_tether/(s.set.segments+1) + mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether + m_kitepoints = (m_kite/4) + mass_halfbridleline + PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] + PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # ----------------------------- + # Equations for Each Segment (Spring Forces, Drag, etc.) + # ----------------------------- + for i in 1:total_segments(s) + eqs = [ + segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], + norm1[i] ~ norm(segment[:, i]), + unit_vector[:, i] ~ -segment[:, i] / norm1[i], + rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], + spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), + spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], + v_apparent[:, i] ~ s.set.v_wind_tether .- (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, + v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], + norm_v_app[i] ~ norm(v_app_perp[:, i]) + ] + if i > 9 # tether segments + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*d_tether/1000) * v_app_perp[:, i]) + elseif i in [1, 3, 4, 7] # bridle lines, try to find Cd_bridlelines later + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*(s.set.d_line/1000)) * v_app_perp[:, i]) + else # kite + push!(eqs, half_drag_force[:, i] ~ zeros(3)) + end + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + end + # ----------------------------- + # Reference Frame and Aerodynamic Coefficients + # ----------------------------- + ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) + eqs2 = vcat(eqs2, ref_frame_eqs) + # only 1 AOA + v_a_kite = s.set.v_wind_tether - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step + eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) + # getting Cl and Cd + Cl = cl_interp(alpha1p) + Cd = cd_interp(alpha1p) + + # ----------------------------- + # Force Balance at Each Point + # ----------------------------- + for i in 1:points(s) + eqs = [] + force = sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - + sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) + v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i] + if i == 1 # KCU drag at bridle point + area_kcu = pi * ((kcu_diameter / 2) ^ 2) + Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) + Dy_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[2, i]*v_app_point[2, i]) + Dz_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[3, i]*v_app_point[3, i]) + D = [Dx_kcu, Dy_kcu, Dz_kcu] + push!(eqs, total_force[:, i] ~ force + D) + elseif i in 2:5 # the kite points that get Aero Forces + + v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 + # Lift calculation + L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) + # Cross product and normalization + cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) + normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) + L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y + # Final lift force vector + L = L_perpoint * L_direction + + # Drag calculation + D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared + # Create drag direction components + D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] + # Final drag force vector components + D = D_perpoint * D_direction + + # Total aerodynamic force + Fa = [L[1]+ D[1], L[2]+ D[2], L[3]+ D[3]] + push!(eqs, total_force[1, i] ~ force[1] + Fa[1]) + push!(eqs, total_force[2, i] ~ force[2] + Fa[2]) + push!(eqs, total_force[3, i] ~ force[3] + Fa[3]) + elseif i != 6 + push!(eqs, total_force[:, i] ~ force) + end + push!(eqs, acc[:, i] ~ s.set.g_earth + total_force[:, i] / PointMasses[i]) + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i]) + end + @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) + simple_sys = structural_simplify(sys) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p +end +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# next step function +# ----------------------------- +function next_step!(s::KPS5; dt=s.set.dt) + s.t_0 = s.integrator.t + steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) + s.iter += 1 + s.integrator.t, steptime +end +# ----------------------------- +# Simulation Function +# ----------------------------- +function simulate(s, logger) + dt = s.set.dt + tol = s.set.tol + tspan = (0.0, dt) + time_range = 0:dt:s.set.sim_time-dt + steps = length(time_range) + iter = 0 + for i in 1:steps + next_step!(s; dt=s.set.dt) + u = s.get_state(s.integrator) + x = u[1][1, :] + y = u[1][2, :] + z = u[1][3, :] + iter += s.iter + sys_state = SysState{points(s)}() + sys_state.X .= x + sys_state.Y .= y + sys_state.Z .= z + println("iter: $iter", " steps: $steps") + log!(logger, sys_state) + println(x[end], ", ", sys_state.X[end]) + end + println("iter: $iter", " steps: $steps") + return nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Generate Getters +# ----------------------------- +function generate_getters!(s) + sys = s.sys + c = collect + get_state = ModelingToolkit.getu(sys, + [c(sys.pos)] + ) + s.get_state = (integ) -> get_state(integ) + return nothing +end +function play(s, lg) + conn = getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + for step in 1:length(0:s.set.dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + for i in 1:points(s) #FIX THIS! + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + # Calculate appropriate limits for the plot + x_min, x_max = 0, 40 + z_min, z_max = 0, 60 + t = s.set.dt * (step-1) + # Plot the kite system at this time step + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Running the Main Function +# ----------------------------- +function main() + global lg, s + set = Settings2() + s = KPS5(set=set) + time_range = 0:set.dt:set.sim_time-set.dt + steps = length(time_range) + logger = Logger(points(s), steps) + init_sim!(s) + generate_getters!(s) + simulate(s, logger) + save_log(logger, "tmp") + lg = load_log("tmp") + play(s, lg) +end + +main() +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \ No newline at end of file From 0a4d73bbbe5208db876f9219a2bf59b812f9da3b Mon Sep 17 00:00:00 2001 From: Friso Broekhuizen Date: Fri, 18 Apr 2025 15:50:14 +0200 Subject: [PATCH 18/25] Exactly this code works in tethers.jl but not here --- examples/KPS5onecode.jl | 453 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 453 insertions(+) create mode 100644 examples/KPS5onecode.jl diff --git a/examples/KPS5onecode.jl b/examples/KPS5onecode.jl new file mode 100644 index 000000000..bf778add5 --- /dev/null +++ b/examples/KPS5onecode.jl @@ -0,0 +1,453 @@ +using ModelingToolkit: t_nounits as t, D_nounits as D +using KiteUtils +using Timers +using Pkg +#using DAEProblemLibrary +# if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) +# using TestEnv; TestEnv.activate() +# end +using ControlPlots +import OrdinaryDiffEqCore.init +import OrdinaryDiffEqCore.step! +import OrdinaryDiffEqCore.solve +import OrdinaryDiffEqCore +using Dierckx, ModelingToolkit, LinearAlgebra, Statistics, Parameters, + OrdinaryDiffEqCore, OrdinaryDiffEqBDF#, OrdinaryDiffEqSDIRK, NonlinearSolve, FiniteDiff, DifferentiationInterface +toc() +@with_kw mutable struct Settings2 @deftype Float64 + g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] + v_wind_tether::Vector{Float64} = [13.5, 0.0, 0.0] # wind velocity [m/s] + rho::Float64 = 1.225 # air density [kg/m³] + sim_time::Float64 = 5 # simulation duration [s] + dt = 0.05 # time step [s] + tol = 0.0006 # tolerance for the solver + save::Bool = false # save animation frames + # kite + mass::Float64 = 6.2 # mass of kite [kg] + Area::Float64 = 10.18 # surface area [m²] + width::Float64 = 5.77 # width of kite [m] + height::Float64 = 2.23 # height of kite [m] + chord_length::Float64 = 2.0 # chord length [m] + # KCU + bridle + h_bridle = 4.9 # height of bridle [m] + d_line::Float64 = 2.5 # bridle line diameter [mm] + l_bridle::Float64 = 33.4 # sum of the lengths of the bridle lines [m] + kcu_cd::Float64 = 0.47 # KCU drag coefficient + kcu_diameter::Float64 = 0.4 # KCU diameter [m] + kcu_mass::Float64 = 8.4 # mass of KCU (at bridle point)[kg] + # spring constants + springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] + springconstant_bridle::Float64 = 614600.0 # BRIDLE unit spring constant [N] + springconstant_kite::Float64 = 614600.0 # KITE unit spring constant [N] + # damping + damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] + rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] + rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] + # tether + l_tether::Float64 = 10.0 # tether length [m] + v_reel_out ::Float64 = 0.5 # reel-out speed [m/s] + rho_tether::Float64 = 724.0 # density of tether [kg/m³] + cd_tether::Float64 = 0.958 # drag coefficient of tether + d_tether::Float64 = 4 # tether diameter [mm] + elevation::Float64 = deg2rad(70.8) # Elevation angle, angle XZ plane, between origin and bridle point [rad] + segments::Int64 = 6 # number of tether segments [-] +end +@with_kw mutable struct KPS5 + "Reference to the settings2 struct" + set::Settings2 = Settings2() + sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + t_0::Float64 = 0.0 + iter::Int64 = 0 + prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing + get_state::Function = () -> nothing +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------ +# deriving a constants for points, total segments, and connections +# --------- +function points(s) + return 5 + s.set.segments +end +function total_segments(s) + return 9 + s.set.segments +end +function getconnections(s) + conn = [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] # connections between KITE points + conn = vcat(conn, [(6+i, 6+i+1) for i in 0:(s.set.segments-2)]...) # connection between tether points + conn = vcat(conn, [(6+s.set.segments-1, 1)]) # connection final tether point to bridle point + return conn +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +# Interpolating polars using Dierckx +# ------------------------------------ +alpha_cl = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] +cl_list = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] +alpha_cd = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] +cd_list = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] +function cl_interp(alpha) + cl_spline = Spline1D(alpha_cl, cl_list) + return cl_spline(alpha) +end +function cd_interp(alpha) + cd_spline = Spline1D(alpha_cd, cd_list) + return cd_spline(alpha) +end +@register_symbolic cl_interp(alpha) +@register_symbolic cd_interp(alpha) +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Initialize the simulation +# ----------------------------------------- +function init_sim!(s::KPS5) + pos, vel = calc_initial_state(s) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) + s.sys = simple_sys + tspan = (0.0, s.set.sim_time) + s.prob = ODEProblem(simple_sys, nothing, tspan) + s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) +end +# ------------------------------ +# Calculate Initial State +# ------------------------------ +function calc_initial_state(s) + p1location = [s.set.l_tether*cos(s.set.elevation) 0 s.set.l_tether*sin(s.set.elevation)] + kitepos0rot = get_kite_points(s) + POS0 = kitepos0rot .+ p1location' + POS0 = hcat(POS0, zeros(3, 1)) + if s.set.segments > 1 + extra_nodes = [POS0[:,6] + (POS0[:,1] - POS0[:,6]) * i / s.set.segments for i in 1:(s.set.segments-1)] + POS0 = hcat(POS0, extra_nodes...) + end + VEL0 = zeros(3, points(s)) + return POS0, VEL0 +end +# ----------------------------------------------------- +# initializing kite points +# ----------------------------------------------------- +function get_kite_points(s) + # Original kite points in local reference frame + kitepos0 = # KITE points + # P1 Bridle P2 P3 P4 P5 + [0.000 s.set.chord_length/2 -s.set.chord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height+s.set.h_bridle s.set.height+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + + beta = s.set.elevation + Y_r = [sin(beta) 0 cos(beta); + 0 1 0; + -cos(beta) 0 sin(beta)] + # Apply rotation to all points + kitepos0rot = Y_r * kitepos0 + + return kitepos0rot +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Define a function to create reference frame update equations +# --------------------------------------------------------------- +function create_reference_frame_equations(pos, e_x, e_y, e_z) + # Calculate vectors for the reference frame + X = pos[:, 2] - pos[:, 3] # Vector from P3 to P2 + Y = pos[:, 5] - pos[:, 4] # Vector from P4 to P5 + Z = cross(X, Y) # Cross product for Z axis + # Normalize these vectors to get unit vectors + norm_X = sqrt(sum(X .^ 2)) + norm_Y = sqrt(sum(Y .^ 2)) + norm_Z = sqrt(sum(Z .^ 2)) + # Create equations to update the reference frame + ref_frame_eqs = [ + e_x[1] ~ X[1] / norm_X, + e_x[2] ~ X[2] / norm_X, + e_x[3] ~ X[3] / norm_X, + e_y[1] ~ Y[1] / norm_Y, + e_y[2] ~ Y[2] / norm_Y, + e_y[3] ~ Y[3] / norm_Y, + e_z[1] ~ Z[1] / norm_Z, + e_z[2] ~ Z[2] / norm_Z, + e_z[3] ~ Z[3] / norm_Z + ] + return ref_frame_eqs +end +# -------------------------------------------------------------------------- +# computing angle of attack +# -------------------------------------------------------------------------- +function compute_alpha1p(v_a, e_z, e_x) + # Calculate the angle Alpha1p + v_a_z = v_a ⋅ e_z # Use the symbolic dot product + v_a_x = v_a ⋅ e_x # Use the symbolic dot product + + # Use atan for the symbolic computation + alpha1p =rad2deg.(atan(v_a_z, v_a_x)) + return alpha1p +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Calculate Rest Lengths +# --------------------------- +function calc_rest_lengths(s) + conn = getconnections(s) + POS0, VEL0 = calc_initial_state(s) + lengths = [norm(POS0[:,conn[i][2]] - POS0[:,conn[i][1]]) for i in 1:9] + l10 = norm(POS0[:,1] - POS0[:,6]) + lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) + return lengths, l10 +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ +#Define the Model +# ----------------------------------------------- +function model(s, pos, vel) + POS0, VEL0 = pos, vel + rest_lengths, l_tether = calc_rest_lengths(s) + @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether + @parameters rho=s.set.rho cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.Area + @parameters kcu_cd=s.set.kcu_cd kcu_diameter=s.set.kcu_diameter + @variables pos(t)[1:3, 1:points(s)] = POS0 + @variables vel(t)[1:3, 1:points(s)] = VEL0 + @variables acc(t)[1:3, 1:points(s)] + @variables v_app_point(t)[1:3, 1:points(s)] + @variables segment(t)[1:3, 1:total_segments(s)] + @variables unit_vector(t)[1:3, 1:total_segments(s)] + @variables norm1(t)[1:total_segments(s)] + @variables rel_vel(t)[1:3, 1:total_segments(s)] + @variables spring_vel(t)[1:total_segments(s)] + @variables k_spring(t)[1:total_segments(s)] + @variables spring_force(t)[1:3, 1:total_segments(s)] + @variables v_apparent(t)[1:3, 1:total_segments(s)] + @variables v_app_perp(t)[1:3, 1:total_segments(s)] + @variables norm_v_app(t)[1:total_segments(s)] + @variables half_drag_force(t)[1:3, 1:total_segments(s)] + @variables drag_force(t)[1:3, 1:total_segments(s)] + @variables total_force(t)[1:3, 1:points(s)] + # local kite reference frame + @variables e_x(t)[1:3] + @variables e_y(t)[1:3] + @variables e_z(t)[1:3] + @variables alpha1p(t)[1:1] + + eqs1 = vcat(D.(pos) .~ vel, + D.(vel) .~ acc) + eqs2 = vcat(eqs1...) + eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded + # ----------------------------- + # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # ----------------------------- + # connections adding segment connections, from origin to bridle + conn = getconnections(s) + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) + k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] + k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) + # unit damping constants (C1 tether, C2 bridle, C3 kite) + c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] + c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) + # masses + mass_bridlelines = ((s.set.d_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint + mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether + mass_tetherpoints = mass_tether/(s.set.segments+1) + mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether + m_kitepoints = (m_kite/4) + mass_halfbridleline + PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] + PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # ----------------------------- + # Equations for Each Segment (Spring Forces, Drag, etc.) + # ----------------------------- + for i in 1:total_segments(s) + eqs = [ + segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], + norm1[i] ~ norm(segment[:, i]), + unit_vector[:, i] ~ -segment[:, i] / norm1[i], + rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], + spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), #define relative compresions stiffnes + spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], + v_apparent[:, i] ~ s.set.v_wind_tether .- (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, + v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], + norm_v_app[i] ~ norm(v_app_perp[:, i]) + ] + if i > 9 # tether segments + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*d_tether/1000) * v_app_perp[:, i]) + elseif i in [1, 3, 4, 7] # bridle lines, try to find Cd_bridlelines later + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*(s.set.d_line/1000)) * v_app_perp[:, i]) + else # kite + push!(eqs, half_drag_force[:, i] ~ zeros(3)) + end + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + end + # ----------------------------- + # Reference Frame and Aerodynamic Coefficients + # ----------------------------- + ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) + eqs2 = vcat(eqs2, ref_frame_eqs) + # only 1 AOA + v_a_kite = s.set.v_wind_tether - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step + eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) + # getting Cl and Cd + Cl = cl_interp(alpha1p) + Cd = cd_interp(alpha1p) + + # ----------------------------- + # Force Balance at Each Point + # ----------------------------- + for i in 1:points(s) + eqs = [] + force = sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - + sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) + v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i] + if i == 1 # KCU drag at bridle point + area_kcu = pi * ((kcu_diameter / 2) ^ 2) + Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) + Dy_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[2, i]*v_app_point[2, i]) + Dz_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[3, i]*v_app_point[3, i]) + D = [Dx_kcu, Dy_kcu, Dz_kcu] + push!(eqs, total_force[:, i] ~ force + D) + elseif i in 2:5 # the kite points that get Aero Forces + + v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 + # Lift calculation + L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) + # Cross product and normalization + cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) + normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) + L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y + # Final lift force vector + L = L_perpoint * L_direction + + # Drag calculation + D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared + # Create drag direction components + D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] + # Final drag force vector components + D = D_perpoint * D_direction + + # Total aerodynamic force + Fa = [L[1]+ D[1], L[2]+ D[2], L[3]+ D[3]] + push!(eqs, total_force[1, i] ~ force[1] + Fa[1]) + push!(eqs, total_force[2, i] ~ force[2] + Fa[2]) + push!(eqs, total_force[3, i] ~ force[3] + Fa[3]) + elseif i != 6 + push!(eqs, total_force[:, i] ~ force) + end + push!(eqs, acc[:, i] ~ s.set.g_earth + total_force[:, i] / PointMasses[i]) + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i]) + end + @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) + simple_sys = structural_simplify(sys) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p +end +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# next step function +# ----------------------------- +function next_step!(s::KPS5; dt=s.set.dt) + s.t_0 = s.integrator.t + steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) + s.iter += 1 + s.integrator.t, steptime +end +# ----------------------------- +# Simulation Function +# ----------------------------- +function simulate(s, logger) + dt = s.set.dt + tol = s.set.tol + tspan = (0.0, dt) + time_range = 0:dt:s.set.sim_time-dt + steps = length(time_range) + iter = 0 + for i in 1:steps + next_step!(s; dt=s.set.dt) + u = s.get_state(s.integrator) + x = u[1][1, :] + y = u[1][2, :] + z = u[1][3, :] + iter += s.iter + sys_state = SysState{points(s)}() + sys_state.X .= x + sys_state.Y .= y + sys_state.Z .= z + println("iter: $iter", " steps: $steps") + log!(logger, sys_state) + println(x[end], ", ", sys_state.X[end]) + end + println("iter: $iter", " steps: $steps") + return nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Generate Getters +# ----------------------------- +function generate_getters!(s) + sys = s.sys + c = collect + get_state = ModelingToolkit.getu(sys, + [c(sys.pos)] + ) + s.get_state = (integ) -> get_state(integ) + return nothing +end +function play(s, lg) + conn = getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + for step in 1:length(0:s.set.dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + for i in 1:points(s) #FIX THIS! + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + # Calculate appropriate limits for the plot + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 + t = s.set.dt * (step-1) + # Plot the kite system at this time step + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Running the Main Function +# ----------------------------- +function main() + global lg, s + set = Settings2() + s = KPS5(set=set) + time_range = 0:set.dt:set.sim_time-set.dt + steps = length(time_range) + logger = Logger(points(s), steps) + init_sim!(s) + generate_getters!(s) + simulate(s, logger) + save_log(logger, "tmp") + lg = load_log("tmp") + play(s, lg) +end + +main() +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \ No newline at end of file From d1b3c88a3f64e22adfdec0f163ad86eff5f8272d Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 00:45:44 +0200 Subject: [PATCH 19/25] Add files to .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 55a910501..8e989afc8 100644 --- a/.gitignore +++ b/.gitignore @@ -48,3 +48,5 @@ data/settings_v9f.yaml Manifest-v1.10.toml.bak Manifest-v1.11.toml.bak Manifest-v1.11.toml.backup +data/prob_dynamic_1.10_3_seg.bin.default +data/prob_dynamic_1.11_3_seg.bin.default From 42b7e3bcde0a9ae225ad8d4dfd5a733914a1d460 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 00:45:56 +0200 Subject: [PATCH 20/25] Bugfix --- src/KPS5.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KPS5.jl b/src/KPS5.jl index e18f2dcf6..271e67ccd 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -30,7 +30,7 @@ are acting on three of the four kite point masses. Four point kite model, included from KiteModels.jl. Scientific background: http://arxiv.org/abs/1406.6218 =# -KiteModels.stiffnea=1000 +# KiteModels.stiffness=1000 """ mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel From f370d0465679251ac19478a8db2e1eb441f6ae09 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 01:07:26 +0200 Subject: [PATCH 21/25] Add example_kps5.jl --- examples/kps5_example_uwe.jl | 106 +++++++++++++++++++++++++++++++++++ examples/plotting_kps5.jl | 95 +++++++++++++++++++++++++++++++ 2 files changed, 201 insertions(+) create mode 100644 examples/kps5_example_uwe.jl create mode 100644 examples/plotting_kps5.jl diff --git a/examples/kps5_example_uwe.jl b/examples/kps5_example_uwe.jl new file mode 100644 index 000000000..4e42f8c2c --- /dev/null +++ b/examples/kps5_example_uwe.jl @@ -0,0 +1,106 @@ +using KiteModels +using Timers +using Pkg +if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end +using ControlPlots + +set = deepcopy(load_settings("system_kps5.yaml")) +kcu::KCU = KCU(set) +s = KPS5(kcu) + +dt = 1/s.set.sample_freq +time_range = 0:dt:set.sim_time-dt +steps = length(time_range) +logger = Logger(KiteModels.points(s), steps) +KiteModels.get_kite_points(s) +KiteModels.calc_initial_state(s) +KiteModels.init_sim!(s) +KiteModels.generate_getters!(s) +KiteModels.simulate(s, logger) +save_log(logger, "tmp") +lg = load_log("tmp") + +function play(s, lg, front_view = false, side_view = true) + dt = 1/s.set.sample_freq + conn = KiteModels.getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + + for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + + # Use comparison operator == instead of assignment = + if front_view == true + for i in 1:KiteModels.points(s) + # For front view: Y on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[y[i], z[i], x[i]]) + end + + # Calculate appropriate limits for front view + y_min, y_max = -10, 10 + z_min, z_max = 0, 50 + t = (dt) * (step-1) + + # Plot the kite system at this time step (front view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (y_min, y_max), + ylim = (z_min, z_max) + ) + elseif side_view == true + for i in 1:KiteModels.points(s) + # For side view: X on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + + # Calculate appropriate limits for side view + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 + t = (dt) * (step-1) + + # Plot the kite system at this time step (side view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + end + + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end +play(s, lg) + +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end \ No newline at end of file diff --git a/examples/plotting_kps5.jl b/examples/plotting_kps5.jl new file mode 100644 index 000000000..a2cd9c999 --- /dev/null +++ b/examples/plotting_kps5.jl @@ -0,0 +1,95 @@ +using KiteModels +using Timers +using Pkg +if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end +using ControlPlots + +include("plotting_kps5.jl") + +set = deepcopy(load_settings("system_kps5.yaml")) +kcu::KCU = KCU(set) +s = KPS5(kcu) + +dt = 1/s.set.sample_freq +time_range = 0:dt:set.sim_time-dt +steps = length(time_range) +logger = Logger(KiteModels.points(s), steps) +KiteModels.get_kite_points(s) +KiteModels.calc_initial_state(s) +init_sim!(s) +KiteModels.generate_getters!(s) +KiteModels.simulate(s, logger) +save_log(logger, "tmp") +lg = load_log("tmp") + +function play(s, lg, front_view = false, side_view = true) + dt = 1/s.set.sample_freq + conn = KiteModels.getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + + for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + + # Use comparison operator == instead of assignment = + if front_view == true + for i in 1:KiteModels.points(s) + # For front view: Y on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[y[i], z[i], x[i]]) + end + + # Calculate appropriate limits for front view + y_min, y_max = -10, 10 + z_min, z_max = 0, 50 + t = (dt) * (step-1) + + # Plot the kite system at this time step (front view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (y_min, y_max), + ylim = (z_min, z_max) + ) + elseif side_view == true + for i in 1:KiteModels.points(s) + # For side view: X on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + + # Calculate appropriate limits for side view + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 + t = (dt) * (step-1) + + # Plot the kite system at this time step (side view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + end + + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end + +play(s, lg) + From 39356d4518f9526196bd1a28583e0782d30172bf Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 01:09:04 +0200 Subject: [PATCH 22/25] Remove two libraries from Project.toml --- Project.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index 4deff4ad7..f11f06ae2 100644 --- a/Project.toml +++ b/Project.toml @@ -6,8 +6,6 @@ version = "0.6.17" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" AtmosphericModels = "c59cac55-771d-4f45-b14d-1c681463a295" -ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" -DAEProblemLibrary = "dfb8ca35-80a1-48ba-a605-84916a45b4f8" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" From 0ba3db8c9e723dea032e37c9ca9bbf87c613d89e Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 01:23:22 +0200 Subject: [PATCH 23/25] Cleanup --- examples/kps5_example_uwe.jl | 22 ++------ examples/plotting_kps5.jl | 98 ++---------------------------------- 2 files changed, 10 insertions(+), 110 deletions(-) diff --git a/examples/kps5_example_uwe.jl b/examples/kps5_example_uwe.jl index 4e42f8c2c..c1c38fef8 100644 --- a/examples/kps5_example_uwe.jl +++ b/examples/kps5_example_uwe.jl @@ -6,6 +6,8 @@ if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) end using ControlPlots +include("plotting_kps5.jl") + set = deepcopy(load_settings("system_kps5.yaml")) kcu::KCU = KCU(set) s = KPS5(kcu) @@ -16,11 +18,10 @@ steps = length(time_range) logger = Logger(KiteModels.points(s), steps) KiteModels.get_kite_points(s) KiteModels.calc_initial_state(s) -KiteModels.init_sim!(s) +init_sim!(s) KiteModels.generate_getters!(s) KiteModels.simulate(s, logger) save_log(logger, "tmp") -lg = load_log("tmp") function play(s, lg, front_view = false, side_view = true) dt = 1/s.set.sample_freq @@ -88,19 +89,6 @@ function play(s, lg, front_view = false, side_view = true) end nothing end -function plot_front_view3(lg) - display(plotxy(lg.y, lg.z; - xlabel="pos_y [m]", - ylabel="height [m]", - fig="front_view")) - nothing -end -play(s, lg) -function plot_front_view3(lg) - display(plotxy(lg.y, lg.z; - xlabel="pos_y [m]", - ylabel="height [m]", - fig="front_view")) - nothing -end \ No newline at end of file +lg = load_log("tmp") +play(s, lg) \ No newline at end of file diff --git a/examples/plotting_kps5.jl b/examples/plotting_kps5.jl index a2cd9c999..facbc3642 100644 --- a/examples/plotting_kps5.jl +++ b/examples/plotting_kps5.jl @@ -1,95 +1,7 @@ -using KiteModels -using Timers -using Pkg -if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end -using ControlPlots - -include("plotting_kps5.jl") - -set = deepcopy(load_settings("system_kps5.yaml")) -kcu::KCU = KCU(set) -s = KPS5(kcu) - -dt = 1/s.set.sample_freq -time_range = 0:dt:set.sim_time-dt -steps = length(time_range) -logger = Logger(KiteModels.points(s), steps) -KiteModels.get_kite_points(s) -KiteModels.calc_initial_state(s) -init_sim!(s) -KiteModels.generate_getters!(s) -KiteModels.simulate(s, logger) -save_log(logger, "tmp") -lg = load_log("tmp") - -function play(s, lg, front_view = false, side_view = true) - dt = 1/s.set.sample_freq - conn = KiteModels.getconnections(s) - sl = lg.syslog - total_segmentsvector = Vector{Int64}[] - for conn_pair in conn - push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) - end - # Add tether segments - for i in 0:(s.set.segments-2) - push!(total_segmentsvector, [6+i, 6+i+1]) - end - # Add final connection from last tether point to bridle point - push!(total_segmentsvector, [6+s.set.segments-1, 1]) - - for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt - # Get positions at this time step - x = sl.X[step] - y = sl.Y[step] - z = sl.Z[step] - - # Create points array for all points in the system - pointsvector = Vector{Float64}[] - - # Use comparison operator == instead of assignment = - if front_view == true - for i in 1:KiteModels.points(s) - # For front view: Y on horizontal axis, Z on vertical axis - push!(pointsvector, Float64[y[i], z[i], x[i]]) - end - - # Calculate appropriate limits for front view - y_min, y_max = -10, 10 - z_min, z_max = 0, 50 - t = (dt) * (step-1) - - # Plot the kite system at this time step (front view) - plot2d(pointsvector, total_segmentsvector, t; - zoom = false, - xlim = (y_min, y_max), - ylim = (z_min, z_max) - ) - elseif side_view == true - for i in 1:KiteModels.points(s) - # For side view: X on horizontal axis, Z on vertical axis - push!(pointsvector, Float64[x[i], y[i], z[i]]) - end - - # Calculate appropriate limits for side view - x_min, x_max = 0, 20 - z_min, z_max = 0, 20 - t = (dt) * (step-1) - - # Plot the kite system at this time step (side view) - plot2d(pointsvector, total_segmentsvector, t; - zoom = false, - xlim = (x_min, x_max), - ylim = (z_min, z_max) - ) - end - - # Add a small delay to control animation speed - sleep(0.05) - end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) nothing end - -play(s, lg) - From d5cd1df244adbec7dd4280dcabc97141cea0b3ac Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 01:27:36 +0200 Subject: [PATCH 24/25] Small fix regarding the packages --- examples/kps5_example.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl index 960fd4ecc..feba1d1e0 100644 --- a/examples/kps5_example.jl +++ b/examples/kps5_example.jl @@ -1,10 +1,9 @@ using KiteModels using Timers using Pkg -using DAEProblemLibrary -# if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) -# using TestEnv; TestEnv.activate() -# end +if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end using ControlPlots import OrdinaryDiffEqCore.init import OrdinaryDiffEqCore.step! From 8e4af33c3eed339fb65efcfdefb3110c84df0efe Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Tue, 22 Apr 2025 01:34:24 +0200 Subject: [PATCH 25/25] Cleanup --- examples/kps5_example.jl | 167 --------------------------------------- src/KPS5.jl | 22 +++--- 2 files changed, 9 insertions(+), 180 deletions(-) delete mode 100644 examples/kps5_example.jl diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl deleted file mode 100644 index feba1d1e0..000000000 --- a/examples/kps5_example.jl +++ /dev/null @@ -1,167 +0,0 @@ -using KiteModels -using Timers -using Pkg -if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end -using ControlPlots -import OrdinaryDiffEqCore.init -import OrdinaryDiffEqCore.step! -import OrdinaryDiffEqCore.solve -import OrdinaryDiffEqCore -using Dierckx, Interpolations, Serialization, StaticArrays, LinearAlgebra, Statistics, Parameters, NLsolve, - DocStringExtensions, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, NonlinearSolve, FiniteDiff, DifferentiationInterface -# using ModelingToolkit: Symbolics, @register_symbolic -# using ModelingToolkit, ControlPlots, Parameters -set = deepcopy(load_settings("system_kps5.yaml")) -kcu::KCU = KCU(set) -s = KPS5(kcu) - -dt = 1/s.set.sample_freq -time_range = 0:dt:set.sim_time-dt -steps = length(time_range) -logger = Logger(KiteModels.points(s), steps) -#solver = DImplicitEuler(autodiff=false)#@AutoFiniteDiff()) -KiteModels.get_kite_points(s) -KiteModels.calc_initial_state(s) -KiteModels.init_sim!(s) #(integrate gen getter in initsim) -KiteModels.generate_getters!(s) -KiteModels.simulate(s, logger) -save_log(logger, "tmp") -lg = load_log("tmp") -function play(s, lg, front_view = false, side_view = true) - dt = 1/s.set.sample_freq - conn = KiteModels.getconnections(s) - sl = lg.syslog - total_segmentsvector = Vector{Int64}[] - for conn_pair in conn - push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) - end - # Add tether segments - for i in 0:(s.set.segments-2) - push!(total_segmentsvector, [6+i, 6+i+1]) - end - # Add final connection from last tether point to bridle point - push!(total_segmentsvector, [6+s.set.segments-1, 1]) - - for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt - # Get positions at this time step - x = sl.X[step] - y = sl.Y[step] - z = sl.Z[step] - - # Create points array for all points in the system - pointsvector = Vector{Float64}[] - - # Use comparison operator == instead of assignment = - if front_view == true - for i in 1:KiteModels.points(s) - # For front view: Y on horizontal axis, Z on vertical axis - push!(pointsvector, Float64[y[i], z[i], x[i]]) - end - - # Calculate appropriate limits for front view - y_min, y_max = -10, 10 - z_min, z_max = 0, 50 - t = (dt) * (step-1) - - # Plot the kite system at this time step (front view) - plot2d(pointsvector, total_segmentsvector, t; - zoom = false, - xlim = (y_min, y_max), - ylim = (z_min, z_max) - ) - elseif side_view == true - for i in 1:KiteModels.points(s) - # For side view: X on horizontal axis, Z on vertical axis - push!(pointsvector, Float64[x[i], y[i], z[i]]) - end - - # Calculate appropriate limits for side view - x_min, x_max = 0, 20 - z_min, z_max = 0, 20 - t = (dt) * (step-1) - - # Plot the kite system at this time step (side view) - plot2d(pointsvector, total_segmentsvector, t; - zoom = false, - xlim = (x_min, x_max), - ylim = (z_min, z_max) - ) - end - - # Add a small delay to control animation speed - sleep(0.05) - end - nothing -end -function plot_front_view3(lg) - display(plotxy(lg.y, lg.z; - xlabel="pos_y [m]", - ylabel="height [m]", - fig="front_view")) - nothing -end -play(s, lg) - - -# @with_kw mutable struct Settings_not_there_yet @deftype Float64 -# g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] -# save::Bool = false # save animation frames -# # spring constants -# springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] -# springconstant_bridle::Float64 = 10000.0 # BRIDLE unit spring constant [N] -# springconstant_kite::Float64 = 60000.0 # KITE unit spring constant [N] -# # damping -# damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] -# rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] -# rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] -# end -# @with_kw mutable struct KPS5 -# "Reference to the settings2 struct" -# #set::Settings2 = Settings2() -# sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing -# t_0::Float64 = 0.0 -# iter::Int64 = 0 -# prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing -# integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing -# get_state::Function = () -> nothing -# end - -# # ----------------------------- -# # Simulation Function -# # ----------------------------- -# function simulate(s, logger) -# dt = 1/s.set.sample_freq -# tol = s.set.tol -# tspan = (0.0, dt) -# time_range = 0:dt:s.set.sim_time-dt -# steps = length(time_range) -# iter = 0 -# for i in 1:steps -# next_step!(s; dt=s.set.dt) -# u = s.get_state(s.integrator) -# x = u[1][1, :] -# y = u[1][2, :] -# z = u[1][3, :] -# iter += s.iter -# sys_state = SysState{points(s)}() -# sys_state.X .= x -# sys_state.Y .= y -# sys_state.Z .= z -# println("iter: $iter", " steps: $steps") -# log!(logger, sys_state) -# println(x[end], ", ", sys_state.X[end]) -# end -# println("iter: $iter", " steps: $steps") -# return nothing -# end - - -function plot_front_view3(lg) - display(plotxy(lg.y, lg.z; - xlabel="pos_y [m]", - ylabel="height [m]", - fig="front_view")) - nothing -end \ No newline at end of file diff --git a/src/KPS5.jl b/src/KPS5.jl index 271e67ccd..06e7d7520 100644 --- a/src/KPS5.jl +++ b/src/KPS5.jl @@ -123,11 +123,6 @@ function init_sim!(s::KPS5) s.sys = simple_sys tspan = (0.0, s.set.sim_time) s.prob = ODEProblem(simple_sys, nothing, tspan) - #differential_vars = ones(Bool, length(y0)) - #prob = DAEProblem{true}(residual!, yd0, y0, tspan, s; differential_vars) - #solver = DFBDF(autodiff=AutoFiniteDiff(), max_order=Val{s.set.max_order}()) - #s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) - #s.integrator = OrdinaryDiffEqCore.init(prob, solver; abstol=abstol, reltol=s.set.rel_tol, save_everystep=false, initializealg=OrdinaryDiffEqCore.NoInit()) s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); dt, abstol=s.set.abs_tol, reltol = s.set.rel_tol, save_on=false) generate_getters!(s) end @@ -221,8 +216,9 @@ function get_wind_vector(s::KPS5, v_wind_height) wind_vector = v_wind_magnitude*[-cos(wind_angle),-sin(wind_angle), 0.0] return wind_vector end + # ------------------------------------------------------------------------------------------------------------------------------------------------------------------ -#Define the Model +# Define the Model # ----------------------------------------------- function model(s::KPS5, pos, vel) POS0, VEL0 = pos, vel @@ -264,12 +260,12 @@ function model(s::KPS5, pos, vel) eqs2 = vcat(eqs1...) eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded # ----------------------------- - # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # Defining the connections and their respective rest lengths, unit spring constants, damping and masses # ----------------------------- - # connections adding segment connections, from origin to bridle + # connections adding segment connections, from origin to bridle conn = getconnections(s) - # final connection last tether point to bridle point - # unit spring constants (K1 tether, K2 bridle, K3 kite) K3 from settings + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) K3 from settings K1 = E_modulus_tether*pi*(d_tether/2000)^2 K2 = E_modulus_tether*pi*(d_bridle_line/2000)^2 k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] @@ -401,7 +397,7 @@ function next_step!(s::KPS5; dt=(1/s.set.sample_freq)) # The wind vector from the state #println("Height at each segment : $(height[:, :])") #println("Time: $(s.integrator.t), Wind Vector at each segment: $(wind_vec[:, :])") - println("Norm of the segments at each segment: $(norm1[:, :])") + #println("Norm of the segments at each segment: $(norm1[:, :])") end s.integrator.t, steptime end @@ -423,9 +419,9 @@ function simulate(s, logger) sys_state.X .= x sys_state.Y .= y sys_state.Z .= z - println("iter: $iter", " steps: $steps") + #println("iter: $iter", " steps: $steps") log!(logger, sys_state) - println(x[end], ", ", sys_state.X[end]) + #println(x[end], ", ", sys_state.X[end]) end println("iter: $iter", " steps: $steps") return nothing