diff --git a/Project.toml b/Project.toml index a5b3ba29..ef4c8267 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,6 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" @@ -36,6 +35,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" +SymbolicAWEModels = "9c9a347c-5289-41db-a9b9-25ccc76c3360" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Timers = "21f18d07-b854-4dab-86f0-c15a3821819a" @@ -70,7 +70,7 @@ KiteUtils = "0.10.15" LaTeXStrings = "1.4.0" LinearAlgebra = "1.10, 1.11" LinearSolve = "~2.39.0, 3" -ModelingToolkit = "~9.78.0" +ModelingToolkit = "10.10.0" NLSolversBase = "~7.8.3" NLsolve = "4.5" NonlinearSolve = "4.8.0" @@ -82,7 +82,6 @@ Parameters = "0.12" Pkg = "1.10, 1.11" PrecompileTools = "1.2, 1.3" REPL = "1.10.0, 1.11.0" -Reexport = "1.1, 1.2" Rotations = "1.7" SHA = "0.7.0" SciMLBase = "<2.115" @@ -91,11 +90,12 @@ StaticArrays = "1.9.7" Statistics = "1" StatsBase = "0.34" Sundials = "4.24" +SymbolicAWEModels = "0.1.2" SymbolicIndexingInterface = "0.3" SymbolicUtils = "3.25" Test = "1" Timers = "0.1.5" -VortexStepMethod = "1.2.6" +VortexStepMethod = "2.0.0" WinchModels = "0.3.6" julia = "1.10, 1.11" diff --git a/data/settings_ram.yaml b/data/settings_ram.yaml index 0a8c52dd..f85c5a33 100644 --- a/data/settings_ram.yaml +++ b/data/settings_ram.yaml @@ -3,14 +3,14 @@ system: # use / as path delimiter, even on Windows time_lapse: 1.0 # relative replay speed sim_time: 100.0 # simulation time [sim only] - segments: 6 # number of tether segments + segments: 3 # 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 initial: - l_tethers: [50.0, 50.0, 50.0] # initial tether length [m] + l_tethers: [50.0, 50.4, 50.4] # initial tether length [m] elevation: 70.8 # initial elevation angle [deg] v_reel_outs: [0.0, 0.0, 0.0] # initial reel out speed [m/s] @@ -69,7 +69,7 @@ environment: 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=F + profile_law: 3 # 1=EXP, 2=LOG, 3=EXPLOG # the following parameters are for calculating the turbulent wind field using the Mann model use_turbulence: 0.0 # turbulence intensity relative to Cabau, NL v_wind_gnds: [3.483, 5.324, 8.163] # wind speeds at ref height for calculating the turbulent wind field [m/s] diff --git a/docs/make.jl b/docs/make.jl index 5a138e2d..2fae365f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,7 +28,6 @@ makedocs(; "Home" => "index.md", "Types" => "types.md", "Functions" => "functions.md", - "SymbolicAWEModel" => "ram_air_kite.md", "Parameters" => "parameters.md", "Examples 1p" => "examples.md", "Examples 4p" => "examples_4p.md", diff --git a/docs/src/functions.md b/docs/src/functions.md index 4be3c824..a2e520a4 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -34,7 +34,6 @@ set_v_wind_ground! unstretched_length tether_length pos_kite -calc_aoa calc_height calc_elevation calc_azimuth diff --git a/docs/src/index.md b/docs/src/index.md index a92cdacf..f388e7c2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -83,7 +83,7 @@ times faster. - the documentation was improved ## Provides -The types [`KPS3`](@ref), [`KPS4`](@ref) and [`SymbolicAWEModel`](@ref), representing the one point, the four point kite model and the ram air kite model, together with the high level simulation interface consisting of the functions [`init!`](@ref) and [`next_step!`](@ref). Other kite models can be added inside or outside of this package by implementing the non-generic methods required for an AbstractKiteModel. +The types [`KPS3`](@ref), [`KPS4`](@ref) and `SymbolicAWEModel`, representing the one point, the four point kite model and the ram air kite model, together with the high level simulation interface consisting of the functions [`init!`](@ref) and [`next_step!`](@ref). Other kite models can be added inside or outside of this package by implementing the non-generic methods required for an AbstractKiteModel. Additional functions to provide inputs and outputs of the model on each time step. In particular the constructor [`SysState`](@ref) can be called once per time step to create a SysState struct for logging or for displaying the state in a viewer. For the KPS3 and KPS4 model, once per time step the [`residual!`](@ref) function is called as many times as needed to find the solution at the end diff --git a/docs/src/ram_air_kite.md b/docs/src/ram_air_kite.md deleted file mode 100644 index de813e6f..00000000 --- a/docs/src/ram_air_kite.md +++ /dev/null @@ -1,53 +0,0 @@ -```@meta -CurrentModule = KiteModels -``` -## Introduction -The [`SystemStructure`](@ref) provides a flexible framework for defining the physical structure of airborne wind energy (AWE) systems using discrete mass-spring-damper models. This structure can represent many different AWE system configurations, from simple single-line kites to complex multi-wing systems with intricate bridle networks. - -The [`SystemStructure`](@ref) serves as input to the [`SymbolicAWEModel`](@ref), which is based on ModelingToolkit and automatically generates symbolic differential algebraic equations from the structural definition. - -## Workflow -1. Define system components ([`Point`](@ref), [`Segment`](@ref), [`Group`](@ref), etc.) -2. Assemble into a [`SystemStructure`](@ref) -3. Pass to [`SymbolicAWEModel`](@ref) for automatic MTK model generation -4. Simulate the resulting symbolic model - -## Public enumerations -```@docs -SegmentType -DynamicsType -``` - -## Public constructors -```@docs -SystemStructure(name, set; points=Point[], groups=Group[], segments=Segment[], - pulleys=Pulley[], tethers=Tether[], winches=Winch[], - wings=Wing[], transforms=Transform[]) -SystemStructure -SymbolicAWEModel(::Settings, ::SystemStructure, ::Vector{<:BodyAerodynamics}, ::Vector{<:VortexStepMethod.Solver}) -SymbolicAWEModel(::Settings) -Point(idx, pos_cad, type) -Point -Group(::Any, ::Any, ::RamAirWing, ::Any, ::Any, ::Any) -Group -Segment(idx, point_idxs, type) -Segment -Pulley(idx, segment_idxs, type) -Pulley -Tether -Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0) -Winch -Wing(idx, group_idxs, R_b_c, pos_cad) -Wing -Transform(idx, elevation, azimuth, heading) -Transform -``` - -## Private functions -```@docs -wing_eqs! -reinit! -scalar_eqs! -linear_vsm_eqs! -force_eqs! -``` \ No newline at end of file diff --git a/docs/src/tutorial_system_structure.md b/docs/src/tutorial_system_structure.md index 73993cd0..1bdeae6e 100644 --- a/docs/src/tutorial_system_structure.md +++ b/docs/src/tutorial_system_structure.md @@ -35,7 +35,7 @@ segments = Segment[] points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0)) ``` -The first point we add is a static point. There are four different [`DynamicsType`](@ref)s to choose from: `STATIC`, `QUASI_STATIC`, `DYNAMIC` and `WING`. `STATIC` just means that the point doesn't move. `DYNAMIC` is a point modeled with acceleration, while `QUASI_STATIC` constrains this acceleration to be zero at all times. A `WING` point is connected to a wing body. +The first point we add is a static point. There are four different `DynamicsType`s to choose from: `STATIC`, `QUASI_STATIC`, `DYNAMIC` and `WING`. `STATIC` just means that the point doesn't move. `DYNAMIC` is a point modeled with acceleration, while `QUASI_STATIC` constrains this acceleration to be zero at all times. A `WING` point is connected to a wing body. Now we can add `DYNAMIC` points and connect them to each other with segments. `BRIDLE` segments don't need to have a tether, because they have a constant unstretched length. ```julia @@ -51,21 +51,21 @@ for i in 1:set.segments end ``` -In order to describe the initial orientation of the structure, we define a [`Transform(idx, elevation, azimuth, heading)`](@ref) with an elevation (-80 degrees), azimuth and heading, and a base position `[0.0, 0.0, 50.0]`. +In order to describe the initial orientation of the structure, we define a `Transform` with an elevation (-80 degrees), azimuth and heading, and a base position `[0.0, 0.0, 50.0]`. ```julia transforms = [Transform(1, deg2rad(-80), 0.0, 0.0; base_pos = [0.0, 0.0, 50.0], base_point_idx=points[1].idx, rot_point_idx=points[end].idx)] ``` -From the points, segments and transform we create a [`SystemStructure(name, set)`](@ref), which can be plotted in 2d to quickly investigate if the model is correct. +From the points, segments and transform we create a `SystemStructure`, which can be plotted in 2d to quickly investigate if the model is correct. ```julia sys_struct = SystemStructure("tether", set; points, segments, transforms) plot(sys_struct, 0.0) ``` ![SystemStructure visualization](tether_sys_struct.png) -If the system looks good, we can easily model it, by first creating a [`SymbolicAWEModel`](@ref), initializing it and stepping through time. +If the system looks good, we can easily model it, by first creating a `SymbolicAWEModel`, initializing it and stepping through time. ```julia sam = SymbolicAWEModel(set, sys_struct) diff --git a/docs/src/types.md b/docs/src/types.md index faffcc9e..ab3c84a1 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -16,7 +16,6 @@ AKM ```@docs KPS3 KPS4 -SymbolicAWEModel ``` These structs store the state of the one point model and four point model. Only in unit tests it is allowed to access the members directly, otherwise use the input and output functions. diff --git a/examples/lin_ram_model.jl b/examples/lin_ram_model.jl index 0b844462..d7ba17cc 100644 --- a/examples/lin_ram_model.jl +++ b/examples/lin_ram_model.jl @@ -1,93 +1,63 @@ # Copyright (c) 2025 Bart van de Lint # SPDX-License-Identifier: MPL-2.0 -#= -This example demonstrates linearized model accuracy by comparing: -1. Nonlinear SymbolicAWEModel model simulation -2. Linearized state-space model simulation - -Both models start from the same operating point and are subjected -to identical steering inputs. The resulting state trajectories are -plotted together to visualize how well the linearized model -approximates the nonlinear dynamics. -=# +# This example demonstrates how to: +# - Load a system configuration from a YAML file, +# - Initialize a SymbolicAWEModel with specified winch torques, +# - Stabilize the system at its operating point, +# - Linearize the system at that point, +# - And plot the Bode plots of the resulting linearized system. using Timers tic() @info "Loading packages " PLOT = true -if PLOT - using Pkg - if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() - end - using ControlPlots, LaTeXStrings, ControlSystemsBase +using Pkg +if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() end +using ControlPlots, LaTeXStrings, ControlSystemsBase -using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore -using ModelingToolkit -using ModelingToolkit: t_nounits +using KiteModels, KiteUtils, LinearAlgebra, Statistics +using KiteModels.SymbolicAWEModels: find_steady_state! +using KiteModels.SymbolicAWEModels.ModelingToolkit +using KiteModels.SymbolicAWEModels.ModelingToolkit: t_nounits toc() -# TODO: use sparse autodiff - -# Simulation parameters -dt = 0.001 -total_time = 1.0 # Increased from 0.1s to 1.0s for better dynamics observation -vsm_interval = 3 -steps = Int(round(total_time / dt)) - -# Steering parameters -steering_freq = 1/2 # Hz - full left-right cycle frequency -steering_magnitude = 5.0 # Magnitude of steering input [Nm] - # Initialize model -set = load_settings("system_ram.yaml") -set.segments = 3 +set = Settings("system_ram.yaml") set_values = [-50.0, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] -set.quasi_static = false -set.physical_model = "simple_ram" -@info "Creating SymbolicAWEModel model..." +@info "Creating SymbolicAWEModel..." s = SymbolicAWEModel(set) -s.set.abs_tol = 1e-2 -s.set.rel_tol = 1e-2 toc() # Define outputs for linearization - heading -lin_outputs = @variables heading(t_nounits)[1] +@variables begin + heading(t_nounits)[1] + angle_of_attack(t_nounits)[1] + tether_len(t_nounits)[1:3] + winch_force(t_nounits)[1:3] +end +lin_outputs = [heading[1], angle_of_attack[1], tether_len[1], winch_force[1]] +@info "Linear outputs: $lin_outputs" # Initialize at elevation with linearization outputs -s.sys_struct.winches[2].tether_length += 0.2 -s.sys_struct.winches[3].tether_length += 0.2 -KiteModels.init!(s; - remake=false, - reload=true, - lin_outputs # Specify which outputs to track in linear model -) +init!(s; lin_outputs) sys = s.sys -@show rad2deg(s.integrator[sys.elevation[1]]) - - @info "System initialized at:" toc() # --- Stabilize system at operating point --- @info "Stabilizing system at operating point..." -s.integrator.ps[sys.stabilize] = true -stabilization_steps = Int(10 ÷ dt) -for i in 1:stabilization_steps - next_step!(s; dt, vsm_interval=0.05÷dt) -end -s.integrator.ps[sys.stabilize] = false +find_steady_state!(s) # --- Linearize at operating point --- @info "Linearizing system at operating point..." -@time (; A, B, C, D) = KiteModels.linearize(s) -@time (; A, B, C, D) = KiteModels.linearize(s) -@show norm(A) +(; A, B, C, D) = linearize!(s) +@time (; A, B, C, D) = linearize!(s) @info "System linearized with matrix dimensions:" A=size(A) B=size(B) C=size(C) D=size(D) sys = ss(A,B,C,D) diff --git a/examples/menu.jl b/examples/menu.jl index f9fa691a..fd377d05 100644 --- a/examples/menu.jl +++ b/examples/menu.jl @@ -18,7 +18,6 @@ options = ["bench = include(\"bench.jl\")", "steering_test_4p = include(\"steering_test_4p.jl\")", "ram_air_kite = SIMPLE=false; include(\"ram_air_kite.jl\")", "simple_ram_air_kite = SIMPLE=true; include(\"ram_air_kite.jl\")", - "lin_ram_model = include(\"lin_ram_model.jl\")", "calc_spectrum = include(\"calc_spectrum.jl\")", "plot_spectrum_ = include(\"plot_spectrum.jl\")", "calculate_rotational_inertia = include(\"calculate_rotational_inertia.jl\")", @@ -39,4 +38,4 @@ function example_menu() end end -example_menu() \ No newline at end of file +example_menu() diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 5ad10857..e7bd8f55 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -5,13 +5,14 @@ using Timers tic() @info "Loading packages " -PLOT = true +PLOT = false using Pkg if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) using TestEnv; TestEnv.activate() end using ControlPlots, LaTeXStrings using KiteModels, LinearAlgebra, Statistics +using KiteModels.SymbolicAWEModels: find_steady_state! if ! @isdefined SIMPLE SIMPLE = false @@ -30,25 +31,16 @@ steering_freq = 1/2 # Hz - full left-right cycle frequency steering_magnitude = 10.0 # Magnitude of steering input [Nm] # Initialize model -set = load_settings("system_ram.yaml") -set.segments = 3 +set = Settings("system_ram.yaml") set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] set.quasi_static = false set.physical_model = SIMPLE ? "simple_ram" : "ram" -@info "Creating wing, aero, vsm_solver, sys_struct and symbolic_awe_model:" +@info "Creating SymbolicAWEModel:" sam = SymbolicAWEModel(set) -sam.set.abs_tol = 1e-2 -sam.set.rel_tol = 1e-2 toc() -# init_Q_b_w, R_b_w = KiteModelsam.initial_orient(sam.set) -# init_kite_pos = init!(sam.sys_struct, sam.set, R_b_w, init_Q_b_w) -# plot(sam.sys_struct, 0.0; zoom=false, front=false) - # Initialize at elevation -set.l_tethers[2] += 0.2 -set.l_tethers[3] += 0.2 init!(sam; remake=false, reload=false) sys = sam.sys @@ -63,7 +55,7 @@ sys_state = SysState(sam) t = 0.0 runtime = 0.0 integ_runtime = 0.0 -bias = set.quasi_static ? 0.45 : 0.40 +bias = set.quasi_static ? 0.45 : 0.35 t0 = sam.integrator.t try @@ -73,14 +65,13 @@ try PLOT && plot(sam, t; zoom=false, front=false) # Calculate steering inputs based on cosine wave - steering = steering_magnitude * cos(2π * steering_freq * t + bias)-0.001*t*t + steering = steering_magnitude * cos(2π * steering_freq * t + bias) set_values = -sam.set.drum_radius .* sam.integrator[sys.winch_force] _vsm_interval = 1 if t > 1.0 set_values .+= [0.0, steering, -steering] # Opposite steering for left/right _vsm_interval = vsm_interval end - # Step simulation steptime = @elapsed next_step!(sam; set_values, dt, vsm_interval=vsm_interval) t_new = sam.integrator.t @@ -149,3 +140,5 @@ display(p) @info "Performance:" times_realtime=(total_time/2)/runtime integrator_times_realtime=(total_time/2)/integ_runtime # 55x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) +# 40-65x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) - commit 6620ed5d0a38e96930615aad9a66e4cd666955f2 +# 40x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) - commit 88a78894038d3cbd50fbff83dfbe5c26266b0637 diff --git a/examples/test_init_4p.jl b/examples/test_init_4p.jl index 5712be1d..d7d4ff90 100644 --- a/examples/test_init_4p.jl +++ b/examples/test_init_4p.jl @@ -6,7 +6,7 @@ using Printf using Pkg using KiteModels, KitePodModels, KiteUtils, LinearAlgebra, Rotations -set = deepcopy(load_settings("system.yaml")) +set = Settings("system.yaml") using Pkg if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) @@ -40,6 +40,8 @@ logger = Logger(set.segments + 5, STEPS) kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) + +set.upwind_dir = UPWIND_DIR integrator = KiteModels.init!(kps4; delta=0.03, stiffness_factor=0.01, prn=STATISTIC) for i=1:80 next_step!(kps4, integrator; set_speed=kps4.set.v_reel_out, dt) diff --git a/examples/tether.jl b/examples/tether.jl deleted file mode 100644 index f74eaef6..00000000 --- a/examples/tether.jl +++ /dev/null @@ -1,37 +0,0 @@ -# SPDX-FileCopyrightText: 2025 Bart van de Lint -# -# SPDX-License-Identifier: MPL-2.0 - -using KiteModels, VortexStepMethod, ControlPlots - -set = se("system_ram.yaml") -set.segments = 20 -dynamics_type = DYNAMIC - -points = Point[] -segments = Segment[] - -points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0)) - -segment_idxs = Int[] -for i in 1:set.segments - global points, segments - point_idx = i+1 - pos = [0.0, 0.0, i * set.l_tether / set.segments] - push!(points, Point(point_idx, pos, dynamics_type; wing_idx=0)) - segment_idx = i - push!(segments, Segment(segment_idx, (point_idx-1, point_idx), BRIDLE)) - push!(segment_idxs, segment_idx) -end - -transforms = [Transform(1, deg2rad(-80), 0.0, 0.0, [0.0, 0.0, 50.0], points[1].idx; rot_point_idx=points[end].idx)] -sys_struct = SystemStructure("tether", set; points, segments, transforms) -plot(sys_struct, 0.0) - -sam = SymbolicAWEModel(set, sys_struct) - -init!(sam; remake=false) -for i in 1:80 - plot(sam, i/set.sample_freq) - next_step!(sam) -end diff --git a/ext/KiteModelsControlPlotsExt.jl b/ext/KiteModelsControlPlotsExt.jl index aa91c2e0..0d9facd1 100644 --- a/ext/KiteModelsControlPlotsExt.jl +++ b/ext/KiteModelsControlPlotsExt.jl @@ -7,41 +7,5 @@ using ControlPlots, KiteModels export plot -function ControlPlots.plot(sys::SystemStructure, reltime; l_tether=50.0, wing_pos=nothing, e_z=zeros(3), zoom=false, front=false, xy=nothing) - pos = [sys.points[i].pos_w for i in eachindex(sys.points)] - !isnothing(wing_pos) && (pos = [pos..., wing_pos...]) - seg = [[sys.segments[i].point_idxs[1], sys.segments[i].point_idxs[2]] for i in eachindex(sys.segments)] - if zoom && !front - xlim = (pos[end][1] - 6, pos[end][1]+6) - ylim = (pos[end][3] - 10, pos[end][3]+2) - elseif zoom && front - xlim = (pos[end][2] - 6, pos[end][2]+6) - ylim = (pos[end][3] - 10, pos[end][3]+2) - elseif !zoom && !front - xlim = (-5, l_tether+10) - ylim = (-5, l_tether+10) - elseif !zoom && front - xlim = (-12.5 - 0.5l_tether, 12.5 + 0.5l_tether) - ylim = (-5, l_tether+10) - end - ControlPlots.plot2d(pos, seg, reltime; zoom, front, xlim, ylim, dz_zoom=0.6, xy) -end - -function ControlPlots.plot(s::SymbolicAWEModel, reltime; kwargs...) - wings = s.sys_struct.wings - pos = s.integrator[s.sys.pos] - if length(wings) > 0 - wing_pos = [s.integrator[s.sys.wing_pos[i, :]] for i in eachindex(wings)] - e_z = [s.integrator[s.sys.e_z[i, :]] for i in eachindex(wings)] - else - wing_pos = nothing - e_z = zeros(3) - end - - for (i, point) in enumerate(s.sys_struct.points) - point.pos_w .= pos[:, i] - end - plot(s.sys_struct, reltime; s.set.l_tether, wing_pos, e_z, kwargs...) -end end diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 9123557c..9a195a10 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -12,43 +12,92 @@ Scientific background: http://arxiv.org/abs/1406.6218 =# module KiteModels +#======================================================================# +# DEPENDENCIES +#======================================================================# + +# --- Julia Standard Library & General Utilities --- using PrecompileTools: @setup_workload, @compile_workload -using Dierckx, Interpolations, Serialization, StaticArrays, LinearAlgebra, Statistics, Parameters, NLsolve, - DocStringExtensions, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqNonlinearSolve, - NonlinearSolve, SHA +using DocStringExtensions +using Pkg +using SHA +using Serialization +using Statistics +using LinearAlgebra + +# --- Numerical & Scientific Computing --- +using Dierckx +using Interpolations +using StaticArrays +using Parameters +using Rotations + +# --- Solvers (Nonlinear, Differential Equations) --- +using NLsolve +using NonlinearSolve +using OrdinaryDiffEqCore +using OrdinaryDiffEqBDF +using OrdinaryDiffEqNonlinearSolve import Sundials -using Reexport, Pkg + +# --- Open Source AWE Packages --- using VortexStepMethod using KiteUtils -import KiteUtils: init!, next_step!, update_sys_state! -import KiteUtils: calc_elevation, calc_heading, calc_course, SysState -@reexport using VortexStepMethod: RamAirWing, BodyAerodynamics, Solver, NONLIN -@reexport using KitePodModels -@reexport using WinchModels -@reexport using AtmosphericModels -using Rotations +using SymbolicAWEModels +using KitePodModels +using WinchModels +using AtmosphericModels + +#======================================================================# +# EXPORTS +# (The Public API of this Module) +#======================================================================# + +# --- Importing functions to be extended or exported --- import Base.zero -import OrdinaryDiffEqCore.init -import OrdinaryDiffEqCore.step! -using ModelingToolkit, SymbolicIndexingInterface -using ModelingToolkit: t_nounits as t, D_nounits as D -using ADTypes: AutoFiniteDiff -import ModelingToolkit.SciMLBase: successful_retcode - -export KPS3, KPS4, SymbolicAWEModel, KVec3, SimFloat, 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!, reinit!, next_step!, init_pos_vel # high level workers -export pos_kite, calc_height, calc_elevation, calc_azimuth, calc_heading, calc_course, calc_orient_quat, calc_aoa # getters -export calc_azimuth_north, calc_azimuth_east -export winch_force, lift_drag, cl_cd, lift_over_drag, unstretched_length, tether_length, v_wind_kite # getters -export calculate_rotational_inertia! -export kite_ref_frame, orient_euler, spring_forces, upwind_dir, copy_model_settings, menu2 -export create_ram_sys_struct, create_simple_ram_sys_struct import LinearAlgebra: norm -export SystemStructure, Point, Group, Segment, Pulley, Tether, Winch, Wing, Transform -export DynamicsType, DYNAMIC, QUASI_STATIC, WING, STATIC -export SegmentType, POWER_LINE, STEERING_LINE, BRIDLE +import OrdinaryDiffEqCore: init, step! + +# --- SymbolicAWEModels --- +import SymbolicAWEModels: Wing +export SymbolicAWEModel, linearize! +export SystemStructure, Point, Segment, Pulley, Winch, Wing, Transform +export DynamicsType, DYNAMIC, STATIC, BRIDLE, QUASI_STATIC + +# --- KiteUtils --- +import KiteUtils: init!, next_step!, update_sys_state!, calc_heading, calc_elevation, + calc_course +export init!, next_step!, update_sys_state!, load_settings, se, calc_heading, + calc_orient_quat +export SysState, Settings + +# --- KitePodModels --- +export KCU + +# --- AtmosphericModels --- +export AtmosphericModel + +# --- Types --- +export SimFloat, KVec3, SVec3, AbstractKiteModel, AKM + +# --- Kite models --- +export KPS3, KPS4 + +# --- Input functions --- +export set_depower_steering!, set_v_wind_ground! + +# --- Output functions --- +export unstretched_length, tether_length, pos_kite, calc_height, calc_elevation, + calc_azimuth, calc_azimuth_east, calc_azimuth_north, calc_heading, calc_course, + cl_cd, winch_force, spring_forces, lift_drag, lift_over_drag, v_wind_kite, + kite_ref_frame, orient_euler, upwind_dir + +# --- Low level simulation interface --- +export clear!, find_steady_state!, residual! + +# --- Helper functions --- +export copy_examples, copy_bin, calc_drag, calculate_rotational_inertia!, calc_set_cl_cd!, + calc_aero_forces!, calc_particle_forces!, inner_loop!, loop! set_zero_subnormals(true) # required to avoid drastic slow down on Intel CPUs when numbers become very small @@ -97,9 +146,6 @@ function __init__() end include("KPS4.jl") # include code, specific for the four point kite model -include("system_structure.jl") -include("symbolic_awe_model.jl") # include code, specific for the ram air kite model -include("mtk_model.jl") include("KPS3.jl") # include code, specific for the one point kite model include("init.jl") # functions to calculate the initial state vector, the initial masses and initial springs @@ -216,7 +262,6 @@ function upwind_dir(v_wind_gnd) wind_dir = atan(v_wind_gnd[2], v_wind_gnd[1]) -(wind_dir + π/2) end -@register_symbolic upwind_dir(v_wind_gnd) """ tether_length(s::AKM) @@ -445,7 +490,7 @@ end # var_16::MyFloat # end -function update_sys_state!(ss::SysState, s::AKM, zoom=1.0) +function KiteUtils.update_sys_state!(ss::SysState, s::AKM, zoom=1.0) ss.time = s.t_0 pos = s.pos P = length(pos) diff --git a/src/init.jl b/src/init.jl index f9b35864..f65dfe79 100644 --- a/src/init.jl +++ b/src/init.jl @@ -144,15 +144,6 @@ function init(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES-1)+1); old=fal MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}(res1), MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}(res2) end -# rotate a 3d vector around the x axis in the yz plane - following the right hand rule -function rotate_around_x(vec, angle::T) where T - result = zeros(T, 3) - result[1] = vec[1] - result[2] = cos(angle) * vec[2] - sin(angle) * vec[3] - result[3] = sin(angle) * vec[2] + cos(angle) * vec[3] - result -end - # rotate a 3d vector around the y axis in the xz plane - following the right hand rule function rotate_around_y(vec, angle::T) where T result = zeros(T, 3) @@ -161,12 +152,3 @@ function rotate_around_y(vec, angle::T) where T result[3] = -sin(angle) * vec[1] + cos(angle) * vec[3] result end - -# rotate a 3d vector around the z axis in the yx plane - following the right hand rule -function rotate_around_z(vec, angle::T) where T - result = zeros(T, 3) - result[1] = cos(angle) * vec[1] - sin(angle) * vec[2] - result[2] = sin(angle) * vec[1] + cos(angle) * vec[2] - result[3] = vec[3] - result -end diff --git a/src/mtk_model.jl b/src/mtk_model.jl deleted file mode 100644 index 6b34c5a9..00000000 --- a/src/mtk_model.jl +++ /dev/null @@ -1,954 +0,0 @@ -# Copyright (c) 2025 Bart van de Lint -# SPDX-License-Identifier: MPL-2.0 - -# ==================== mtk model functions ================================================ -# Implementation of the ram air wing model using ModelingToolkit.jl - -function calc_speed_acc(winch::AsyncMachine, tether_vel, norm_, set_speed) - calc_acceleration(winch, tether_vel, norm_; set_speed, set_torque=nothing, use_brake=false) # TODO: add brake setting -end -function calc_moment_acc(winch::TorqueControlledMachine, tether_vel, norm_, set_torque) - calc_acceleration(winch, tether_vel, norm_; set_speed=nothing, set_torque, use_brake=false) -end - -function calc_heading(e_x) - return atan(e_x[2], e_x[1]) -end - -function calc_angle_of_attack(va_wing_b) - return atan(va_wing_b[3], va_wing_b[1]) -end - -function sym_normalize(vec) - return vec / norm(vec) -end - -function quaternion_to_rotation_matrix(q) - w, x, y, z = q[1], q[2], q[3], q[4] - - return [ - 1 - 2*(y*y + z*z) 2*(x*y - z*w) 2*(x*z + y*w); - 2*(x*y + z*w) 1 - 2*(x*x + z*z) 2*(y*z - x*w); - 2*(x*z - y*w) 2*(y*z + x*w) 1 - 2*(x*x + y*y) - ] -end - -function rotation_matrix_to_quaternion(R) - tr_ = tr(R) - - if tr_ > 0 - S = sqrt(tr_ + 1.0) * 2 - w = 0.25 * S - x = (R[3,2] - R[2,3]) / S - y = (R[1,3] - R[3,1]) / S - z = (R[2,1] - R[1,2]) / S - elseif (R[1,1] > R[2,2]) && (R[1,1] > R[3,3]) - S = sqrt(1.0 + R[1,1] - R[2,2] - R[3,3]) * 2 - w = (R[3,2] - R[2,3]) / S - x = 0.25 * S - y = (R[1,2] + R[2,1]) / S - z = (R[1,3] + R[3,1]) / S - elseif R[2,2] > R[3,3] - S = sqrt(1.0 + R[2,2] - R[1,1] - R[3,3]) * 2 - w = (R[1,3] - R[3,1]) / S - x = (R[1,2] + R[2,1]) / S - y = 0.25 * S - z = (R[2,3] + R[3,2]) / S - else - S = sqrt(1.0 + R[3,3] - R[1,1] - R[2,2]) * 2 - w = (R[2,1] - R[1,2]) / S - x = (R[1,3] + R[3,1]) / S - y = (R[2,3] + R[3,2]) / S - z = 0.25 * S - end - - return [w, x, y, z] -end - -function get_pos_w(sys_struct::SystemStructure, idx::Int16) - return sys_struct.points[idx].pos_w -end -@register_array_symbolic get_pos_w(sys::SystemStructure, point_idx::Int16) begin - size=size(KVec3) - eltype=SimFloat -end - -""" - force_eqs!(s, system, eqs, defaults, guesses; kwargs...) - -Generate the force equations for the wing system including spring forces, drag forces, -pulley dynamics and winch forces. - -# Arguments -- `s::SymbolicAWEModel`: The wing system state -- `system::SystemStructure`: The point mass representation -- `eqs`: Current system equations -- `defaults`: Default values for variables -- `guesses`: Initial guesses for variables -- `R_b_w`: Body to world rotation matrix -- `wing_pos`: Kite position vector -- `wing_vel`: Kite velocity vector -- `wind_vec_gnd`: Ground wind vector -- `group_aero_moment`: Aerodynamic moments per group -- `twist_angle`: Twist angles per group -- `stabilize`: Whether in stabilize mode - -# Returns -Tuple containing: -- Updated equations -- Updated defaults -- Updated guesses -- Tether forces on wing -- Tether moments on wing -""" -function force_eqs!(s, system, eqs, defaults, guesses; - R_b_w, wing_pos, wing_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff) - - @parameters acc_multiplier = 1 - - @unpack points, groups, segments, pulleys, tethers, winches, wings = system - - # ==================== POINTS ==================== # - tether_wing_force = zeros(Num, length(wings), 3) - tether_wing_moment = zeros(Num, length(wings), 3) - @parameters psys::SystemStructure = system - @variables begin - pos(t)[1:3, eachindex(points)] - vel(t)[1:3, eachindex(points)] - acc(t)[1:3, eachindex(points)] - point_force(t)[1:3, eachindex(points)] - tether_r(t)[1:3, eachindex(points)] - point_mass(t)[eachindex(points)] - - spring_force_vec(t)[1:3, eachindex(segments)] - drag_force(t)[1:3, eachindex(segments)] - l0(t)[eachindex(segments)] - end - for point in points - F::Vector{Num} = zeros(Num, 3) - mass = 0.0 - in_bridle = false - for segment in segments - if point.idx in segment.point_idxs - mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2 - inverted = segment.point_idxs[2] == point.idx - if inverted - F .-= spring_force_vec[:, segment.idx] - else - F .+= spring_force_vec[:, segment.idx] - end - mass += mass_per_meter * l0[segment.idx] / 2 - F .+= 0.5drag_force[:, segment.idx] - - if segment.type == BRIDLE - in_bridle = true - end - end - end - @assert !iszero(mass) - - eqs = [ - eqs - point_mass[point.idx] ~ mass - point_force[:, point.idx] ~ F - ] - - if point.type == WING - found = 0 - group = nothing - for group_ in groups - if point.idx in group_.point_idxs - group = group_ - found += 1 - end - end - !(found in [0,1]) && error("Kite point number $(point.idx) is part of $found groups, - and should be part of exactly 0 or 1 groups.") - - if found == 1 - found = 0 - wing = nothing - for wing_ in wings - if group.idx in wing_.group_idxs - wing = wing_ - found += 1 - end - end - !(found == 1) && error("Kite group number $(group.idx) is part of $found wings, - and should be part of exactly 1 wing.") - - fixed_pos = group.le_pos - chord_b = point.pos_b - fixed_pos - normal = chord_b × group.y_airf - pos_b = fixed_pos + cos(twist_angle[group.idx]) * chord_b - sin(twist_angle[group.idx]) * normal - elseif found == 0 - pos_b = point.pos_b - eqs = [ - eqs - tether_r[:, point.idx] ~ pos[:, point.idx] - wing_pos[point.wing_idx, :] - ] - tether_wing_moment[point.wing_idx, :] .+= tether_r[:, point.idx] × point_force[:, point.idx] - end - tether_wing_force[point.wing_idx, :] .+= point_force[:, point.idx] - - eqs = [ - eqs - pos[:, point.idx] ~ wing_pos[point.wing_idx, :] + R_b_w[point.wing_idx, :, :] * pos_b - vel[:, point.idx] ~ zeros(3) - acc[:, point.idx] ~ zeros(3) - ] - elseif point.type == STATIC - eqs = [ - eqs - pos[:, point.idx] ~ get_pos_w(psys, point.idx) - vel[:, point.idx] ~ zeros(3) - acc[:, point.idx] ~ zeros(3) - ] - elseif point.type == DYNAMIC - # p = pos[:, point.idx] - # n = sym_normalize(wing_pos) - # n = n * (p ⋅ n) - # r = (p - n) # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation - @parameters bridle_damp = 1.0 - @parameters measured_ω_z = 0.6 - if in_bridle && length(wings) > 0 - bridle_damp_vec = bridle_damp * (vel[:, point.idx] - wing_vel[point.wing_idx, :]) - else - bridle_damp_vec = zeros(Num, 3) - end - eqs = [ - eqs - D(pos[:, point.idx]) ~ vel[:, point.idx] - D(vel[:, point.idx]) ~ acc_multiplier * acc[:, point.idx] - bridle_damp_vec - acc[:, point.idx] ~ point_force[:, point.idx] / mass + [0, 0, -G_EARTH] - # ifelse.(stabilize==true, r * norm(measured_ω_z)^2, zeros(3)) - ] - defaults = [ - defaults - [pos[j, point.idx] => get_pos_w(psys, point.idx)[j] for j in 1:3] - [vel[j, point.idx] => 0 for j in 1:3] - ] - elseif point.type == QUASI_STATIC - eqs = [ - eqs - vel[:, point.idx] ~ zeros(3) - acc[:, point.idx] ~ zeros(3) - acc[:, point.idx] ~ point_force[:, point.idx] / mass + [0, 0, -G_EARTH] - ] - guesses = [ - guesses - [acc[j, point.idx] => 0 for j in 1:3] - [pos[j, point.idx] => get_pos_w(psys, point.idx)[j] for j in 1:3] - [point_force[j, point.idx] => 0 for j in 1:3] - ] - else - error("Unknown point type: $(typeof(point))") - end - end - - # ==================== GROUPS ==================== # - if length(groups) > 0 - @variables begin - trailing_edge_angle(t)[eachindex(groups)] # angle left / right - trailing_edge_ω(t)[eachindex(groups)] # angular rate - trailing_edge_α(t)[eachindex(groups)] # angular acc - free_twist_angle(t)[eachindex(groups)] - twist_α(t)[eachindex(groups)] # angular acc - group_tether_moment(t)[eachindex(groups)] - tether_force(t)[eachindex(groups), eachindex(groups[1].point_idxs)] - tether_moment(t)[eachindex(groups), eachindex(groups[1].point_idxs)] - end - end - - for group in groups - found = 0 - wing = nothing - for wing_ in wings - if group.idx in wing_.group_idxs - wing = wing_ - found += 1 - end - end - !(found == 1) && error("Kite group number $(group.idx) is part of $found wings, - and should be part of exactly 1 wing.") - - all(iszero.(tether_wing_moment[wing.idx, :])) && - error("Tether wing moment is zero. At least one of the wing connection points should not be part of a deforming group.") - - x_airf = normalize(group.chord) - init_z_airf = x_airf × group.y_airf - z_airf = x_airf * sin(twist_angle[group.idx]) + init_z_airf * cos(twist_angle[group.idx]) - for (i, point_idx) in enumerate(group.point_idxs) - r = (points[point_idx].pos_b - (group.le_pos + group.moment_frac*group.chord)) ⋅ normalize(group.chord) - eqs = [ - eqs - tether_force[group.idx, i] ~ (point_force[:, point_idx] ⋅ (R_b_w[wing.idx, :, :] * -z_airf)) - tether_moment[group.idx, i] ~ r * tether_force[group.idx, i] - ] - end - - inertia = 1/3 * (s.set.mass/length(groups)) * (norm(group.chord))^2 # plate inertia around leading edge - @assert !(inertia ≈ 0.0) - @parameters twist_damp = 50 - @parameters max_twist = deg2rad(90) - - eqs = [ - eqs - group_tether_moment[group.idx] ~ sum(tether_moment[group.idx, :]) - twist_α[group.idx] ~ (group_aero_moment[group.idx] + group_tether_moment[group.idx]) / inertia - twist_angle[group.idx] ~ clamp(free_twist_angle[group.idx], -max_twist, max_twist) - ] - if group.type == DYNAMIC - eqs = [ - eqs - D(free_twist_angle[group.idx]) ~ ifelse(fix_nonstiff==true, 0, twist_ω[group.idx]) - D(twist_ω[group.idx]) ~ ifelse(fix_nonstiff==true, 0, twist_α[group.idx] - twist_damp * twist_ω[group.idx]) - ] - defaults = [ - defaults - free_twist_angle[group.idx] => 0 - twist_ω[group.idx] => 0 - ] - elseif group.type == QUASI_STATIC - eqs = [ - eqs - twist_ω[group.idx] ~ 0 - twist_α[group.idx] ~ 0 - ] - guesses = [ - guesses - free_twist_angle[group.idx] => 0 - twist_angle[group.idx] => 0 - ] - else - error("Wrong group type.") - end - end - - # ==================== SEGMENTS ==================== # - @variables begin - segment_vec(t)[1:3, eachindex(segments)] - unit_vec(t)[1:3, eachindex(segments)] - len(t)[eachindex(segments)] - rel_vel(t)[1:3, eachindex(segments)] - spring_vel(t)[eachindex(segments)] - spring_force(t)[eachindex(segments)] - stiffness(t)[eachindex(segments)] - damping(t)[eachindex(segments)] - - height(t)[eachindex(segments)] - segment_vel(t)[1:3, eachindex(segments)] - segment_rho(t)[eachindex(segments)] - wind_vel(t)[1:3, eachindex(segments)] - va(t)[1:3, eachindex(segments)] - area(t)[eachindex(segments)] - app_perp_vel(t)[1:3, eachindex(segments)] - drag_force(t)[1:3, eachindex(segments)] - - pulley_l0(t)[eachindex(pulleys)] - - tether_length(t)[eachindex(winches)] - end - for segment in segments - p1, p2 = segment.point_idxs[1], segment.point_idxs[2] - # if s.set.quasi_static - guesses = [ - guesses - [segment_vec[i, segment.idx] => points[p2].pos_w[i] - points[p1].pos_w[i] for i in 1:3] - ] - # end - - if segment.type == BRIDLE - in_pulley = 0 - for pulley in pulleys - if segment.idx == pulley.segment_idxs[1] # each bridle segment has to be part of no pulley or one pulley - eqs = [ - eqs - l0[segment.idx] ~ pulley_l0[pulley.idx] - ] - in_pulley += 1 - end - if segment.idx == pulley.segment_idxs[2] - eqs = [ - eqs - l0[segment.idx] ~ pulley.sum_length - pulley_l0[pulley.idx] - ] - in_pulley += 1 - end - end - if in_pulley == 0 - eqs = [ - eqs - l0[segment.idx] ~ segment.l0 - ] - end - (in_pulley > 1) && error("Bridle segment number $(segment.idx) is part of - $in_pulley pulleys, and should be part of either 0 or 1 pulleys.") - elseif segment.type == POWER_LINE || segment.type == STEERING_LINE - in_tether = 0 - for tether in tethers - if segment.idx in tether.segment_idxs # each tether segment has to be part of exactly one tether - in_winch = 0 - winch_idx = 0 - for winch in winches - if tether.idx in winch.tether_idxs - winch_idx = winch.idx - in_winch += 1 - end - end - !(in_winch in [0,1]) && error("Tether number $(tether.idx) is part of - $(in_winch) winches, and should be part of exactly 0 or 1 winch.") - - if in_winch == 1 - eqs = [ - eqs - l0[segment.idx] ~ tether_length[winch_idx] / length(tether.segment_idxs) - ] - elseif in_winch == 0 - eqs = [ - eqs - l0[segment.idx] ~ segment.l0 - ] - end - in_tether += 1 - end - end - (in_tether != 1) && error("Segment number $(segment.idx) is part of - $in_tether tethers, and should be part of exactly 1 tether.") - else - error("Unknown segment type: $(segment.type)") - end - - stiffness_m = s.set.e_tether * (segment.diameter/2)^2 * pi - @parameters stiffness_frac = 0.01 - (segment.type == BRIDLE) && (stiffness_m = stiffness_frac * stiffness_m) - - damping_m = (s.set.damping / s.set.c_spring) * stiffness_m - - if segment.compression_frac ≈ 1.0 - eqs = [eqs; stiffness[segment.idx] ~ stiffness_m / len[segment.idx]] - else - eqs = [eqs; stiffness[segment.idx] ~ ifelse(len[segment.idx] > l0[segment.idx], - stiffness_m / len[segment.idx], - segment.compression_frac * stiffness_m / len[segment.idx])] - end - eqs = [ - eqs - # spring force equations - segment_vec[:, segment.idx] ~ pos[:, p2] - pos[:, p1] - len[segment.idx] ~ norm(segment_vec[:, segment.idx]) - unit_vec[:, segment.idx] ~ segment_vec[:, segment.idx]/len[segment.idx] - rel_vel[:, segment.idx] ~ vel[:, p1] - vel[:, p2] - spring_vel[segment.idx] ~ rel_vel[:, segment.idx] ⋅ unit_vec[:, segment.idx] - damping[segment.idx] ~ damping_m / len[segment.idx] - spring_force[segment.idx] ~ (stiffness[segment.idx] * (len[segment.idx] - l0[segment.idx]) - - damping[segment.idx] * spring_vel[segment.idx]) - spring_force_vec[:, segment.idx] ~ spring_force[segment.idx] * unit_vec[:, segment.idx] - - # drag force equations - height[segment.idx] ~ max(0.0, 0.5(pos[:, p1][3] + pos[:, p2][3])) - segment_vel[:, segment.idx] ~ 0.5(vel[:, p1] + vel[:, p2]) - segment_rho[segment.idx] ~ calc_rho(s.am, height[segment.idx]) - wind_vel[:, segment.idx] ~ AtmosphericModels.calc_wind_factor(s.am, height[segment.idx], s.set.profile_law) * wind_vec_gnd - va[:, segment.idx] ~ wind_vel[:, segment.idx] - segment_vel[:, segment.idx] - area[segment.idx] ~ len[segment.idx] * segment.diameter - app_perp_vel[:, segment.idx] ~ va[:, segment.idx] - - (va[:, segment.idx] ⋅ unit_vec[:, segment.idx]) * unit_vec[:, segment.idx] - drag_force[:, segment.idx] ~ (0.5 * segment_rho[segment.idx] * s.set.cd_tether * norm(va[:, segment.idx]) * - area[segment.idx]) * app_perp_vel[:, segment.idx] - ] - end - - # ==================== PULLEYS ==================== # - @variables begin - pulley_l0(t)[eachindex(pulleys)] - pulley_vel(t)[eachindex(pulleys)] - pulley_force(t)[eachindex(pulleys)] - pulley_acc(t)[eachindex(pulleys)] - end - @parameters pulley_damp = 1.0 - for pulley in pulleys - segment = segments[pulley.segment_idxs[1]] - mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2 - mass = pulley.sum_length * mass_per_meter - @assert !(mass ≈ 0.0) - eqs = [ - eqs - pulley_force[pulley.idx] ~ spring_force[pulley.segment_idxs[1]] - spring_force[pulley.segment_idxs[2]] - pulley_acc[pulley.idx] ~ pulley_force[pulley.idx] / mass - ] - if pulley.type == DYNAMIC - eqs = [ - eqs - D(pulley_l0[pulley.idx]) ~ pulley_vel[pulley.idx] - D(pulley_vel[pulley.idx]) ~ acc_multiplier * pulley_acc[pulley.idx] - pulley_damp * pulley_vel[pulley.idx] - ] - defaults = [ - defaults - pulley_l0[pulley.idx] => segments[pulley.segment_idxs[1]].l0 - pulley_vel[pulley.idx] => 0 - ] - elseif pulley.type == QUASI_STATIC - eqs = [ - eqs - pulley_vel[pulley.idx] ~ 0 - pulley_acc[pulley.idx] ~ 0 - ] - guesses = [ - guesses - pulley_l0[pulley.idx] => segments[pulley.segment_idxs[1]].l0 - ] - else - error("Wrong pulley type") - end - end - - # ==================== WINCHES ==================== # - @variables begin - tether_vel(t)[eachindex(winches)] - tether_acc(t)[eachindex(winches)] - winch_force(t)[eachindex(winches)] - end - for winch in winches - F = zero(Num) - for tether_idx in winch.tether_idxs - found = 0 - point_idx = 0 - for segment_idx in tethers[tether_idx].segment_idxs - segment = segments[segment_idx] - for point_idx_ in segment.point_idxs - if points[point_idx_].type == STATIC && point_idx != point_idx_ - found += 1 - point_idx = point_idx_ - end - end - end - (found != 1) && error("Tether number $tether_idx has $found static points and should have exactly 1 static point.") - F += norm(point_force[:, point_idx]) - end - eqs = [ - eqs - D(tether_length[winch.idx]) ~ ifelse(fix_nonstiff==true, 0, ifelse(stabilize==true, 0, tether_vel[winch.idx])) - D(tether_vel[winch.idx]) ~ ifelse(fix_nonstiff==true, 0, ifelse(stabilize==true, 0, tether_acc[winch.idx])) - - tether_acc[winch.idx] ~ calc_moment_acc( # TODO: moment and speed control - winch.model, tether_vel[winch.idx], - winch_force[winch.idx], - set_values[winch.idx] - ) - winch_force[winch.idx] ~ F - ] - defaults = [ - defaults - tether_length[winch.idx] => winch.tether_length - tether_vel[winch.idx] => 0 - ] - end - - # ==================== TETHERS ==================== # - @variables begin - stretched_length(t)[eachindex(tethers)] - end - for tether in tethers - slen = zero(Num) - for segment_idx in tether.segment_idxs - slen += len[segment_idx] - end - eqs = [ - eqs - stretched_length[tether.idx] ~ slen - ] - end - - return eqs, defaults, guesses, tether_wing_force, tether_wing_moment -end - -""" - wing_eqs!(s, eqs, defaults; kwargs...) - -Generate the differential equations for wing dynamics including quaternion kinematics, -angular velocities and accelerations, and forces/moments. - -# Arguments -- `s::SymbolicAWEModel`: The wing system state -- `eqs`: Current system equations -- `defaults`: Default values for variables -- `tether_wing_force`: Forces from tethers on wing -- `tether_wing_moment`: Moments from tethers on wing -- `aero_force_b`: Aerodynamic forces in body frame -- `aero_moment_b`: Aerodynamic moments in body frame -- `ω_b`: Angular velocity in body frame -- `R_b_w`: Body to world rotation matrix -- `wing_pos`: Kite position vector -- `wing_vel`: Kite velocity vector -- `wing_acc`: Kite acceleration vector -- `stabilize`: Whether in stabilize mode - -# Returns -Tuple of updated equations and defaults -""" -function wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero_force_b, - aero_moment_b, ω_b, α_b, R_b_w, wing_pos, wing_vel, wing_acc, stabilize, fix_nonstiff -) - wings = s.sys_struct.wings - @variables begin - # potential differential variables - wing_acc_b(t)[eachindex(wings), 1:3] - α_b_damped(t)[eachindex(wings), 1:3] - ω_b_stable(t)[eachindex(wings), 1:3] - - # rotations and frames - Q_b_w(t)[eachindex(wings), 1:4] # quaternion orientation of the wing body frame relative to the world frame - Q_vel(t)[eachindex(wings), 1:4] # quaternion rate of change - - # rest: forces, moments, vectors and scalar values - moment_b(t)[eachindex(wings), 1:3] # moment in principal frame - end - @parameters ω_damp = 150 - - Ω(ω) = [0 -ω[1] -ω[2] -ω[3]; - ω[1] 0 ω[3] -ω[2]; - ω[2] -ω[3] 0 ω[1]; - ω[3] ω[2] -ω[1] 0] - - for wing in wings - vsm_wing = s.vsm_wings[wing.idx] - I_b = [vsm_wing.inertia_tensor[1,1], vsm_wing.inertia_tensor[2,2], vsm_wing.inertia_tensor[3,3]] - @assert !(s.set.mass ≈ 0) - axis = sym_normalize(wing_pos[wing.idx, :]) - axis_b = R_b_w[wing.idx, :, :]' * axis - eqs = [ - eqs - [D(Q_b_w[wing.idx, i]) ~ Q_vel[wing.idx, i] for i in 1:4] - [Q_vel[wing.idx, i] ~ 0.5 * sum(Ω(ω_b_stable[wing.idx, :])[i, j] * Q_b_w[wing.idx, j] for j in 1:4) for i in 1:4] - ω_b_stable[wing.idx, :] ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - ω_b[wing.idx, :] - ω_b[wing.idx, :] ⋅ axis_b * axis_b, - ω_b[wing.idx, :] - ) - ) - D(ω_b[wing.idx, :]) ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - α_b_damped[wing.idx, :] - α_b_damped[wing.idx, :] ⋅ axis_b * axis_b, - α_b_damped[wing.idx, :] - ) - ) - α_b_damped[wing.idx, :] ~ [α_b[wing.idx, 1], α_b[wing.idx, 2] - ω_damp*ω_b[wing.idx, 2], α_b[wing.idx, 3]] - - [R_b_w[wing.idx, :, i] ~ quaternion_to_rotation_matrix(Q_b_w[wing.idx, :])[:, i] for i in 1:3] - - α_b[wing.idx, 1] ~ (moment_b[wing.idx, 1] + (I_b[2] - I_b[3]) * ω_b[wing.idx, 2] * ω_b[wing.idx, 3]) / I_b[1] - α_b[wing.idx, 2] ~ (moment_b[wing.idx, 2] + (I_b[3] - I_b[1]) * ω_b[wing.idx, 3] * ω_b[wing.idx, 1]) / I_b[2] - α_b[wing.idx, 3] ~ (moment_b[wing.idx, 3] + (I_b[1] - I_b[2]) * ω_b[wing.idx, 1] * ω_b[wing.idx, 2]) / I_b[3] - moment_b[wing.idx, :] ~ aero_moment_b[wing.idx, :] + R_b_w[wing.idx, :, :]' * tether_wing_moment[wing.idx, :] - - D(wing_pos[wing.idx, :]) ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - wing_vel[wing.idx, :] ⋅ axis * axis, - wing_vel[wing.idx, :] - ) - ) - D(wing_vel[wing.idx, :]) ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - wing_acc[wing.idx, :] ⋅ axis * axis, - wing_acc[wing.idx, :] - ) - ) - wing_acc[wing.idx, :] ~ (tether_wing_force[wing.idx, :] + R_b_w[wing.idx, :, :] * aero_force_b[wing.idx, :]) / s.set.mass - ] - defaults = [ - defaults - [Q_b_w[wing.idx, i] => wing.orient[i] for i in 1:4] - [ω_b[wing.idx, i] => wing.angular_vel[i] for i in 1:3] - [wing_pos[wing.idx, i] => wing.pos_w[i] for i in 1:3] - [wing_vel[wing.idx, i] => wing.vel_w[i] for i in 1:3] - ] - end - - return eqs, defaults -end - -function rotate_v_around_k(v, k, θ) - k = sym_normalize(k) - v_rot = v * cos(θ) + (k × v) * sin(θ) + k * (k ⋅ v) * (1 - cos(θ)) - return v_rot -end - -function calc_R_v_w(wing_pos, e_x) - z = sym_normalize(wing_pos) - y = sym_normalize(z × e_x) - x = y × z - return [x y z] -end - -function calc_R_t_w(elevation, azimuth) - x = rotate_around_z(rotate_around_y([0, 0, -1], -elevation), azimuth) - z = rotate_around_z(rotate_around_y([1, 0, 0], -elevation), azimuth) - y = z × x - return [x y z] -end - -""" - scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω) - -Generate equations for scalar quantities like elevation, azimuth, heading and course angles. - - # Arguments - - `s::SymbolicAWEModel`: The wing system state - - `eqs`: Current system equations - - `R_b_w`: Body to world rotation matrix - - `wind_vec_gnd`: Ground wind vector - - `va_wing_b`: Apparent wind velocity in body frame - - `wing_pos`: Kite position vector - - `wing_vel`: Kite velocity vector - - `wing_acc`: Kite acceleration vector - - # Returns - - Updated system equations including: - - Heading angle from x-axis - - Elevation angle - - Azimuth angle - - Course angle - - Angular velocities and accelerations - """ -function scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b) - @unpack wings = s.sys_struct - @parameters wind_scale_gnd = s.set.v_wind - @parameters upwind_dir = deg2rad(s.set.upwind_dir) - @variables begin - e_x(t)[eachindex(wings), 1:3] - e_y(t)[eachindex(wings), 1:3] - e_z(t)[eachindex(wings), 1:3] - wind_vel_wing(t)[eachindex(wings), 1:3] - va_wing(t)[eachindex(wings), 1:3] - end - eqs = [ - eqs - wind_vec_gnd ~ wind_scale_gnd * rotate_around_z([0, -1, 0], -upwind_dir) - ] - for wing in wings - eqs = [ - eqs - e_x[wing.idx, :] ~ R_b_w[wing.idx, :,1] - e_y[wing.idx, :] ~ R_b_w[wing.idx, :,2] - e_z[wing.idx, :] ~ R_b_w[wing.idx, :,3] - wind_vel_wing[wing.idx, :] ~ AtmosphericModels.calc_wind_factor(s.am, wing_pos[wing.idx, 3], s.set.profile_law) * wind_vec_gnd - va_wing[wing.idx, :] ~ wind_vel_wing[wing.idx, :] - wing_vel[wing.idx, :] - va_wing_b[wing.idx, :] ~ R_b_w[wing.idx, :, :]' * va_wing[wing.idx, :] - ] - end - @variables begin - heading(t)[eachindex(wings)] - turn_rate(t)[eachindex(wings), 1:3] - turn_acc(t)[eachindex(wings), 1:3] - azimuth(t)[eachindex(wings)] - azimuth_vel(t)[eachindex(wings)] - azimuth_acc(t)[eachindex(wings)] - elevation(t)[eachindex(wings)] - elevation_vel(t)[eachindex(wings)] - elevation_acc(t)[eachindex(wings)] - course(t)[eachindex(wings)] - x_acc(t)[eachindex(wings)] - y_acc(t)[eachindex(wings)] - sphere_pos(t)[eachindex(wings), 1:2, 1:2] - sphere_vel(t)[eachindex(wings), 1:2, 1:2] - sphere_acc(t)[eachindex(wings), 1:2, 1:2] - angle_of_attack(t)[eachindex(wings)] - R_v_w(t)[eachindex(wings), 1:3, 1:3] - R_t_w(t)[eachindex(wings), 1:3, 1:3] - distance(t)[eachindex(wings)] - distance_vel(t)[eachindex(wings)] - distance_acc(t)[eachindex(wings)] - end - - for wing in wings - x, y, z = wing_pos[wing.idx, :] - x´, y´, z´ = wing_vel[wing.idx, :] - x´´, y´´, z´´ = wing_acc[wing.idx, :] - - half_len = wing.group_idxs[1] + length(wing.group_idxs)÷2 - 1 - heading_vec = R_t_w[wing.idx, :, :]' * R_v_w[wing.idx, :, 1] - - eqs = [ - eqs - vec(R_v_w[wing.idx, :, :]) .~ vec(calc_R_v_w(wing_pos[wing.idx, :], e_x[wing.idx, :])) - vec(R_t_w[wing.idx, :, :]) .~ vec(calc_R_t_w(elevation[wing.idx], azimuth[wing.idx])) - heading[wing.idx] ~ atan(heading_vec[2], heading_vec[1]) - turn_rate[wing.idx, :] ~ R_v_w[wing.idx, :, :]' * (R_b_w[wing.idx, :, :] * ω_b[wing.idx, :]) # Project angular velocity onto view frame - turn_acc[wing.idx, :] ~ R_v_w[wing.idx, :, :]' * (R_b_w[wing.idx, :, :] * α_b[wing.idx, :]) - distance[wing.idx] ~ norm(wing_pos[wing.idx, :]) - distance_vel[wing.idx] ~ wing_vel[wing.idx, :] ⋅ R_v_w[wing.idx, :, 3] - distance_acc[wing.idx] ~ wing_acc[wing.idx, :] ⋅ R_v_w[wing.idx, :, 3] - - elevation[wing.idx] ~ KiteUtils.calc_elevation(wing_pos[wing.idx, :]) - # elevation_vel = d/dt(atan(z/x)) = (x*ż' - z*ẋ')/(x^2 + z^2) according to wolframalpha - elevation_vel[wing.idx] ~ (x*z´ - z*x´) / - (x^2 + z^2) - elevation_acc[wing.idx] ~ ((x^2 + z^2)*(x*z´´ - z*x´´) + 2(z*x´ - x*z´)*(x*x´ + z*z´))/(x^2 + z^2)^2 - azimuth[wing.idx] ~ -KiteUtils.azimuth_east(wing_pos[wing.idx, :]) - # azimuth_vel = d/dt(atan(y/x)) = (-y*x´ + x*y´)/(x^2 + y^2) # TODO: check if correct - azimuth_vel[wing.idx] ~ (-y*x´ + x*y´) / - (x^2 + y^2) - azimuth_acc[wing.idx] ~ ((x^2 + y^2)*(-y*x´´ + x*y´´) + 2(y*x´ - x*y´)*(x*x´ + y*y´))/(x^2 + y^2)^2 - course[wing.idx] ~ atan(-azimuth_vel[wing.idx], elevation_vel[wing.idx]) - x_acc[wing.idx] ~ wing_acc ⋅ e_x - y_acc[wing.idx] ~ wing_acc ⋅ e_y - - angle_of_attack[wing.idx] ~ calc_angle_of_attack(va_wing_b[wing.idx, :]) + - 0.5twist_angle[half_len] + 0.5twist_angle[half_len+1] - ] - end - return eqs -end - -function Base.getindex(x::ModelingToolkit.Symbolics.SymArray, idxs::Vector{Int16}) - Num[Base.getindex(x, idx) for idx in idxs] -end - -""" -linear_vsm_eqs!(s, eqs; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_wing_b, ω_b) - -Generate linearized aerodynamic equations using the Vortex Step Method (VSM). -Uses linearization around current operating point to approximate aerodynamic forces -and moments. The Jacobian is computed using the VSM solver. - -# Arguments -- `s::SymbolicAWEModel`: The wing system state -- `eqs`: Current system equations -- `aero_force_b`: Aerodynamic forces in body frame -- `aero_moment_b`: Aerodynamic moments in body frame -- `group_aero_moment`: Aerodynamic moments per group -- `init_va_b`: Initial apparent wind velocity -- `twist_angle`: Twist angles per group -- `va_wing_b`: Apparent wind velocity in body frame -- `ω_b`: Angular velocity in body frame - -# Returns -- Updated system equations including linearized aerodynamics: -- Force and moment calculations -- Group moment distributions -- Jacobian matrix for state derivatives -""" -function linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_wing_b, ω_b) - @unpack groups, wings = s.sys_struct - if length(wings) == 0 - return eqs, guesses - end - - ny = 3+length(wings[1].group_idxs)+3 - nx = 3+3+length(wings[1].group_idxs) - y_ = zeros(length(wings), ny) - x_ = zeros(length(wings), nx) - jac_ = zeros(length(wings), nx, ny) - for wing in wings - y_[wing.idx, :] .= [init_va_b[wing.idx, :]; zeros(length(wing.group_idxs)); zeros(3)] - end - @parameters begin - last_y[eachindex(wings), 1:ny] = y_ - last_x[eachindex(wings), 1:nx] = x_ - vsm_jac[eachindex(wings), 1:nx, 1:ny] = jac_ - end - - @variables begin - y(t)[eachindex(wings), 1:ny] - dy(t)[eachindex(wings), 1:ny] - end - - for wing in wings - last_y_ = [last_y[wing.idx, i] for i in 1:ny] # https://github.com/SciML/ModelingToolkit.jl/issues/3730 - last_x_ = [last_x[wing.idx, i] for i in 1:nx] - vsm_jac_ = [vsm_jac[wing.idx, i, j] for i in 1:nx, j in 1:ny] - eqs = [ - eqs - y[wing.idx, :] ~ [va_wing_b[wing.idx, :]; ω_b[wing.idx, :]; twist_angle[wing.group_idxs]] - dy[wing.idx, :] ~ y[wing.idx, :] - last_y_ - [aero_force_b[wing.idx, :]; aero_moment_b[wing.idx, :]; group_aero_moment[wing.group_idxs]] ~ last_x_ + vsm_jac_ * dy[wing.idx, :] - ] - - if s.set.quasi_static - guesses = [guesses; [y[wing.idx, i] => y_[wing.idx, i] for i in 1:ny]] - end - end - return eqs, guesses -end - -function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b) - eqs = [] - defaults = Pair{Num, Any}[] - guesses = Pair{Num, Any}[] - - @unpack wings, groups, winches = system - - @parameters begin - stabilize = false - fix_nonstiff = false - end - @variables begin - # potential differential variables - set_values(t)[eachindex(winches)] = zeros(length(winches)) - wing_pos(t)[eachindex(wings), 1:3] # xyz pos of wing in world frame - wing_vel(t)[eachindex(wings), 1:3] - wing_acc(t)[eachindex(wings), 1:3] - ω_b(t)[eachindex(wings), 1:3] # turn rate in principal frame - α_b(t)[eachindex(wings), 1:3] - - # rotations and frames - R_b_w(t)[eachindex(wings), 1:3, 1:3] # rotation of the wing body frame relative to the world frame - - # rest: forces, moments, vectors and scalar values - aero_force_b(t)[eachindex(wings), 1:3] - aero_moment_b(t)[eachindex(wings), 1:3] - twist_angle(t)[eachindex(groups)] - twist_ω(t)[eachindex(groups)] - group_aero_moment(t)[eachindex(groups)] - wind_vec_gnd(t)[1:3] - va_wing_b(t)[eachindex(wings), 1:3] - end - - eqs, defaults, guesses, tether_wing_force, tether_wing_moment = - force_eqs!(s, system, eqs, defaults, guesses; - R_b_w, wing_pos, wing_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff) - eqs, guesses = linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_wing_b, ω_b) - eqs, defaults = wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero_force_b, aero_moment_b, - ω_b, α_b, R_b_w, wing_pos, wing_vel, wing_acc, stabilize, fix_nonstiff) - eqs = scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b) - - # te_I = (1/3 * (s.set.mass/8) * te_length^2) - # # -damping / I * ω = α_damping - # # solve for c: (c * (k*m/s^2) / (k*m^2)) * (m/s)=m/s^2 in wolframalpha - # # damping should be in N*m*s - # rot_damping = 0.1s.damping * te_length - - # eqs = [ - # eqs - # trailing_edge_α[1] ~ (force[:, s.i_A]) ⋅ e_te_A * te_length / te_I - (rot_damping[1] / te_I) * trailing_edge_ω[1] # TODO: add trailing edge - # trailing_edge_α[2] ~ (force[:, s.i_B]) ⋅ e_te_B * te_length / te_I - (rot_damping[2] / te_I) * trailing_edge_ω[2] - # ] - - eqs = Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs))) - - # discrete_events = [ - # true => [ - # [Q_b_w[i] ~ normalize(Q_b_w)[i] for i in 1:4] - # [twist_angle[i] ~ clamp(twist_angle[i], -π/2, π/2) for i in eachindex(s.point_groups)] - # ] - # ] - - @info "Creating ODESystem" - # @named sys = ODESystem(eqs, t; discrete_events) - @time @named sys = ODESystem(eqs, t) - - defaults = [ - defaults - [set_values[i] => [-50.0, -1.0, -1.0][i] for i in eachindex(winches)] - ] - - s.defaults = defaults - s.guesses = guesses - s.full_sys = sys - return set_values -end diff --git a/src/precompile.jl b/src/precompile.jl index e716be33..77a6191e 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -51,50 +51,5 @@ end # they belong to your package or not (on Julia 1.8 and higher) integrator = KiteModels.init!(kps3_; stiffness_factor=0.035, prn=false) integrator = KiteModels.init!(kps4_; delta=0.03, stiffness_factor=0.05, prn=false) - - sam_set = load_settings("system_ram.yaml") - sam_set.segments = 3 - set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] - sam_set.quasi_static = false - sam_set.physical_model = "ram" - model_name = get_model_name(sam_set) - model_file = normpath(joinpath(path, "..", "data", model_name)) - output_file = normpath(joinpath(path, "..", "data", model_name * ".default")) - input_file = normpath(joinpath(path, "..", "data", model_name * ".default.xz")) - if isfile(input_file) && ! isfile(output_file) - using CodecXz - decompress_binary(input_file, output_file) - @info "Decompressed $input_file to $output_file" - elseif isfile(output_file) - @info "Output file $output_file already exists, skipping decompression." - else - @error "Input file $input_file does not exist, skipping decompression." - end - if VERSION.minor == 11 - m1 = "Manifest-v1.11.toml" - m2 = "Manifest-v1.11.toml.default" - else - m1 = "Manifest-v1.10.toml" - m2 = "Manifest-v1.10.toml.default" - end - if filecmp(m1, m2) - @info "Manifest files match, using the default xz files will work!" - else - @warn "Manifest files differ, no precompilation will be done." - end - # Check if the output file exists and is the same as the input file - if isfile(output_file) && filecmp(m1, m2) - s = SymbolicAWEModel(sam_set) - - # Initialize at elevation - KiteModels.init!(s; prn=false, precompile=true) - @info "Copying $output_file to $model_file !" - cp(output_file, model_file; force=true) - find_steady_state!(s) - steps = Int(round(10 / 0.05)) - logger = Logger(length(s.sys_struct.points), steps) - sys_state = SysState(s) - end - nothing end -end \ No newline at end of file +end diff --git a/src/symbolic_awe_model.jl b/src/symbolic_awe_model.jl deleted file mode 100644 index 0adb6c22..00000000 --- a/src/symbolic_awe_model.jl +++ /dev/null @@ -1,886 +0,0 @@ -# Copyright (c) 2024, 2025 Bart van de Lint and Uwe Fechner -# SPDX-License-Identifier: MIT - -@with_kw mutable struct SerializedModel - set_hash::Vector{UInt8} - "Reference to the geometric wing model" - vsm_wings::Vector{VortexStepMethod.RamAirWing} - "Reference to the aerodynamic wing model" - vsm_aeros::Vector{VortexStepMethod.BodyAerodynamics} - "Reference to the VSM aerodynamics solver" - vsm_solvers::Vector{VortexStepMethod.Solver} - sys_struct_hash::Vector{UInt8} - "Simplified system of the mtk model" - sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing - "Unsimplified system of the mtk model" - full_sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing - "Linearization function of the mtk model" - lin_prob::Union{ModelingToolkit.LinearizationProblem, Nothing} = nothing - "ODE function of the mtk model" - prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing - - unknowns_vec::Vector{SimFloat} = zeros(SimFloat, 3) - defaults::Vector{Pair} = Pair[] - guesses::Vector{Pair} = Pair[] - - set_psys::Function = (_, _) -> nothing - set_set_values::Function = (_, _) -> nothing - set_wind::Function = (_, _) -> nothing - set_vsm::Function = (_, _) -> nothing - set_unknowns::Function = (_, _) -> nothing - set_nonstiff::Function = (_, _) -> nothing - set_lin_vsm::Function = (_, _) -> nothing - set_lin_set_values::Function = (_, _) -> nothing - set_lin_unknowns::Function = (_, _) -> nothing - set_stabilize::Function = (_, _) -> nothing - - get_vsm::Function = (_) -> nothing - get_set_values::Function = (_) -> nothing - get_unknowns::Function = (_) -> nothing - get_wing_state::Function = (_) -> nothing - get_winch_state::Function = (_) -> nothing - get_point_state::Function = (_) -> nothing - get_y::Function = (_) -> nothing - get_unstretched_length::Function = (_) -> nothing - get_tether_length::Function = (_) -> nothing - get_wing_pos::Function = (_) -> nothing - get_winch_force::Function = (_) -> nothing - get_spring_force::Function = (_) -> nothing - get_stabilize::Function = (_) -> nothing - get_pos::Function = (_) -> nothing -end - -""" - mutable struct SymbolicAWEModel{S, V, P} <: AbstractKiteModel - -State of the kite power system, using a quaternion kite model and three steering lines to the ground. 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. -- V: Vector type, e.g. KVec3 -- P: number of tether points of the system, 3segments+3 -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 SymbolicAWEModel <: AbstractKiteModel - "Reference to the settings struct" - set::Settings - "Reference to the point mass system with points, segments, pulleys and tethers" - sys_struct::SystemStructure - serialized_model::SerializedModel - "Reference to the atmospheric model as implemented in the package AtmosphericModels" - am::AtmosphericModel = AtmosphericModel(set) - integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing - "relative start time of the current time interval" - t_0::SimFloat = 0.0 - "Number of solve! calls" - iter::Int64 = 0 - t_vsm::SimFloat = zero(SimFloat) - t_step::SimFloat = zero(SimFloat) - set_tether_length::Vector{SimFloat} = zeros(SimFloat, 3) -end - -function Base.getproperty(sam::SymbolicAWEModel, sym::Symbol) - if hasfield(SymbolicAWEModel, sym) - getfield(sam, sym) - else - getproperty(getfield(sam, :serialized_model), sym) - end -end -function Base.setproperty!(sam::SymbolicAWEModel, sym::Symbol, val) - if hasfield(SymbolicAWEModel, sym) - setfield!(sam, sym, val) - else - serialized_model = getfield(sam, :serialized_model) - setproperty!(serialized_model, sym, val) - end -end - -""" - SymbolicAWEModel(set::Settings, sys_struct::SystemStructure, - vsm_aeros::Vector{<:BodyAerodynamics}=BodyAerodynamics[], - vsm_solvers::Vector{<:VortexStepMethod.Solver}=VortexStepMethod.Solver[]) - -Constructs a SymbolicAWEModel that can generate ModelingToolkit equations -from the discrete mass-spring-damper representation defined in the [`SystemStructure`](@ref). -The aerodynamic models provide forces and moments acting on wing components. - -# Arguments -- `set::Settings`: Configuration parameters (see [KiteUtils.Settings](https://ufechner7.github.io/KiteUtils.jl/stable/types/#KiteUtils.Settings)) -- `sys_struct::SystemStructure`: Physical system definition with points, segments, groups, etc. -- `vsm_aeros::Vector{<:BodyAerodynamics}=BodyAerodynamics[]`: Aerodynamic models for each wing -- `vsm_solvers::Vector{<:VortexStepMethod.Solver}=VortexStepMethod.Solver[]`: VSM solvers for aerodynamic calculations - -# Returns -- `SymbolicAWEModel`: Model ready for symbolic equation generation via [`init!`](@ref) - -# Example -```julia -# Create wing geometry and aerodynamics -set = se() -wing = RamAirWing(set) -aero = BodyAerodynamics([wing]) -solver = Solver(aero; solver_type=NONLIN) - -# Create system structure -sys_struct = SystemStructure(set, wing) - -# Create symbolic model -model = SymbolicAWEModel(set, sys_struct, [aero], [solver]) -``` -""" -function SymbolicAWEModel( - set::Settings, - sys_struct::SystemStructure, - vsm_aeros::Vector{<:BodyAerodynamics}=BodyAerodynamics[], - vsm_solvers::Vector{<:VortexStepMethod.Solver}=VortexStepMethod.Solver[] -) - vsm_wings = [aero.wings[1] for aero in vsm_aeros] - set_hash = get_set_hash(set) - sys_struct_hash = get_sys_struct_hash(sys_struct) - serialized_model = SerializedModel(; set_hash, sys_struct_hash, vsm_wings, vsm_aeros, vsm_solvers) - return SymbolicAWEModel(; set, sys_struct, serialized_model) -end - -""" - SymbolicAWEModel(set::Settings) - -Constructs a default SymbolicAWEModel with automatically generated components. - -This convenience constructor creates a complete AWE model using default configurations: -- Generates a ram-air wing from settings -- Creates aerodynamic model and VSM solver -- Builds system structure based on the wing geometry -- Assembles everything into a ready-to-use symbolic model - -# Arguments -- `set::Settings`: Configuration parameters (see [KiteUtils.Settings](https://ufechner7.github.io/KiteUtils.jl/stable/types/#KiteUtils.Settings)) - -# Returns -- `SymbolicAWEModel`: Model ready for symbolic equation generation via [`init!`](@ref) - -# Example -```julia -set = se() # Load default settings -model = SymbolicAWEModel(set) - -# Initialize and run simulation -init!(model) -for i in 1:1000 - next_step!(model) -end -``` -""" -function SymbolicAWEModel(set::Settings) - wing = RamAirWing(set; prn=false) - aero = BodyAerodynamics([wing]) - vsm_solver = Solver(aero; solver_type=NONLIN, atol=2e-8, rtol=2e-8) - sys_struct = SystemStructure(set, wing) - return SymbolicAWEModel(set, sys_struct, [aero], [vsm_solver]) -end - -function update_sys_state!(ss::SysState, s::SymbolicAWEModel, zoom=1.0) - isnothing(s.integrator) && error("run init!(s) first") - ss.time = s.integrator.t # Use integrator time - - # Get the state vectors from the integrator - set_values, tether_length, tether_vel, winch_force = s.get_winch_state(s.integrator) - Q_b_w, elevation, azimuth, course, heading, twist_angle, wing_vel, aero_force_b, - aero_moment_b, turn_rate, va_wing_b, wind_vel_wing = s.get_wing_state(s.integrator) - pos, acc, wind_vec_gnd = s.get_point_state(s.integrator) - - P = length(s.sys_struct.points) - for i in 1:P - ss.X[i] = pos[1, i] * zoom - ss.Y[i] = pos[2, i] * zoom - ss.Z[i] = pos[3, i] * zoom - end - - # --- Populate SysState fields --- - ss.acc = norm(acc) # Use the norm of the wing's acceleration vector - ss.orient .= Q_b_w[1, :] - ss.turn_rates .= turn_rate[1, :] - ss.elevation = elevation[1] - ss.azimuth = azimuth[1] - - # Handle potential size mismatch for tether/winch related arrays - num_winches = length(s.sys_struct.winches) - ss.l_tether[1:num_winches] .= tether_length - ss.v_reelout[1:num_winches] .= tether_vel - ss.force[1:num_winches] .= winch_force - - # Depower and Steering from twist angles - num_groups = length(s.sys_struct.wings[1].group_idxs) - ss.twist_angles[1:num_groups] .= twist_angle[1:num_groups] - ss.depower = rad2deg(mean(ss.twist_angles)) # Average twist for depower - ss.steering = rad2deg(ss.twist_angles[num_groups] - ss.twist_angles[1]) - ss.heading = heading[1] # Use heading from MTK model - ss.course = course[1] - # Apparent Wind and Aerodynamics - ss.v_app = norm(va_wing_b[1, :]) - ss.v_wind_gnd .= wind_vec_gnd - ss.v_wind_kite .= wind_vel_wing[1, :] - # Calculate AoA and Side Slip from apparent wind in body frame - # AoA: angle between v_app projected onto xz-plane and x-axis - # Side Slip: angle between v_app and the xz-plane - if ss.v_app > 1e-6 # Avoid division by zero - ss.AoA = atan(va_wing_b[1, 3], va_wing_b[1, 1]) - ss.side_slip = asin(va_wing_b[1, 2] / ss.v_app) - else - ss.AoA = 0.0 - ss.side_slip = 0.0 - end - ss.aero_force_b .= aero_force_b[1, :] - ss.aero_moment_b .= aero_moment_b[1, :] - ss.vel_kite .= wing_vel[1, :] - # Calculate Roll, Pitch, Yaw from Quaternion - q = Q_b_w[1, :] - # roll (x-axis rotation) - sinr_cosp = 2 * (q[1] * q[2] + q[3] * q[4]) - cosr_cosp = 1 - 2 * (q[2] * q[2] + q[3] * q[3]) - ss.roll = atan(sinr_cosp, cosr_cosp) - # pitch (y-axis rotation) - sinp = 2 * (q[1] * q[3] - q[4] * q[2]) - if abs(sinp) >= 1 - ss.pitch = copysign(pi / 2, sinp) # use 90 degrees if out of range - else - ss.pitch = asin(sinp) - end - # yaw (z-axis rotation) - siny_cosp = 2 * (q[1] * q[4] + q[2] * q[3]) - cosy_cosp = 1 - 2 * (q[3] * q[3] + q[4] * q[4]) - ss.yaw = atan(siny_cosp, cosy_cosp) - ss.set_torque[1:3] .= set_values - nothing -end - -function SysState(s::SymbolicAWEModel, zoom=1.0) - ss = SysState{length(s.sys_struct.points)}() - update_sys_state!(ss, s, zoom) - ss -end - -""" - init!(s::SymbolicAWEModel; solver=nothing, adaptive=true, prn=true, - precompile=false, remake=false, reload=false, - lin_outputs=Num[]) -> OrdinaryDiffEqCore.ODEIntegrator - -Initialize a kite power system model. - -If a serialized model exists for the current configuration, it will load that model -and only update the state variables. Otherwise, it will create a new model from scratch. - -# Fast path (serialized model exists): -1. Loads existing ODEProblem from disk -2. Calls `reinit!` to update state variables -3. Sets up integrator with initial settings - -# Slow path (no serialized model): -1. Creates symbolic MTK system with all equations -2. Simplifies system equations -3. Creates ODEProblem and serializes to disk -4. Proceeds with fast path - -# Arguments -- `s::SymbolicAWEModel`: The kite system state object - -# Keyword arguments -- `solver`: Solver algorithm to use. If `nothing`, defaults to `FBDF()` or `QNDF()` based on `s.set.solver`. -- `adaptive::Bool=true`: Whether to use adaptive time stepping. -- `stiffness_factor=nothing`: ignored, for backwards compatibility -- `delta=nothing`: ignored, for backwards compatibility -- `prn::Bool=true`: Whether to print progress information. -- `precompile::Bool=false`: Whether to build problem for precompilation. -- `remake::Bool=false`: If true, forces the system to be rebuilt, even if a serialized model exists. -- `reload::Bool=false`: If true, forces the system to reload the serialized model from disk. -- `lin_outputs::Vector{Num}=Num[]`: List of symbolic variables for which to generate a linearization function. - -# Returns -- `integrator::OrdinaryDiffEqCore.ODEIntegrator`: The initialized ODE integrator. -""" -function init!(s::SymbolicAWEModel; - solver=nothing, stiffness_factor = nothing, delta = nothing, adaptive=true, prn=true, - precompile=false, remake=false, reload=false, - lin_outputs=Num[] -) - if isnothing(solver) - solver = if s.set.solver == "FBDF" - if s.set.quasi_static - FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=s.set.relaxation)) - else - FBDF() - end - elseif s.set.solver == "QNDF" - @warn "This solver is not tested." - QNDF() - else - error("Unavailable solver for SymbolicAWEModel: $(s.set.solver).") - end - end - function init(s) - init_Q_b_w, R_b_w, init_va_b = initial_orient(s) - init!(s.sys_struct, s.set) - - inputs = create_sys!(s, s.sys_struct; init_va_b) - prn && @info "Simplifying the system" - prn ? (@time (sys, _) = structural_simplify(s.full_sys, (inputs, []))) : - ((sys, _) = structural_simplify(s.full_sys, (inputs, []))) - s.sys = sys - dt = SimFloat(1/s.set.sample_freq) - if prn - @info "Creating ODEProblem" - @time s.prob = ODEProblem(s.sys, s.defaults, (0.0, dt); s.guesses) - else - s.prob = ODEProblem(s.sys, s.defaults, (0.0, dt); s.guesses) - end - if length(lin_outputs) > 0 - lin_fun, _ = linearization_function(s.full_sys, [inputs...], lin_outputs; op=s.defaults, guesses=s.guesses) - s.lin_prob = LinearizationProblem(lin_fun, 0.0) - end - sym_vec = get_unknowns(s.sys_struct, s.sys) - s.unknowns_vec = zeros(SimFloat, length(sym_vec)) - generate_getters!(s, sym_vec) - s.set_hash = get_set_hash(s.set) - s.sys_struct_hash = get_sys_struct_hash(s.sys_struct) - serialize(model_path, s.serialized_model) - s.integrator = nothing - return nothing - end - model_path = joinpath(KiteUtils.get_data_path(), get_model_name(s.set; precompile)) - if !ispath(model_path) || remake - init(s) - end - _, success = reinit!(s, solver; adaptive, precompile, reload, lin_outputs, prn) - if !success - rm(model_path) - @info "Rebuilding the system. This can take some minutes..." - init(s) - reinit!(s, solver; adaptive, precompile, lin_outputs, prn, reload=true) - end - return s.integrator -end - -function linearize(s::SymbolicAWEModel; set_values=s.get_set_values(s.integrator)) - isnothing(s.lin_prob) && error("Run init! with remake=true and lin_outputs=...") - s.set_lin_vsm(s.lin_prob, s.get_vsm(s.integrator)) - s.set_lin_set_values(s.lin_prob, set_values) - s.set_lin_unknowns(s.lin_prob, s.get_unknowns(s.integrator)) - return solve(s.lin_prob) -end - -""" - reinit!(s::SymbolicAWEModel, solver; prn=true, precompile=false) -> Nothing - -Reinitialize an existing kite power system model with new state values. -The new state is coming from the init section of the settings, stored -in the struct `s.set`. - -This function performs the following operations: -1. If no integrator exists yet: - - Loads a serialized ODEProblem from disk - - Initializes a new ODE integrator - - Generates getter/setter functions for the system -2. Converts initial settings to quaternion orientation -3. Initializes the point mass system with new positions -4. Sets initial values for all state variables -5. Reinitializes the ODE integrator -6. Updates the linearized aerodynamic model - -This is more efficient than `init!` as it reuses the existing model structure -and only updates the state variables to match the current initial settings. - -# Arguments -- `s::SymbolicAWEModel`: The kite power system state object -- `solver`: The solver to be used -- `prn::Bool=true`: Whether to print progress information - -# Returns -- `Nothing` - -# Throws -- `ArgumentError`: If no serialized problem exists (run `init!` first) -""" -function reinit!( - s::SymbolicAWEModel, - solver; - adaptive=true, - prn=true, - reload=true, - precompile=false, - lin_outputs=Num[] -) - isnothing(s.sys_struct) && error("SystemStructure not defined") - - # init_Q_b_w, R_b_w, init_va_b = initial_orient(s) - init!(s.sys_struct, s.set) - - if isnothing(s.prob) || reload - model_path = joinpath(KiteUtils.get_data_path(), get_model_name(s.set; precompile)) - !ispath(model_path) && error("$model_path not found. Run init!(s::SymbolicAWEModel) first.") - try - s.serialized_model = deserialize(model_path) - catch e - @warn "Failure to deserialize $model_path !" - return s.integrator, false - end - if length(lin_outputs) > 0 && isnothing(s.lin_prob) - @warn "lin_prob is nothing." - return s.integrator, false - elseif (get_set_hash(s.set) != s.serialized_model.set_hash) - @warn "The Settings have changed." - return s.integrator, false - elseif (get_sys_struct_hash(s.sys_struct) != s.serialized_model.sys_struct_hash) - @warn "The SystemStructure has changed." - return s.integrator, false - end - end - if isnothing(s.integrator) || !successful_retcode(s.integrator.sol) || reload - t = @elapsed begin - dt = SimFloat(1/s.set.sample_freq) - s.sys = s.prob.f.sys - s.integrator = OrdinaryDiffEqCore.init(s.prob, solver; - adaptive, dt, abstol=s.set.abs_tol, reltol=s.set.rel_tol, save_on=false, save_everystep=false) - !s.set.quasi_static && (length(s.unknowns_vec) != length(s.integrator.u)) && - error("sam.integrator unknowns of length $(length(s.integrator.u)) should equal sam.unknowns_vec of length $(length(s.unknowns_vec)). - Maybe you forgot to run init!(model; remake=true)?") - end - prn && @info "Initialized integrator in $t seconds" - end - - init_unknowns_vec!(s, s.sys_struct, s.unknowns_vec) - s.set_unknowns(s.integrator, s.unknowns_vec) - s.set_psys(s.integrator, s.sys_struct) - OrdinaryDiffEqCore.reinit!(s.integrator, s.integrator.u; reinit_dae=true) - linearize_vsm!(s) - return s.integrator, true -end - -function generate_getters!(s, sym_vec) - sys = s.sys - c = collect - @unpack wings, winches, tethers = s.sys_struct - - if length(wings) > 0 - vsm_sym = c.([ - sys.last_x, - sys.last_y, - sys.vsm_jac, - ]) - get_vsm = getp(sys, vsm_sym) - s.get_vsm = (integ) -> get_vsm(integ) - get_y = getu(sys, sys.y) - s.get_y = (integ) -> get_y(integ) - get_wing_state = getu(sys, c.([ - sys.Q_b_w, # Orientation quaternion - sys.elevation, # Elevation angle - sys.azimuth, # Azimuth angle - sys.course, # Course angle - sys.heading, # Heading angle (based on body x-axis projection) - sys.twist_angle, # Twist angle per group - sys.wing_vel, # Kite center velocity vector (world frame) - sys.aero_force_b, # Aerodynamic force (body frame) - sys.aero_moment_b, # Aerodynamic moment (body frame) - sys.turn_rate, # Angular velocity (body frame) - sys.va_wing_b, # Apparent wind velocity (body frame) - sys.wind_vel_wing # Wind vector at wing height (world frame) - ])) - s.get_wing_state = (integ) -> get_wing_state(integ) - get_wing_pos = getu(sys, sys.wing_pos) - s.get_wing_pos = (integ) -> get_wing_pos(integ) - - set_vsm = setp(sys, vsm_sym) - s.set_vsm = (integ, val) -> set_vsm(integ, val) - if !isnothing(s.lin_prob) - set_lin_vsm = setp(s.lin_prob, vsm_sym) - s.set_lin_vsm = (lin_prob, val) -> set_lin_vsm(lin_prob, val) - end - end - - if length(winches) > 0 - get_winch_state = getu(sys, c.([ - sys.set_values, - sys.tether_length, # Unstretched length per winch - sys.tether_vel, # Reeling velocity per winch - sys.winch_force, # Force at winch connection point per winch - ])) - s.get_winch_state = (integ) -> get_winch_state(integ) - get_set_values = getp(sys, sys.set_values) - s.get_set_values = (integ) -> get_set_values(integ) - get_tether_length = getu(sys, sys.tether_length) - s.get_tether_length = (integ) -> get_tether_length(integ) - get_winch_force = getu(sys, sys.winch_force) - s.get_winch_force = (integ) -> get_winch_force(integ) - - set_set_values = setp(sys, sys.set_values) - s.set_set_values = (integ, val) -> set_set_values(integ, val) - if !isnothing(s.lin_prob) - set_lin_set_values = setp(s.lin_prob, sys.set_values) - s.set_lin_set_values = (lin_prob, val) -> set_lin_set_values(lin_prob, val) - end - end - - if length(tethers) > 0 - get_unstretched_length = getu(sys, sys.tether_length) - s.get_unstretched_length = (integ) -> get_unstretched_length(integ) - end - - if length(winches) + length(wings) > 0 - set_stabilize = setp(sys, sys.stabilize) - s.set_stabilize = (integ, val) -> set_stabilize(integ, val) - - get_stabilize = getp(sys, sys.stabilize) - s.get_stabilize = (integ) -> get_stabilize(integ) - end - - set_psys = setp(sys, sys.psys) - s.set_psys = (integ, val) -> set_psys(integ, val) - set_wind = setp(sys, [sys.wind_scale_gnd, sys.upwind_dir]) - s.set_wind = (integ, val) -> set_wind(integ, val) - set_unknowns = setu(sys, sym_vec) - s.set_unknowns = (integ, val) -> set_unknowns(integ, val) - set_nonstiff = setu(sys, get_nonstiff_unknowns(s.sys_struct, s.sys)) - s.set_nonstiff = (integ, val) -> set_nonstiff(integ, val) - - get_unknowns = getu(sys, sym_vec) - s.get_unknowns = (integ) -> get_unknowns(integ) - get_point_state = getu(sys, c.([ - sys.pos, # Particle positions - sys.acc, # Kite center acceleration vector (world frame) - sys.wind_vec_gnd, # Ground wind vector (world frame) - ])) - s.get_point_state = (integ) -> get_point_state(integ) - get_spring_force = getu(sys, sys.spring_force) - s.get_spring_force = (integ) -> get_spring_force(integ) - get_pos = getu(sys, sys.pos) - s.get_pos = (integ) -> get_pos(integ) - - if !isnothing(s.lin_prob) - set_lin_unknowns = setu(s.lin_prob, Initial.(sym_vec)) - s.set_lin_unknowns = (lin_prob, val) -> set_lin_unknowns(lin_prob, val) - end - nothing -end - -function linearize_vsm!(s::SymbolicAWEModel, integ=s.integrator) - @unpack wings, y, x, jac = s.sys_struct - if length(wings) > 0 - y .= s.get_y(integ) - for wing in wings - res = VortexStepMethod.linearize( - s.vsm_solvers[wing.idx], - s.vsm_aeros[wing.idx], - y[wing.idx, :]; - va_idxs=1:3, - omega_idxs=4:6, - theta_idxs=7:6+length(s.sys_struct.groups), - moment_frac=s.sys_struct.groups[1].moment_frac - ) - jac[wing.idx, :, :] .= res[1] - x[wing.idx, :] .= res[2] - end - s.set_vsm(integ, [x, y, jac]) - end - nothing -end - -""" - next_step!(s::SymbolicAWEModel, integrator::ODEIntegrator; set_values=nothing, upwind_dir=nothing, dt=1/s.set.sample_freq, vsm_interval=1) - -Take a simulation step, using the internal integrator. - -This function performs the following steps: -1. Optionally update the set values (control inputs) -2. Optionally update the upwind direction -3. Optionally linearize the VSM (Vortex Step Method) model -4. Step the ODE integrator forward by `dt` seconds -5. Check for a successful return code from the integrator -6. Increment the iteration counter - -# Arguments -- `s::SymbolicAWEModel`: The kite power system state object -- `integrator::ODEIntegrator`: The ODE integrator to use - -# Keyword Arguments -- `set_values=nothing`: New values for the set variables (control inputs). If `nothing`, the current values are used. -- `dt=1/s.set.sample_freq`: Time step size in seconds. Defaults to the inverse of the sample frequency. -- `vsm_interval=1`: Interval (in number of steps) at which to linearize the VSM model. If 0, the VSM model is not linearized. - -# Returns -- `Nothing` -""" -function next_step!(s::SymbolicAWEModel, integrator::OrdinaryDiffEqCore.ODEIntegrator; set_values=nothing, dt=1/s.set.sample_freq, vsm_interval=1) - !(s.integrator === integrator) && error("The ODEIntegrator doesn't belong to the SymbolicAWEModel") - next_step!(s; set_values, upwind_dir, dt, vsm_interval) -end - -function next_step!(s::SymbolicAWEModel; set_values=nothing, dt=1/s.set.sample_freq, vsm_interval=1) - if (!isnothing(set_values)) - s.set_set_values(s.integrator, set_values) - end - if vsm_interval != 0 && s.iter % vsm_interval == 0 - s.t_vsm = @elapsed linearize_vsm!(s) - end - - s.t_0 = s.integrator.t - s.t_step = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) - if !successful_retcode(s.integrator.sol) - @warn "Return code for solution: $(s.integrator.sol.retcode)" - end - @assert successful_retcode(s.integrator.sol) - s.iter += 1 - return nothing -end - -function get_model_name(set::Settings; precompile=false) - suffix = "" - ver = "$(VERSION.major).$(VERSION.minor)" - if precompile - suffix = ".default" - end - dynamics_type = ifelse(set.quasi_static, "static", "dynamic") - return "model_$(ver)_$(set.physical_model)_$(dynamics_type)_$(set.segments)_seg.bin$suffix" -end - -""" -Calculate and return the angle of attack in rad -""" -function calc_aoa(s::SymbolicAWEModel) - alpha_array = s.vsm_solvers[1].sol.alpha_array - middle = length(alpha_array) ÷ 2 - if iseven(length(alpha_array)) - return 0.5alpha_array[middle] + 0.5alpha_array[middle+1] - else - return alpha_array[middle+1] - end -end - -function init_unknowns_vec!( - s::SymbolicAWEModel, - system::SystemStructure, - vec::Vector{SimFloat} -) - !s.set.quasi_static && (length(vec) != length(s.integrator.u)) && - error("Unknowns of length $(length(s.integrator.u)) but vector provided of length $(length(vec))") - - @unpack points, groups, segments, pulleys, winches, wings = system - vec_idx = 1 - - for point in points - if point.type == DYNAMIC - for i in 1:3 - vec[vec_idx] = point.pos_w[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = point.vel_w[i] - vec_idx += 1 - end - end - end - for pulley in pulleys - if pulley.type == DYNAMIC - vec[vec_idx] = pulley.length - vec_idx += 1 - vec[vec_idx] = pulley.vel - vec_idx += 1 - end - end - for group in groups - if group.type == DYNAMIC - vec[vec_idx] = group.twist - vec_idx += 1 - vec[vec_idx] = group.twist_vel - vec_idx += 1 - end - end - for winch in winches - vec[vec_idx] = winch.tether_length - vec_idx += 1 - vec[vec_idx] = winch.tether_vel - vec_idx += 1 - end - - for wing in wings - for i in 1:4 - vec[vec_idx] = wing.orient[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = wing.angular_vel[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = wing.pos_w[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = wing.vel_w[i] - vec_idx += 1 - end - end - - (vec_idx-1 != length(vec)) && - error("Unknowns vec is of length $(length(vec)) but the last index is $(vec_idx-1)") - nothing -end - -function get_unknowns(sys_struct::SystemStructure, sys::ODESystem) - vec = Num[] - @unpack points, groups, segments, pulleys, winches, wings = sys_struct - for point in points - for i in 1:3 - point.type == DYNAMIC && push!(vec, sys.pos[i, point.idx]) - end - for i in 1:3 - point.type == DYNAMIC && push!(vec, sys.vel[i, point.idx]) - end - end - for pulley in pulleys - pulley.type == DYNAMIC && push!(vec, sys.pulley_l0[pulley.idx]) - pulley.type == DYNAMIC && push!(vec, sys.pulley_vel[pulley.idx]) - end - vec = get_nonstiff_unknowns(sys_struct, sys, vec) - return vec -end - -function get_nonstiff_unknowns(sys_struct::SystemStructure, sys::ODESystem, vec=Num[]) - @unpack points, groups, segments, pulleys, winches, wings = sys_struct - for group in groups - group.type == DYNAMIC && push!(vec, sys.free_twist_angle[group.idx]) - group.type == DYNAMIC && push!(vec, sys.twist_ω[group.idx]) - end - for winch in winches - push!(vec, sys.tether_length[winch.idx]) - push!(vec, sys.tether_vel[winch.idx]) - end - for wing in wings - [push!(vec, sys.Q_b_w[wing.idx, i]) for i in 1:4] - [push!(vec, sys.ω_b[wing.idx, i]) for i in 1:3] - [push!(vec, sys.wing_pos[wing.idx, i]) for i in 1:3] - [push!(vec, sys.wing_vel[wing.idx, i]) for i in 1:3] - end - return vec -end - -function find_steady_state!(s::SymbolicAWEModel; dt=1/s.set.sample_freq) - old_state = s.get_stabilize(s.integrator) - s.set_stabilize(s.integrator, true) - for _ in 1:1÷dt - next_step!(s; dt, vsm_interval=1) - end - s.set_stabilize(s.integrator, old_state) - return nothing -end - -function initial_orient(s::SymbolicAWEModel) - set = s.set - wings = s.sys_struct.wings - R_b_w = zeros(length(wings), 3, 3) - Q_b_w = zeros(length(wings), 4) - init_va_b = zeros(length(wings), 3) - for wing in wings - R_cad_body = s.vsm_wings[wing.idx].R_cad_body - x = [0, 0, -1] # laying flat along x axis - z = [1, 0, 0] # laying flat along x axis - x = rotate_around_y(x, -deg2rad(set.elevation)) - z = rotate_around_y(z, -deg2rad(set.elevation)) - x = rotate_around_z(x, deg2rad(set.azimuth)) - z = rotate_around_z(z, deg2rad(set.azimuth)) - R_b_w[wing.idx, :, :] .= R_cad_body' * hcat(x, z × x, z) - Q_b_w[wing.idx, :] .= rotation_matrix_to_quaternion(R_b_w[wing.idx, :, :]) - init_va_b[wing.idx, :] .= R_b_w[wing.idx, :, :]' * [s.set.v_wind, 0., 0.] - end - return Q_b_w, R_b_w, init_va_b -end - -unstretched_length(s::SymbolicAWEModel) = s.get_unstretched_length(s.integrator) -tether_length(s::SymbolicAWEModel) = s.get_tether_length(s.integrator) -calc_height(s::SymbolicAWEModel) = s.get_wing_pos(s.integrator)[3] -winch_force(s::SymbolicAWEModel) = s.get_winch_force(s.integrator) -spring_forces(s::SymbolicAWEModel) = s.get_spring_force(s.integrator) -function pos(s::SymbolicAWEModel) - pos = s.get_pos(s.integrator) - return [pos[:,i] for i in eachindex(pos[1,:])] -end - -function min_chord_length(s::SymbolicAWEModel) - min_len = Inf - for wing in s.vsm_wings - le_pos = [wing.le_interp[i](wing.gamma_tip) for i in 1:3] - te_pos = [wing.te_interp[i](wing.gamma_tip) for i in 1:3] - min_len = min(norm(le_pos - te_pos), min_len) - end - return min_len -end - -""" - set_depower_steering!(s::SymbolicAWEModel, depower, steering) -> Nothing - -Set kite depower and steering by adjusting tether lengths. Depower controls angle of attack, -steering controls left/right differential. Values are scaled by minimum chord length. -""" -function set_depower_steering!(s::SymbolicAWEModel, depower, steering) - len = s.set_tether_length - len .= tether_length(s) - depower *= min_chord_length(s) - steering *= min_chord_length(s) - len[2] = 0.5 * (2*depower + 2*len[1] + steering) - len[3] = 0.5 * (2*depower + 2*len[1] - steering) - return nothing -end - -""" - set_v_wind_ground!(s::SymbolicAWEModel, v_wind_gnd=s.set.v_wind, upwind_dir=-π/2) -> Nothing - -Set ground wind speed (m/s) and upwind direction (radians). Direction: 0=north, π/2=east, -π=zouth, -π/2=west (default). -""" -function set_v_wind_ground!(s::SymbolicAWEModel, v_wind_gnd=s.set.v_wind, upwind_dir=-pi/2) - s.set_wind(s.integrator, [v_wind_gnd, upwind_dir]) - return nothing -end - -function get_set_hash(set::Settings; - fields=[:segments, :model, :foil_file, :physical_model, :quasi_static, :winch_model] - ) - h = zeros(UInt8, 1) - for field in fields - value = getfield(set, field) - h = sha1(string((value, h))) - end - return h -end - -function get_sys_struct_hash(sys_struct::SystemStructure) - @unpack points, groups, segments, pulleys, tethers, winches, wings, transforms = sys_struct - data_parts = [] - for point in points - push!(data_parts, ("point", point.idx, point.wing_idx, Int(point.type))) - end - for segment in segments - push!(data_parts, ("segment", segment.idx, segment.point_idxs, Int(segment.type))) - end - for group in groups - push!(data_parts, ("group", group.idx, group.point_idxs, Int(group.type))) - end - for pulley in pulleys - push!(data_parts, ("pulley", pulley.idx, pulley.segment_idxs, Int(pulley.type))) - end - for tether in tethers - push!(data_parts, ("tether", tether.idx, tether.segment_idxs)) - end - for winch in winches - model_type = winch.model isa TorqueControlledMachine - push!(data_parts, ("winch", winch.idx, model_type, winch.tether_idxs)) - end - for wing in wings - push!(data_parts, ("wing", wing.idx, wing.group_idxs)) - end - for transform in transforms - push!(data_parts, ("transform", transform.idx, transform.wing_idx, transform.rot_point_idx, - transform.base_point_idx, transform.base_transform_idx)) - end - content = string(data_parts) - return sha1(content) -end diff --git a/src/system_structure.jl b/src/system_structure.jl deleted file mode 100644 index df855531..00000000 --- a/src/system_structure.jl +++ /dev/null @@ -1,1093 +0,0 @@ - -# Copyright (c) 2025 Bart van de Lint -# SPDX-License-Identifier: MPL-2.0 - -function VortexStepMethod.RamAirWing(set::Settings; prn=true, kwargs...) - obj_path = joinpath(dirname(get_data_path()), set.model) - dat_path = joinpath(dirname(get_data_path()), set.foil_file) - return RamAirWing(obj_path, dat_path; - mass=set.mass, crease_frac=set.crease_frac, n_groups=4, - align_to_principal=true, prn, kwargs... - ) -end - -""" - SegmentType `POWER_LINE` `STEERING_LINE` `BRIDLE` - -Type of segment. - -# Elements -- POWER_LINE: Belongs to a power line -- STEERING_LINE: Belongs to a steering line -- BRIDLE: Belongs to the bridle -""" -@enum SegmentType begin - POWER_LINE - STEERING_LINE - BRIDLE -end - -""" - DynamicsType `DYNAMIC` `QUASI_STATIC` `WING` `STATIC` - -Enumeration of the models that are attached to a point. - -# Elements -- DYNAMIC: Belongs to a dynamic tether model -- QUASI_STATIC: Belongs to a quasi static tether model -- WING: Connected to the rigid wing body -- STATIC: Does not change position -""" -@enum DynamicsType begin - DYNAMIC - QUASI_STATIC - WING - STATIC -end - -""" - mutable struct Point - -A point mass. - -$(TYPEDFIELDS) -""" -struct Point - idx::Int16 - transform_idx::Int16 # idx of wing used for initial orientation - wing_idx::Int16 - pos_cad::KVec3 - pos_b::KVec3 # pos relative to wing COM in body frame - pos_w::KVec3 # pos in world frame - vel_w::KVec3 # vel in world frame - type::DynamicsType -end - -""" - Point(idx, pos_cad, type; wing_idx=1, vel_w=zeros(KVec3), transform_idx=1) - -Constructs a Point object. A point can be of four different [`DynamicsType`](@ref)s: -- `STATIC`: the point doesn't move. ``\\ddot{\\mathbf{r}} = \\mathbf{0}`` -- `DYNAMIC`: the point moves according to Newton's second law. ``\\ddot{\\mathbf{r}} = \\mathbf{F}/m`` -- `QUASI_STATIC`: the acceleration is constrained to be zero, by solving a nonlinear problem. ``\\mathbf{F}/m = \\mathbf{0}`` -- `WING`: the point has a static position in the rigid body wing frame. ``\\mathbf{r}_w = \\mathbf{r}_{wing} + \\mathbf{R}_{b\\rightarrow w} \\mathbf{r}_b`` - -where: -- ``\\mathbf{r}`` is the point position vector -- ``\\mathbf{F}`` is the net force acting on the point -- ``m`` is the point mass -- ``\\mathbf{r}_w`` is the position in world frame -- ``\\mathbf{r}_{wing}`` is the wing center position -- ``\\mathbf{R}_{b\\rightarrow w}`` is the rotation matrix from body to world frame -- ``\\mathbf{r}_b`` is the position in body frame - -# Arguments -- `idx::Int16`: Unique identifier for the point. -- `pos_cad::KVec3`: Position of the point in the CAD frame. -- `type::DynamicsType`: Dynamics type of the point (STATIC, DYNAMIC, etc.). - -# Keyword Arguments -- `wing_idx::Int16=1`: Index of the wing this point is attached to. -- `vel_w::KVec3=zeros(KVec3)`: Initial velocity of the point in world frame. -- `transform_idx::Int16=1`: Index of the transform used for initial positioning. - -# Returns -- `Point`: A new Point object. - -# Example -To create a Point: -```julia - point = Point(1, [1.0, 2.0, 3.0], DYNAMIC; wing_idx=1) -``` -""" -function Point(idx, pos_cad, type; wing_idx=1, vel_w=zeros(KVec3), transform_idx=1) - Point(idx, transform_idx, wing_idx, pos_cad, zeros(KVec3), zeros(KVec3), vel_w, type) -end - -""" - struct Group - -Set of bridle lines that share the same twist angle and trailing edge angle. - -$(TYPEDFIELDS) -""" -mutable struct Group - const idx::Int16 - const point_idxs::Vector{Int16} - const le_pos::KVec3 # point which the group rotates around under wing deformation - const chord::KVec3 # chord vector in body frame which the group rotates around under wing deformation - const y_airf::KVec3 # spanwise vector in local panel frame which the group rotates around under wing deformation - const type::DynamicsType - moment_frac::SimFloat - twist::SimFloat - twist_vel::SimFloat -end - -""" - Group(idx, point_idxs, vsm_wing::RamAirWing, gamma, type, moment_frac) - -Constructs a Group object representing a collection of points on a kite body that share -a common twist deformation. - -A Group models the local deformation of a kite wing section through twist dynamics. -All points within a group undergo the same twist rotation about the chord vector. - -The governing equation is: -```math -\\begin{aligned} -\\tau = \\underbrace{\\sum_{i=1}^{4} r_{b,i} \\times (\\mathbf{F}_{b,i} \\cdot \\hat{\\mathbf{z}})}_{\\text{bridles}} + \\underbrace{r_a \\times (\\mathbf{F}_a \\cdot \\hat{\\mathbf{z}})}_{\\text{aero}} -\\end{aligned} -``` - -![System Overview](group-slice.svg) - -where: -- ``\\tau`` is the total torque about the twist axis -- ``r_{b,i}`` is the position vector of bridle point ``i`` relative to the twist center -- ``\\mathbf{F}_{b,i}`` is the force at bridle point ``i`` -- ``\\hat{\\mathbf{z}}`` is the unit vector along the twist axis (chord direction) -- ``r_a`` is the position vector of the aerodynamic center relative to the twist center -- ``\\mathbf{F}_a`` is the aerodynamic force at the group's aerodynamic center - -The group can have two [`DynamicsType`](@ref)s: -- `DYNAMIC`: the group rotates according to Newton's second law: ``I\\ddot{\\theta} = \\tau`` -- `QUASI_STATIC`: the rotational acceleration is zero: ``\\tau = 0`` - -# Arguments -- `idx::Int16`: Unique identifier for the group -- `point_idxs::Vector{Int16}`: Indices of points that move together with this group's twist -- `vsm_wing::RamAirWing`: Wing geometry object used to extract local chord and spanwise vectors -- `gamma`: Spanwise parameter (typically -1 to 1) defining the group's location along the wing -- `type::DynamicsType`: Dynamics type (DYNAMIC for time-varying twist, QUASI_STATIC for equilibrium) -- `moment_frac::SimFloat`: Chordwise position (0=leading edge, 1=trailing edge) about which the group rotates - -# Returns -- `Group`: A new Group object with twist dynamics capability - -# Example -Create a group at mid-span with quarter of the wing moment: -```julia - group = Group(1, [1, 2, 3], vsm_wing, 0.0, DYNAMIC, 0.25) -``` -""" -function Group(idx, point_idxs, vsm_wing::RamAirWing, gamma, type, moment_frac) - le_pos = [vsm_wing.le_interp[i](gamma) for i in 1:3] - chord = [vsm_wing.te_interp[i](gamma) for i in 1:3] .- le_pos - y_airf = normalize([vsm_wing.le_interp[i](gamma-0.01) for i in 1:3] - le_pos) - Group(idx, point_idxs, le_pos, chord, y_airf, type, moment_frac, 0.0, 0.0) -end - -""" - Group(idx, point_idxs, le_pos, chord, y_airf, type, moment_frac) - -Constructs a Group object representing a collection of points on a kite body that share -a common twist deformation. See: [`Group(::Any, ::Any, ::RamAirWing, ::Any, ::Any, ::Any)`](@ref). - -""" -function Group(idx, point_idxs, le_pos, chord, y_airf, type, moment_frac) - Group(idx, point_idxs, le_pos, chord, y_airf, type, moment_frac, 0.0, 0.0) -end - -""" - mutable struct Segment - -A segment from one point index to another point index. - -$(TYPEDFIELDS) -""" -mutable struct Segment - const idx::Int16 - const point_idxs::Tuple{Int16, Int16} - const type::SegmentType - l0::SimFloat - compression_frac::SimFloat - diameter::SimFloat -end - -""" - Segment(idx, point_idxs, type; l0=zero(SimFloat), compression_frac=0.1) - -Constructs a Segment object representing an elastic spring-damper connection between two points. - -The segment follows Hooke's law with damping and aerodynamic drag: - -**Spring-Damper Force:** -```math -\\mathbf{F}_{spring} = \\left[k(l - l_0) - c\\dot{l}\\right]\\hat{\\mathbf{u}} -``` - -**Aerodynamic Drag:** -```math -\\mathbf{F}_{drag} = \\frac{1}{2}\\rho C_d A |\\mathbf{v}_a| \\mathbf{v}_{a,\\perp} -``` - -**Total Force:** -```math -\\mathbf{F}_{total} = \\mathbf{F}_{spring} + \\mathbf{F}_{drag} -``` - -where: -- ``k = \\frac{E \\pi d^2/4}{l}`` is the axial stiffness -- ``l`` is current length, ``l_0`` is unstretched length -- ``c = \\frac{\\xi}{c_{spring}} k`` is damping coefficient -- ``\\hat{\\mathbf{u}} = \\frac{\\mathbf{r}_2 - \\mathbf{r}_1}{l}`` is unit vector along segment -- ``\\dot{l} = (\\mathbf{v}_1 - \\mathbf{v}_2) \\cdot \\hat{\\mathbf{u}}`` is extension rate -- ``\\mathbf{v}_{a,\\perp}`` is apparent wind velocity perpendicular to segment - -# Arguments -- `idx::Int16`: Unique identifier for the segment. -- `point_idxs::Tuple{Int16, Int16}`: Tuple containing the indices of the two points connected by this segment. -- `type::SegmentType`: Type of the segment (POWER_LINE, STEERING_LINE, BRIDLE). - -# Keyword Arguments -- `l0::SimFloat=zero(SimFloat)`: Unstretched length of the segment. Calculated from point positions if zero. -- `compression_frac::SimFloat=0.1`: Compression fraction of stiffness for compression behavior. - -# Returns -- `Segment`: A new Segment object. - -# Example -To create a Segment: -```julia - segment = Segment(1, (1, 2), BRIDLE; l0=10.0) -``` -""" -function Segment(idx, point_idxs, type; l0=zero(SimFloat), compression_frac=0.1) - Segment(idx, point_idxs, type, l0, compression_frac, zero(SimFloat)) -end - -""" - mutable struct Pulley - -A pulley described by two segments with the common point of the segments being the pulley. - -$(TYPEDFIELDS) -""" -mutable struct Pulley - const idx::Int16 - const segment_idxs::Tuple{Int16, Int16} - const type::DynamicsType - sum_length::SimFloat - length::SimFloat - vel::SimFloat -end - -""" - Pulley(idx, segment_idxs, type) - -Constructs a Pulley object that enforces length redistribution between two segments. - -The pulley constraint maintains constant total length while allowing force transmission: - -**Constraint Equations:** -```math -l_1 + l_2 = l_{total} = \\text{constant} -``` - -**Force Balance:** -```math -F_{pulley} = F_1 - F_2 -``` - -**Dynamics:** -```math -m\\ddot{l}_1 = F_{pulley} = F_1 - F_2 -``` - -where: -- ``l_1, l_2`` are the lengths of connected segments -- ``F_1, F_2`` are the spring forces in the segments -- ``m = \\rho_{tether} \\pi (d/2)^2 l_{total}`` is the total mass of both segments -- ``\\dot{l}_1 + \\dot{l}_2 = 0`` (velocity constraint) - -The pulley can have two [`DynamicsType`](@ref)s: -- `DYNAMIC`: the length redistribution follows Newton's second law: ``m\\ddot{l}_1 = F_1 - F_2`` -- `QUASI_STATIC`: the forces are balanced instantaneously: ``F_1 = F_2`` - -# Arguments -- `idx::Int16`: Unique identifier for the pulley. -- `segment_idxs::Tuple{Int16, Int16}`: Tuple containing the indices of the two segments connected by this pulley. -- `type::DynamicsType`: Dynamics type of the pulley (DYNAMIC or QUASI_STATIC). - -# Returns -- `Pulley`: A new Pulley object. - -# Example -To create a Pulley: -```julia - pulley = Pulley(1, (1, 2), DYNAMIC) -``` -""" -function Pulley(idx, segment_idxs, type) - return Pulley(idx, segment_idxs, type, 0.0, 0.0, 0.0) -end - -""" - Tether(idx, segment_idxs) - -Constructs a Tether object representing a flexible line composed of multiple segments. - -A tether enforces a shared unstretched length constraint across all its constituent segments: - -**Length Constraint:** -```math -\\sum_{i \\in \\text{segments}} l_{0,i} = L -``` - -**Winch Control:** -The unstretched tether length is controlled by winch acceleration: -```math -\\ddot L = \\alpha(v, F, u) -``` - -where: -- ``L`` is the tether length -- ``l_{0,i}`` is the segment unstretched length -- ``\\alpha(v, F, u)`` is the winch acceleration function depending on model type - -# Arguments -- `idx::Int16`: Unique identifier for the tether -- `segment_idxs::Vector{Int16}`: Indices of segments that form this tether - -# Returns -- `Tether`: A new Tether object - -# Example -Create a tether from segments 1, 2, and 3: -```julia - tether = Tether(1, [1, 2, 3]) -``` -""" -struct Tether - idx::Int16 - segment_idxs::Vector{Int16} -end - -""" - mutable struct Winch - -A set of tethers or just one tether connected to a winch. - -$(TYPEDFIELDS) -""" -mutable struct Winch - const idx::Int16 - const model::AbstractWinchModel - const tether_idxs::Vector{Int16} - tether_length::SimFloat - tether_vel::SimFloat -end - -""" - Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0) - -Constructs a Winch object that controls tether length through torque or speed regulation. - -**Tether Length Control:** -```math -\\ddot{L} = \\alpha(v, F, u) -``` -where: -- ``L`` is the tether length -- ``v`` is the reel out velocity (tether extension rate) -- ``F`` is the tether force -- ``u`` is the applied torque or speed setpoint -- ``\\alpha(v, F, u)`` is the winch acceleration function depending on model type - -where the winch acceleration function `f_winch` depends on the winch model type: -- **Torque-controlled**: Direct torque input with motor dynamics -- **Speed-controlled**: Velocity regulation with internal control loops - -For detailed mathematical models of winch dynamics, motor characteristics, and control algorithms, -see the [WinchModels.jl documentation](https://github.com/aenarete/WinchModels.jl/blob/main/docs/winch.md). - -# Arguments -- `idx::Int16`: Unique identifier for the winch. -- `model::AbstractWinchModel`: The winch model (TorqueControlledMachine, AsyncMachine, etc.). -- `tether_idxs::Vector{Int16}`: Vector containing the indices of the tethers connected to this winch. -- `tether_length::SimFloat`: Initial tether length. - -# Keyword Arguments -- `tether_vel::SimFloat=0.0`: Initial tether velocity (reel-out rate). - -# Returns -- `Winch`: A new Winch object. - -# Example -To create a Winch: -```julia - winch = Winch(1, TorqueControlledMachine(set), [1, 2], 100.0) -``` -""" -function Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0) - return Winch(idx, model, tether_idxs, tether_length, tether_vel) -end - -""" - struct Wing - -A rigid wing body that can have multiple groups of points attached to it. - -# Fields -- `idx::Int16`: Unique identifier for the wing -- `group_idxs::Vector{Int16}`: Indices of groups attached to this wing -- `transform_idx::Int16`: Transform used for initial positioning and orientation -- `R_b_c::Matrix{SimFloat}`: Rotation matrix from body frame to CAD frame -- `angular_vel::KVec3`: Angular velocity of the wing in world frame -- `pos_w::KVec3`: Position of wing center of mass in world frame -- `pos_cad::KVec3`: Position of wing center of mass in CAD frame -- `vel_w::KVec3`: Velocity of wing center of mass in world frame - -The wing provides a rigid body reference frame for attached points and groups. -Points with `type == WING` move rigidly with the wing body according to the -wing's orientation matrix `R_b_c` and position `pos_w`. - -# Extended help -The wing's orientation can be accessed as a quaternion through the `orient` property: -```julia -wing = Wing(1, [1,2], I(3), zeros(3)) -quat = wing.orient # Returns quaternion representation of R_b_c -``` -""" -struct Wing - idx::Int16 - group_idxs::Vector{Int16} - transform_idx::Int16 - R_b_c::Matrix{SimFloat} - angular_vel::KVec3 - pos_w::KVec3 - pos_cad::KVec3 - vel_w::KVec3 -end -function Base.getproperty(wing::Wing, sym::Symbol) - if sym == :orient - return rotation_matrix_to_quaternion(wing.R_b_c) - else - return getfield(wing, sym) - end -end - -""" - Wing(idx, group_idxs, R_b_c, pos_cad; transform_idx=1, angular_vel=zeros(KVec3), - pos_w=zeros(KVec3), vel_w=zeros(KVec3)) - -Constructs a Wing object representing a rigid body that serves as a reference frame for attached points and groups. - -A Wing provides a rigid body coordinate system for kite components. Points with `type == WING` move rigidly -with the wing body according to the wing's orientation matrix and position. Groups attached to the wing -undergo local deformation (twist) relative to the rigid wing body frame. - -**Rigid Body Dynamics:** -The wing follows standard rigid body equations of motion: - -```math -\\begin{aligned} -\\frac{\\delta \\mathbf{q}_b^w}{\\delta t} &= \\frac{1}{2} \\Omega(\\boldsymbol{\\omega}_b) \\mathbf{q}_b^w \\\\ -\\boldsymbol{\\tau}_w &= \\mathbf{I} \\frac{\\delta \\boldsymbol{\\omega}}{\\delta t} + \\boldsymbol{\\omega}_b \\times (\\mathbf{I}\\boldsymbol{\\omega}_b) -\\end{aligned} -``` - -where: -- ``\\mathbf{q}_b^w`` is the quaternion from world to body frame -- ``\\boldsymbol{\\omega}_b`` is the angular velocity in body frame -- ``\\Omega(\\boldsymbol{\\omega}_b)`` is the quaternion multiplication matrix -- ``\\mathbf{I}`` is the inertia tensor in body frame -- ``\\boldsymbol{\\tau}_w`` is the total applied torque to the rigid wing body (aerodynamic + tether forces) - -**Coordinate Transformations:** -Points attached to the wing transform as: -```math -\\mathbf{r}_w = \\mathbf{r}_{w} + \\mathbf{R}_{b \\rightarrow w} \\mathbf{r}_b -``` - -where: -- ``\\mathbf{r}_w`` is the position in world frame -- ``\\mathbf{r}_{w}`` is the wing position in world frame -- ``\\mathbf{R}_{b \\rightarrow w}`` is the rotation from body to world frame -- ``\\mathbf{r}_b`` is the point position in body frame - -# Arguments -- `idx::Int16`: Unique identifier for the wing -- `group_idxs::Vector{Int16}`: Indices of groups attached to this wing that can deform relative to the body -- `R_b_c::Matrix{SimFloat}`: Rotation matrix from body frame to CAD frame (3×3 orthogonal matrix) -- `pos_cad::KVec3`: Position of wing center of mass in CAD frame - -# Keyword Arguments -- `transform_idx::Int16=1`: Transform used for initial positioning and orientation -- `angular_vel::KVec3=zeros(KVec3)`: Initial angular velocity of the wing in world frame -- `pos_w::KVec3=zeros(KVec3)`: Initial position of wing center of mass in world frame -- `vel_w::KVec3=zeros(KVec3)`: Initial velocity of wing center of mass in world frame - -# Special Properties -The wing orientation can be accessed as a quaternion: -```julia - quat = wing.orient # Returns quaternion representation of R_b_c -``` - -# Returns -- `Wing`: A new Wing object providing a rigid body reference frame - -# Example -Create a wing with identity orientation and two attached groups: -```julia - R_b_c = I(3) # identity matrix - pos_cad = [0.0, 0.0, 0.0] - wing = Wing(1, [1, 2], R_b_c, pos_cad) -``` -""" -function Wing(idx, group_idxs, R_b_c, pos_cad; transform_idx=1, angular_vel=zeros(KVec3), - pos_w=zeros(KVec3), vel_w=zeros(KVec3)) - return Wing(idx, group_idxs, transform_idx, R_b_c, angular_vel, pos_w, pos_cad, vel_w) -end - -""" - mutable struct Transform - -Describes the spatial transformation (position and orientation) of system components -relative to a base reference point. - -$(TYPEDFIELDS) -""" -mutable struct Transform - const idx::Int16 - const wing_idx::Union{Int16, Nothing} - const rot_point_idx::Union{Int16, Nothing} - const base_point_idx::Union{Int16, Nothing} - const base_transform_idx::Union{Int16, Nothing} - elevation::SimFloat # The elevation of the rotating point or kite as seen from the base point - azimuth::SimFloat # The azimuth of the rotating point or kite as seen from the base point - heading::SimFloat - base_pos::Union{KVec3, Nothing} -end - -""" - Transform(idx, elevation, azimuth, heading; - base_point_idx=nothing, base_pos=nothing, base_transform_idx=nothing, - wing_idx=nothing, rot_point_idx=nothing) - -Constructs a Transform object that orients system components using spherical coordinates. - -**All points and wings with matching `transform_idx` are transformed together as a rigid body:** -1. **Translation**: Translate such that base is at specified base pos -1. **Rotation 1**: Rotate so target is at (elevation, azimuth) relative to base -2. **Rotation 2**: Rotate all components by `heading` around the base-target vector - -```math -\\mathbf{r}_{transformed} = \\mathbf{r}_{base} + \\mathbf{R}_{heading} \\circ \\mathbf{R}_{elevation,azimuth}(\\mathbf{r} - \\mathbf{r}_{base}) -``` - -# Arguments -- `idx::Int16`: Unique identifier for the transform -- `elevation::SimFloat`: Target elevation angle from base (radians) -- `azimuth::SimFloat`: Target azimuth angle from base (radians) -- `heading::SimFloat`: Rotation around base-target vector (radians) - -# Keyword Arguments -**Base Reference (choose one):** -- `base_pos + base_point_idx`: Fixed position and reference point -- `base_transform_idx`: Chain to another transform's position - -**Target Object (choose one):** -- `wing_idx`: Wing to position at (elevation, azimuth) -- `rot_point_idx`: Point to position at (elevation, azimuth) - -# Returns -- `Transform`: Transform affecting all components with matching `transform_idx` - -# Examples -```julia -# Position wing and all associated points at 45° elevation, 30° azimuth -transform = Transform(1, deg2rad(45), deg2rad(30), 0.0; - base_pos=[0,0,0], base_point_idx=1, wing_idx=1) - -# Chain transforms for multi-kite systems -transform2 = Transform(2, deg2rad(30), deg2rad(45), deg2rad(10); - base_transform_idx=1, wing_idx=2) -``` -""" -function Transform(idx, elevation, azimuth, heading; - base_point_idx=nothing, base_pos=nothing, base_transform_idx=nothing, - wing_idx=nothing, rot_point_idx=nothing) - (isnothing(wing_idx) == isnothing(rot_point_idx)) && error("Either provide a wing_idx or a rot_point_idx, not both or none.") - (isnothing(base_pos) == isnothing(base_transform_idx)) && error("Either provide the base_pos or the base_transform_idx, not both or none.") - (isnothing(base_pos) !== isnothing(base_point_idx)) && error("When providing a base_pos, also provide a base_point_idx.") - Transform(idx, wing_idx, rot_point_idx, base_point_idx, base_transform_idx, elevation, azimuth, heading, base_pos) -end -function Transform(idx, set, base_point_idx; kwargs...) - Transform(idx, set.elevations[idx], set.azimuths[idx], set.headings[idx], base_point_idx; kwargs...) -end - -function get_rot_pos(transform::Transform, wings, points) - if !isnothing(transform.wing_idx) - return wings[transform.wing_idx].pos_w - elseif !isnothing(transform.rot_point_idx) - return points[transform.rot_point_idx].pos_w - end -end - -function get_base_pos(transform::Transform, wings, points) - curr_base_pos = points[transform.base_point_idx].pos_cad - if !isnothing(transform.base_pos) - return transform.base_pos, curr_base_pos - elseif !isnothing(transform.base_transform_idx) - base_transform = transforms[transform.base_transform_idx] - return get_rot_pos(base_transform, wings, points), curr_base_pos - end -end - -""" - struct SystemStructure - -A discrete mass-spring-damper representation of a kite system, where point masses -connected by elastic segments model the kite and tether dynamics. - -# Components -- [`Point`](@ref): Point masses representing wing attachment points, dynamic bridle/tether points, and fixed ground anchors -- [`Group`](@ref): Collections of points that move together according to wing deformation (twist and trailing edge deflection) -- [`Segment`](@ref): Spring-damper elements connecting points -- [`Pulley`](@ref): Elements that redistribute line lengths between segments -- [`Tether`](@ref): Groups of segments with a common unstretched length -- [`Winch`](@ref): Ground-based winches that control tether lengths -- [`Wing`](@ref): Rigid wing bodies that serve as reference frames -- [`Transform`](@ref): Spatial transformations for initial positioning and orientation - -See the individual component documentation for detailed mathematical models and governing equations. -""" -struct SystemStructure - name::String - set::Settings - points::Vector{Point} - groups::Vector{Group} - segments::Vector{Segment} - pulleys::Vector{Pulley} - tethers::Vector{Tether} - winches::Vector{Winch} - wings::Vector{Wing} - transforms::Vector{Transform} - y::Array{Float64, 2} - x::Array{Float64, 2} - jac::Array{Float64, 3} -end - -""" - SystemStructure(name, set; points=Point[], groups=Group[], segments=Segment[], - pulleys=Pulley[], tethers=Tether[], winches=Winch[], - wings=Wing[], transforms=Transform[]) - -Constructs a SystemStructure object representing a complete kite system using a discrete mass-spring-damper model. - -## Components - -- **Points** - See [`Point`](@ref) for discrete mass dynamics -- **Segments** - See [`Segment`](@ref) for elastic spring-damper connections -- **Groups** - See [`Group`](@ref) for wing twist deformation modeling -- **Wings** - See [`Wing`](@ref) for rigid body dynamics -- **Pulleys** - See [`Pulley`](@ref) for length redistribution between segments -- **Tethers** - See [`Tether`](@ref) for segment groups with shared unstretched length -- **Winches** - See [`Winch`](@ref) for ground-based tether length control -- **Transforms** - See [`Transform`](@ref) for initial positioning and orientation - -## Physical Models -- **"ram"**: 4 deformable wing groups, complex pulley bridle system -- **"simple_ram"**: 4 deformable wing groups, direct bridle connections - -# Arguments -- `name::String`: Model identifier. "ram" and "simple_ram" are defined inside KiteModels.jl, provide a different name for a custom model. -- `set::Settings`: Configuration parameters (see [KiteUtils.Settings](https://ufechner7.github.io/KiteUtils.jl/stable/types/#KiteUtils.Settings)) - -# Returns -- `SystemStructure`: Complete system ready for building a [`SymbolicAWEModel`](@ref) - -# Examples -```julia -# Auto-generate from wing geometry -wing = RamAirWing(set) -sys_struct = SystemStructure(set, wing) - -# Manual construction -points = [Point(1, [0,0,0], STATIC), Point(2, [0,0,10], DYNAMIC)] -segments = [Segment(1, (1,2), BRIDLE)] -sys_struct = SystemStructure("custom", set; points, segments) -``` -""" -function SystemStructure(name, set; - points=Point[], - groups=Group[], - segments=Segment[], - pulleys=Pulley[], - tethers=Tether[], - winches=Winch[], - wings=Wing[], - transforms=Transform[], - ) - for (i, point) in enumerate(points) - @assert point.idx == i - @assert point.transform_idx <= length(transforms) - end - for (i, group) in enumerate(groups) - @assert group.idx == i - end - for (i, segment) in enumerate(segments) - @assert segment.idx == i - end - for (i, pulley) in enumerate(pulleys) - @assert pulley.idx == i - end - for (i, tether) in enumerate(tethers) - @assert tether.idx == i - end - for (i, winch) in enumerate(winches) - @assert winch.idx == i - set.l_tethers[i] = winch.tether_length - set.v_reel_outs[i] = winch.tether_vel - end - for (i, wing) in enumerate(wings) - @assert wing.idx == i - end - for (i, transform) in enumerate(transforms) - @assert transform.idx == i - set.elevations[i] = rad2deg(transform.elevation) - set.azimuths[i] = rad2deg(transform.azimuth) - set.headings[i] = rad2deg(transform.heading) - end - if length(wings) > 0 - ny = 3+length(wings[1].group_idxs)+3 - nx = 3+3+length(wings[1].group_idxs) - else - ny = 0 - nx = 0 - end - y = zeros(length(wings), ny) - x = zeros(length(wings), nx) - jac = zeros(length(wings), nx, ny) - set.physical_model = name - sys_struct = SystemStructure(name, set, points, groups, segments, pulleys, tethers, winches, wings, transforms, y, x, jac) - init!(sys_struct, set) - return sys_struct -end - -function SystemStructure(set::Settings, wing::RamAirWing) - length(set.bridle_fracs) != 4 && throw(ArgumentError("4 bridle fracs should be provided for all models.")) - - if set.physical_model == "ram" - return create_ram_sys_struct(set, wing) - elseif set.physical_model == "simple_ram" - return create_simple_ram_sys_struct(set, wing) - else - throw(ArgumentError("Undefined physical model")) - end -end - -function init!(transforms::Vector{Transform}, sys_struct::SystemStructure) - @unpack points, wings = sys_struct - for transform in transforms - # ==================== TRANSLATE ==================== # - base_pos, curr_base_pos = get_base_pos(transform, wings, points) - T = base_pos - curr_base_pos - for point in points - if point.transform_idx == transform.idx - point.pos_w .= point.pos_cad .+ T - end - end - for wing in wings - if wing.transform_idx == transform.idx - wing.pos_w .= wing.pos_cad .+ T - end - end - - # ==================== ROTATE ==================== # - curr_rot_pos = get_rot_pos(transform, wings, points) - curr_elevation = KiteUtils.calc_elevation(curr_rot_pos - base_pos) - curr_azimuth = -KiteUtils.azimuth_east(curr_rot_pos - base_pos) - curr_R_t_w = calc_R_t_w(curr_elevation, curr_azimuth) - R_t_w = calc_R_t_w(transform.elevation, transform.azimuth) - - for point in points - if point.transform_idx == transform.idx - vec = point.pos_w - base_pos - vec_along_z = rotate_around_z(curr_R_t_w' * vec, transform.heading) - point.pos_w .= base_pos + R_t_w * vec_along_z - end - if point.type == WING - wing = wings[point.wing_idx] - point.pos_b .= wing.R_b_c' * (point.pos_cad - wing.pos_cad) # TODO: test this - end - end - for wing in wings - if wing.transform_idx == transform.idx - vec = wing.pos_w - base_pos - vec_along_z = rotate_around_x(curr_R_t_w' * vec, transform.heading) - wing.pos_w .= base_pos + R_t_w * vec_along_z - for i in 1:3 - wing.R_b_c[:, i] .= R_t_w * rotate_around_x(curr_R_t_w' * wing.R_b_c[:, i], transform.heading) - end - end - end - - end -end - -function calc_pos(wing::RamAirWing, gamma, frac) - le_pos = [wing.le_interp[i](gamma) for i in 1:3] - chord = [wing.te_interp[i](gamma) for i in 1:3] .- le_pos - pos = le_pos .+ chord .* frac - return pos -end - -function create_tether(tether_idx, set, points, segments, tethers, attach_point, type, dynamics_type, z=[0,0,1]) - winch_pos = find_axis_point(attach_point.pos_cad, set.l_tether, z) - dir = winch_pos - attach_point.pos_cad - segment_idxs = Int16[] - for i in 1:set.segments - frac = i / set.segments - pos = attach_point.pos_cad + frac * dir - point_idx = length(points)+1 # last point idx - segment_idx = length(segments)+1 # last segment idx - if i == 1 - points = [points; Point(point_idx, pos, dynamics_type)] - segments = [segments; Segment(segment_idx, (attach_point.idx, point_idx), type)] - elseif i == set.segments - points = [points; Point(point_idx, pos, STATIC)] - segments = [segments; Segment(segment_idx, (point_idx-1, point_idx), type)] - else - points = [points; Point(point_idx, pos, dynamics_type)] - segments = [segments; Segment(segment_idx, (point_idx-1, point_idx), type)] - end - push!(segment_idxs, segment_idx) - end - tethers = [tethers; Tether(tether_idx, segment_idxs)] - return points, segments, tethers, tethers[end].idx -end - -function cad_to_body_frame(wing::RamAirWing, pos) - return wing.R_cad_body * (pos + wing.T_cad_body) -end - -# Find the point on the z-axis with distance l from P in the negative direction -# TODO: rename P to pos -function find_axis_point(P, l, v=[0,0,1]) - # Compute dot product v · P - v ⋅ P = v[1] * P[1] + v[2] * P[2] + v[3] * P[3] - # Compute discriminant - D = (v ⋅ P)^2 - norm(v)^2 * (norm(P)^2 - l^2) - D < 0 && error("No real solution: l is too small or parameters invalid") - # Solve quadratic for t, choose solution for negative direction - t = (v ⋅ P - √D) / norm(v)^2 - # Compute point Q = t * v - return [t * v[1], t * v[2], t * v[3]] -end - -function create_ram_sys_struct(set::Settings, vsm_wing::RamAirWing) - points = Point[] - groups = Group[] - segments = Segment[] - pulleys = Pulley[] - tethers = Tether[] - winches = Winch[] - wings = Wing[] - - attach_points = Point[] - - bridle_top_left = [cad_to_body_frame(vsm_wing, set.top_bridle_points[i]) for i in eachindex(set.top_bridle_points)] - bridle_top_right = [bridle_top_left[i] .* [1, -1, 1] for i in eachindex(set.top_bridle_points)] - - dynamics_type = set.quasi_static ? QUASI_STATIC : DYNAMIC - z = vsm_wing.R_cad_body[:,3] - - function create_bridle(bridle_top, gammas) - i_pnt = length(points) # last point idx - i_seg = length(segments) # last segment idx - i_pul = length(pulleys) # last pulley idx - i_grp = length(groups) # last group idx - - # ==================== CREATE DEFORMING WING GROUPS ==================== # - points = [ - points - Point(1+i_pnt, calc_pos(vsm_wing, gammas[1], set.bridle_fracs[1]), WING) - Point(2+i_pnt, calc_pos(vsm_wing, gammas[1], set.bridle_fracs[3]), WING) - Point(3+i_pnt, calc_pos(vsm_wing, gammas[1], set.bridle_fracs[4]), WING) - - Point(4+i_pnt, calc_pos(vsm_wing, gammas[2], set.bridle_fracs[1]), WING) - Point(5+i_pnt, calc_pos(vsm_wing, gammas[2], set.bridle_fracs[3]), WING) - Point(6+i_pnt, calc_pos(vsm_wing, gammas[2], set.bridle_fracs[4]), WING) - ] - groups = [ - groups - Group(1+i_grp, [1+i_pnt, 2+i_pnt, 3+i_pnt], vsm_wing, gammas[1], DYNAMIC, set.bridle_fracs[2]) - Group(2+i_grp, [4+i_pnt, 5+i_pnt, 6+i_pnt], vsm_wing, gammas[2], DYNAMIC, set.bridle_fracs[2]) - ] - - # ==================== CREATE PULLEY BRIDLE SYSTEM ==================== # - points = [ - points - Point(7+i_pnt, bridle_top[1], dynamics_type) - Point(8+i_pnt, bridle_top[2], WING) - Point(9+i_pnt, bridle_top[3], dynamics_type) - Point(10+i_pnt, bridle_top[4], dynamics_type) - - Point(11+i_pnt, bridle_top[2] - 1z, dynamics_type) - - Point(12+i_pnt, bridle_top[1] - 2z, dynamics_type) - Point(13+i_pnt, bridle_top[3] - 2z, dynamics_type) - - Point(14+i_pnt, bridle_top[1] - 4z, dynamics_type) - Point(15+i_pnt, bridle_top[3] - 4z, dynamics_type) - ] - segments = [ - segments - Segment(1+i_seg, (1+i_pnt, 7+i_pnt), BRIDLE) - Segment(2+i_seg, (2+i_pnt, 9+i_pnt), BRIDLE) - Segment(3+i_seg, (3+i_pnt, 10+i_pnt), BRIDLE) - - Segment(4+i_seg, (4+i_pnt, 7+i_pnt), BRIDLE) - Segment(5+i_seg, (5+i_pnt, 9+i_pnt), BRIDLE) - Segment(6+i_seg, (6+i_pnt, 10+i_pnt), BRIDLE) - - Segment(7+i_seg, (7+i_pnt, 12+i_pnt), BRIDLE; l0=2) - Segment(8+i_seg, (8+i_pnt, 11+i_pnt), BRIDLE; l0=1) - Segment(9+i_seg, (9+i_pnt, 13+i_pnt), BRIDLE; l0=2) - Segment(10+i_seg, (10+i_pnt, 15+i_pnt), BRIDLE; l0=4) - - Segment(11+i_seg, (11+i_pnt, 12+i_pnt), BRIDLE; l0=1) - Segment(12+i_seg, (11+i_pnt, 13+i_pnt), BRIDLE; l0=1) - - Segment(13+i_seg, (12+i_pnt, 14+i_pnt), BRIDLE; l0=2) - Segment(14+i_seg, (13+i_pnt, 14+i_pnt), BRIDLE; l0=2) - Segment(15+i_seg, (13+i_pnt, 15+i_pnt), BRIDLE; l0=2) - ] - pulleys = [ - pulleys - Pulley(1+i_pul, (11+i_seg, 12+i_seg), dynamics_type) - Pulley(2+i_pul, (14+i_seg, 15+i_seg), dynamics_type) - ] - push!(attach_points, points[end-1]) - push!(attach_points, points[end]) - return nothing - end - - gammas = [-3/4, -1/4, 1/4, 3/4] * vsm_wing.gamma_tip - create_bridle(bridle_top_left, gammas[[1,2]]) - create_bridle(bridle_top_right, gammas[[3,4]]) - - points, segments, tethers, left_power_idx = create_tether(1, set, points, segments, tethers, attach_points[1], POWER_LINE, dynamics_type, z) - points, segments, tethers, right_power_idx = create_tether(2, set, points, segments, tethers, attach_points[3], POWER_LINE, dynamics_type, z) - points, segments, tethers, left_steering_idx = create_tether(3, set, points, segments, tethers, attach_points[2], STEERING_LINE, dynamics_type, z) - points, segments, tethers, right_steering_idx = create_tether(4, set, points, segments, tethers, attach_points[4], STEERING_LINE, dynamics_type, z) - - winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx], set.l_tether)] - winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)] - winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)] - - wings = [Wing(1, [1,2,3,4], I(3), zeros(3))] - transforms = [Transform(1, deg2rad(set.elevation), deg2rad(set.azimuth), deg2rad(set.heading); - base_pos= zeros(3), base_point_idx=points[end].idx, wing_idx=1)] - - return SystemStructure(set.physical_model, set; points, groups, segments, pulleys, tethers, winches, wings, transforms) -end - -function create_simple_ram_sys_struct(set::Settings, wing::RamAirWing) - points = Point[] - groups = Group[] - segments = Segment[] - pulleys = Pulley[] - tethers = Tether[] - winches = Winch[] - - dynamics_type = set.quasi_static ? QUASI_STATIC : DYNAMIC - gammas = [-3/4, -1/4, 1/4, 3/4] * wing.gamma_tip - - bridle_top_left = [wing.R_cad_body * (set.top_bridle_points[i] + wing.T_cad_body) for i in eachindex(set.top_bridle_points)] # cad to kite frame - bridle_top_right = [bridle_top_left[i] .* [1, -1, 1] for i in eachindex(set.top_bridle_points)] - - # ==================== CREATE DEFORMING WING GROUPS ==================== # - points = [ - points - Point(1, calc_pos(wing, gammas[1], set.bridle_fracs[4]), WING) - Point(2, calc_pos(wing, gammas[2], set.bridle_fracs[4]), WING) - Point(3, calc_pos(wing, gammas[3], set.bridle_fracs[4]), WING) - Point(4, calc_pos(wing, gammas[4], set.bridle_fracs[4]), WING) - ] - groups = [ - groups - Group(1, [1], wing, gammas[1], DYNAMIC, set.bridle_fracs[2]) - Group(2, [2], wing, gammas[2], DYNAMIC, set.bridle_fracs[2]) - Group(3, [3], wing, gammas[3], DYNAMIC, set.bridle_fracs[2]) - Group(4, [4], wing, gammas[4], DYNAMIC, set.bridle_fracs[2]) - ] - # ==================== CREATE PULLEY BRIDLE SYSTEM ==================== # - points = [ - points - Point(5, bridle_top_left[2], WING) - Point(6, bridle_top_left[4], dynamics_type) - Point(7, bridle_top_right[2], WING) - Point(8, bridle_top_right[4], dynamics_type) - ] - - segments = [ - segments - Segment(1, (1, 6), BRIDLE) - Segment(2, (2, 6), BRIDLE) - Segment(3, (3, 8), BRIDLE) - Segment(4, (4, 8), BRIDLE) - ] - - points, segments, tethers, left_power_idx = create_tether(1, set, points, segments, tethers, points[5], POWER_LINE, dynamics_type) - points, segments, tethers, right_power_idx = create_tether(2, set, points, segments, tethers, points[7], POWER_LINE, dynamics_type) - points, segments, tethers, left_steering_idx = create_tether(3, set, points, segments, tethers, points[6], STEERING_LINE, dynamics_type) - points, segments, tethers, right_steering_idx = create_tether(4, set, points, segments, tethers, points[8], STEERING_LINE, dynamics_type) - - winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx], set.l_tether)] - winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)] - winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)] - - wings = [Wing(1, [1,2,3,4], I(3), zeros(3))] - transforms = [Transform(1, deg2rad(set.elevation), deg2rad(set.azimuth), deg2rad(set.heading); - base_pos= zeros(3), base_point_idx=points[end].idx, wing_idx=1)] - - return SystemStructure(set.physical_model, set; points, groups, segments, pulleys, tethers, winches, wings, transforms) -end - -function init!(sys_struct::SystemStructure, set::Settings) - @unpack points, groups, segments, pulleys, tethers, winches, wings, transforms = sys_struct - - for segment in segments - (segment.type === BRIDLE) && (segment.diameter = 0.001set.bridle_tether_diameter) - (segment.type === POWER_LINE) && (segment.diameter = 0.001set.power_tether_diameter) - (segment.type === STEERING_LINE) && (segment.diameter = 0.001set.steering_tether_diameter) - @assert (0 < segment.diameter < 1) - end - - for winch in winches - winch.tether_length = set.l_tethers[winch.idx] - winch.tether_vel = set.v_reel_outs[winch.idx] - end - - (length(groups) > 0) && (first_moment_frac = groups[1].moment_frac) - for group in groups - group.twist = 0.0 - group.twist_vel = 0.0 - @assert group.moment_frac ≈ first_moment_frac "All group.moment_frac must be the same." - end - - for transform in transforms - transform.elevation = deg2rad(set.elevations[transform.idx]) - transform.azimuth = deg2rad(set.azimuths[transform.idx]) - transform.heading = deg2rad(set.headings[transform.idx]) - end - - for segment in segments - (segment.l0 ≈ 0) && (segment.l0 = norm(points[segment.point_idxs[1]].pos_cad - points[segment.point_idxs[2]].pos_cad)) - @assert (segment.l0 > 0) - end - - for pulley in pulleys - segment1, segment2 = segments[pulley.segment_idxs[1]], segments[pulley.segment_idxs[2]] - pulley.sum_length = segment1.l0 + segment2.l0 - pulley.length = segment1.l0 - pulley.vel = 0.0 - @assert !(pulley.sum_length ≈ 0) - end - - init!(transforms, sys_struct) - return nothing -end diff --git a/test/aqua.jl b/test/aqua.jl index 7522a0df..f592486a 100644 --- a/test/aqua.jl +++ b/test/aqua.jl @@ -1,7 +1,12 @@ # SPDX-FileCopyrightText: 2024 Uwe Fechner # SPDX-License-Identifier: MIT -using Aqua +using Pkg +if ! ("Aqua" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end + +using KiteModels, Aqua, Test @testset "Aqua.jl" begin Aqua.test_all( KiteModels; diff --git a/test/bench3.jl b/test/bench3.jl index 226e6e67..4a61b677 100644 --- a/test/bench3.jl +++ b/test/bench3.jl @@ -2,7 +2,7 @@ # SPDX-License-Identifier: MIT using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils -using KiteModels, KitePodModels +using KiteModels, KitePodModels, AtmosphericModels if ! @isdefined SEGMENTS const SEGMENTS = se().segments @@ -105,4 +105,4 @@ println(msg, "\n") # ███▆███▆▆▆▅▁▅▄▄▅▄▅▃▄▃▃▃▃▁▄▇▇▇▆▇▇▇▆▇▆▅▄▃▅▄▃▁▃▁▄▃▄▁▄▄▁▁▃▁▁▃█▇▅▅ █ # 523 ns Histogram: log(frequency) by time 676 ns < -# Memory estimate: 0 bytes, allocs estimate: 0. \ No newline at end of file +# Memory estimate: 0 bytes, allocs estimate: 0. diff --git a/test/bench4.jl b/test/bench4.jl index 0ea9c5d1..0b50a433 100644 --- a/test/bench4.jl +++ b/test/bench4.jl @@ -7,7 +7,7 @@ if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies)) Pkg.update() end using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils -using KiteModels, KitePodModels +using KiteModels, KitePodModels, AtmosphericModels set_data_path(joinpath(dirname(dirname(pathof(KiteModels))), "data")) set = load_settings("system.yaml") diff --git a/test/runtests.jl b/test/runtests.jl index 92ef7faf..0a542ceb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,9 +30,6 @@ KiteUtils.set_data_path("") include("bench3.jl") include("bench4.jl") end - if ! haskey(ENV, "NO_MTK") - include("test_ram_air_kite.jl") - end include("test_helpers.jl") println("--> 4") include("test_inertia_calculation.jl") @@ -40,4 +37,4 @@ KiteUtils.set_data_path("") include("test_interface.jl") println("--> 6") include("aqua.jl") -end \ No newline at end of file +end diff --git a/test/test_for_precompile.jl b/test/test_for_precompile.jl index b71269a9..12cc096d 100644 --- a/test/test_for_precompile.jl +++ b/test/test_for_precompile.jl @@ -113,7 +113,7 @@ if ! haskey(ENV, "NO_MTK") # Initialize at elevation KiteModels.init!(s; prn=false, precompile=true) - find_steady_state!(s) + KiteModels.SymbolicAWEModels.find_steady_state!(s) steps = Int(round(10 / 0.05)) logger = Logger(length(s.sys_struct.points), steps) sys_state = SysState(s) diff --git a/test/test_interface.jl b/test/test_interface.jl index bf78d62d..4e0127a3 100644 --- a/test/test_interface.jl +++ b/test/test_interface.jl @@ -9,8 +9,6 @@ using KiteModels: init! set = deepcopy(load_settings("system.yaml")) kps3::KPS3 = KPS3(set) kps4::KPS4 = KPS4(set) -set_ram = deepcopy(load_settings("system_ram.yaml")) -sam::SymbolicAWEModel = SymbolicAWEModel(set_ram) @testset "KPS3 constructor interface" begin @test kps3 isa KiteUtils.AbstractKiteModel @@ -26,13 +24,6 @@ end @test kps4.am isa AtmosphericModel @test kps4.iter isa Int64 end -@testset "SymbolicAWEModel constructor interface" begin - @test sam isa KiteUtils.AbstractKiteModel - @test sam.set isa KiteUtils.Settings - @test isnothing(sam.integrator) - @test sam.am isa AtmosphericModel - @test sam.iter isa Int64 -end @testset "KPS3 init! interface" begin kps4.set.upwind_dir = rad2deg(-pi/2) @@ -48,10 +39,4 @@ end @test kps4.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator end -@testset "SymbolicAWEModel init! interface" begin - sam.set.upwind_dir = rad2deg(-pi/2) - integ = init!(sam; stiffness_factor=0.5, delta=0.0001, prn=false) - @test integ isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator - @test kps4.integrator isa KiteModels.OrdinaryDiffEqCore.ODEIntegrator -end -nothing \ No newline at end of file +nothing diff --git a/test/test_kps3.jl b/test/test_kps3.jl index 0013c4e6..0135e65f 100644 --- a/test/test_kps3.jl +++ b/test/test_kps3.jl @@ -3,7 +3,7 @@ using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils -using KiteModels, KitePodModels +using KiteModels, KitePodModels, AtmosphericModels const SEGMENTS = se().segments kcu::KCU = KCU(se()) diff --git a/test/test_kps3_fails.jl b/test/test_kps3_fails.jl index 8da92236..77199e7d 100644 --- a/test/test_kps3_fails.jl +++ b/test/test_kps3_fails.jl @@ -8,7 +8,7 @@ end using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils -using KiteModels, KitePodModels +using KiteModels, KitePodModels, AtmosphericModels const SEGMENTS = se().segments if ! @isdefined kcu diff --git a/test/test_kps4.jl b/test/test_kps4.jl index 61003637..1f1a0758 100644 --- a/test/test_kps4.jl +++ b/test/test_kps4.jl @@ -7,7 +7,7 @@ if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies)) Pkg.update() end using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils -using KiteModels, KitePodModels +using KiteModels, KitePodModels, AtmosphericModels set_data_path(joinpath(dirname(dirname(pathof(KiteModels))), "data")) set = deepcopy(load_settings("system.yaml")) @@ -591,4 +591,4 @@ end end end -nothing \ No newline at end of file +nothing diff --git a/test/test_orientation.jl b/test/test_orientation.jl index 1c1d458c..8e4c02e6 100644 --- a/test/test_orientation.jl +++ b/test/test_orientation.jl @@ -4,6 +4,7 @@ using KiteUtils using KiteModels using KitePodModels +using AtmosphericModels using LinearAlgebra using StaticArrays @@ -315,4 +316,4 @@ end end nothing -calc_azimuth \ No newline at end of file +calc_azimuth diff --git a/test/test_ram_air_kite.jl b/test/test_ram_air_kite.jl deleted file mode 100644 index a9bca9d8..00000000 --- a/test/test_ram_air_kite.jl +++ /dev/null @@ -1,330 +0,0 @@ -# SPDX-FileCopyrightText: 2025 Bart van de Lint -# SPDX-License-Identifier: MIT - -using Test, LinearAlgebra, KiteUtils, VortexStepMethod -using ControlPlots -using KiteModels -using Statistics - -old_path = get_data_path() -package_data_path = joinpath(dirname(dirname(pathof(KiteModels))), "data") -temp_data_path = joinpath(tempdir(), "data") -Base.Filesystem.cptree(package_data_path, temp_data_path; force=true) -set_data_path(temp_data_path) - -# Testing tolerance -const TOL = 1e-5 -const BUILD_SYS = true - -@testset verbose = true "SymbolicAWEModel MTK Model Tests" begin - # Initialize model - set = se("system_ram.yaml") - set.segments = 3 - set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] - set.quasi_static = false - set.physical_model = "ram" - - @info "Creating s:" - @time s = SymbolicAWEModel(set) - - s.set.abs_tol = 1e-4 - s.set.rel_tol = 1e-4 - - # Initialize at elevation - set.elevation = 80.0 - - @testset "Model Initialization Chain" begin - # if BUILD_SYS - # Delete existing problem file to force init! - @info "Data path: $(get_data_path())" - model_path = joinpath(get_data_path(), KiteModels.get_model_name(s.set)) - # if isfile(model_path) - # @info "Removing existing serialized problem from $model_path to test full initialization" - # rm(model_path) - # end - - # 1. First time initialization - should create new model - @info "Testing initial init! (should create new model if it doesn't exist yet)..." - @time KiteModels.init!(s; prn=true) - - # Check that serialization worked - @test isfile(model_path) - - # Check initialization results - @test !isnothing(s.integrator) - @test !isnothing(s.sys) - @test !isnothing(s.sys_struct) - # end - s.integrator = nothing - s.sys = nothing - - # Keep references to first integrator and point system - first_integrator_ptr = objectid(s.integrator) - first_sys_struct_ptr = objectid(s.sys_struct) - - # 2. First init! - should load from serialized file - @info "Testing first init! (should load serialized file)..." - @time KiteModels.init!(s; prn=true, reload=false) - next_step!(s) - - # Check that it's a new integrator - second_integrator_ptr = objectid(s.integrator) - second_sys_struct_ptr = objectid(s.sys_struct) - @test first_integrator_ptr != second_integrator_ptr - @test first_sys_struct_ptr == second_sys_struct_ptr - - # 3. Second init! - should reuse existing integrator - @info "Testing second init! (should reuse integrator)..." - @time KiteModels.init!(s; prn=true, reload=false) - - # This should create a new point system but reuse the existing integrator - third_integrator_ptr = objectid(s.integrator) - third_sys_struct_ptr = objectid(s.sys_struct) - @test second_integrator_ptr == third_integrator_ptr # Should be the same - @test second_sys_struct_ptr == third_sys_struct_ptr - - # Get positions using SysState - sys_state = KiteModels.SysState(s) - - # Check dimension consistency - # Note: pos_integrator is no longer directly fetched, comparing SysState to sys_struct - @test length(sys_state.X) == length(s.sys_struct.points) - - # Compare positions in different representations - for (i, point) in enumerate(s.sys_struct.points) - # Points' world positions should match SysState positions - point_pos = point.pos_w - sys_state_pos = [sys_state.X[i], sys_state.Y[i], sys_state.Z[i]] - - # Use norm for comparison as exact vector match might be too strict due to float precision - @test isapprox(norm(point_pos), norm(sys_state_pos), rtol=1e-2) - - # Positions should not be zero (except ground points) - if point.type != KiteModels.STATIC # Skip ground points which might be at origin - @test norm(point_pos) > 0.1 - @test norm(sys_state_pos) > 0.1 - end - end - end - - @testset "State Consistency" begin - KiteModels.init!(s, prn=true, reload=false) - sys_state_before = KiteModels.SysState(s) - @test isapprox(norm(s.integrator[s.sys.Q_b_w]), 1.0, atol=TOL) - @test isapprox(sys_state_before.elevation, deg2rad(set.elevation), atol=1e-2) - - # Change measurement and reinitialize - old_elevation = set.elevation - set.elevation = 85.0 - KiteModels.init!(s, prn=true, reload=false) - - # Get new state using SysState - sys_state_after = KiteModels.SysState(s) - - # Verify state changed according to measurement - @test !isapprox(sys_state_after.elevation, old_elevation, atol=1e-2) - @test isapprox(sys_state_after.elevation, deg2rad(85.0), atol=1e-2) - - @testset "set_depower_steering!" begin - initial_tether_lengths = s.get_tether_length(s.integrator) - depower = 0.1 - steering = 0.05 - KiteModels.set_depower_steering!(s, depower, steering) - new_tether_lengths = s.set_tether_length - @test !isapprox(new_tether_lengths, initial_tether_lengths) - # Verify the changes based on the equations - len = s.set_tether_length - len1 = initial_tether_lengths[1] - len2 = 0.5 * (2*depower*KiteModels.min_chord_length(s) + 2*len1 + steering*KiteModels.min_chord_length(s)) - len3 = 0.5 * (2*depower*KiteModels.min_chord_length(s) + 2*len1 - steering*KiteModels.min_chord_length(s)) - @test isapprox(len[2], len2) - @test isapprox(len[3], len3) - @test isapprox(KiteModels.min_chord_length(s), 0.434108) - end - - @testset "set_v_wind_ground!" begin - initial_wind_speed = s.set.v_wind - initial_upwind_dir = deg2rad(s.set.upwind_dir) - @test initial_upwind_dir == -π/2 - @test s.integrator[s.sys.wind_vec_gnd[1]] == s.set.v_wind - - # Set new wind speed and direction - new_wind_speed = 10.0 - new_upwind_dir = -pi/4 - KiteModels.set_v_wind_ground!(s, new_wind_speed, new_upwind_dir) - - # Check if wind speed and direction have been updated correctly - @test norm(s.integrator[s.sys.wind_vec_gnd]) ≈ new_wind_speed - @test s.integrator.ps[s.sys.upwind_dir] ≈ -pi/4 # Default upwind_dir - @test s.integrator[s.sys.wind_vec_gnd[1]] ≈ -s.integrator[s.sys.wind_vec_gnd[2]] - - KiteModels.set_v_wind_ground!(s, initial_wind_speed, initial_upwind_dir) - @test s.integrator[s.sys.wind_vec_gnd[1]] ≈ initial_wind_speed - @test s.integrator.ps[s.sys.upwind_dir] ≈ -pi/2 # Default upwind_dir - end - end - - function test_step(s, d_set_values=zeros(3); dt=0.05, steps=5) - s.integrator.ps[s.sys.stabilize] = true - for i in 1:1÷0.1 - next_step!(s; dt, vsm_interval=1) - end - s.integrator.ps[s.sys.stabilize] = false - @info "Stepping" - for _ in 1:steps - set_values = -s.set.drum_radius * s.integrator[s.sys.winch_force] + d_set_values - next_step!(s; set_values, dt) - # Use SysState to get heading if needed, or directly from integrator if simpler - # sys_state_step = KiteModels.SysState(s) - # @show sys_state_step.heading # Example if heading is in SysState - @show s.integrator[s.sys.heading] # Keep direct access if simpler for this specific value - end - end - - function test_plot(s) - @testset "Plotting of SymbolicAWEModel" begin - function plot_(zoom, front) - plt.figure("Kite") - lines, sc, txt = plot(s, 0.0; zoom, front) - plt.show(block=false) - sleep(1) - @test !isnothing(lines) - @test length(lines) ≥ 1 # Should have at least one line - @test !isnothing(sc) # Should have scatter points - @test !isnothing(txt) # Should have time text - end - plot_(false, false) - plot_(false, true) - plot_(true, false) - plot_(true, true) - end - end - - @testset "Simulation Step with SysState" begin - # Basic step and time advancement test - KiteModels.init!(s; prn=true, reload=false) - sys_state_before = KiteModels.SysState(s) - - # Run a simulation step with zero set values - set_values = [0.0, 0.0, 0.0] - dt = 1/s.set.sample_freq - next_step!(s; set_values, dt=dt) - # Update sys_state_before *after* the step to compare with the state *before* the loop - KiteModels.update_sys_state!(sys_state_before, s) - @test isapprox(s.integrator.t, dt, atol=TOL) - - # Run multiple steps - num_steps = 10 - for _ in 1:num_steps - next_step!(s; set_values, dt=dt) - end - sys_state_after = KiteModels.SysState(s) # Get state after the loop - # Compare state after loop with state after first step (stored in sys_state_before) - @test any(abs.(sys_state_after.X .- sys_state_before.X) .> 0.01) # Reduced tolerance slightly - - @testset "Course Direction at 60 Degrees Elevation" begin # Corrected description - # Initialize at 60 degrees elevation - set.elevation = 60.0 - - KiteModels.init!(s; prn=true) - - # Verify initial conditions using SysState - sys_state_init = KiteModels.SysState(s) - @test sys_state_init.elevation ≈ deg2rad(set.elevation) atol=1e-2 # Use relaxed tolerance consistent with other tests - @test sys_state_init.azimuth ≈ deg2rad(set.azimuth) atol=1e-2 - - # Run simulation steps - test_step(s) - - # Check course direction using SysState - sys_state = KiteModels.SysState(s) - @info "Course at 60 deg elevation:" sys_state.course - - # At 60 degrees elevation, course should be roughly forward - @show sys_state.course - @test sys_state.course ≈ 0.0 atol=deg2rad(45.0) - end - - # Utility function to calculate the signed angle difference between two angles - function angle_diff(angle1, angle2) - diff = angle1 - angle2 - # Normalize to [-π, π] - while diff > π; diff -= 2π; end - while diff < -π; diff += 2π; end - return diff - end - - if build_is_production_build - @testset "Steering Response Using SysState" begin - # Initialize model at moderate elevation - set.elevation = 70.0 - KiteModels.init!(s; prn=true, reload=false) - test_step(s) - sys_state_initial = KiteModels.SysState(s) - - # steering right - KiteModels.init!(s; prn=true, reload=false) - test_step(s, [0, 10, -10]; steps=20) - sys_state_right = KiteModels.SysState(s) - - # steering left - KiteModels.init!(s; prn=true, reload=false) - test_step(s, [0, -10, 10]; steps=20) - sys_state_left = KiteModels.SysState(s) - - # Check steering values - @info "Steering:" sys_state_right.steering sys_state_left.steering - @test sys_state_right.steering > 3.0 - @test sys_state_left.steering < -3.0 - - # Check heading changes - right_heading_diff = angle_diff(sys_state_right.heading, sys_state_initial.heading) - @test right_heading_diff ≈ 0.9 atol=0.2 - left_heading_diff = angle_diff(sys_state_left.heading, sys_state_initial.heading) - @test left_heading_diff ≈ -0.9 atol=0.2 - end - test_plot(s) - end - end - - @testset "Just a tether, without winch or kite" begin - set.segments = 20 - dynamics_type = DYNAMIC - - points = Point[] - segments = Segment[] - - points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0)) - - segment_idxs = Int[] - for i in 1:set.segments - point_idx = i+1 - pos = [0.0, 0.0, i * set.l_tether / set.segments] - push!(points, Point(point_idx, pos, dynamics_type; wing_idx=0)) - segment_idx = i - push!(segments, Segment(segment_idx, (point_idx-1, point_idx), BRIDLE)) - push!(segment_idxs, segment_idx) - end - - transforms = [Transform(1, deg2rad(-80), 0.0, 0.0; - base_pos=[0.0, 0.0, 50.0], base_point_idx=points[1].idx, rot_point_idx=points[end].idx)] - sys_struct = SystemStructure("tether", set; points, segments, transforms) - - sam = SymbolicAWEModel(set, sys_struct) - sys = sam.sys - KiteModels.init!(sam; remake=false) - @test isapprox(sam.integrator[sam.sys.pos[:, end]], [8.682408883346524, 0.0, 0.7596123493895988], atol=1e-2) - for i in 1:100 - next_step!(sam) - end - @test sam.integrator[sam.sys.pos[1, end]] > 0.8set.l_tether - @test isapprox(sam.integrator[sam.sys.pos[2, end]], 0.0, atol=1.0) - test_plot(s) - set.elevation = 80 - end -end - -# Restore original data path -set_data_path(old_path) - -nothing diff --git a/test/test_steady_state.jl b/test/test_steady_state.jl index c1d55a9c..c57e6065 100644 --- a/test/test_steady_state.jl +++ b/test/test_steady_state.jl @@ -2,7 +2,7 @@ # SPDX-License-Identifier: MIT using StaticArrays, LinearAlgebra, KiteUtils -using KiteModels, KitePodModels +using KiteModels, KitePodModels, AtmosphericModels set_data_path(joinpath(dirname(dirname(pathof(KiteModels))), "data")) set = deepcopy(load_settings("system.yaml"))