Revert "Remove SymbolicAWEModel code from the package (#247)"

This reverts commit 0e5fb65.

Author: ufechner7

Project.toml

@@ -27,6 +27,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -35,7 +36,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
-SymbolicAWEModels = "9c9a347c-5289-41db-a9b9-25ccc76c3360"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Timers = "21f18d07-b854-4dab-86f0-c15a3821819a"
@@ -70,7 +70,7 @@ KiteUtils = "0.10.15"
 LaTeXStrings = "1.4.0"
 LinearAlgebra = "1.10, 1.11"
 LinearSolve = "~2.39.0, 3"
-ModelingToolkit = "10.10.0"
+ModelingToolkit = "~9.78.0"
 NLSolversBase = "~7.8.3"
 NLsolve = "4.5"
 NonlinearSolve = "4.8.0"
@@ -82,6 +82,7 @@ Parameters = "0.12"
 Pkg = "1.10, 1.11"
 PrecompileTools = "1.2, 1.3"
 REPL = "1.10.0, 1.11.0"
+Reexport = "1.1, 1.2"
 Rotations = "1.7"
 SHA = "0.7.0"
 SciMLBase = "<2.115"
@@ -90,12 +91,11 @@ StaticArrays = "1.9.7"
 Statistics = "1"
 StatsBase = "0.34"
 Sundials = "4.24"
-SymbolicAWEModels = "0.1.2"
 SymbolicIndexingInterface = "0.3"
 SymbolicUtils = "3.25"
 Test = "1"
 Timers = "0.1.5"
-VortexStepMethod = "2.0.0"
+VortexStepMethod = "1.2.6"
 WinchModels = "0.3.6"
 julia = "1.10, 1.11"
 

data/settings_ram.yaml

@@ -3,14 +3,14 @@ system:
                                    #   use / as path delimiter, even on Windows 
     time_lapse:   1.0              # relative replay speed
     sim_time:   100.0              # simulation time             [sim only]
-    segments:       3              # number of tether segments
+    segments:       6              # number of tether segments
     sample_freq:   20              # sample frequency in Hz
     zoom:        0.03              # zoom factor for the system view
     kite_scale:   3.0              # relative zoom factor for the 4 point kite
     fixed_font: ""                 # name or filepath+filename of alternative fixed pitch font
 
 initial:
-    l_tethers: [50.0, 50.4, 50.4]  # initial tether length       [m]
+    l_tethers: [50.0, 50.0, 50.0]  # initial tether length       [m]
     elevation: 70.8                # initial elevation angle   [deg]
     v_reel_outs: [0.0, 0.0, 0.0]   # initial reel out speed    [m/s]
 
@@ -69,7 +69,7 @@ environment:
     rho_0:  1.225            # air density at zero height and 15 °C    [kg/m³]
     alpha:  0.08163          # exponent of the wind profile law
     z0:     0.0002           # surface roughness                       [m]
-    profile_law: 3           # 1=EXP, 2=LOG, 3=EXPLOG
+    profile_law: 3           # 1=EXP, 2=LOG, 3=EXPLOG, 4=F
     # the following parameters are for calculating the turbulent wind field using the Mann model
     use_turbulence: 0.0      # turbulence intensity relative to Cabau, NL
     v_wind_gnds: [3.483, 5.324, 8.163] # wind speeds at ref height for calculating the turbulent wind field [m/s]

docs/make.jl

@@ -28,6 +28,7 @@ makedocs(;
         "Home" => "index.md",
         "Types" => "types.md",
         "Functions" => "functions.md",
+        "SymbolicAWEModel" => "ram_air_kite.md",
         "Parameters" => "parameters.md",
         "Examples 1p" => "examples.md",
         "Examples 4p" => "examples_4p.md",

docs/src/functions.md

@@ -34,6 +34,7 @@ set_v_wind_ground!
 unstretched_length
 tether_length
 pos_kite
+calc_aoa
 calc_height
 calc_elevation
 calc_azimuth

docs/src/index.md

@@ -83,7 +83,7 @@ times faster.
 - the documentation was improved
 
 ## Provides
-The types [`KPS3`](@ref), [`KPS4`](@ref) and `SymbolicAWEModel`, representing the one point, the four point kite model and the ram air kite model, together with the high level simulation interface consisting of the functions [`init!`](@ref) and [`next_step!`](@ref). Other kite models can be added inside or outside of this package by implementing the non-generic methods required for an AbstractKiteModel.
+The types [`KPS3`](@ref), [`KPS4`](@ref) and [`SymbolicAWEModel`](@ref), representing the one point, the four point kite model and the ram air kite model, together with the high level simulation interface consisting of the functions [`init!`](@ref) and [`next_step!`](@ref). Other kite models can be added inside or outside of this package by implementing the non-generic methods required for an AbstractKiteModel.
 
 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

docs/src/ram_air_kite.md

@@ -0,0 +1,53 @@
+```@meta
+CurrentModule = KiteModels
+```
+## Introduction
+The [`SystemStructure`](@ref) provides a flexible framework for defining the physical structure of airborne wind energy (AWE) systems using discrete mass-spring-damper models. This structure can represent many different AWE system configurations, from simple single-line kites to complex multi-wing systems with intricate bridle networks.
+
+The [`SystemStructure`](@ref) serves as input to the [`SymbolicAWEModel`](@ref), which is based on ModelingToolkit and automatically generates symbolic differential algebraic equations from the structural definition.
+
+## Workflow
+1. Define system components ([`Point`](@ref), [`Segment`](@ref), [`Group`](@ref), etc.) 
+2. Assemble into a [`SystemStructure`](@ref)
+3. Pass to [`SymbolicAWEModel`](@ref) for automatic MTK model generation
+4. Simulate the resulting symbolic model
+
+## Public enumerations
+```@docs
+SegmentType
+DynamicsType
+```
+
+## Public constructors
+```@docs
+SystemStructure(name, set; points=Point[], groups=Group[], segments=Segment[], 
+                   pulleys=Pulley[], tethers=Tether[], winches=Winch[], 
+                   wings=Wing[], transforms=Transform[])
+SystemStructure
+SymbolicAWEModel(::Settings, ::SystemStructure, ::Vector{<:BodyAerodynamics}, ::Vector{<:VortexStepMethod.Solver})
+SymbolicAWEModel(::Settings)
+Point(idx, pos_cad, type)
+Point
+Group(::Any, ::Any, ::RamAirWing, ::Any, ::Any, ::Any)
+Group
+Segment(idx, point_idxs, type)
+Segment
+Pulley(idx, segment_idxs, type)
+Pulley
+Tether
+Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0)
+Winch
+Wing(idx, group_idxs, R_b_c, pos_cad)
+Wing
+Transform(idx, elevation, azimuth, heading)
+Transform
+```
+
+## Private functions
+```@docs
+wing_eqs!
+reinit!
+scalar_eqs!
+linear_vsm_eqs!
+force_eqs!
+```

docs/src/tutorial_system_structure.md

@@ -35,7 +35,7 @@ segments = Segment[]
 points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0))
 ```
 
-The first point we add is a static point. There are four different `DynamicsType`s to choose from: `STATIC`, `QUASI_STATIC`, `DYNAMIC` and `WING`. `STATIC` just means that the point doesn't move. `DYNAMIC` is a point modeled with acceleration, while `QUASI_STATIC` constrains this acceleration to be zero at all times. A `WING` point is connected to a wing body.
+The first point we add is a static point. There are four different [`DynamicsType`](@ref)s to choose from: `STATIC`, `QUASI_STATIC`, `DYNAMIC` and `WING`. `STATIC` just means that the point doesn't move. `DYNAMIC` is a point modeled with acceleration, while `QUASI_STATIC` constrains this acceleration to be zero at all times. A `WING` point is connected to a wing body.
 
 Now we can add `DYNAMIC` points and connect them to each other with segments. `BRIDLE` segments don't need to have a tether, because they have a constant unstretched length.
 ```julia
@@ -51,21 +51,21 @@ for i in 1:set.segments
 end
 ```
 
-In order to describe the initial orientation of the structure, we define a `Transform` with an elevation (-80 degrees), azimuth and heading, and a base position `[0.0, 0.0, 50.0]`.
+In order to describe the initial orientation of the structure, we define a [`Transform(idx, elevation, azimuth, heading)`](@ref) with an elevation (-80 degrees), azimuth and heading, and a base position `[0.0, 0.0, 50.0]`.
 ```julia
 transforms = [Transform(1, deg2rad(-80), 0.0, 0.0; 
               base_pos = [0.0, 0.0, 50.0], base_point_idx=points[1].idx,
               rot_point_idx=points[end].idx)]
 ```
 
-From the points, segments and transform we create a `SystemStructure`, which can be plotted in 2d to quickly investigate if the model is correct.
+From the points, segments and transform we create a [`SystemStructure(name, set)`](@ref), which can be plotted in 2d to quickly investigate if the model is correct.
 ```julia
 sys_struct = SystemStructure("tether", set; points, segments, transforms)
 plot(sys_struct, 0.0)
 ```
 ![SystemStructure visualization](tether_sys_struct.png)
 
-If the system looks good, we can easily model it, by first creating a `SymbolicAWEModel`, initializing it and stepping through time.
+If the system looks good, we can easily model it, by first creating a [`SymbolicAWEModel`](@ref), initializing it and stepping through time.
 ```julia
 sam = SymbolicAWEModel(set, sys_struct)
 

docs/src/types.md

@@ -16,6 +16,7 @@ AKM
 ```@docs
 KPS3
 KPS4
+SymbolicAWEModel
 ```
 These structs store the state of the one point model and four point model. Only in unit tests
 it is allowed to access the members directly, otherwise use the input and output functions.

examples/lin_ram_model.jl

@@ -1,63 +1,93 @@
 # Copyright (c) 2025 Bart van de Lint
 # SPDX-License-Identifier: MPL-2.0
 
-# This example demonstrates how to:
-# - Load a system configuration from a YAML file,
-# - Initialize a SymbolicAWEModel with specified winch torques,
-# - Stabilize the system at its operating point,
-# - Linearize the system at that point,
-# - And plot the Bode plots of the resulting linearized system.
+#=
+This example demonstrates linearized model accuracy by comparing:
+1. Nonlinear SymbolicAWEModel model simulation 
+2. Linearized state-space model simulation
+
+Both models start from the same operating point and are subjected
+to identical steering inputs. The resulting state trajectories are
+plotted together to visualize how well the linearized model
+approximates the nonlinear dynamics.
+=#
 
 using Timers
 tic()
 @info "Loading packages "
 
 PLOT = true
-using Pkg
-if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
-    using TestEnv; TestEnv.activate()
+if PLOT
+    using Pkg
+    if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
+        using TestEnv; TestEnv.activate()
+    end
+    using ControlPlots, LaTeXStrings, ControlSystemsBase
 end
-using ControlPlots, LaTeXStrings, ControlSystemsBase
 
-using KiteModels, KiteUtils, LinearAlgebra, Statistics
-using KiteModels.SymbolicAWEModels: find_steady_state!
-using KiteModels.SymbolicAWEModels.ModelingToolkit
-using KiteModels.SymbolicAWEModels.ModelingToolkit: t_nounits
+using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore 
+using ModelingToolkit
+using ModelingToolkit: t_nounits
 toc()
 
+# TODO: use sparse autodiff
+
+# Simulation parameters
+dt = 0.001
+total_time = 1.0  # Increased from 0.1s to 1.0s for better dynamics observation
+vsm_interval = 3
+steps = Int(round(total_time / dt))
+
+# Steering parameters
+steering_freq = 1/2  # Hz - full left-right cycle frequency
+steering_magnitude = 5.0      # Magnitude of steering input [Nm]
+
 # Initialize model
-set = Settings("system_ram.yaml")
+set = load_settings("system_ram.yaml")
+set.segments = 3
 set_values = [-50.0, 0.0, 0.0]  # Set values of the torques of the three winches. [Nm]
+set.quasi_static = false
+set.physical_model = "simple_ram"
 
-@info "Creating SymbolicAWEModel..."
+@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
-@variables begin
-    heading(t_nounits)[1]
-    angle_of_attack(t_nounits)[1]
-    tether_len(t_nounits)[1:3]
-    winch_force(t_nounits)[1:3]
-end
-lin_outputs = [heading[1], angle_of_attack[1], tether_len[1], winch_force[1]]
-@info "Linear outputs: $lin_outputs"
+lin_outputs = @variables heading(t_nounits)[1]
 
 # Initialize at elevation with linearization outputs
-init!(s; lin_outputs)
+s.sys_struct.winches[2].tether_length += 0.2
+s.sys_struct.winches[3].tether_length += 0.2
+KiteModels.init!(s; 
+    remake=false,
+    reload=true,
+    lin_outputs  # Specify which outputs to track in linear model
+)
 sys = s.sys
 
+@show rad2deg(s.integrator[sys.elevation[1]])
+
+
 @info "System initialized at:"
 toc()
 
 # --- Stabilize system at operating point ---
 @info "Stabilizing system at operating point..."
-find_steady_state!(s)
+s.integrator.ps[sys.stabilize] = true
+stabilization_steps = Int(10 ÷ dt)
+for i in 1:stabilization_steps
+    next_step!(s; dt, vsm_interval=0.05÷dt)
+end
+s.integrator.ps[sys.stabilize] = false
 
 # --- Linearize at operating point ---
 @info "Linearizing system at operating point..."
-(; A, B, C, D) = linearize!(s)
-@time (; A, B, C, D) = linearize!(s)
+@time (; A, B, C, D) = KiteModels.linearize(s)
+@time (; A, B, C, D) = KiteModels.linearize(s)
+@show norm(A)
 @info "System linearized with matrix dimensions:" A=size(A) B=size(B) C=size(C) D=size(D)
 
 sys = ss(A,B,C,D)

examples/menu.jl

@@ -18,6 +18,7 @@ options = ["bench = include(\"bench.jl\")",
            "steering_test_4p = include(\"steering_test_4p.jl\")",
            "ram_air_kite = SIMPLE=false; include(\"ram_air_kite.jl\")",
            "simple_ram_air_kite = SIMPLE=true; include(\"ram_air_kite.jl\")",
+           "lin_ram_model = include(\"lin_ram_model.jl\")",
            "calc_spectrum = include(\"calc_spectrum.jl\")",
            "plot_spectrum_ = include(\"plot_spectrum.jl\")",
            "calculate_rotational_inertia = include(\"calculate_rotational_inertia.jl\")",
@@ -38,4 +39,4 @@ function example_menu()
     end
 end
 
-example_menu()
+example_menu()