Benchmark simplify operation

src/mtk_model.jl

@@ -876,7 +876,7 @@ function linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aer
     return eqs, guesses
 end
 
-function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b)
+function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b, bench=false)
     eqs = []
     defaults = Pair{Num, Any}[]
     guesses = Pair{Num, Any}[]
@@ -916,31 +916,15 @@ function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b)
     eqs, defaults = wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero_force_b, aero_moment_b, 
         ω_b, α_b, R_b_w, wing_pos, wing_vel, wing_acc, stabilize, fix_nonstiff)
     eqs = scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b)
-
-    # te_I = (1/3 * (s.set.mass/8) * te_length^2)
-    # # -damping / I * ω = α_damping
-    # # solve for c: (c * (k*m/s^2) / (k*m^2)) * (m/s)=m/s^2 in wolframalpha
-    # # damping should be in N*m*s
-    # rot_damping = 0.1s.damping * te_length
-
-    # eqs = [
-    #     eqs
-    #     trailing_edge_α[1] ~ (force[:, s.i_A]) ⋅ e_te_A * te_length / te_I - (rot_damping[1] / te_I) * trailing_edge_ω[1] # TODO: add trailing edge
-    #     trailing_edge_α[2] ~ (force[:, s.i_B]) ⋅ e_te_B * te_length / te_I - (rot_damping[2] / te_I) * trailing_edge_ω[2]
-    # ]
-
+
     eqs = Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs)))
 
-    # discrete_events = [
-    #     true => [
-    #         [Q_b_w[i] ~ normalize(Q_b_w)[i] for i in 1:4]
-    #         [twist_angle[i] ~ clamp(twist_angle[i], -π/2, π/2) for i in eachindex(s.point_groups)]
-    #         ]
-    #     ]
-
-    @info "Creating ODESystem"
-    # @named sys = ODESystem(eqs, t; discrete_events)
-    @time @named sys = ODESystem(eqs, t)
+    ! bench && @info "Creating ODESystem"
+    if bench
+        @named sys = ODESystem(eqs, t)
+    else
+        @time @named sys = ODESystem(eqs, t)
+    end
 
     defaults = [
         defaults

src/symbolic_awe_model.jl

@@ -301,7 +301,7 @@ and only update the state variables. Otherwise, it will create a new model from
 """
 function init!(s::SymbolicAWEModel; 
     solver=nothing, stiffness_factor = nothing, delta = nothing, adaptive=true, prn=true, 
-    precompile=false, remake=false, reload=false, 
+    precompile=false, remake=false, reload=false, bench=false,
     lin_outputs=Num[]
 )
     if isnothing(solver)
@@ -322,11 +322,16 @@ function init!(s::SymbolicAWEModel;
         init_Q_b_w, R_b_w, init_va_b = initial_orient(s)
         init!(s.sys_struct, s.set)
 
-        inputs = create_sys!(s, s.sys_struct; init_va_b)
+        inputs = create_sys!(s, s.sys_struct; init_va_b, bench)
         prn && @info "Simplifying the system"
         @suppress_err begin
-            prn ? (@time (sys, _) = structural_simplify(s.full_sys, (inputs, []))) :
-                ((sys, _) = structural_simplify(s.full_sys, (inputs, [])))
+            if prn && !bench
+                @time (sys, _) = structural_simplify(s.full_sys, (inputs, []))
+            elseif bench
+                return @elapsed structural_simplify(s.full_sys, (inputs, []))
+            else
+                (sys, _) = structural_simplify(s.full_sys, (inputs, []))
+            end
             s.sys = sys
         end
         dt = SimFloat(1/s.set.sample_freq)
@@ -353,13 +358,19 @@ function init!(s::SymbolicAWEModel;
     end
     model_path = joinpath(KiteUtils.get_data_path(), get_model_name(s.set; precompile))
     if !ispath(model_path) || remake
-        init(s)
+        res = init(s)
+        if bench
+            return res
+        end
     end
     _, success = reinit!(s, solver; adaptive, precompile, reload, lin_outputs, prn)
     if !success
         rm(model_path)
         @info "Rebuilding the system. This can take some minutes..."
-        init(s)
+        res = init(s)
+        if bench
+            return res
+        end
         reinit!(s, solver; adaptive, precompile, lin_outputs, prn, reload=true)
     end
     return s.integrator

test/bench_simplify.jl

@@ -0,0 +1,48 @@
+# Copyright (c) 2024, 2025 Bart van de Lint, Uwe Fechner
+# SPDX-License-Identifier: MIT
+
+SIMPLE = false
+
+using Pkg
+if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
+    using TestEnv; TestEnv.activate()
+end
+using KiteModels, LinearAlgebra, Statistics
+
+# Simulation parameters
+dt = 0.05
+total_time = 10.0  # Longer simulation to see oscillations
+vsm_interval = 3
+steps = Int(round(total_time / dt))
+
+# Steering parameters
+steering_freq = 1/2  # Hz - full left-right cycle frequency
+steering_magnitude = 10.0      # Magnitude of steering input [Nm]
+
+# Initialize model
+set = load_settings("system_ram.yaml")
+set.segments = 3
+set_values = [-50, 0.0, 0.0]  # Set values of the torques of the three winches. [Nm]
+set.quasi_static = false
+set.physical_model = SIMPLE ? "simple_ram" : "ram"
+
+sam = SymbolicAWEModel(set)
+sam.set.abs_tol = 1e-2
+sam.set.rel_tol = 1e-2
+rm("data/model_1.11_ram_dynamic_3_seg.bin"; force=true)
+
+# Initialize at elevation
+set.l_tethers[2] += 0.2
+set.l_tethers[3] += 0.2
+time = init!(sam; remake=false, reload=true, bench=true)
+@info "Simplify took $time seconds"
+sys = sam.sys
+nothing
+
+# Desktop, AMD Ryzen 9 7950X, Julia 1.11:
+# - first  run 34.5 seconds
+# - second run 21.1 seconds
+
+# Laptop, AMD Ryzen 7 7840U, Julia 1.11:
+# - first  run 36.5 seconds
+# - second run 24.3 seconds