diff --git a/AUTHORS.md b/AUTHORS.md index e15262068..a144bee7e 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -11,7 +11,7 @@ Bart van de Lint joined the project in 2024. His major contribution is the first ## Developers - Uwe Fechner, Delft/ Den Haag, The Netherlands - original author of the `KPS3` and `KPS4` kite models -- Bart van de Lint, Trondheim, Norway - author of the `RamAirKite` model +- Bart van de Lint, Trondheim, Norway - author of the `SymbolicAWEModel` model ## Contributors - Daan van Wolffelaar, Delft, The Netherlands - contributed the scripts diff --git a/CHANGELOG.md b/CHANGELOG.md index 7d720fa79..f31a0954d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,16 +8,16 @@ SPDX-License-Identifier: MIT - added licenses to each file, the command `pipx run reuse lint` succeeds now - add the command above to the CI scripts - the script `create_xz_file` -- the option to linearize the RamAirKite system using ModelingToolkit +- the option to linearize the SymbolicAWEModel system using ModelingToolkit - a simplified ram air kite model for faster development and testing - the example `lin_ram_model.jl` to show how to linearize a model -- add the page `Examples RamAirKite` do the documentation +- add the page `Examples SymbolicAWEModel` do the documentation #### Changed - the example `ram_air_kite.jl` can now be run like this `SIMPLE=true; include("examples/ram_air_kite.jl")` - the package `Rotations` is no longer re-exported - improved documentation #### Fixed -- small fixes of the RamAirKite model +- small fixes of the SymbolicAWEModel model ### KiteModels v0.7.3 2025-05-05 #### Fixed @@ -32,7 +32,7 @@ SPDX-License-Identifier: MIT - fixed or documented issues found by `Aqua.jl` - made `DSP` a test dependency - remove package `OrdinaryDiffEqSDIRK` -- improve documentation for `RamAirKite` +- improve documentation for `SymbolicAWEModel` #### Added - the examples `calc_spectrum.jl` and `plot_spectrum.jl` to the menu - the quality insurance package `Aqua.jl` @@ -45,10 +45,10 @@ SPDX-License-Identifier: MIT #### Added - added `mwe_26.jl` for debugging the initial state solver - the example `ram_air_kite.jl` -- the struct `PointMassSystem` for easy definition of the kite power system +- the struct `SystemStructure` for easy definition of the kite power system #### Changed -- BREAKING: the model KPS_3L was renamed to RamAirKite -- the RamAirKite model is using the **VortexStepMethod** with a deforming wing now +- BREAKING: the model KPS_3L was renamed to SymbolicAWEModel +- the SymbolicAWEModel model is using the **VortexStepMethod** with a deforming wing now - bump KiteUtils to `v0.10` - bump ModellingToolkit to `9.72` - bump VortexStepMethod to `1.2.5` @@ -128,7 +128,7 @@ SPDX-License-Identifier: MIT - added tests for calc_azimuth(s::AKM), the azimuth in wind reference frame - re-enable logging of the angles of attack of the three plates - `steering_test_4p.jl` now calculates both `c1` and `c2` of the turn-rate law -- the environment variable `NO_MTK` disables the pre-compilation of the `RamAirKite` model +- the environment variable `NO_MTK` disables the pre-compilation of the `SymbolicAWEModel` model to save time during development - the script `menu2.jl` for model verification was added @@ -168,7 +168,7 @@ SPDX-License-Identifier: MIT - always specify the `system.yaml` file to use in the examples, always use `load_settings` instead of `se`. This ensures that the settings are always freshly loaded from the file when the script is launched, so any changes to the settings become immediately effective. -- the RamAirKite model was replaced by the pure ModelingToolkit (MTK) based version. This allows not only a much faster simulation, but the results are also much more accurate. +- the SymbolicAWEModel model was replaced by the pure ModelingToolkit (MTK) based version. This allows not only a much faster simulation, but the results are also much more accurate. ### KiteModels v0.6.5 - 2024-08-12 #### Changed @@ -247,7 +247,7 @@ to the settings become immediately effective. - update Documenter to v1.0 #### Added -- add type `RamAirKite`, which is now only a copy of `KPS4`, but shall implement a kite with the steering +- add type `SymbolicAWEModel`, which is now only a copy of `KPS4`, but shall implement a kite with the steering lines going to the ground ### KiteModels v0.5.11 - 2024-04-04 diff --git a/README.md b/README.md index 6fbb58fc4..25c5100ea 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,7 @@ This package is part of Julia Kite Power Tools, which consists of the following - a new 5-point model based on ModelingToolkit (MTK) is in development; this will allow to create linearized models around any operation point and to do analysis in the frequency domain. #### April 2025 -- a new model `RamAirKite` was contributed, based on the package [VortexStepMethod](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl) +- a new model `SymbolicAWEModel` was contributed, based on the package [VortexStepMethod](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl) #### November 2024 - the four point kite model KPS4 was extended to include aerodynamic damping of pitch oscillations; for this purpose, the parameters `cmq` and `cord_length` must be defined in `settings.yaml` diff --git a/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz b/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz deleted file mode 100644 index 30824cac9..000000000 Binary files a/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz and /dev/null differ diff --git a/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz b/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz deleted file mode 100644 index 39c8e8bac..000000000 Binary files a/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz and /dev/null differ diff --git a/data/settings_ram.yaml b/data/settings_ram.yaml index e03d9d805..c2b113464 100644 --- a/data/settings_ram.yaml +++ b/data/settings_ram.yaml @@ -23,7 +23,7 @@ solver: kite: model: "data/ram_air_kite_body.obj" # 3D model of the kite foil_file: "data/ram_air_kite_foil.dat" # filename for the foil shape - physical_model: "ram" # name of the kite model to use (KPS3, KPS4 or RamAirKite) + physical_model: "ram" # name of the kite model to use (KPS3, KPS4 or SymbolicAWEModel) top_bridle_points: # top bridle points that are not on the kite body in CAD frame - [0.290199, 0.784697, -2.61305] - [0.392683, 0.785271, -2.61201] diff --git a/docs/Project.toml b/docs/Project.toml index 7d3f57cfb..feb2e7f2a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,11 +5,12 @@ KiteModels = "b94af626-7959-4878-9336-2adc27959007" KitePodModels = "9de5dc81-f971-414a-927b-652b2f41c539" KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533" OpenSSL_jll = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" [compat] -OpenSSL_jll = "~3.0.0" ControlPlots = "0.2.7" Documenter = "1.11" KitePodModels = "0.3" KiteUtils = "0.7, 0.8, 0.9, 0.10" +OpenSSL_jll = "~3.0.0" julia = "1.10, 1.11" diff --git a/docs/make.jl b/docs/make.jl index 7568bf3ca..cc7ea67c5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,17 +8,14 @@ if ("TestEnv" ∈ keys(Pkg.project().dependencies)) using TestEnv; TestEnv.activate() end end -using ControlPlots +using ControlPlots, VortexStepMethod using KiteModels using Documenter DocMeta.setdocmeta!(KiteModels, :DocTestSetup, :(using KitePodSimulator); recursive=true) makedocs(; - modules=[KiteModels, - isdefined(Base, :get_extension) ? - Base.get_extension(KiteModels, :KiteModelsControlPlotsExt) : - KiteModels.KiteModelsControlPlotsExt], + modules=[KiteModels], authors="Uwe Fechner , Bart van de Lint and contributors", repo="https://github.com/ufechner7/KiteModels.jl/blob/{commit}{path}#{line}", sitename="KiteModels.jl", @@ -32,11 +29,12 @@ makedocs(; "Home" => "index.md", "Types" => "types.md", "Functions" => "functions.md", - "RamAirKite" => "ram_air_kite.md", + "SymbolicAWEModel" => "ram_air_kite.md", "Parameters" => "parameters.md", "Examples 1p" => "examples.md", "Examples 4p" => "examples_4p.md", - "Examples RamAirKite" => "examples_ram_air.md", + "Examples SymbolicAWEModel" => "examples_ram_air.md", + "SystemStructure for custom models" => "tutorial_system_structure.md", "Quickstart" => "quickstart.md", "Advanced usage" => "advanced.md", ], diff --git a/docs/src/examples_ram_air.md b/docs/src/examples_ram_air.md index bbc5bb416..dc54e32c0 100644 --- a/docs/src/examples_ram_air.md +++ b/docs/src/examples_ram_air.md @@ -26,7 +26,7 @@ Expected output for first run: ``` [ Info: Loading packages Time elapsed: 7.483472342 s -[ Info: Creating wing, aero, vsm_solver, point_system and s: +[ Info: Creating wing, aero, vsm_solver, system_structure and s: Time elapsed: 15.341197455 s [ Info: Creating ODESystem 4.316010 seconds (8.72 M allocations: 222.606 MiB, 1.42% gc time, 25.46% compilation time: 14% of which was recompilation) @@ -48,7 +48,7 @@ file in the `data` folder: ``` [ Info: Loading packages Time elapsed: 7.396961284 s -[ Info: Creating wing, aero, vsm_solver, point_system and s: +[ Info: Creating wing, aero, vsm_solver, system_structure and s: Time elapsed: 15.387790726 s [ Info: Initialized integrator in 29.545349428 seconds [ Info: System initialized at: @@ -83,7 +83,7 @@ include("examples/lin_ram_model.jl") ``` See: [`lin_ram_model.jl`](https://github.com/ufechner7/KiteModels.jl/blob/main/examples/lin_ram_model.jl) -## How to create a RamAirKite +## How to create a SymbolicAWEModel The following code is a minimal example that shows how to create a ram air kite struct: ```julia using KiteModels @@ -91,5 +91,5 @@ using KiteModels # Initialize model set = load_settings("system_ram.yaml") -rak = RamAirKite(set) +sam = SymbolicAWEModel(set) ``` \ No newline at end of file diff --git a/docs/src/functions.md b/docs/src/functions.md index 7c0d712e3..0d34cda90 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -18,7 +18,7 @@ s = KPS4(KCU(set)) Or, if you want to use the ram-air kite model: ```julia set = load_settings("system_ram.yaml") -s = RamAirKite(set) +s = SymbolicAWEModel(set) ``` Functions with an "!" as last character of the function name modify one of more of their parameters, in this context mostly the variable s. diff --git a/docs/src/index.md b/docs/src/index.md index ac01bbcac..7ea5617bd 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -54,7 +54,7 @@ include("examples/menu.jl") - a new 5-point model based on ModelingToolkit (MTK) is in development; this will allow to create linearized models around any operation point and to do analysis in the frequency domain. #### April 2025 -- a new model `RamAirKite` was contributed, based on the package [VortexStepMethod](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl) +- a new model `SymbolicAWEModel` was contributed, based on the package [VortexStepMethod](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl) #### November 2024 - the four point kite model KPS4 was extended to include aerodynamic damping of pitch oscillations; for this purpose, the parameters `cmq` and `cord_length` must be defined in `settings.yaml` @@ -69,7 +69,7 @@ include("examples/menu.jl") - the documentation was improved ## Provides -The type [`AbstractKiteModel`](@ref) with the implementation [`KPS3`](@ref), [`KPS4`](@ref) and [`RamAirKite`](@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_sim!`](@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 type [`AbstractKiteModel`](@ref) with the implementation [`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_sim!`](@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 index 01a6e35e6..039f77b78 100644 --- a/docs/src/ram_air_kite.md +++ b/docs/src/ram_air_kite.md @@ -2,7 +2,7 @@ CurrentModule = KiteModels ``` ## Introduction -The [`RamAirKite`](@ref) is based on ModelingToolkit, which allows to define the differential algebraic equations in symbolic form. It does not use a KCU. Instead, the kite is controlled from the ground (e.g. a ship) using three winches and four tethers. +The [`SymbolicAWEModel`](@ref) is based on ModelingToolkit, which allows to define the differential algebraic equations in symbolic form. It does not use a KCU. Instead, the kite is controlled from the ground (e.g. a ship) using three winches and four tethers. ## Private enumerations ```@docs @@ -17,13 +17,13 @@ Pulley Segment Tether Winch -KitePointGroup -PointMassSystem +Group +SystemStructure ``` ## Private functions ```@docs -diff_eqs! +wing_eqs! reinit! scalar_eqs! linear_vsm_eqs! diff --git a/docs/src/test_plan.md b/docs/src/test_plan.md index c4caa748f..0cb41f4b9 100644 --- a/docs/src/test_plan.md +++ b/docs/src/test_plan.md @@ -31,7 +31,7 @@ Expected output: julia> include("examples/ram_air_kite.jl") [ Info: Loading packages Time elapsed: 2.863292943 s -[ Info: Creating wing, aero, vsm_solver, point_system and s: +[ Info: Creating wing, aero, vsm_solver, system_structure and s: Time elapsed: 4.51880941 s [ Info: Initialized integrator in 6.000656365 seconds [ Info: System initialized at: diff --git a/docs/src/tutorial_system_structure.md b/docs/src/tutorial_system_structure.md new file mode 100644 index 000000000..557699022 --- /dev/null +++ b/docs/src/tutorial_system_structure.md @@ -0,0 +1,62 @@ +```@meta +CurrentModule = KiteModels +``` +# Custom SystemStructure and SymbolicAWESystem + +A custom `SystemStructure` can be used to create models of kite power systems of almost any configuration. +- custom amount of tethers +- custom bridle configurations +- quasi-static or dynamic point masses +- different amounts of stiffness, damping and diameter on different tether segments + +## Creating a simple tether + +We start by loading the necessary packages and defining settings and parameters. + +```@example 1 +using KiteModels, VortexStepMethod, ControlPlots + +set = se("system_ram.yaml") +set.segments = 20 +dynamics_type = DYNAMIC +``` +Then, we define vectors of the system structure types we are going to use. For this simple example we only need points and segments. +```@example 1 +points = Point[] +segments = Segment[] + +points = push!(points, Point(1, [0.0, 0.0, set.l_tether], 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 rigid wing body. + +Now we can add `DYNAMIC` points and connect them to eachother with segments. `BRIDLE` segments don't need to have a tether, because they have a constant unstretched length. +```@example 1 +segment_idxs = Int[] +for i in 1:set.segments + global points, segments + point_idx = i+1 + pos = [0.0, 0.0, set.l_tether] - set.l_tether / set.segments * [0.0, 0.0, i] + 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 +``` +From these arrays of points and segments we can create a [`SystemStructure`](@ref), which can be plotted in 2d to quickly investigate if the model is correct. +```@example 1 +system_structure = SystemStructure("tether"; points, segments) +plot(system_structure, 0.0) +plt.gcf() +``` + +If the system looks good, we can easily model it, by first creating a [`SymbolicAWEModel`](@ref), initializing it and stepping through time. +```@example 1 +model = SymbolicAWEModel(set, system_structure) + +init_sim!(model; remake=false) +for i in 1:100 + plot(model, i/set.sample_freq) + next_step!(model; dt=1/set.sample_freq) +end +plt.gcf() +``` diff --git a/docs/src/types.md b/docs/src/types.md index e2c097fb0..2cfa94078 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -13,11 +13,11 @@ AbstractKiteModel AKM ``` -## Struct KPS3 and KPS4 and RamAirKite +## Struct KPS3 and KPS4 and SymbolicAWEModel ```@docs KPS3 KPS4 -RamAirKite +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 20f3b1b48..f32468bd7 100644 --- a/examples/lin_ram_model.jl +++ b/examples/lin_ram_model.jl @@ -3,7 +3,7 @@ #= This example demonstrates linearized model accuracy by comparing: -1. Nonlinear RamAirKite model simulation +1. Nonlinear SymbolicAWEModel model simulation 2. Linearized state-space model simulation Both models start from the same operating point and are subjected @@ -49,18 +49,18 @@ set_values = [-50.0, 0.0, 0.0] # Set values of the torques of the three winches set.quasi_static = false set.physical_model = "simple_ram" -@info "Creating RamAirKite model..." -s = RamAirKite(set) +@info "Creating SymbolicAWEModel model..." +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) +lin_outputs = @variables heading(t)[1] # Initialize at elevation with linearization outputs -s.point_system.winches[2].tether_length += 0.2 -s.point_system.winches[3].tether_length += 0.2 +s.system_structure.winches[2].tether_length += 0.2 +s.system_structure.winches[3].tether_length += 0.2 KiteModels.init_sim!(s; remake=false, reload=true, @@ -68,7 +68,7 @@ KiteModels.init_sim!(s; ) sys = s.sys -@show rad2deg(s.integrator[sys.elevation]) +@show rad2deg(s.integrator[sys.elevation[1]]) @info "System initialized at:" diff --git a/examples/plotting.jl b/examples/plotting.jl deleted file mode 100644 index 0b1d356f2..000000000 --- a/examples/plotting.jl +++ /dev/null @@ -1,32 +0,0 @@ -# Copyright (c) 2022, 2024 Uwe Fechner -# SPDX-License-Identifier: MIT - -function plot(sys::PointMassSystem, reltime; kite_pos=nothing, e_z=zeros(3), zoom=false, front=false) - pos = [sys.points[i].pos_w for i in eachindex(sys.points)] - !isnothing(kite_pos) && (pos = [pos..., kite_pos]) - seg = [[sys.segments[i].points[1], sys.segments[i].points[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 = (0, 60) - ylim = (0, 60) - elseif !zoom && front - xlim = (-30, 30) - ylim = (0, 60) - end - ControlPlots.plot2d(pos, seg, reltime; zoom, front, xlim, ylim, dz_zoom=0.6) -end - -function plot(s::RamAirKite, reltime; kwargs...) - pos = s.integrator[s.sys.pos] - kite_pos = s.integrator[s.sys.kite_pos] - for (i, point) in enumerate(s.point_system.points) - point.pos_w .= pos[:, i] - end - e_z = s.integrator[s.sys.e_z] - plot(s.point_system, reltime; kite_pos, e_z, kwargs...) -end \ No newline at end of file diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 89aca7928..ddaadc0e7 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -36,45 +36,45 @@ set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. set.quasi_static = false set.physical_model = SIMPLE ? "simple_ram" : "ram" -@info "Creating wing, aero, vsm_solver, point_system and s:" -s = RamAirKite(set) -s.set.abs_tol = 1e-2 -s.set.rel_tol = 1e-2 +@info "Creating wing, aero, vsm_solver, system_structure and s:" +sam = SymbolicAWEModel(set) +sam.set.abs_tol = 1e-2 +sam.set.rel_tol = 1e-2 toc() -# init_Q_b_w, R_b_w = initial_orient(s.set) -# init_kite_pos = init!(s.point_system, s.set, R_b_w, init_Q_b_w) -# plot(s.point_system, 0.0; zoom=false, front=false) +# init_Q_b_w, R_b_w = KiteModelsam.initial_orient(sam.set) +# init_kite_pos = init!(sam.system_structure, sam.set, R_b_w, init_Q_b_w) +# plot(sam.system_structure, 0.0; zoom=false, front=false) # Initialize at elevation -s.point_system.winches[2].tether_length += 0.2 -s.point_system.winches[3].tether_length += 0.2 -init_sim!(s; remake=false, reload=true) -sys = s.sys +sam.system_structure.winches[2].tether_length += 0.2 +sam.system_structure.winches[3].tether_length += 0.2 +init_sim!(sam; remake=false, reload=true) +sys = sam.sys @info "System initialized at:" toc() # Stabilize system -find_steady_state!(s) +find_steady_state!(sam) -logger = Logger(length(s.point_system.points), steps) -sys_state = SysState(s) +logger = Logger(length(sam.system_structure.points), steps) +sys_state = SysState(sam) t = 0.0 runtime = 0.0 integ_runtime = 0.0 bias = set.quasi_static ? 0.45 : 0.35 -t0 = s.integrator.t +t0 = sam.integrator.t try while t < total_time local steering global t, set_values, runtime, integ_runtime - PLOT && plot(s, t; zoom=false, front=false) + PLOT && plot(sam, t; zoom=false, front=false) # Calculate steering inputs based on cosine wave steering = steering_magnitude * cos(2π * steering_freq * t + bias) - set_values = -s.set.drum_radius .* s.integrator[sys.winch_force] + 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 @@ -82,18 +82,18 @@ try end # Step simulation - steptime = @elapsed (t_new, integ_steptime) = next_step!(s; set_values, dt, vsm_interval=vsm_interval) + steptime = @elapsed (t_new, integ_steptime) = next_step!(sam; set_values, dt, vsm_interval=vsm_interval) t = t_new - t0 # Adjust for initial stabilization time # Track performance after initial transient if (t > total_time/2) runtime += steptime integ_runtime += integ_steptime - s.integrator.ps[sys.twist_damp] = 10 + sam.integrator.ps[sys.twist_damp] = 10 end # Log state variables - update_sys_state!(sys_state, s) + update_sys_state!(sys_state, sam) sys_state.time = t log!(logger, sys_state) end diff --git a/examples/tether.jl b/examples/tether.jl new file mode 100644 index 000000000..f230cb81c --- /dev/null +++ b/examples/tether.jl @@ -0,0 +1,36 @@ +# 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, [0.0, 0.0, set.l_tether], 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, set.l_tether] - set.l_tether / set.segments * [0.0, 0.0, i] + 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 + +system_structure = SystemStructure("tether"; points, segments) +plot(system_structure, 0.0) + +model = SymbolicAWEModel(set, system_structure) + +init_sim!(model; remake=false) +for i in 1:100 + plot(model, i/set.sample_freq) + next_step!(model) +end \ No newline at end of file diff --git a/ext/KiteModelsControlPlotsExt.jl b/ext/KiteModelsControlPlotsExt.jl index cc935a810..0e82a0a5e 100644 --- a/ext/KiteModelsControlPlotsExt.jl +++ b/ext/KiteModelsControlPlotsExt.jl @@ -7,10 +7,10 @@ using ControlPlots, KiteModels export plot -function ControlPlots.plot(sys::PointMassSystem, reltime; kite_pos=nothing, e_z=zeros(3), zoom=false, front=false) +function ControlPlots.plot(sys::SystemStructure, reltime; l_tether=50.0, wing_pos=nothing, e_z=zeros(3), zoom=false, front=false) pos = [sys.points[i].pos_w for i in eachindex(sys.points)] - !isnothing(kite_pos) && (pos = [pos..., kite_pos]) - seg = [[sys.segments[i].points[1], sys.segments[i].points[2]] for i in eachindex(sys.segments)] + !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) @@ -18,23 +18,30 @@ function ControlPlots.plot(sys::PointMassSystem, reltime; kite_pos=nothing, e_z= xlim = (pos[end][2] - 6, pos[end][2]+6) ylim = (pos[end][3] - 10, pos[end][3]+2) elseif !zoom && !front - xlim = (0, 60) - ylim = (0, 60) + xlim = (-5, l_tether+10) + ylim = (-5, l_tether+10) elseif !zoom && front - xlim = (-30, 30) - ylim = (0, 60) + 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) end -function ControlPlots.plot(s::RamAirKite, reltime; kwargs...) +function ControlPlots.plot(s::SymbolicAWEModel, reltime; kwargs...) + wings = s.system_structure.wings pos = s.integrator[s.sys.pos] - kite_pos = s.integrator[s.sys.kite_pos] - for (i, point) in enumerate(s.point_system.points) + 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.system_structure.points) point.pos_w .= pos[:, i] end - e_z = s.integrator[s.sys.e_z] - plot(s.point_system, reltime; kite_pos, e_z, kwargs...) + plot(s.system_structure, reltime; s.set.l_tether, wing_pos, e_z, kwargs...) end end \ No newline at end of file diff --git a/mwes/mwe_10.jl b/mwes/mwe_10.jl index 172be07af..df16f645e 100644 --- a/mwes/mwe_10.jl +++ b/mwes/mwe_10.jl @@ -4,7 +4,7 @@ using KiteModels -kps4_3L::RamAirKite = RamAirKite(KCU(se())) +kps4_3L::SymbolicAWEModel = SymbolicAWEModel(KCU(se())) integrator = KiteModels.init_sim!(kps4_3L; stiffness_factor=0.035, prn=false) nothing \ No newline at end of file diff --git a/mwes/mwe_27.jl b/mwes/mwe_27.jl index 864957012..399e4d4e1 100644 --- a/mwes/mwe_27.jl +++ b/mwes/mwe_27.jl @@ -4,5 +4,5 @@ using KiteModels set = load_settings("system_ram.yaml") -rak = RamAirKite(set) -KiteModels.init_sim!(rak; remake=false, reload=false) \ No newline at end of file +sam = SymbolicAWEModel(set) +KiteModels.init_sim!(sam; remake=false, reload=false) \ No newline at end of file diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 62d4dac8a..ab0d0669f 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -36,7 +36,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D using ADTypes: AutoFiniteDiff import ModelingToolkit.SciMLBase: successful_retcode -export KPS3, KPS4, RamAirKite, KVec3, SimFloat, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG # constants and types +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_sim!, init!, reinit!, next_step!, init_pos_vel # high level workers @@ -45,15 +45,17 @@ 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_point_system, create_simple_ram_point_system +export create_ram_system_structure, create_simple_ram_system_structure import LinearAlgebra: norm +export SystemStructure, Point, Group, Segment, Pulley, Tether, Winch, Wing +export DynamicsType, DYNAMIC, QUASI_STATIC, WING, STATIC +export SegmentType, POWER, STEERING, BRIDLE set_zero_subnormals(true) # required to avoid drastic slow down on Intel CPUs when numbers become very small # Constants const G_EARTH = 9.81 # gravitational acceleration const BRIDLE_DRAG = 1.1 # should probably be removed -const SYS_3L = "system_3l.yaml" # default system project for the 3L model # Type definitions """ diff --git a/src/mtk_model.jl b/src/mtk_model.jl index c3c34e660..a68717330 100644 --- a/src/mtk_model.jl +++ b/src/mtk_model.jl @@ -2,7 +2,7 @@ # SPDX-License-Identifier: MPL-2.0 # ==================== mtk model functions ================================================ -# Implementation of the ram air kite model using ModelingToolkit.jl +# 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 @@ -15,8 +15,8 @@ function calc_heading(e_x) return atan(e_x[2], e_x[1]) end -function calc_angle_of_attack(va_kite_b) - return atan(va_kite_b[3], va_kite_b[1]) +function calc_angle_of_attack(va_wing_b) + return atan(va_wing_b[3], va_wing_b[1]) end function sym_normalize(vec) @@ -68,18 +68,18 @@ end """ force_eqs!(s, system, eqs, defaults, guesses; kwargs...) -Generate the force equations for the kite system including spring forces, drag forces, +Generate the force equations for the wing system including spring forces, drag forces, pulley dynamics and winch forces. # Arguments -- `s::RamAirKite`: The kite system state -- `system::PointMassSystem`: The point mass representation +- `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 -- `kite_pos`: Kite position vector -- `kite_vel`: Kite velocity vector +- `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 @@ -90,25 +90,24 @@ Tuple containing: - Updated equations - Updated defaults - Updated guesses -- Tether forces on kite -- Tether moments on kite +- Tether forces on wing +- Tether moments on wing """ function force_eqs!(s, system, eqs, defaults, guesses; - R_b_w, kite_pos, kite_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff) + 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 - points, groups, segments, pulleys, tethers, winches = - system.points, system.groups, system.segments, system.pulleys, system.tethers, system.winches + @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) @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_kite_force(t)[1:3, eachindex(points)] - tether_kite_moment(t)[1:3, eachindex(points)] tether_r(t)[1:3, eachindex(points)] point_mass(t)[eachindex(points)] @@ -121,9 +120,9 @@ function force_eqs!(s, system, eqs, defaults, guesses; mass = 0.0 in_bridle = false for segment in segments - if point.idx in segment.points + if point.idx in segment.point_idxs mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2 - inverted = segment.points[2] == point.idx + inverted = segment.point_idxs[2] == point.idx if inverted F .-= spring_force_vec[:, segment.idx] else @@ -143,61 +142,51 @@ function force_eqs!(s, system, eqs, defaults, guesses; point_force[:, point.idx] ~ F ] - winch_point = false - for tether in tethers - if point.idx == tether.winch_point - winch_point = true - break - end - end - - if point.type !== KITE - eqs = [ - eqs - tether_kite_force[:, point.idx] ~ zeros(3) - tether_kite_moment[:, point.idx] ~ zeros(3) - ] - end - - if point.type == KITE + if point.type == WING found = 0 - group_idx = 0 - for group in groups - if point.idx in group.points - group_idx = group.idx + group = nothing + for group_ in groups + if point.idx in group_.point_idxs + group = group_ found += 1 end end - !(found in [0,1]) && throw(ArgumentError("Kite point number $(point.idx) is part of $found groups, - and should be part of exactly 0 or 1 groups.")) + !(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 - group = groups[group_idx] + 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 - eqs = [ - eqs - tether_kite_moment[:, point.idx] ~ zeros(3) - ] + pos_b = fixed_pos + cos(twist_angle[group.idx]) * chord_b - sin(twist_angle[group.idx]) * normal else pos_b = point.pos_b eqs = [ eqs - tether_r[:, point.idx] ~ pos[:, point.idx] - kite_pos - tether_kite_moment[:, point.idx] ~ tether_r[:, point.idx] × tether_kite_force[:, point.idx] + 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 - tether_kite_force[:, point.idx] ~ point_force[:, point.idx] - pos[:, point.idx] ~ kite_pos + R_b_w * pos_b + 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 == WINCH + elseif point.type == STATIC eqs = [ eqs pos[:, point.idx] ~ point.pos_w @@ -205,26 +194,30 @@ function force_eqs!(s, system, eqs, defaults, guesses; acc[:, point.idx] ~ zeros(3) ] elseif point.type == DYNAMIC - p = pos[:, point.idx] - n = sym_normalize(kite_pos) - n = n * (p ⋅ n) - r = (p - n) # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation + # 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] - in_bridle * bridle_damp * (vel[:, point.idx] - kite_vel) - acc[:, point.idx] ~ point_force[:, point.idx] / mass + - [0, 0, -G_EARTH] + - ifelse.(stabilize==true, r * norm(measured_ω_z)^2, zeros(3)) # TODO: add other stabilize accelerations + 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] => point.pos_w[j] for j in 1:3] [vel[j, point.idx] => 0 for j in 1:3] ] - elseif point.type == STATIC + elseif point.type == QUASI_STATIC eqs = [ eqs vel[:, point.idx] ~ zeros(3) @@ -238,32 +231,47 @@ function force_eqs!(s, system, eqs, defaults, guesses; [point_force[j, point.idx] => 0 for j in 1:3] ] else - throw(ArgumentError("Unknown point type: $(typeof(point))")) + error("Unknown point type: $(typeof(point))") end end # ==================== GROUPS ==================== # - @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].points)] - tether_moment(t)[eachindex(groups), eachindex(groups[1].points)] + 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.points) + 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 * -z_airf)) + 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 @@ -290,7 +298,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; free_twist_angle[group.idx] => 0 twist_ω[group.idx] => 0 ] - elseif group.type == STATIC + elseif group.type == QUASI_STATIC eqs = [ eqs twist_ω[group.idx] ~ 0 @@ -302,7 +310,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; twist_angle[group.idx] => 0 ] else - throw(ArgumentError("Wrong group type.")) + error("Wrong group type.") end end @@ -331,7 +339,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; tether_length(t)[eachindex(winches)] end for segment in segments - p1, p2 = segment.points[1], segment.points[2] + p1, p2 = segment.point_idxs[1], segment.point_idxs[2] if s.set.quasi_static guesses = [ guesses @@ -342,14 +350,14 @@ function force_eqs!(s, system, eqs, defaults, guesses; if segment.type == BRIDLE in_pulley = 0 for pulley in pulleys - if segment.idx == pulley.segments[1] # each bridle segment has to be part of no pulley or one pulley + 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.segments[2] + if segment.idx == pulley.segment_idxs[2] eqs = [ eqs l0[segment.idx] ~ pulley.sum_length - pulley_l0[pulley.idx] @@ -363,34 +371,41 @@ function force_eqs!(s, system, eqs, defaults, guesses; l0[segment.idx] ~ segment.l0 ] end - (in_pulley > 1) && throw(ArgumentError("Bridle segment number $(segment.idx) is part of - $in_pulley pulleys, and should be part of either 0 or 1 pulleys.")) + (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 || segment.type == STEERING in_tether = 0 for tether in tethers - if segment.idx in tether.segments # each tether segment has to be part of exactly one tether + 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.tethers + if tether.idx in winch.tether_idxs winch_idx = winch.idx in_winch += 1 end end - (in_winch != 1) && throw(ArgumentError("Tether number $(tether.idx) is part of - $(in_winch) winches, and should be part of exactly 1 winch.")) - - eqs = [ - eqs - l0[segment.idx] ~ tether_length[winch_idx] / length(tether.segments) - ] + !(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) && throw(ArgumentError("Segment number $(segment.idx) is part of - $in_tether tethers, and should be part of exactly 1 tether.")) + (in_tether != 1) && error("Segment number $(segment.idx) is part of + $in_tether tethers, and should be part of exactly 1 tether.") else - throw(ArgumentError("Unknown segment type: $(segment.type)")) + error("Unknown segment type: $(segment.type)") end stiffness_m = s.set.e_tether * (segment.diameter/2)^2 * pi @@ -442,13 +457,13 @@ function force_eqs!(s, system, eqs, defaults, guesses; end @parameters pulley_damp = 1.0 for pulley in pulleys - segment = segments[pulley.segments[1]] + 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.segments[1]] - spring_force[pulley.segments[2]] + 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 @@ -459,10 +474,10 @@ function force_eqs!(s, system, eqs, defaults, guesses; ] defaults = [ defaults - pulley_l0[pulley.idx] => segments[pulley.segments[1]].l0 + pulley_l0[pulley.idx] => segments[pulley.segment_idxs[1]].l0 pulley_vel[pulley.idx] => 0 ] - elseif pulley.type == STATIC + elseif pulley.type == QUASI_STATIC eqs = [ eqs pulley_vel[pulley.idx] ~ 0 @@ -470,10 +485,10 @@ function force_eqs!(s, system, eqs, defaults, guesses; ] guesses = [ guesses - pulley_l0[pulley.idx] => segments[pulley.segments[1]].l0 + pulley_l0[pulley.idx] => segments[pulley.segment_idxs[1]].l0 ] else - throw(ArgumentError("Wrong pulley type")) + error("Wrong pulley type") end end @@ -485,8 +500,19 @@ function force_eqs!(s, system, eqs, defaults, guesses; end for winch in winches F = zero(Num) - for tether_idx in winch.tethers - point_idx = tethers[tether_idx].winch_point + 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 = [ @@ -497,8 +523,8 @@ function force_eqs!(s, system, eqs, defaults, guesses; 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] # TODO: SISO mode to test linear model accuracy, MHE to check if it is better, use stabilize as state estimator - ) # TODO: add brakes and P. This makes SISO system easy + set_values[winch.idx] + ) winch_force[winch.idx] ~ F ] defaults = [ @@ -514,7 +540,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; end for tether in tethers ulen = zero(Num) - for segment_idx in tether.segments + for segment_idx in tether.segment_idxs ulen += len[segment_idx] end eqs = [ @@ -523,51 +549,49 @@ function force_eqs!(s, system, eqs, defaults, guesses; ] end - return eqs, defaults, guesses, tether_kite_force, tether_kite_moment + return eqs, defaults, guesses, tether_wing_force, tether_wing_moment end """ - diff_eqs!(s, eqs, defaults; kwargs...) + wing_eqs!(s, eqs, defaults; kwargs...) -Generate the differential equations for kite dynamics including quaternion kinematics, +Generate the differential equations for wing dynamics including quaternion kinematics, angular velocities and accelerations, and forces/moments. # Arguments -- `s::RamAirKite`: The kite system state +- `s::SymbolicAWEModel`: The wing system state - `eqs`: Current system equations - `defaults`: Default values for variables -- `tether_kite_force`: Forces from tethers on kite -- `tether_kite_moment`: Moments from tethers on kite +- `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 -- `kite_pos`: Kite position vector -- `kite_vel`: Kite velocity vector -- `kite_acc`: Kite acceleration vector +- `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 diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero_force_b, - aero_moment_b, ω_b, α_b, R_b_w, kite_pos, kite_vel, kite_acc, stabilize, fix_nonstiff +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 ) - kite = s.point_system.kite + wings = s.system_structure.wings @variables begin # potential differential variables - kite_acc_b(t)[1:3] - α_b_damped(t)[1:3] - ω_b_stable(t)[1:3] + 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)[1:4] # quaternion orientation of the kite body frame relative to the world frame - Q_vel(t)[1:4] # quaternion rate of change + 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)[1:3] # moment in principal frame - total_tether_kite_force(t)[1:3] - total_tether_kite_moment(t)[1:3] + moment_b(t)[eachindex(wings), 1:3] # moment in principal frame end @parameters ω_damp = 150 @@ -576,67 +600,69 @@ function diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero ω[2] -ω[3] 0 ω[1]; ω[3] ω[2] -ω[1] 0] - I_b = [s.wing.inertia_tensor[1,1], s.wing.inertia_tensor[2,2], s.wing.inertia_tensor[3,3]] - @assert !(s.set.mass ≈ 0) - axis = sym_normalize(kite_pos) - axis_b = R_b_w' * axis - eqs = [ - eqs - [D(Q_b_w[i]) ~ Q_vel[i] for i in 1:4] - [Q_vel[i] ~ 0.5 * sum(Ω(ω_b_stable)[i, j] * Q_b_w[j] for j in 1:4) for i in 1:4] - ω_b_stable ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - ω_b - ω_b ⋅ axis_b * axis_b, - ω_b + 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) ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - α_b_damped - α_b_damped ⋅ axis_b * axis_b, - α_b_damped + 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 ~ [α_b[1], α_b[2] - ω_damp*ω_b[2], α_b[3]] - - [R_b_w[:, i] ~ quaternion_to_rotation_matrix(Q_b_w)[:, i] for i in 1:3] - - α_b[1] ~ (moment_b[1] + (I_b[2] - I_b[3]) * ω_b[2] * ω_b[3]) / I_b[1] - α_b[2] ~ (moment_b[2] + (I_b[3] - I_b[1]) * ω_b[3] * ω_b[1]) / I_b[2] - α_b[3] ~ (moment_b[3] + (I_b[1] - I_b[2]) * ω_b[1] * ω_b[2]) / I_b[3] - total_tether_kite_moment ~ [sum(tether_kite_moment[i, :]) for i in 1:3] - moment_b ~ aero_moment_b + R_b_w' * total_tether_kite_moment - - D(kite_pos) ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - kite_vel ⋅ axis * axis, - kite_vel + α_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(kite_vel) ~ ifelse.(fix_nonstiff==true, - zeros(3), - ifelse.(stabilize==true, - kite_acc ⋅ axis * axis, - kite_acc + D(wing_vel[wing.idx, :]) ~ ifelse.(fix_nonstiff==true, + zeros(3), + ifelse.(stabilize==true, + wing_acc[wing.idx, :] ⋅ axis * axis, + wing_acc[wing.idx, :] + ) ) - ) - kite_acc ~ (total_tether_kite_force + R_b_w * aero_force_b) / s.set.mass - total_tether_kite_force ~ [sum(tether_kite_force[i, :]) for i in 1:3] - ] - defaults = [ - defaults - [Q_b_w[i] => kite.orient[i] for i in 1:4] - [ω_b[i] => kite.angular_vel[i] for i in 1:3] - [kite_pos[i] => kite.pos[i] for i in 1:3] - [kite_vel[i] => kite.vel[i] for i in 1:3] - ] + 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[i] for i in 1:3] + [wing_vel[wing.idx, i] => wing.vel[i] for i in 1:3] + ] + end + return eqs, defaults end -function calc_R_v_w(kite_pos, e_x) - z = sym_normalize(kite_pos) +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] @@ -650,19 +676,19 @@ function calc_R_t_w(elevation, azimuth) end """ - scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω) + 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::RamAirKite`: The kite system state + - `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_kite_b`: Apparent wind velocity in body frame - - `kite_pos`: Kite position vector - - `kite_vel`: Kite velocity vector - - `kite_acc`: Kite acceleration 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: @@ -672,109 +698,116 @@ Generate equations for scalar quantities like elevation, azimuth, heading and co - Course angle - Angular velocities and accelerations """ -function scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b) +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.system_structure @parameters wind_scale_gnd = s.set.v_wind @parameters upwind_dir = deg2rad(s.set.upwind_dir) @variables begin - e_x(t)[1:3] - e_y(t)[1:3] - e_z(t)[1:3] - wind_vel_kite(t)[1:3] - va_kite(t)[1:3] + 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 - e_x ~ R_b_w[:,1] - e_y ~ R_b_w[:,2] - e_z ~ R_b_w[:,3] wind_vec_gnd ~ wind_scale_gnd * rotate_around_z([0, -1, 0], -upwind_dir) - wind_vel_kite ~ AtmosphericModels.calc_wind_factor(s.am, kite_pos[3], s.set.profile_law) * wind_vec_gnd - va_kite ~ wind_vel_kite - kite_vel - va_kite_b ~ R_b_w' * va_kite ] + 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) - turn_rate(t)[1:3] - turn_acc(t)[1:3] - azimuth(t) - azimuth_vel(t) - azimuth_acc(t) - elevation(t) - elevation_vel(t) - elevation_acc(t) - course(t) - x_acc(t) - y_acc(t) - sphere_pos(t)[1:2, 1:2] - sphere_vel(t)[1:2, 1:2] - sphere_acc(t)[1:2, 1:2] - angle_of_attack(t) - simple_twist_angle(t)[1:2] - simple_twist_ω(t)[1:2] - R_v_w(t)[1:3, 1:3] - R_t_w(t)[1:3, 1:3] - distance(t) - distance_vel(t) - distance_acc(t) + 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 - x, y, z = kite_pos - x´, y´, z´ = kite_vel - x´´, y´´, z´´ = kite_acc + 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 = length(twist_angle)÷2 - heading_vec = R_t_w' * R_v_w[:,1] + 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) .~ vec(calc_R_v_w(kite_pos, e_x)) - vec(R_t_w) .~ vec(calc_R_t_w(elevation, azimuth)) - heading ~ atan(heading_vec[2], heading_vec[1]) - turn_rate ~ R_v_w' * (R_b_w * ω_b) # Project angular velocity onto view frame - turn_acc ~ R_v_w' * (R_b_w * α_b) - distance ~ norm(kite_pos) - distance_vel ~ kite_vel ⋅ R_v_w[:,3] - distance_acc ~ kite_acc ⋅ R_v_w[:,3] - - elevation ~ atan(z / x) - # elevation_vel = d/dt(atan(z/x)) = (x*ż' - z*ẋ')/(x^2 + z^2) according to wolframalpha - elevation_vel ~ (x*z´ - z*x´) / - (x^2 + z^2) - elevation_acc ~ ((x^2 + z^2)*(x*z´´ - z*x´´) + 2(z*x´ - x*z´)*(x*x´ + z*z´))/(x^2 + z^2)^2 - azimuth ~ atan(y / x) - # azimuth_vel = d/dt(atan(y/x)) = (-y*x´ + x*y´)/(x^2 + y^2) # TODO: check if correct - azimuth_vel ~ (-y*x´ + x*y´) / - (x^2 + y^2) - azimuth_acc ~ ((x^2 + y^2)*(-y*x´´ + x*y´´) + 2(y*x´ - x*y´)*(x*x´ + y*y´))/(x^2 + y^2)^2 - course ~ atan(-azimuth_vel, elevation_vel) - x_acc ~ kite_acc ⋅ e_x - y_acc ~ kite_acc ⋅ e_y - - angle_of_attack ~ calc_angle_of_attack(va_kite_b) + 0.5twist_angle[half_len] + 0.5twist_angle[half_len+1] - simple_twist_angle[1] ~ sum(twist_angle[1:half_len]) / half_len - simple_twist_angle[2] ~ sum(twist_angle[half_len+1:end]) / half_len - simple_twist_ω[1] ~ sum(twist_ω[1:half_len]) / half_len - simple_twist_ω[2] ~ sum(twist_ω[half_len+1:end]) / half_len - ] + 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] ~ atan(z / x) + # 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] ~ atan(y / x) + # 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_kite_b, ω_b) +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::RamAirKite`: The kite system state +- `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_kite_b`: Apparent wind velocity in body frame +- `va_wing_b`: Apparent wind velocity in body frame - `ω_b`: Angular velocity in body frame # Returns @@ -783,77 +816,89 @@ and moments. The Jacobian is computed using the VSM solver. - 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_kite_b, ω_b) - sol = s.vsm_solver.sol - @assert length(s.point_system.groups) == length(sol.group_moment_dist) - - y_ = [init_va_b; zeros(length(s.point_system.groups)); zeros(3)] - x_ = zeros(3+3+length(s.point_system.groups)) - jac_ = zeros(length(x_), length(y_)) +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.system_structure + 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(y_)] = y_ - last_x[eachindex(x_)] = x_ - vsm_jac[eachindex(x_), eachindex(y_)] = jac_ + 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(y_)] - dy(t)[eachindex(y_)] + y(t)[eachindex(wings), 1:ny] + dy(t)[eachindex(wings), 1:ny] end - eqs = [ - eqs - y ~ [va_kite_b; ω_b; twist_angle] - dy ~ y - last_y - [aero_force_b; aero_moment_b; group_aero_moment] ~ last_x + vsm_jac * dy - ] - - if s.set.quasi_static - guesses = [guesses; [y[i] => y_[i] for i in eachindex(y_)]] + 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::RamAirKite, system::PointMassSystem; init_va_b) +function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b) eqs = [] defaults = Pair{Num, Real}[] guesses = Pair{Num, Real}[] + @unpack wings, groups, winches = system + @parameters begin stabilize = false fix_nonstiff = false end @variables begin # potential differential variables - set_values(t)[eachindex(system.winches)] = zeros(length(system.winches)) - kite_pos(t)[1:3] # xyz pos of kite in world frame - kite_vel(t)[1:3] - kite_acc(t)[1:3] - ω_b(t)[1:3] # turn rate in principal frame - α_b(t)[1:3] + 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)[1:3, 1:3] # rotation of the kite body frame relative to the world frame + 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)[1:3] - aero_moment_b(t)[1:3] - twist_angle(t)[eachindex(system.groups)] - twist_ω(t)[eachindex(system.groups)] - group_aero_moment(t)[eachindex(system.groups)] + 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_kite_b(t)[1:3] + va_wing_b(t)[eachindex(wings), 1:3] end - eqs, defaults, guesses, tether_kite_force, tether_kite_moment = + eqs, defaults, guesses, tether_wing_force, tether_wing_moment = force_eqs!(s, system, eqs, defaults, guesses; - R_b_w, kite_pos, kite_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_kite_b, ω_b) - eqs, defaults = diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero_force_b, aero_moment_b, - ω_b, α_b, R_b_w, kite_pos, kite_vel, kite_acc, stabilize, fix_nonstiff) - eqs = scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b) + 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 @@ -872,7 +917,7 @@ function create_sys!(s::RamAirKite, system::PointMassSystem; init_va_b) # 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_system.groups)] + # [twist_angle[i] ~ clamp(twist_angle[i], -π/2, π/2) for i in eachindex(s.point_groups)] # ] # ] @@ -882,7 +927,7 @@ function create_sys!(s::RamAirKite, system::PointMassSystem; init_va_b) defaults = [ defaults - [set_values[i] => [-50.0, -1.0, -1.0][i] for i in 1:3] + [set_values[i] => [-50.0, -1.0, -1.0][i] for i in eachindex(winches)] ] s.defaults = defaults diff --git a/src/point_mass_system.jl b/src/point_mass_system.jl index e391f29ae..9aeb3f7ee 100644 --- a/src/point_mass_system.jl +++ b/src/point_mass_system.jl @@ -28,21 +28,21 @@ Type of segment. end """ - DynamicsType `DYNAMIC` `STATIC` `KITE` `WINCH` + DynamicsType `DYNAMIC` `QUASI_STATIC` `WING` `STATIC` Enumeration of the models that are attached to a point. # Elements - DYNAMIC: Belongs to a dynamic tether model -- STATIC: Belongs to a static tether model -- KITE: Rigid body -- WINCH: Winch +- 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 - KITE - WINCH end """ @@ -54,38 +54,42 @@ $(TYPEDFIELDS) """ mutable struct Point idx::Int16 - pos_b::KVec3 # pos relative to kite COM in body frame + wing_idx::Int16 # idx of wing used for initial orientation + 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 -function Point(idx, pos_b, type, vel_w=zeros(KVec3)) - Point(idx, pos_b, copy(pos_b), vel_w, type) +function Point(idx, pos_b, type; vel_w=zeros(KVec3), wing_idx=1) + Point(idx, wing_idx, pos_b, copy(pos_b), vel_w, type) end """ - struct KitePointGroup + struct Group Set of bridle lines that share the same twist angle and trailing edge angle. $(TYPEDFIELDS) """ -mutable struct KitePointGroup +mutable struct Group idx::Int16 - points::Vector{Int16} - le_pos::KVec3 # point which the group rotates around under kite deformation - chord::KVec3 # chord vector in body frame which the group rotates around under kite deformation - y_airf::KVec3 # spanwise vector in local panel frame which the group rotates around under kite deformation + point_idxs::Vector{Int16} + le_pos::KVec3 # point which the group rotates around under wing deformation + chord::KVec3 # chord vector in body frame which the group rotates around under wing deformation + y_airf::KVec3 # spanwise vector in local panel frame which the group rotates around under wing deformation type::DynamicsType moment_frac::SimFloat twist::SimFloat twist_vel::SimFloat end -function KitePointGroup(idx, points, wing, gamma, type, moment_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 - y_airf = normalize([wing.le_interp[i](gamma-0.01) for i in 1:3] - le_pos) - KitePointGroup(idx, points, le_pos, chord, y_airf, type, moment_frac, 0.0, 0.0) +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 +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 """ @@ -97,14 +101,14 @@ $(TYPEDFIELDS) """ mutable struct Segment idx::Int16 - points::Tuple{Int16, Int16} + point_idxs::Tuple{Int16, Int16} type::SegmentType l0::SimFloat compression_frac::SimFloat diameter::SimFloat end -function Segment(idx, points, type, l0=zero(SimFloat), compression_frac=0.1) - Segment(idx, points, type, l0, compression_frac, zero(SimFloat)) +function Segment(idx, point_idxs, type; l0=zero(SimFloat), compression_frac=0.1) + Segment(idx, point_idxs, type, l0, compression_frac, zero(SimFloat)) end """ @@ -116,13 +120,13 @@ $(TYPEDFIELDS) """ mutable struct Pulley idx::Int16 - segments::Tuple{Int16, Int16} + segment_idxs::Tuple{Int16, Int16} type::DynamicsType sum_length::SimFloat length::SimFloat vel::SimFloat - function Pulley(idx, segments, type) - new(idx, segments, type, 0.0, 0.0, 0.0) + function Pulley(idx, segment_idxs, type) + new(idx, segment_idxs, type, 0.0, 0.0, 0.0) end end @@ -135,8 +139,7 @@ $(TYPEDFIELDS) """ struct Tether idx::Int16 - segments::Vector{Int16} - winch_point::Int16 + segment_idxs::Vector{Int16} end """ @@ -149,50 +152,66 @@ $(TYPEDFIELDS) mutable struct Winch idx::Int16 model::AbstractWinchModel - tethers::Vector{Int16} + tether_idxs::Vector{Int16} tether_length::SimFloat tether_vel::SimFloat - function Winch(idx, model, tethers, tether_length, tether_vel=0.0) - new(idx, model, tethers, tether_length, tether_vel) + function Winch(idx, model, tether_idxs, tether_length, tether_vel=0.0) + new(idx, model, tether_idxs, tether_length, tether_vel) end end -@with_kw struct Kite +@with_kw struct Wing + idx::Int16 + group_idxs::Vector{Int16} orient::KVec4 = zeros(KVec4) angular_vel::KVec3 = zeros(KVec3) pos::KVec3 = zeros(KVec3) vel::KVec3 = zeros(KVec3) end +function Wing(idx, group_idxs) + Wing(; idx, group_idxs) +end """ - struct PointMassSystem + 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: - `points::Vector{Point}`: Point masses representing: - - Kite attachment points + - Wing attachment points - Dynamic bridle/tether points - Fixed ground anchor points -- `groups::Vector{KitePointGroup}`: Collections of points that move together, - according to kite deformation (twist and trailing edge deflection) +- `groups::Vector{Group}`: Collections of points that move together, + according to wing deformation (twist and trailing edge deflection) - `segments::Vector{Segment}`: Spring-damper elements between points - `pulleys::Vector{Pulley}`: Elements that redistribute line lengths - `tethers::Vector{Tether}`: Groups of segments with a common unstretched length - `winches::Vector{Winch}`: Ground-based winches that control the tether lengths -See also: [`Point`](@ref), [`Segment`](@ref), [`KitePointGroup`](@ref), [`Pulley`](@ref) +See also: [`Point`](@ref), [`Segment`](@ref), [`Group`](@ref), [`Pulley`](@ref) """ -struct PointMassSystem +struct SystemStructure name::String points::Vector{Point} - groups::Vector{KitePointGroup} + groups::Vector{Group} segments::Vector{Segment} pulleys::Vector{Pulley} tethers::Vector{Tether} winches::Vector{Winch} - kite::Kite - function PointMassSystem(name, points, groups, segments, pulleys, tethers, winches, kite=Kite()) + wings::Vector{Wing} + y::Array{Float64, 2} + x::Array{Float64, 2} + jac::Array{Float64, 3} + function SystemStructure(name; + points=Point[], + groups=Group[], + segments=Segment[], + pulleys=Pulley[], + tethers=Tether[], + winches=Winch[], + wings=Wing[] + ) for (i, point) in enumerate(points) @assert point.idx == i end @@ -211,18 +230,28 @@ struct PointMassSystem for (i, winch) in enumerate(winches) @assert winch.idx == i end - return new(name, points, groups, segments, pulleys, tethers, winches, kite) + 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) + return new(name, points, groups, segments, pulleys, tethers, winches, wings, y, x, jac) end end -function PointMassSystem(set::Settings, wing::RamAirWing) +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_point_system(set, wing) + return create_ram_system_structure(set, wing) elseif set.physical_model == "simple_ram" - return create_simple_ram_point_system(set, wing) + return create_simple_ram_system_structure(set, wing) else throw(ArgumentError("Undefined physical model")) end @@ -235,30 +264,28 @@ function calc_pos(wing::RamAirWing, gamma, frac) return pos end -function create_tether(idx, set, points, segments, tethers, attach_point, type, dynamics_type, z=[0,0,1]) - segment_idxs = Int16[] +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_b, set.l_tether, z) dir = winch_pos - attach_point.pos_b + segment_idxs = Int16[] for i in 1:set.segments frac = i / set.segments pos = attach_point.pos_b + frac * dir - i_pnt = length(points) # last point idx - i_seg = length(segments) # last segment idx + point_idx = length(points)+1 # last point idx + segment_idx = length(segments)+1 # last segment idx if i == 1 - points = [points; Point(1+i_pnt, pos, dynamics_type)] - segments = [segments; Segment(1+i_seg, (attach_point.idx, 1+i_pnt), type)] + 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(1+i_pnt, pos, WINCH)] - segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] + points = [points; Point(point_idx, pos, STATIC)] + segments = [segments; Segment(segment_idx, (point_idx-1, point_idx), type)] else - points = [points; Point(1+i_pnt, pos, dynamics_type)] - segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] + points = [points; Point(point_idx, pos, dynamics_type)] + segments = [segments; Segment(segment_idx, (point_idx-1, point_idx), type)] end - push!(segment_idxs, 1+i_seg) - i_pnt = length(points) + push!(segment_idxs, segment_idx) end - winch_point_idx = points[end].idx - tethers = [tethers; Tether(idx, segment_idxs, winch_point_idx)] + tethers = [tethers; Tether(tether_idx, segment_idxs)] return points, segments, tethers, tethers[end].idx end @@ -280,21 +307,22 @@ function find_axis_point(P, l, v=[0,0,1]) return [t * v[1], t * v[2], t * v[3]] end -function create_ram_point_system(set::Settings, wing::RamAirWing) +function create_ram_system_structure(set::Settings, vsm_wing::RamAirWing) points = Point[] - groups = KitePointGroup[] + groups = Group[] segments = Segment[] pulleys = Pulley[] tethers = Tether[] winches = Winch[] + wings = Wing[] attach_points = Point[] - bridle_top_left = [cad_to_body_frame(wing, set.top_bridle_points[i]) for i in eachindex(set.top_bridle_points)] + 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 ? STATIC : DYNAMIC - z = wing.R_cad_body[:,3] + 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 @@ -302,21 +330,21 @@ function create_ram_point_system(set::Settings, wing::RamAirWing) i_pul = length(pulleys) # last pulley idx i_grp = length(groups) # last group idx - # ==================== CREATE DEFORMING KITE GROUPS ==================== # + # ==================== CREATE DEFORMING WING GROUPS ==================== # points = [ points - Point(1+i_pnt, calc_pos(wing, gammas[1], set.bridle_fracs[1]), KITE) - Point(2+i_pnt, calc_pos(wing, gammas[1], set.bridle_fracs[3]), KITE) - Point(3+i_pnt, calc_pos(wing, gammas[1], set.bridle_fracs[4]), KITE) + 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(wing, gammas[2], set.bridle_fracs[1]), KITE) - Point(5+i_pnt, calc_pos(wing, gammas[2], set.bridle_fracs[3]), KITE) - Point(6+i_pnt, calc_pos(wing, gammas[2], set.bridle_fracs[4]), KITE) + 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 - KitePointGroup(1+i_grp, [1+i_pnt, 2+i_pnt, 3+i_pnt], wing, gammas[1], DYNAMIC, set.bridle_fracs[2]) - KitePointGroup(2+i_grp, [4+i_pnt, 5+i_pnt, 6+i_pnt], wing, gammas[2], DYNAMIC, set.bridle_fracs[2]) + 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 ==================== # @@ -324,7 +352,7 @@ function create_ram_point_system(set::Settings, wing::RamAirWing) points = [ points Point(7+i_pnt, bridle_top[1], dynamics_type) - Point(8+i_pnt, bridle_top[2], KITE) + 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) @@ -346,17 +374,17 @@ function create_ram_point_system(set::Settings, wing::RamAirWing) 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, 2) - Segment(8+i_seg, (8+i_pnt, 11+i_pnt), BRIDLE, 1) - Segment(9+i_seg, (9+i_pnt, 13+i_pnt), BRIDLE, 2) - Segment(10+i_seg, (10+i_pnt, 15+i_pnt), BRIDLE, 4) + 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, 1) - Segment(12+i_seg, (11+i_pnt, 13+i_pnt), BRIDLE, 1) + 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, 2) - Segment(14+i_seg, (13+i_pnt, 14+i_pnt), BRIDLE, 2) - Segment(15+i_seg, (13+i_pnt, 15+i_pnt), BRIDLE, 2) + 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 @@ -368,7 +396,7 @@ function create_ram_point_system(set::Settings, wing::RamAirWing) return nothing end - gammas = [-3/4, -1/4, 1/4, 3/4] * wing.gamma_tip + 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]]) @@ -381,44 +409,46 @@ function create_ram_point_system(set::Settings, wing::RamAirWing) winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)] winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)] - return PointMassSystem(set.physical_model, points, groups, segments, pulleys, tethers, winches) + wings = [Wing(1, [1,2,3,4])] + + return SystemStructure(set.physical_model; points, groups, segments, pulleys, tethers, winches, wings) end -function create_simple_ram_point_system(set::Settings, wing::RamAirWing) +function create_simple_ram_system_structure(set::Settings, wing::RamAirWing) points = Point[] - groups = KitePointGroup[] + groups = Group[] segments = Segment[] pulleys = Pulley[] tethers = Tether[] winches = Winch[] - dynamics_type = set.quasi_static ? STATIC : DYNAMIC + 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 KITE GROUPS ==================== # + # ==================== CREATE DEFORMING WING GROUPS ==================== # points = [ points - Point(1, calc_pos(wing, gammas[1], set.bridle_fracs[4]), KITE) - Point(2, calc_pos(wing, gammas[2], set.bridle_fracs[4]), KITE) - Point(3, calc_pos(wing, gammas[3], set.bridle_fracs[4]), KITE) - Point(4, calc_pos(wing, gammas[4], set.bridle_fracs[4]), KITE) + 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 - KitePointGroup(1, [1], wing, gammas[1], DYNAMIC, set.bridle_fracs[2]) - KitePointGroup(2, [2], wing, gammas[2], DYNAMIC, set.bridle_fracs[2]) - KitePointGroup(3, [3], wing, gammas[3], DYNAMIC, set.bridle_fracs[2]) - KitePointGroup(4, [4], wing, gammas[4], DYNAMIC, set.bridle_fracs[2]) + 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], KITE) + Point(5, bridle_top_left[2], WING) Point(6, bridle_top_left[4], dynamics_type) - Point(7, bridle_top_right[2], KITE) + Point(7, bridle_top_right[2], WING) Point(8, bridle_top_right[4], dynamics_type) ] @@ -439,24 +469,26 @@ function create_simple_ram_point_system(set::Settings, wing::RamAirWing) winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)] winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)] - return PointMassSystem(set.physical_model, points, groups, segments, pulleys, tethers, winches) + wings = [Wing(1, [1,2,3,4])] + + return SystemStructure(set.physical_model; points, groups, segments, pulleys, tethers, winches, wings) end -function init!(system::PointMassSystem, set::Settings, R_b_w, Q_b_w) - @unpack points, groups, segments, pulleys, tethers, winches, kite = system +function init!(system::SystemStructure, set::Settings, R_b_w, Q_b_w) + @unpack points, groups, segments, pulleys, tethers, winches, wings = system for segment in segments (segment.type === BRIDLE) && (segment.diameter = 0.001set.bridle_tether_diameter) (segment.type === POWER) && (segment.diameter = 0.001set.power_tether_diameter) (segment.type === STEERING) && (segment.diameter = 0.001set.steering_tether_diameter) - (segment.l0 ≈ 0) && (segment.l0 = norm(points[segment.points[1]].pos_b - points[segment.points[2]].pos_b) * 0.9999) + (segment.l0 ≈ 0) && (segment.l0 = norm(points[segment.point_idxs[1]].pos_b - points[segment.point_idxs[2]].pos_b) * 0.9999) @assert (0 < segment.diameter < 1) @assert (segment.l0 > 0) end for pulley in pulleys - segment1, segment2 = segments[pulley.segments[1]], segments[pulley.segments[2]] + 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 @@ -469,7 +501,7 @@ function init!(system::PointMassSystem, set::Settings, R_b_w, Q_b_w) @assert !(winch.tether_length ≈ 0) end - first_moment_frac = groups[1].moment_frac + (length(groups) > 0) && (first_moment_frac = groups[1].moment_frac) for group in groups group.twist = 0.0 group.twist_vel = 0.0 @@ -483,24 +515,19 @@ function init!(system::PointMassSystem, set::Settings, R_b_w, Q_b_w) end end for point in points - point.pos_w .= R_b_w * (point.pos_b .- min_point) + if point.wing_idx == 0 + R = I(3) + else + R = R_b_w[point.wing_idx, :, :] + end + point.pos_w .= R * (point.pos_b .- min_point) point.vel_w .= 0.0 end - kite.pos .= R_b_w * -min_point - kite.orient .= Q_b_w - kite.vel .= 0.0 - kite.angular_vel .= 0.0 + for wing in wings + wing.pos .= R_b_w[wing.idx, :, :] * -min_point + wing.orient .= Q_b_w[wing.idx, :] + wing.vel .= 0.0 + wing.angular_vel .= 0.0 + end return nothing end - -function initial_orient(set::Settings, R_cad_body=I(3)) - 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 = R_cad_body' * hcat(x, z × x, z) - Q_b_w = rotation_matrix_to_quaternion(R_b_w) - return Q_b_w, R_b_w -end diff --git a/src/precompile.jl b/src/precompile.jl index 70904d35a..1be397137 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -86,7 +86,7 @@ end 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" - s = RamAirKite(set) + s = SymbolicAWEModel(set) # Initialize at elevation KiteModels.init_sim!(s; prn=false, precompile=true) @@ -95,4 +95,4 @@ end end nothing end -end \ No newline at end of file +end \ No newline at end of file diff --git a/src/ram_air_kite.jl b/src/ram_air_kite.jl index 7b4748628..0f0be130c 100644 --- a/src/ram_air_kite.jl +++ b/src/ram_air_kite.jl @@ -2,7 +2,7 @@ # SPDX-License-Identifier: MIT """ - mutable struct RamAirKite{S, V, P} <: AbstractKiteModel + 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 @@ -14,17 +14,17 @@ use the input and output functions instead. $(TYPEDFIELDS) """ -@with_kw mutable struct RamAirKite{S, V, P} <: AbstractKiteModel +@with_kw mutable struct SymbolicAWEModel{S, V, P} <: AbstractKiteModel "Reference to the settings struct" set::Settings "Reference to the geometric wing model" - wing::VortexStepMethod.RamAirWing + vsm_wings::Vector{VortexStepMethod.RamAirWing} "Reference to the aerodynamic wing model" - aero::VortexStepMethod.BodyAerodynamics + vsm_aeros::Vector{VortexStepMethod.BodyAerodynamics} "Reference to the VSM aerodynamics solver" - vsm_solver::VortexStepMethod.Solver + vsm_solvers::Vector{VortexStepMethod.Solver} "Reference to the point mass system with points, segments, pulleys and tethers" - point_system::PointMassSystem + system_structure::SystemStructure "Reference to the atmospheric model as implemented in the package AtmosphericModels" am::AtmosphericModel = AtmosphericModel() "relative start time of the current time interval" @@ -44,67 +44,75 @@ $(TYPEDFIELDS) defaults::Vector{Pair{Num, Real}} = Pair{Num, Real}[] guesses::Vector{Pair{Num, Real}} = Pair{Num, Real}[] - set_set_values::Function = () -> nothing - set_wind_dir::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 + set_set_values::Function = (_, _) -> nothing + set_wind_dir::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_state::Function = () -> nothing - get_y::Function = () -> nothing - get_unstretched_length::Function = () -> nothing - get_tether_length::Function = () -> nothing - get_kite_pos::Function = () -> nothing - get_winch_force::Function = () -> nothing - get_spring_force::Function = () -> nothing - get_stabilize::Function = () -> nothing - get_pos::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 prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing end -function RamAirKite(set::Settings, aero::BodyAerodynamics, vsm_solver::VortexStepMethod.Solver, point_system::PointMassSystem) - length(aero.wings) > 1 && throw(ArgumentError("Just one wing allowed in BodyAerodynamics object")) - wing = aero.wings[1] +function SymbolicAWEModel( + set::Settings, + system_structure::SystemStructure, + vsm_aeros::Vector{<:BodyAerodynamics}=BodyAerodynamics[], + vsm_solvers::Vector{<:VortexStepMethod.Solver}=VortexStepMethod.Solver[] +) + vsm_wings = [aero.wings[1] for aero in vsm_aeros] if set.winch_model == "TorqueControlledMachine" - s = RamAirKite{SimFloat, Vector{SimFloat}, 3*(set.segments + 1)}( - ; set, wing, aero, vsm_solver, point_system + s = SymbolicAWEModel{SimFloat, Vector{SimFloat}, 3*(set.segments + 1)}( + ; set, system_structure, vsm_wings, vsm_aeros, vsm_solvers ) s.torque_control = true else - s = RamAirKite{SimFloat, Vector{SimFloat}, 3*(set.segments + 1)}( - ; set, wing, aero, vsm_solver, point_system + s = SymbolicAWEModel{SimFloat, Vector{SimFloat}, 3*(set.segments + 1)}( + ; set, system_structure, vsm_wings, vsm_aeros, vsm_solvers ) s.torque_control = false end return s end -function RamAirKite(set::Settings) +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) - point_system = PointMassSystem(set, wing) - return RamAirKite(set, aero, vsm_solver, point_system) + system_structure = SystemStructure(set, wing) + return SymbolicAWEModel(set, system_structure, [aero], [vsm_solver]) end -function update_sys_state!(ss::SysState, s::RamAirKite, zoom=1.0) - isnothing(s.integrator) && throw(ArgumentError("run init_sim!(s) first")) +function update_sys_state!(ss::SysState, s::SymbolicAWEModel, zoom=1.0) + isnothing(s.integrator) && error("run init_sim!(s) first") ss.time = s.integrator.t # Use integrator time - # Get the extended state vector from the integrator - set_values, pos, acc_vec, Q_b_w, elevation, azimuth, course, heading, tether_length, tether_vel, winch_force, - twist_angle, kite_vel, aero_force_b, aero_moment_b, turn_rate, va_kite_b, wind_vec_gnd, wind_vel_kite = s.get_state(s.integrator) + # 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.point_system.points) + P = length(s.system_structure.points) for i in 1:P ss.X[i] = pos[1, i] * zoom ss.Y[i] = pos[2, i] * zoom @@ -112,44 +120,44 @@ function update_sys_state!(ss::SysState, s::RamAirKite, zoom=1.0) end # --- Populate SysState fields --- - ss.acc = norm(acc_vec) # Use the norm of the kite's acceleration vector - ss.orient .= Q_b_w - ss.turn_rates .= turn_rate - ss.elevation = elevation - ss.azimuth = azimuth + 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.point_system.winches) + num_winches = length(s.system_structure.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.point_system.groups) - ss.twist_angles[1:num_groups] .= twist_angle - ss.depower = rad2deg(mean(twist_angle)) # Average twist for depower - ss.steering = rad2deg(twist_angle[end] - twist_angle[1]) - ss.heading = heading # Use heading from MTK model - ss.course = course + num_groups = length(s.system_structure.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_kite_b) + ss.v_app = norm(va_wing_b[1, :]) ss.v_wind_gnd .= wind_vec_gnd - ss.v_wind_kite .= wind_vel_kite + 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_kite_b[3], va_kite_b[1]) - ss.side_slip = asin(va_kite_b[2] / ss.v_app) + 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 - ss.aero_moment_b .= aero_moment_b - ss.vel_kite .= kite_vel + 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 + 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]) @@ -169,14 +177,14 @@ function update_sys_state!(ss::SysState, s::RamAirKite, zoom=1.0) nothing end -function SysState(s::RamAirKite, zoom=1.0) - ss = SysState{length(s.point_system.points)}() +function SysState(s::SymbolicAWEModel, zoom=1.0) + ss = SysState{length(s.system_structure.points)}() update_sys_state!(ss, s, zoom) ss end """ - init_sim!(s::RamAirKite; solver=nothing, adaptive=true, prn=true, + init_sim!(s::SymbolicAWEModel; solver=nothing, adaptive=true, prn=true, precompile=false, remake=false, reload=false, lin_outputs=Num[]) -> OrdinaryDiffEqCore.ODEIntegrator @@ -197,7 +205,7 @@ and only update the state variables. Otherwise, it will create a new model from 4. Proceeds with fast path # Arguments -- `s::RamAirKite`: The kite system state object +- `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`. @@ -211,7 +219,7 @@ and only update the state variables. Otherwise, it will create a new model from # Returns - `integrator::OrdinaryDiffEqCore.ODEIntegrator`: The initialized ODE integrator. """ -function init_sim!(s::RamAirKite; +function init_sim!(s::SymbolicAWEModel; solver=nothing, adaptive=true, prn=true, precompile=false, remake=false, reload=false, lin_outputs=Num[] @@ -227,16 +235,14 @@ function init_sim!(s::RamAirKite; @warn "This solver is not tested." QNDF() else - error("Unavailable solver for RamAirKite: $(s.set.solver).") + error("Unavailable solver for SymbolicAWEModel: $(s.set.solver).") end end function init(s) - init_Q_b_w, R_b_w = initial_orient(s.set, s.wing.R_cad_body) - init!(s.point_system, s.set, R_b_w, init_Q_b_w) - - init_va_b = R_b_w' * [s.set.v_wind, 0., 0.] + init_Q_b_w, R_b_w, init_va_b = initial_orient(s) + init!(s.system_structure, s.set, R_b_w, init_Q_b_w) - inputs = create_sys!(s, s.point_system; init_va_b) + inputs = create_sys!(s, s.system_structure; init_va_b) prn && @info "Simplifying the system" prn ? (@time (sys, _) = structural_simplify(s.full_sys, (inputs, []))) : ((sys, _) = structural_simplify(sys, (inputs, []))) @@ -270,8 +276,8 @@ function init_sim!(s::RamAirKite; return s.integrator end -function linearize(s::RamAirKite; set_values=s.get_set_values(s.integrator)) - isnothing(s.lin_prob) && throw(ArgumentError("Run init_sim! with remake=true and lin_outputs=...")) +function linearize(s::SymbolicAWEModel; set_values=s.get_set_values(s.integrator)) + isnothing(s.lin_prob) && error("Run init_sim! 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)) @@ -279,7 +285,7 @@ function linearize(s::RamAirKite; set_values=s.get_set_values(s.integrator)) end """ - reinit!(s::RamAirKite, solver; prn=true, precompile=false) -> Nothing + 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 @@ -300,7 +306,7 @@ 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::RamAirKite`: The kite power system state object +- `s::SymbolicAWEModel`: The kite power system state object - `solver`: The solver to be used - `prn::Bool=true`: Whether to print progress information @@ -311,7 +317,7 @@ and only updates the state variables to match the current initial settings. - `ArgumentError`: If no serialized problem exists (run `init_sim!` first) """ function reinit!( - s::RamAirKite, + s::SymbolicAWEModel, solver; adaptive=true, prn=true, @@ -319,17 +325,17 @@ function reinit!( precompile=false, lin_outputs=Num[] ) - isnothing(s.point_system) && throw(ArgumentError("PointMassSystem not defined")) + isnothing(s.system_structure) && error("SystemStructure not defined") - init_Q_b_w, R_b_w = initial_orient(s.set, s.wing.R_cad_body) - init!(s.point_system, s.set, R_b_w, init_Q_b_w) + init_Q_b_w, R_b_w = initial_orient(s) + init!(s.system_structure, s.set, R_b_w, init_Q_b_w) if isnothing(s.prob) || reload prob_path = joinpath(KiteUtils.get_data_path(), get_prob_name(s.set; precompile)) - !ispath(prob_path) && throw(ArgumentError("$prob_path not found. Run init_sim!(s::RamAirKite) first.")) + !ispath(prob_path) && error("$prob_path not found. Run init_sim!(s::SymbolicAWEModel) first.") try (s.prob, s.full_sys, s.lin_prob, s.defaults, s.guesses) = deserialize(prob_path) - length(lin_outputs) > 0 && isnothing(s.lin_prob) && throw(ArgumentError("lin_prob is nothing.")) + length(lin_outputs) > 0 && isnothing(s.lin_prob) && error("lin_prob is nothing.") catch e @warn "Failure to deserialize $prob_path !" return s.integrator, false @@ -348,7 +354,7 @@ function reinit!( prn && @info "Initialized integrator in $t seconds" end - init_unknowns_vec!(s, s.point_system, s.unknowns_vec) + init_unknowns_vec!(s, s.system_structure, s.unknowns_vec) s.set_unknowns(s.integrator, s.unknowns_vec) OrdinaryDiffEqCore.reinit!(s.integrator, s.integrator.u; reinit_dae=true) linearize_vsm!(s) @@ -358,101 +364,131 @@ end function generate_getters!(s, sym_vec) sys = s.sys c = collect - vsm_sym = c.([ - sys.last_x, - sys.last_y, - sys.vsm_jac, - ]) + @unpack wings, winches, tethers = s.system_structure + + 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 - set_set_values = setp(sys, sys.set_values) - set_wind_dir = setp(sys, sys.upwind_dir) - set_vsm = setp(sys, vsm_sym) - set_unknowns = setu(sys, sym_vec) - set_nonstiff = setu(sys, get_nonstiff_unknowns(s)) - set_stabilize = setp(sys, sys.stabilize) - - get_vsm = getp(sys, vsm_sym) - get_set_values = getp(sys, sys.set_values) - get_unknowns = getu(sys, sym_vec) - get_state = getu(sys, - [c(sys.set_values), - c(sys.pos), # Particle positions - c(sys.acc), # Kite center acceleration vector (world frame) - 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) - c(sys.tether_length), # Unstretched length per winch - c(sys.tether_vel), # Reeling velocity per winch - c(sys.winch_force), # Force at winch connection point per winch - c(sys.twist_angle), # Twist angle per group - c(sys.kite_vel), # Kite center velocity vector (world frame) - c(sys.aero_force_b), # Aerodynamic force (body frame) - c(sys.aero_moment_b), # Aerodynamic moment (body frame) - c(sys.turn_rate), # Angular velocity (body frame) - c(sys.va_kite_b), # Apparent wind velocity (body frame) - c(sys.wind_vec_gnd), # Ground wind vector (world frame) - c(sys.wind_vel_kite) # Wind vector at kite height (world frame) - ] - ) - get_y = getu(sys, sys.y) - get_unstretched_length = getu(sys, sys.unstretched_length) - get_tether_length = getu(sys, sys.tether_length) - get_kite_pos = getu(sys, sys.kite_pos) - get_winch_force = getu(sys, sys.winch_force) - get_spring_force = getu(sys, sys.spring_force) - get_stabilize = getp(sys, sys.stabilize) - get_pos = getu(sys, sys.pos) + 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 - s.set_set_values = (integ, val) -> set_set_values(integ, val) + if length(tethers) > 0 + get_unstretched_length = getu(sys, sys.unstretched_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_wind_dir = setp(sys, sys.upwind_dir) s.set_wind_dir = (integ, val) -> set_wind_dir(integ, val) - s.set_vsm = (integ, val) -> set_vsm(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)) s.set_nonstiff = (integ, val) -> set_nonstiff(integ, val) - s.set_stabilize = (integ, val) -> set_stabilize(integ, val) - s.get_vsm = (integ) -> get_vsm(integ) - s.get_set_values = (integ) -> get_set_values(integ) + get_unknowns = getu(sys, sym_vec) s.get_unknowns = (integ) -> get_unknowns(integ) - s.get_state = (integ) -> get_state(integ) - s.get_y = (integ) -> get_y(integ) - s.get_unstretched_length = (integ) -> get_unstretched_length(integ) - s.get_tether_length = (integ) -> get_tether_length(integ) - s.get_kite_pos = (integ) -> get_kite_pos(integ) - s.get_winch_force = (integ) -> get_winch_force(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) - s.get_stabilize = (integ) -> get_stabilize(integ) + get_pos = getu(sys, sys.pos) s.get_pos = (integ) -> get_pos(integ) if !isnothing(s.lin_prob) - set_lin_set_values = setp(s.lin_prob, sys.set_values) set_lin_unknowns = setu(s.lin_prob, Initial.(sym_vec)) - set_lin_vsm = setp(s.lin_prob, vsm_sym) - - s.set_lin_set_values = (lin_prob, val) -> set_lin_set_values(lin_prob, val) s.set_lin_unknowns = (lin_prob, val) -> set_lin_unknowns(lin_prob, val) - s.set_lin_vsm = (lin_prob, val) -> set_lin_vsm(lin_prob, val) end nothing end -function linearize_vsm!(s::RamAirKite, integ=s.integrator) - y = s.get_y(integ) - jac, x = VortexStepMethod.linearize( - s.vsm_solver, - s.aero, - y; - va_idxs=1:3, - omega_idxs=4:6, - theta_idxs=7:6+length(s.point_system.groups), - moment_frac=s.point_system.groups[1].moment_frac) - s.set_vsm(integ, [x, y, jac]) +function linearize_vsm!(s::SymbolicAWEModel, integ=s.integrator) + @unpack wings, y, x, jac = s.system_structure + 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.system_structure.groups), + moment_frac=s.system_structure.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::RamAirKite, integrator::ODEIntegrator; set_values=nothing, upwind_dir=nothing, dt=1/s.set.sample_freq, vsm_interval=1) + 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. @@ -465,7 +501,7 @@ This function performs the following steps: 6. Increment the iteration counter # Arguments -- `s::RamAirKite`: The kite power system state object +- `s::SymbolicAWEModel`: The kite power system state object - `integrator::ODEIntegrator`: The ODE integrator to use # Keyword Arguments @@ -477,12 +513,12 @@ This function performs the following steps: # Returns - `Tuple{SimFloat, Float64}`: A tuple containing the current simulation time and the time taken for the step. """ -function next_step!(s::RamAirKite, integrator::OrdinaryDiffEqCore.ODEIntegrator; set_values=nothing, upwind_dir=nothing, dt=1/s.set.sample_freq, vsm_interval=1) - !(s.integrator === integrator) && error("The ODEIntegrator doesn't belong to the RamAirKite") +function next_step!(s::SymbolicAWEModel, integrator::OrdinaryDiffEqCore.ODEIntegrator; set_values=nothing, upwind_dir=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::RamAirKite; set_values=nothing, upwind_dir=nothing, dt=1/s.set.sample_freq, vsm_interval=1) +function next_step!(s::SymbolicAWEModel; set_values=nothing, upwind_dir=nothing, dt=1/s.set.sample_freq, vsm_interval=1) if (!isnothing(set_values)) s.set_set_values(s.integrator, set_values) end @@ -516,8 +552,8 @@ end """ Calculate and return the angle of attack in rad """ -function calc_aoa(s::RamAirKite) - alpha_array = s.vsm_solver.sol.alpha_array +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] @@ -527,14 +563,14 @@ function calc_aoa(s::RamAirKite) end function init_unknowns_vec!( - s::RamAirKite, - system::PointMassSystem, + s::SymbolicAWEModel, + system::SystemStructure, vec::Vector{SimFloat} ) !s.set.quasi_static && (length(vec) != length(s.integrator.u)) && - throw(ArgumentError("Unknowns of length $(length(s.integrator.u)) but vector provided of length $(length(vec))")) + error("Unknowns of length $(length(s.integrator.u)) but vector provided of length $(length(vec))") - @unpack points, groups, segments, pulleys, winches, kite = system + @unpack points, groups, segments, pulleys, winches, wings = system vec_idx = 1 for point in points @@ -571,31 +607,34 @@ function init_unknowns_vec!( vec[vec_idx] = winch.tether_vel vec_idx += 1 end - for i in 1:4 - vec[vec_idx] = kite.orient[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = kite.angular_vel[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = kite.pos[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = kite.vel[i] - vec_idx += 1 + + 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[i] + vec_idx += 1 + end + for i in 1:3 + vec[vec_idx] = wing.vel[i] + vec_idx += 1 + end end (vec_idx-1 != length(vec)) && - throw(ArgumentError("Unknowns vec is of length $(length(vec)) but the last index is $(vec_idx-1)")) + error("Unknowns vec is of length $(length(vec)) but the last index is $(vec_idx-1)") nothing end -function get_unknowns(s::RamAirKite) +function get_unknowns(s::SymbolicAWEModel) vec = Num[] - @unpack points, groups, segments, pulleys, winches, kite = s.point_system + @unpack points, groups, segments, pulleys, winches, wings = s.system_structure sys = s.sys for point in points for i in 1:3 @@ -611,12 +650,13 @@ function get_unknowns(s::RamAirKite) end vec = get_nonstiff_unknowns(s, vec) !s.set.quasi_static && (length(vec) != length(s.integrator.u)) && - throw(ArgumentError("Integrator unknowns of length $(length(s.integrator.u)) should equal vec of length $(length(vec))")) + error("Integrator unknowns of length $(length(s.integrator.u)) should equal vec of length $(length(vec)). + Maybe you forgot to run init_sim!(model; remake=true)?") return vec end -function get_nonstiff_unknowns(s::RamAirKite, vec=Num[]) - @unpack points, groups, segments, pulleys, winches, kite = s.point_system +function get_nonstiff_unknowns(s::SymbolicAWEModel, vec=Num[]) + @unpack points, groups, segments, pulleys, winches, wings = s.system_structure sys = s.sys for group in groups @@ -627,14 +667,16 @@ function get_nonstiff_unknowns(s::RamAirKite, vec=Num[]) push!(vec, sys.tether_length[winch.idx]) push!(vec, sys.tether_vel[winch.idx]) end - [push!(vec, sys.Q_b_w[i]) for i in 1:4] - [push!(vec, sys.ω_b[i]) for i in 1:3] - [push!(vec, sys.kite_pos[i]) for i in 1:3] - [push!(vec, sys.kite_vel[i]) for i in 1:3] + 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::RamAirKite; dt=1/s.set.sample_freq) +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 @@ -644,12 +686,33 @@ function find_steady_state!(s::RamAirKite; dt=1/s.set.sample_freq) return nothing end -unstretched_length(s::RamAirKite) = s.get_unstretched_length(s.integrator) -tether_length(s::RamAirKite) = s.get_tether_length(s.integrator) -calc_height(s::RamAirKite) = s.get_kite_pos(s.integrator)[3] -winch_force(s::RamAirKite) = s.get_winch_force(s.integrator) -spring_forces(s::RamAirKite) = s.get_spring_force(s.integrator) -function pos(s::RamAirKite) +function initial_orient(s::SymbolicAWEModel) + set = s.set + wings = s.system_structure.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 +end diff --git a/test/create_xz_file.jl b/test/create_xz_file.jl index c6407c85c..a8bb227ea 100644 --- a/test/create_xz_file.jl +++ b/test/create_xz_file.jl @@ -26,15 +26,15 @@ set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. set.quasi_static = false set.physical_model = "ram" -@info "Creating wing, aero, vsm_solver, point_system and s:" -s = RamAirKite(set) +@info "Creating wing, aero, vsm_solver, system_structure and s:" +s = SymbolicAWEModel(set) s.set.abs_tol = 1e-2 s.set.rel_tol = 1e-2 toc() # Initialize at elevation -s.point_system.winches[2].tether_length += 0.2 -s.point_system.winches[3].tether_length += 0.2 +s.system_structure.winches[2].tether_length += 0.2 +s.system_structure.winches[3].tether_length += 0.2 KiteModels.init_sim!(s; remake=false, reload=true) sys = s.sys diff --git a/test/test_for_precompile.jl b/test/test_for_precompile.jl index 3bf58a029..d1c1a72e2 100644 --- a/test/test_for_precompile.jl +++ b/test/test_for_precompile.jl @@ -110,9 +110,9 @@ if ! haskey(ENV, "NO_MTK") set = se("system_ram.yaml") set.segments = 2 set_values = [-50, -1.1, -1.1] - mtk_kite = RamAirKite(set) + mtk_kite = SymbolicAWEModel(set) KiteModels.init_sim!(mtk_kite) - logger = Logger(length(mtk_kite.point_system.points), 5) + logger = Logger(length(mtk_kite.system_structure.points), 5) for i in 1:5 next_step!(mtk_kite, set_values) diff --git a/test/test_ram_air_kite.jl b/test/test_ram_air_kite.jl index 1be8501f8..b844a1af5 100644 --- a/test/test_ram_air_kite.jl +++ b/test/test_ram_air_kite.jl @@ -16,7 +16,7 @@ set_data_path(temp_data_path) const TOL = 1e-5 const BUILD_SYS = true -@testset verbose = true "RamAirKite MTK Model Tests" begin +@testset verbose = true "SymbolicAWEModel MTK Model Tests" begin # Initialize model set = se("system_ram.yaml") set.segments = 3 @@ -25,7 +25,7 @@ const BUILD_SYS = true set.physical_model = "ram" @info "Creating s:" - @time s = RamAirKite(set) + @time s = SymbolicAWEModel(set) s.set.abs_tol = 1e-2 s.set.rel_tol = 1e-2 @@ -53,14 +53,14 @@ const BUILD_SYS = true # Check initialization results @test !isnothing(s.integrator) @test !isnothing(s.sys) - @test !isnothing(s.point_system) + @test !isnothing(s.system_structure) end s.integrator = nothing s.sys = nothing # Keep references to first integrator and point system first_integrator_ptr = objectid(s.integrator) - first_point_system_ptr = objectid(s.point_system) + first_system_structure_ptr = objectid(s.system_structure) # 2. First init_sim! - should load from serialized file @info "Testing first init_sim! (should load serialized file)..." @@ -69,9 +69,9 @@ const BUILD_SYS = true # Check that it's a new integrator second_integrator_ptr = objectid(s.integrator) - second_point_system_ptr = objectid(s.point_system) + second_system_structure_ptr = objectid(s.system_structure) @test first_integrator_ptr != second_integrator_ptr - @test first_point_system_ptr == second_point_system_ptr + @test first_system_structure_ptr == second_system_structure_ptr # 3. Second init_sim! - should reuse existing integrator @info "Testing second init_sim! (should reuse integrator)..." @@ -79,19 +79,19 @@ const BUILD_SYS = true # This should create a new point system but reuse the existing integrator third_integrator_ptr = objectid(s.integrator) - third_point_system_ptr = objectid(s.point_system) + third_system_structure_ptr = objectid(s.system_structure) @test second_integrator_ptr == third_integrator_ptr # Should be the same - @test second_point_system_ptr == third_point_system_ptr + @test second_system_structure_ptr == third_system_structure_ptr # Get positions using SysState sys_state = KiteModels.SysState(s) # Check dimension consistency - # Note: pos_integrator is no longer directly fetched, comparing SysState to point_system - @test length(sys_state.X) == length(s.point_system.points) + # Note: pos_integrator is no longer directly fetched, comparing SysState to system_structure + @test length(sys_state.X) == length(s.system_structure.points) # Compare positions in different representations - for (i, point) in enumerate(s.point_system.points) + for (i, point) in enumerate(s.system_structure.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]] @@ -100,7 +100,7 @@ const BUILD_SYS = true @test isapprox(norm(point_pos), norm(sys_state_pos), rtol=1e-2) # Positions should not be zero (except ground points) - if point.type != KiteModels.WINCH # Skip ground points which might be at origin + 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 @@ -233,7 +233,7 @@ const BUILD_SYS = true end end - @testset "Plotting of RamAirKite" begin + @testset "Plotting of SymbolicAWEModel" begin plt.figure("Kite") lines, sc, txt = plot(s, 0.0) plt.show(block=false)