Add an option for a quasi-static model

data/settings_ram.yaml

@@ -35,6 +35,7 @@ kite:
     bridle_fracs: [0.088, 0.31, 0.58, 0.93] # distances along normalized foil chord for bridle attachment points
     fixed_index: 1                          # the bridle frac index around which the kite twists
     mass: 0.9                               # kite mass [kg]
+    quasi_static: false                     # whether to use quasi static kite points or not
 
 tether:
     bridle_tether_diameter: 2 # [mm]

examples/ram_air_kite_loop.jl

@@ -8,13 +8,14 @@ end
 include("./plotting.jl")
 
 dt = 0.05
-total_time = 4.5
+total_time = 1.8
 vsm_interval = 5
 steps = Int(round(total_time / dt))
 
 set = se("system_ram.yaml")
-set.segments = 2
+set.segments = 7
 set_values = [-50, -1.1, -1.1]
+set.quasi_static = true
 
 if !@isdefined s
     wing = RamAirWing(set)
@@ -28,16 +29,19 @@ end
 s.set.abs_tol = 1e-5
 s.set.rel_tol = 1e-3
 
+measure.sphere_pos .= deg2rad.([50.0 50.0; 1.0 -1.0])
 if !ispath(joinpath(get_data_path(), "prob.bin"))
     KiteModels.init_sim!(s, measure)
 end
-measure.sphere_pos .= deg2rad.([50.0 50.0; 1.0 -1.0])
-@time KiteModels.reinit!(s, measure)
+@time KiteModels.reinit!(s, measure; reload=true)
+sys = s.sys
+s.integrator.ps[sys.steady] = true
+next_step!(s; dt=10.0, vsm_interval=1)
+s.integrator.ps[sys.steady] = false
 
 logger = Logger(length(s.point_system.points), steps)
 
 sys_state = KiteModels.SysState(s)
-sys = s.sys
 l = s.set.l_tether + 10
 t = 0.
 runtime = 0.
@@ -47,8 +51,9 @@ try
         global t, runtime, integ_runtime
         PLOT && plot(s, t; zoom=false, front=false)
         global set_values = -s.set.drum_radius .* s.integrator[sys.winch_force]
-        if (t < 3.0); set_values -= [0, 0, 5]; end
+        if (t < 1.5); set_values += [0, 10, 5]; end
         steptime = @elapsed (t, integ_steptime) = next_step!(s, set_values; dt, vsm_interval)
+        t -= 10.0
         if (t > total_time/2); runtime += steptime; end
         if (t > total_time/2); integ_runtime += integ_steptime; end
         KiteModels.update_sys_state!(sys_state, s)
@@ -63,10 +68,17 @@ try
         sys_state.var_07 = s.integrator[sys.aero_moment_b[2]]
         sys_state.var_08 = s.integrator[sys.group_aero_moment[1]]
 
-        sys_state.var_09 = s.integrator[sys.twist_angle[1]]
-        sys_state.var_10 = s.integrator[sys.twist_angle[2]]
-        sys_state.var_11 = s.integrator[sys.twist_angle[3]]
-        sys_state.var_12 = s.integrator[sys.twist_angle[4]]
+        sys_state.var_09 = s.integrator[sys.twist_α[1]]
+        sys_state.var_10 = s.integrator[sys.twist_α[2]]
+        sys_state.var_11 = s.integrator[sys.twist_α[3]]
+        sys_state.var_12 = s.integrator[sys.twist_α[4]]
+
+        sys_state.var_13 = clamp(s.integrator[sys.pulley_acc[1]], -100, 100)
+        sys_state.var_14 = clamp(s.integrator[sys.pulley_acc[2]], -100, 100)
+
+        va_kite_b = s.integrator[sys.va_kite_b]
+        e_x = s.integrator[sys.e_x]
+        sys_state.var_15 = rad2deg(acos(dot(normalize(va_kite_b), e_x)))
 
         log!(logger, sys_state)
     end
@@ -83,14 +95,18 @@ p=plotx(logger.time_vec,
         [logger.var_04_vec, logger.var_05_vec],
         [logger.var_06_vec, logger.var_07_vec, logger.var_08_vec],
         [logger.var_09_vec, logger.var_10_vec, logger.var_11_vec, logger.var_12_vec],
+        [logger.var_13_vec, logger.var_14_vec],
+        [logger.var_15_vec],
         [logger.heading_vec]
         ;
-    ylabels=["kite", "tether", "vsm", "twist", "heading"], 
+    ylabels=["kite", "tether", "vsm", "twist", "pulley", "wind", "heading"], 
     labels=[
         ["ω_b[1]", "ω_b[2]", "ω_b[3]"],
         ["vel[1]", "vel[2]"],
         ["force[3]", "kite moment[2]", "group moment[1]"],
-        ["twist_angle[1]", "twist_angle[2]", "twist_angle[3]", "twist_angle[4]"],
+        ["twist_alpha[1]", "twist_alpha[2]", "twist_alpha[3]", "twist_alpha[4]"],
+        ["pulley acc[1]", "pulley acc[2]"],
+        ["alpha"],
         ["heading"]
         ],
     fig="Steering and heading MTK model")
@@ -99,4 +115,13 @@ display(p)
 println("Times realtime: ", (total_time/2) / runtime)
 println("Times realtime, just integrator: ", (total_time/2) / integ_runtime)
 
+@show norm(logger.var_13_vec[steps÷2:end])
+
+# diffs = [norm(diff([u[i] for u in s.integrator.sol.u[1000:end]])) for i in eachindex(s.integrator.u)]
+# names = unknowns(s.prob.f.sys)
+# for (n, d) in zip(names, diffs)
+#     println(round(d; digits=2), "\t", n)
+# end
+
+
 nothing

src/mtk_model.jl

@@ -189,7 +189,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             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
-            @parameters bridle_damp = 20
+            @parameters bridle_damp = 100
             @parameters measured_ω_z = 0.6
             eqs = [
                 eqs
@@ -250,7 +250,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
 
         inertia = 1/3 * (s.set.mass/length(groups)) * (norm(group.chord))^2 # plate inertia around leading edge
         @assert !(inertia ≈ 0.0)
-        @parameters twist_damp = 20
+        @parameters twist_damp = s.set.quasi_static ? 200 : 100
         eqs = [
             eqs
             group_tether_moment[group.idx] ~ sum(tether_moment[group.idx, :])
@@ -261,7 +261,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             eqs = [
                 eqs
                 D(free_twist_angle[group.idx]) ~ twist_ω[group.idx]
-                D(twist_ω[group.idx]) ~ acc_multiplier * twist_α[group.idx] - twist_damp * twist_ω[group.idx]
+                D(twist_ω[group.idx]) ~ twist_α[group.idx] - twist_damp * twist_ω[group.idx]
             ]
             defaults = [
                 defaults
@@ -377,8 +377,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         @parameters stiffness_frac = 0.1
         (segment.type == BRIDLE) && (stiffness_m = stiffness_frac * stiffness_m)
 
-        @parameters set_damping = s.set.damping
-        damping_m = (set_damping / s.set.c_spring) * stiffness_m
+        damping_m = (s.set.damping / s.set.c_spring) * stiffness_m
 
         eqs = [
             eqs
@@ -418,7 +417,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         pulley_force(t)[eachindex(pulleys)]
         pulley_acc(t)[eachindex(pulleys)]
     end
-    @parameters pulley_damping = 10
+    @parameters pulley_damp = 20
     for pulley in pulleys
         segment = segments[pulley.segments[1]]
         mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2
@@ -431,9 +430,9 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         ]
         if pulley.type == DYNAMIC
             eqs = [
-                eqs 
+                eqs
                 D(pulley_l0[pulley.idx])  ~ pulley_vel[pulley.idx]
-                D(pulley_vel[pulley.idx]) ~ acc_multiplier * pulley_acc[pulley.idx] - pulley_damping * pulley_vel[pulley.idx]
+                D(pulley_vel[pulley.idx]) ~ acc_multiplier * pulley_acc[pulley.idx] - pulley_damp * pulley_vel[pulley.idx]
             ]
             defaults = [
                 defaults
@@ -506,6 +505,7 @@ function diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero
         total_tether_kite_force(t)[1:3]
         total_tether_kite_moment(t)[1:3]
     end
+    @parameters ω_damp = 150
 
     Ω = [0       -ω_b[1]  -ω_b[2]  -ω_b[3];
         ω_b[1]    0        ω_b[3]  -ω_b[2];
@@ -521,7 +521,7 @@ function diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero
         [R_b_w[:, i] ~ quaternion_to_rotation_matrix(Q_b_w)[:, i] for i in 1:3]
 
         D(ω_b[1]) ~ α_b[1]
-        D(ω_b[2]) ~ α_b[2]
+        D(ω_b[2]) ~ α_b[2] - ω_damp * ω_b[2]
         D(ω_b[3]) ~ ifelse(steady==true, 0, α_b[3])
 
         α_b[1] ~ (moment_b[1] + (I_b[2] - I_b[3]) * ω_b[2] * ω_b[3]) / I_b[1]
@@ -653,7 +653,7 @@ function init_unknowns_vec!(
     init_kite_pos;
     non_observed=true
 )
-    non_observed && (length(vec) != length(s.integrator.u)) && 
+    !s.set.quasi_static && non_observed && (length(vec) != length(s.integrator.u)) && 
         throw(ArgumentError("Unknowns of length $(length(s.integrator.u)) but vector provided of length $(length(vec))"))
 
     @unpack points, groups, segments, pulleys, winches = system
@@ -724,34 +724,28 @@ function get_unknowns(s::RamAirKite)
     vec = Num[]
     vec = get_stiff_unknowns(s, vec)
     vec = get_nonstiff_unknowns(s, vec)
-    (length(vec) != length(s.integrator.u)) &&
+    !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))"))
     return vec
 end
 
 function get_stiff_unknowns(s, vec=Num[])
     @unpack points, groups, segments, pulleys, winches = s.point_system
     for point in points
-        if point.type == DYNAMIC
-            for i in 1:3
-                push!(vec, s.sys.pos[i, point.idx])
-            end
-            for i in 1:3 # TODO: add speed to init
-                push!(vec, s.sys.vel[i, point.idx])
-            end
+        for i in 1:3
+            point.type == DYNAMIC && push!(vec, s.sys.pos[i, point.idx])
+        end
+        for i in 1:3 # TODO: add speed to init
+            point.type == DYNAMIC && push!(vec, s.sys.vel[i, point.idx])
         end
     end
     for group in groups
-        if group.type == DYNAMIC
-            push!(vec, s.sys.free_twist_angle[group.idx])
-            push!(vec, s.sys.twist_ω[group.idx])
-        end
+        group.type == DYNAMIC && push!(vec, s.sys.free_twist_angle[group.idx])
+        group.type == DYNAMIC && push!(vec, s.sys.twist_ω[group.idx])
     end
     for pulley in pulleys
-        if pulley.type == DYNAMIC
-            push!(vec, s.sys.pulley_l0[pulley.idx])
-            push!(vec, s.sys.pulley_vel[pulley.idx])
-        end
+        pulley.type == DYNAMIC && push!(vec, s.sys.pulley_l0[pulley.idx])
+        pulley.type == DYNAMIC && push!(vec, s.sys.pulley_vel[pulley.idx])
     end
     return vec
 end

src/point_mass_system.jl

@@ -22,6 +22,8 @@ function PointMassSystem(set::Settings, wing::RamAirWing)
     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)]
 
+    dynamics_type = set.quasi_static ? STATIC : DYNAMIC
+
     function create_bridle(bridle_top, gammas)
         i_pnt = length(points) # last point idx
         i_seg = length(segments) # last segment idx
@@ -47,18 +49,18 @@ function PointMassSystem(set::Settings, wing::RamAirWing)
 
         points = [
             points
-            Point(9+i_pnt, bridle_top[1], DYNAMIC)
-            Point(10+i_pnt, bridle_top[2], DYNAMIC)
-            Point(11+i_pnt, bridle_top[3], DYNAMIC)
-            Point(12+i_pnt, bridle_top[4], DYNAMIC)
+            Point(9+i_pnt, bridle_top[1], dynamics_type)
+            Point(10+i_pnt, bridle_top[2], dynamics_type)
+            Point(11+i_pnt, bridle_top[3], dynamics_type)
+            Point(12+i_pnt, bridle_top[4], dynamics_type)
 
-            Point(13+i_pnt, bridle_top[2] .+ [0, 0, -1], DYNAMIC)
+            Point(13+i_pnt, bridle_top[2] .+ [0, 0, -1], dynamics_type)
 
-            Point(14+i_pnt, bridle_top[1] .+ [0, 0, -2], DYNAMIC)
-            Point(15+i_pnt, bridle_top[3] .+ [0, 0, -2], DYNAMIC)
+            Point(14+i_pnt, bridle_top[1] .+ [0, 0, -2], dynamics_type)
+            Point(15+i_pnt, bridle_top[3] .+ [0, 0, -2], dynamics_type)
 
-            Point(16+i_pnt, bridle_top[1] .+ [0, 0, -3], DYNAMIC)
-            Point(17+i_pnt, bridle_top[3] .+ [0, 0, -3], DYNAMIC)
+            Point(16+i_pnt, bridle_top[1] .+ [0, 0, -3], dynamics_type)
+            Point(17+i_pnt, bridle_top[3] .+ [0, 0, -3], dynamics_type)
         ]
         l1 = norm(points[9+i_pnt].pos_b - points[1+i_pnt].pos_b)
         l2 = norm(points[9+i_pnt].pos_b - points[5+i_pnt].pos_b)
@@ -88,8 +90,8 @@ function PointMassSystem(set::Settings, wing::RamAirWing)
         ]
         pulleys = [
             pulleys
-            Pulley(1+i_pul, (13+i_seg, 14+i_seg), DYNAMIC)
-            Pulley(2+i_pul, (16+i_seg, 17+i_seg), DYNAMIC)
+            Pulley(1+i_pul, (13+i_seg, 14+i_seg), dynamics_type)
+            Pulley(2+i_pul, (16+i_seg, 17+i_seg), dynamics_type)
         ]
         push!(attach_points, points[end-1])
         push!(attach_points, points[end])
@@ -107,13 +109,13 @@ function PointMassSystem(set::Settings, wing::RamAirWing)
             i_pnt = length(points) # last point idx
             i_seg = length(segments) # last segment idx
             if i == 1
-                points = [points; Point(1+i_pnt, pos, DYNAMIC)]
+                points = [points; Point(1+i_pnt, pos, dynamics_type)]
                 segments = [segments; Segment(1+i_seg, (attach_point.idx, 1+i_pnt), 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)]
             else
-                points = [points; Point(1+i_pnt, pos, DYNAMIC)]
+                points = [points; Point(1+i_pnt, pos, dynamics_type)]
                 segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)]
             end
             push!(segment_idxs, 1+i_seg)

src/ram_air_kite.jl

@@ -400,7 +400,7 @@
     return nothing
 end
 
-function init_sim!(s::RamAirKite; prn=true)
+function init_sim!(::RamAirKite; prn=true)
     throw(ArgumentError("Use the function init_sim!(s::RamAirKite, measure::Measurement) instead."))
 end
 
@@ -448,9 +448,8 @@
             s.prob = deserialize(prob_path)
             s.sys = s.prob.f.sys
             s.integrator = OrdinaryDiffEqCore.init(s.prob, solver; dt, abstol=s.set.abs_tol, reltol=s.set.rel_tol, save_on=false, save_everystep=false)
-            sym_vec = zeros(Num, length(s.integrator.u))
-            s.unknowns_vec = zeros(SimFloat, length(s.integrator.u))
             sym_vec = get_unknowns(s)
+            s.unknowns_vec = zeros(SimFloat, length(sym_vec))
             generate_getters!(s, sym_vec)
         end
         prn && @info "Loaded problem from $prob_path and initialized integrator in $t seconds"

test/test_ram_air_kite.jl

@@ -225,7 +225,7 @@ const BUILD_SYS = true
             right_heading_diff = angle_diff(sys_state_right.heading, sys_state_initial.heading)
             @test right_heading_diff > 0.6
             left_heading_diff = angle_diff(sys_state_left.heading, sys_state_initial.heading)
-            @test left_heading_diff < -0.6
+            @test left_heading_diff < -0.5
             @test abs(right_heading_diff) ≈ abs(left_heading_diff) atol=0.3
         end
     end