Skip to content

Commit 3d05733

Browse files
1-Bart-11-Bart-1
committed
New package
1 parent 6f9ec98 commit 3d05733

23 files changed

+4916
-4
lines changed

‎examples/lin_ram_model.jl‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
# Copyright (c) 2025 Bart van de Lint
2+
# SPDX-License-Identifier: MPL-2.0
3+
4+
#=
5+
This example demonstrates linearized model accuracy by comparing:
6+
1. Nonlinear SymbolicAWEModel model simulation
7+
2. Linearized state-space model simulation
8+
9+
Both models start from the same operating point and are subjected
10+
to identical steering inputs. The resulting state trajectories are
11+
plotted together to visualize how well the linearized model
12+
approximates the nonlinear dynamics.
13+
=#
14+
15+
using Timers
16+
tic()
17+
@info "Loading packages "
18+
19+
PLOT = true
20+
if PLOT
21+
using Pkg
22+
if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
23+
using TestEnv; TestEnv.activate()
24+
end
25+
using ControlPlots, LaTeXStrings, ControlSystemsBase
26+
end
27+
28+
using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore
29+
using ModelingToolkit
30+
using ModelingToolkit: t_nounits
31+
toc()
32+
33+
# TODO: use sparse autodiff
34+
35+
# Simulation parameters
36+
dt = 0.001
37+
total_time = 1.0 # Increased from 0.1s to 1.0s for better dynamics observation
38+
vsm_interval = 3
39+
steps = Int(round(total_time / dt))
40+
41+
# Steering parameters
42+
steering_freq = 1/2 # Hz - full left-right cycle frequency
43+
steering_magnitude = 5.0 # Magnitude of steering input [Nm]
44+
45+
# Initialize model
46+
set = load_settings("system_ram.yaml")
47+
set.segments = 3
48+
set_values = [-50.0, 0.0, 0.0] # Set values of the torques of the three winches. [Nm]
49+
set.quasi_static = false
50+
set.physical_model = "simple_ram"
51+
52+
@info "Creating SymbolicAWEModel model..."
53+
s = SymbolicAWEModel(set)
54+
s.set.abs_tol = 1e-2
55+
s.set.rel_tol = 1e-2
56+
toc()
57+
58+
# Define outputs for linearization - heading
59+
lin_outputs = @variables heading(t_nounits)[1]
60+
61+
# Initialize at elevation with linearization outputs
62+
s.sys_struct.winches[2].tether_length += 0.2
63+
s.sys_struct.winches[3].tether_length += 0.2
64+
KiteModels.init_sim!(s;
65+
remake=false,
66+
reload=true,
67+
lin_outputs # Specify which outputs to track in linear model
68+
)
69+
sys = s.sys
70+
71+
@show rad2deg(s.integrator[sys.elevation[1]])
72+
73+
74+
@info "System initialized at:"
75+
toc()
76+
77+
# --- Stabilize system at operating point ---
78+
@info "Stabilizing system at operating point..."
79+
s.integrator.ps[sys.stabilize] = true
80+
stabilization_steps = Int(10 ÷ dt)
81+
for i in 1:stabilization_steps
82+
next_step!(s; dt, vsm_interval=0.05÷dt)
83+
end
84+
s.integrator.ps[sys.stabilize] = false
85+
86+
# --- Linearize at operating point ---
87+
@info "Linearizing system at operating point..."
88+
@time (; A, B, C, D) = KiteModels.linearize(s)
89+
@time (; A, B, C, D) = KiteModels.linearize(s)
90+
@show norm(A)
91+
@info "System linearized with matrix dimensions:" A=size(A) B=size(B) C=size(C) D=size(D)
92+
93+
sys = ss(A,B,C,D)
94+
bode_plot(sys[1,1]; from=1e-4)
95+
bode_plot(sys[1,2]; from=1e-4)
96+
bode_plot(sys[1,3]; from=1e-4)

‎examples/menu.jl‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
# Copyright (c) 2022, 2024 Uwe Fechner
2+
# SPDX-License-Identifier: MIT
3+
using REPL.TerminalMenus
4+
5+
options = ["bench = include(\"bench.jl\")",
6+
"bench_4p = include(\"bench_4p.jl\")",
7+
"test_init_1p = include(\"test_init_1p.jl\")",
8+
"test_init_4p = include(\"test_init_4p.jl\")",
9+
"plot_cl_cd_plate = include(\"plot_cl_cd_plate.jl\")",
10+
"plot_side_cl = include(\"plot_side_cl.jl\")",
11+
"compare_kps3_kps4 = include(\"compare_kps3_kps4.jl\")",
12+
"reel_out_1p = include(\"reel_out_1p.jl\")",
13+
"reel_out_4p = include(\"reel_out_4p.jl\")",
14+
"reel_out_4p_torque_control = include(\"reel_out_4p_torque_control.jl\")",
15+
"simulate_simple = include(\"simulate_simple.jl\")",
16+
"simulate_steering = include(\"simulate_steering.jl\")",
17+
"steering_test_1p = include(\"steering_test_1p.jl\")",
18+
"steering_test_4p = include(\"steering_test_4p.jl\")",
19+
"ram_air_kite = SIMPLE=false; include(\"ram_air_kite.jl\")",
20+
"simple_ram_air_kite = SIMPLE=true; include(\"ram_air_kite.jl\")",
21+
"lin_ram_model = include(\"lin_ram_model.jl\")",
22+
"calc_spectrum = include(\"calc_spectrum.jl\")",
23+
"plot_spectrum_ = include(\"plot_spectrum.jl\")",
24+
"calculate_rotational_inertia = include(\"calculate_rotational_inertia.jl\")",
25+
"quit"]
26+
27+
function example_menu()
28+
active = true
29+
while active
30+
menu = RadioMenu(options, pagesize=8)
31+
choice = request("\nChoose function to execute or `q` to quit: ", menu)
32+
33+
if choice != -1 && choice != length(options)
34+
eval(Meta.parse(options[choice]))
35+
else
36+
println("Left menu. Press <ctrl><d> to quit Julia!")
37+
active = false
38+
end
39+
end
40+
end
41+
42+
example_menu()

‎examples/pulley.jl‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
# SPDX-FileCopyrightText: 2025 Bart van de Lint
2+
#
3+
# SPDX-License-Identifier: MPL-2.0
4+
5+
using KiteModels, VortexStepMethod, ControlPlots
6+
7+
set = se("system_ram.yaml")
8+
set.v_wind = 10.0
9+
set.l_tether = 5.0
10+
set.abs_tol = 1e-4
11+
set.rel_tol = 1e-4
12+
dynamics_type = DYNAMIC
13+
14+
points = Point[]
15+
segments = Segment[]
16+
pulleys = Pulley[]
17+
18+
push!(points, Point(1, [0.0, 0.0, 2.0], STATIC))
19+
push!(points, Point(2, [2.0, 0.0, 2.0], STATIC))
20+
push!(points, Point(3, [0.1, 0.0, 1.0], DYNAMIC))
21+
push!(points, Point(4, [0.1, 0.0, 0.0], DYNAMIC; mass=0.1))
22+
23+
push!(segments, Segment(1, (3,1), BRIDLE))
24+
push!(segments, Segment(2, (3,2), BRIDLE))
25+
push!(segments, Segment(3, (3,4), BRIDLE))
26+
27+
push!(pulleys, Pulley(1, (1,2), DYNAMIC))
28+
29+
transforms = [Transform(1, -deg2rad(0.0), 0.0, 0.0; base_pos=[1.0, 0.0, 4.0], base_point_idx=1, rot_point_idx=2)]
30+
sys_struct = KiteModels.SystemStructure("pulley", set; points, segments, pulleys, transforms)
31+
plot(sys_struct, 0.0; zoom=false, l_tether=set.l_tether)
32+
33+
sam = SymbolicAWEModel(set, sys_struct)
34+
35+
init_sim!(sam; remake=false)
36+
for i in 1:100
37+
plot(sam, i/set.sample_freq; zoom=false)
38+
next_step!(sam)
39+
end

‎examples/ram_air_kite.jl‎

Lines changed: 149 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,149 @@
1+
# Copyright (c) 2024, 2025 Bart van de Lint, Uwe Fechner
2+
# SPDX-License-Identifier: MIT
3+
4+
using Timers
5+
tic()
6+
@info "Loading packages "
7+
8+
PLOT = false
9+
using Pkg
10+
if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
11+
using TestEnv; TestEnv.activate()
12+
end
13+
using ControlPlots, LaTeXStrings
14+
using KiteModels, LinearAlgebra, Statistics
15+
16+
if ! @isdefined SIMPLE
17+
SIMPLE = false
18+
end
19+
20+
toc()
21+
22+
# Simulation parameters
23+
dt = 0.05
24+
total_time = 10.0 # Longer simulation to see oscillations
25+
vsm_interval = 3
26+
steps = Int(round(total_time / dt))
27+
28+
# Steering parameters
29+
steering_freq = 1/2 # Hz - full left-right cycle frequency
30+
steering_magnitude = 10.0 # Magnitude of steering input [Nm]
31+
32+
# Initialize model
33+
set = Settings("system_ram.yaml")
34+
set.segments = 3
35+
set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm]
36+
set.quasi_static = false
37+
set.physical_model = SIMPLE ? "simple_ram" : "ram"
38+
39+
@info "Creating wing, aero, vsm_solver, sys_struct and symbolic_awe_model:"
40+
sam = SymbolicAWEModel(set)
41+
sam.set.abs_tol = 1e-2
42+
sam.set.rel_tol = 1e-2
43+
toc()
44+
45+
# Initialize at elevation
46+
set.l_tethers[2] += 0.4
47+
set.l_tethers[3] += 0.4
48+
init_sim!(sam; remake=false, reload=false)
49+
sys = sam.sys
50+
51+
@info "System initialized at:"
52+
toc()
53+
54+
# Stabilize system
55+
find_steady_state!(sam)
56+
57+
logger = Logger(length(sam.sys_struct.points), steps)
58+
sys_state = SysState(sam)
59+
t = 0.0
60+
runtime = 0.0
61+
integ_runtime = 0.0
62+
bias = set.quasi_static ? 0.45 : 0.35
63+
t0 = sam.integrator.t
64+
65+
try
66+
while t < total_time
67+
local steering
68+
global t, set_values, runtime, integ_runtime
69+
PLOT && plot(sam, t; zoom=false, front=false)
70+
71+
# Calculate steering inputs based on cosine wave
72+
steering = steering_magnitude * cos(2π * steering_freq * t + bias)
73+
set_values = -sam.set.drum_radius .* sam.integrator[sys.winch_force]
74+
_vsm_interval = 1
75+
if t > 1.0
76+
set_values .+= [0.0, steering, -steering] # Opposite steering for left/right
77+
_vsm_interval = vsm_interval
78+
end
79+
80+
# Step simulation
81+
steptime = @elapsed next_step!(sam; set_values, dt, vsm_interval=vsm_interval)
82+
t_new = sam.integrator.t
83+
integ_steptime = sam.t_step
84+
t = t_new - t0 # Adjust for initial stabilization time
85+
86+
# Track performance after initial transient
87+
if (t > total_time/2)
88+
runtime += steptime
89+
integ_runtime += integ_steptime
90+
sam.integrator.ps[sys.twist_damp] = 10
91+
end
92+
93+
# Log state variables
94+
update_sys_state!(sys_state, sam)
95+
sys_state.time = t
96+
log!(logger, sys_state)
97+
end
98+
catch e
99+
if isa(e, AssertionError)
100+
@show t
101+
println(e)
102+
else
103+
rethrow(e)
104+
end
105+
end
106+
@info "Total time without plotting:"
107+
toc()
108+
109+
# Plot results
110+
c = collect
111+
save_log(logger, "tmp")
112+
lg =load_log("tmp")
113+
sl = lg.syslog
114+
115+
# --- Updated Plotting ---
116+
# Extract necessary data using meaningful names
117+
turn_rates_deg = rad2deg.(hcat(sl.turn_rates...))
118+
v_reelout_23 = [sl.v_reelout[i][2] for i in eachindex(sl.v_reelout)], [sl.v_reelout[i][3] for i in eachindex(sl.v_reelout)] # Winch 2 and 3
119+
aero_force_z = [sl.aero_force_b[i][3] for i in eachindex(sl.aero_force_b)]
120+
aero_moment_z = [sl.aero_moment_b[i][3] for i in eachindex(sl.aero_moment_b)]
121+
twist_angles_deg = rad2deg.(hcat(sl.twist_angles...))
122+
AoA_deg = rad2deg.(sl.AoA)
123+
heading_deg = rad2deg.(sl.heading)
124+
125+
p = plotx(sl.time,
126+
[turn_rates_deg[1,:], turn_rates_deg[2,:], turn_rates_deg[3,:]],
127+
v_reelout_23,
128+
[aero_force_z, aero_moment_z],
129+
[twist_angles_deg[1,:], twist_angles_deg[2,:], twist_angles_deg[3,:], twist_angles_deg[4,:]],
130+
[AoA_deg],
131+
[heading_deg];
132+
ylabels=["turn rates [°/s]", L"v_{ro}~[m/s]", "aero F/M", "twist [°]", "AoA [°]", "heading [°]"],
133+
ysize=10,
134+
labels=[
135+
[L"\omega_x", L"\omega_y", L"\omega_z"],
136+
["v_ro[2]", "v_ro[3]"],
137+
[L"F_{aero,z}", L"M_{aero,z}"],
138+
["twist[1]", "twist[2]", "twist[3]", "twist[4]"],
139+
["AoA"],
140+
["heading"]
141+
],
142+
fig="Oscillating Steering Input Response")
143+
display(p)
144+
145+
@info "Performance:" times_realtime=(total_time/2)/runtime integrator_times_realtime=(total_time/2)/integ_runtime
146+
147+
# 55x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz)
148+
# 40-65x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) - commit 6620ed5d0a38e96930615aad9a66e4cd666955f2
149+
# 40x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) - commit 88a78894038d3cbd50fbff83dfbe5c26266b0637

0 commit comments

Comments
 (0)