StatProfilerHTML.jl report
Generated on Mon, 01 Apr 2024 21:01:18
File source code
Line Exclusive Inclusive Code
1 # SPDX-License-Identifier: MIT
2
3 #= MIT License
4
5 Copyright (c) 2020, 2021, 2022 Uwe Fechner
6
7 Permission is hereby granted, free of charge, to any person obtaining a copy
8 of this software and associated documentation files (the "Software"), to deal
9 in the Software without restriction, including without limitation the rights
10 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 copies of the Software, and to permit persons to whom the Software is
12 furnished to do so, subject to the following conditions:
13
14 The above copyright notice and this permission notice shall be included in all
15 copies or substantial portions of the Software.
16
17 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23 SOFTWARE. =#
24
25 #= Models of a kite-power system in implicit form: residual = f(y, yd)
26
27 This model implements a 3D mass-spring system with reel-out. It uses six tether segments (the number can be
28 configured in the file data/settings.yaml). Two kite models are provided, the one point and the four point
29 kite model. The spring constant and the damping decrease with the segment length. The aerodynamic kite forces are
30 calculated, depending on reel-out speed, depower and steering settings.
31
32 Scientific background: http://arxiv.org/abs/1406.6218 =#
33
34 module KiteModels
35
36 using PrecompileTools: @setup_workload, @compile_workload
37 using Dierckx, StaticArrays, Rotations, LinearAlgebra, Parameters, NLsolve, DocStringExtensions, Sundials
38 using Reexport
39 @reexport using KitePodModels
40 @reexport using WinchModels
41 @reexport using AtmosphericModels
42 import Base.zero
43 import KiteUtils.calc_elevation
44 import KiteUtils.calc_azimuth
45 import KiteUtils.calc_heading
46 import KiteUtils.calc_course
47 import KiteUtils.SysState
48 import Sundials.init
49 import Sundials.step!
50
51 export KPS3, KPS4, KVec3, SimFloat, ProfileLaw, EXP, LOG, EXPLOG # constants and types
52 export calc_set_cl_cd!, copy_examples, copy_bin, update_sys_state! # helper functions
53 export clear!, find_steady_state!, residual! # low level workers
54 export init_sim!, next_step! # high level workers
55 export pos_kite, calc_height, calc_elevation, calc_azimuth, calc_heading, calc_course # getters
56 export winch_force, lift_drag, lift_over_drag, unstretched_length, tether_length, v_wind_kite # getters
57 export kite_ref_frame, orient_euler, spring_forces
58 import LinearAlgebra: norm
59
60 set_zero_subnormals(true) # required to avoid drastic slow down on Intel CPUs when numbers become very small
61
62 # Constants
63 const G_EARTH = 9.81 # gravitational acceleration
64 const BRIDLE_DRAG = 1.1 # should probably be removed
65
66 # Type definitions
67 """
68 const SimFloat = Float64
69
70 This type is used for all real variables, used in the Simulation. Possible alternatives: Float32, Double64, Dual
71 Other types than Float64 or Float32 do require support of Julia types by the solver.
72 """
73 const SimFloat = Float64
74
75 """
76 const KVec3 = MVector{3, SimFloat}
77
78 Basic 3-dimensional vector, stack allocated, mutable.
79 """
80 const KVec3 = MVector{3, SimFloat}
81
82 """
83 const SVec3 = SVector{3, SimFloat}
84
85 Basic 3-dimensional vector, stack allocated, immutable.
86 """
87 const SVec3 = SVector{3, SimFloat}
88
89 # the following two definitions speed up the function residual! from 940ns to 540ns
90 # disadvantage: changing the cl and cd curves requires a restart of the program
91 const calc_cl = Spline1D(se().alpha_cl, se().cl_list)
92 const calc_cd = Spline1D(se().alpha_cd, se().cd_list)
93
94 """
95 abstract type AbstractKiteModel
96
97 All kite models must inherit from this type. All methods that are defined on this type must work
98 with all kite models. All exported methods must work on this type.
99 """
100 abstract type AbstractKiteModel end
101
102 """
103 const AKM = AbstractKiteModel
104
105 Short alias for the AbstractKiteModel.
106 """
107 const AKM = AbstractKiteModel
108
109 function __init__()
110 if isdir(joinpath(pwd(), "data")) && isfile(joinpath(pwd(), "data", "system.yaml"))
111 set_data_path(joinpath(pwd(), "data"))
112 end
113 end
114
115 include("KPS4.jl") # include code, specific for the four point kite model
116 include("KPS3.jl") # include code, specific for the one point kite model
117 include("init.jl") # functions to calculate the inital state vector, the inital masses and initial springs
118
119 # Calculate the lift and drag coefficient as a function of the angle of attack alpha.
120 function set_cl_cd!(s::AKM, alpha)
121 angle = alpha * 180.0 / π
122 if angle > 180.0
123 angle -= 360.0
124 end
125 if angle < -180.0
126 angle += 360.0
127 end
128 s.param_cl = calc_cl(angle)
129 s.param_cd = calc_cd(angle)
130 nothing
131 end
132
133 # Calculate the angle of attack alpha from the apparend wind velocity vector
134 # v_app and the z unit vector of the kite reference frame.
135 function calc_alpha(v_app, vec_z)
136 π/2.0 - acos(-(v_app ⋅ vec_z) / norm(v_app))
137 end
138
139 """
140 set_depower_steering!(s::AKM, depower, steering)
141
142 Setter for the depower and steering model inputs.
143
144 Parameters:
145 - depower: Relative depower, must be between 0 .. 1.0
146 - steering: Relative steering, must be between -1.0 .. 1.0.
147
148 This function sets the variables s.depower, s.steering and s.alpha_depower.
149
150 It takes the depower offset c0 and the dependency of the steering sensitivity from
151 the depower settings into account.
152 """
153 function set_depower_steering!(s::AKM, depower, steering)
154 s.depower = depower
155 s.alpha_depower = calc_alpha_depower(s.kcu, depower)
156 s.steering = (steering - s.set.c0) / (1.0 + s.set.k_ds * (s.alpha_depower / deg2rad(s.set.alpha_d_max)))
157 nothing
158 end
159
160
161 """
162 set_v_reel_out!(s::AKM, v_reel_out, t_0, period_time = 1.0 / s.set.sample_freq)
163
164 Setter for the reel-out speed. Must be called on every timestep (before each simulation).
165 It also updates the tether length, therefore it must be called even if `v_reel_out` has
166 not changed.
167
168 - t_0 the start time of the next timestep relative to the start of the simulation [s]
169 """
170 function set_v_reel_out!(s::AKM, v_reel_out, t_0, period_time = 1.0 / s.set.sample_freq)
171 s.sync_speed = v_reel_out
172 s.last_v_reel_out = s.v_reel_out
173 s.t_0 = t_0
174 end
175
176 """
177 unstretched_length(s::AKM)
178
179 Getter for the unstretched tether reel-out lenght (at zero force).
180 """
181 function unstretched_length(s::AKM) s.l_tether end
182
183 """
184 lift_drag(s::AKM)
185
186 Return a tuple of the scalar lift and drag forces.
187
188 **Example:**
189
190 lift, drag = lift_drag(s)
191 """
192 function lift_drag(s::AKM) return (norm(s.lift_force), norm(s.drag_force)) end
193
194 """
195 lift_over_drag(s::AKM)
196
197 Return the lift-over-drag ratio.
198 """
199 function lift_over_drag(s::AKM)
200 lift, drag = lift_drag(s)
201 return lift / drag
202 end
203
204 """
205 v_wind_kite(s::AKM)
206
207 Return the vector of the wind speed at the height of the kite.
208 """
209 function v_wind_kite(s::AKM) s.v_wind end
210
211 """
212 set_v_wind_ground!(s::AKM, height, v_wind_gnd=s.set.v_wind, wind_dir=0.0)
213
214 Set the vector of the wind-velocity at the height of the kite. As parameter the height,
215 the ground wind speed [m/s] and the wind direction [radians] are needed.
216 Must be called every at each timestep.
217 """
218 function set_v_wind_ground!(s::AKM, height, v_wind_gnd=s.set.v_wind, wind_dir=0.0)
219 if height < 6.0
220 height = 6.0
221 end
222 s.v_wind .= v_wind_gnd * calc_wind_factor(s.am, height) .* [cos(wind_dir), sin(wind_dir), 0]
223 s.v_wind_gnd .= [v_wind_gnd * cos(wind_dir), v_wind_gnd * sin(wind_dir), 0.0]
224 s.v_wind_tether .= v_wind_gnd * calc_wind_factor(s.am, height / 2.0) .* [cos(wind_dir), sin(wind_dir), 0]
225 s.rho = calc_rho(s.am, height)
226 nothing
227 end
228
229 """
230 tether_length(s::AKM)
231
232 Calculate and return the real, stretched tether lenght.
233 """
234 function tether_length(s::AKM)
235 length = 0.0
236 for i in 1:s.set.segments
237 length += norm(s.pos[i+1] - s.pos[i])
238 end
239 return length
240 end
241
242 """
243 orient_euler(s::AKM)
244
245 Calculate and return the orientation of the kite in euler angles (roll, pitch, yaw)
246 as SVector.
247 """
248 function orient_euler(s::AKM)
249 x, y, z = kite_ref_frame(s)
250 roll = atan(y[3], z[3]) - π/2
251 if roll < -π/2
252 roll += 2π
253 end
254 pitch = asin(-x[3])
255 yaw = -atan(x[2], x[1]) - π/2
256 if yaw < -π/2
257 yaw += 2π
258 end
259 SVector(roll, pitch, yaw)
260 end
261
262 """
263 calc_elevation(s::AKM)
264
265 Determine the elevation angle of the kite in radian.
266 """
267 function calc_elevation(s::AKM)
268 KiteUtils.calc_elevation(pos_kite(s))
269 end
270
271 """
272 calc_azimuth(s::AKM)
273
274 Determine the azimuth angle of the kite in radian.
275 """
276 function calc_azimuth(s::AKM)
277 KiteUtils.azimuth_east(pos_kite(s))
278 end
279
280 """
281 calc_heading(s::AKM)
282
283 Determine the heading angle of the kite in radian.
284 """
285 function calc_heading(s::AKM)
286 orientation = orient_euler(s)
287 elevation = calc_elevation(s)
288 azimuth = calc_azimuth(s)
289 KiteUtils.calc_heading(orientation, elevation, azimuth)
290 end
291
292 """
293 calc_course(s::AKM)
294
295 Determine the course angle of the kite in radian.
296 Undefined if the velocity of the kite is near zero.
297 """
298 function calc_course(s::AKM)
299 elevation = calc_elevation(s)
300 azimuth = calc_azimuth(s)
301 KiteUtils.calc_course(s.vel_kite, elevation, azimuth)
302 end
303
304 # mutable struct SysState{P}
305 # "time since start of simulation in seconds"
306 # time::Float64
307 # "time needed for one simulation timestep"
308 # t_sim::Float64
309 # "state of system state control"
310 # sys_state::Int16
311 # "mechanical energy [Wh]"
312 # e_mech::Float64
313 # "orientation of the kite (quaternion, order w,x,y,z)"
314 # orient::MVector{4, Float32}
315 # "elevation angle in radians"
316 # elevation::MyFloat
317 # "azimuth angle in radians"
318 # azimuth::MyFloat
319 # "tether length [m]"
320 # l_tether::MyFloat
321 # "reel out velocity [m/s]"
322 # v_reelout::MyFloat
323 # "tether force [N]"
324 # force::MyFloat
325 # "depower settings [0..1]"
326 # depower::MyFloat
327 # "steering settings [-1..1]"
328 # steering::MyFloat
329 # "heading angle in radian"
330 # heading::MyFloat
331 # "course angle in radian"
332 # course::MyFloat
333 # "norm of apparent wind speed [m/s]"
334 # v_app::MyFloat
335 # "velocity vector of the kite"
336 # vel_kite::MVector{3, MyFloat}
337 # "vector of particle positions in x"
338 # X::MVector{P, MyFloat}
339 # "vector of particle positions in y"
340 # Y::MVector{P, MyFloat}
341 # "vector of particle positions in z"
342 # Z::MVector{P, MyFloat}
343 # var_01::MyFloat
344 # var_02::MyFloat
345 # var_03::MyFloat
346 # var_04::MyFloat
347 # var_05::MyFloat
348 # end
349
350 function update_sys_state!(ss::SysState, s::AKM, zoom=1.0)
351 ss.time = s.t_0
352 pos = s.pos
353 P = length(pos)
354 for i in 1:P
355 ss.X[i] = pos[i][1] * zoom
356 ss.Y[i] = pos[i][2] * zoom
357 ss.Z[i] = pos[i][3] * zoom
358 end
359 x, y, z = kite_ref_frame(s)
360 pos_kite_ = pos_kite(s)
361 pos_before = pos_kite_ + z
362
363 rotation = rot(pos_kite_, pos_before, -x)
364 q = QuatRotation(rotation)
365 ss.orient .= Rotations.params(q)
366 ss.elevation = calc_elevation(s)
367 ss.azimuth = calc_azimuth(s)
368 ss.force = winch_force(s)
369 ss.heading = calc_heading(s)
370 ss.course = calc_course(s)
371 ss.v_app = norm(s.v_apparent)
372 ss.l_tether = s.l_tether
373 ss.v_reelout = s.v_reel_out
374 ss.depower = s.depower
375 ss.steering = s.steering
376 ss.vel_kite .= s.vel_kite
377 nothing
378 end
379
380 """
381 SysState(s::AKM, zoom=1.0)
382
383 Constructor for creating a SysState object from a kite model (KPS3 or KPS4).
384 The SysState object can be used either for logging or for displaying the
385 system state in a viewer. Optionally the position arrays can be zoomed
386 according to the requirements of the viewer.
387 """
388 function SysState(s::AKM, zoom=1.0)
389 pos = s.pos
390 P = length(pos)
391 X = zeros(MVector{P, MyFloat})
392 Y = zeros(MVector{P, MyFloat})
393 Z = zeros(MVector{P, MyFloat})
394 for i in 1:P
395 X[i] = pos[i][1] * zoom
396 Y[i] = pos[i][2] * zoom
397 Z[i] = pos[i][3] * zoom
398 end
399
400 x, y, z = kite_ref_frame(s)
401 pos_kite_ = pos_kite(s)
402 pos_before = pos_kite_ + z
403
404 rotation = rot(pos_kite_, pos_before, -x)
405 q = QuatRotation(rotation)
406 orient = MVector{4, Float32}(Rotations.params(q))
407
408 elevation = calc_elevation(s)
409 azimuth = calc_azimuth(s)
410 force = winch_force(s)
411 heading = calc_heading(s)
412 course = calc_course(s)
413 v_app_norm = norm(s.v_apparent)
414 t_sim = 0
415 KiteUtils.SysState{P}(s.t_0, t_sim, 0, 0, orient, elevation, azimuth, s.l_tether, s.v_reel_out, force, s.depower, s.steering,
416 heading, course, v_app_norm, s.vel_kite, X, Y, Z, 0, 0, 0, 0, 0)
417 end
418
419 function calc_pre_tension(s::AKM)
420 forces = spring_forces(s)
421 av_force = 0.0
422 for i in 1:s.set.segments
423 av_force += forces[i]
424 end
425 av_force /= s.set.segments
426 res = av_force/s.set.c_spring
427 if res < 0.0 res = 0.0 end
428 if isnan(res) res = 0.0 end
429 return res + 1.0
430 end
431
432 """
433 init_sim!(s; t_end=1.0, stiffness_factor=0.035, prn=false)
434
435 Initialises the integrator of the model.
436
437 Parameters:
438 - s: an instance of an abstract kite model
439 - t_end: end time of the simulation; normally not needed
440 - stiffness_factor: factor applied to the tether stiffness during initialisation
441 - prn: if set to true, print the detailed solver results
442
443 Returns:
444 An instance of a DAE integrator.
445 """
446 function init_sim!(s::AKM; t_end=1.0, stiffness_factor=0.035, prn=false)
447 clear!(s)
448 s.stiffness_factor = stiffness_factor
449 y0, yd0 = KiteModels.find_steady_state!(s; stiffness_factor=stiffness_factor, prn=prn)
450 y0 = Vector{Float64}(y0)
451 yd0 = Vector{Float64}(yd0)
452 differential_vars = ones(Bool, length(y0))
453 solver = IDA(linear_solver=Symbol(s.set.linear_solver), max_order = s.set.max_order)
454 tspan = (0.0, t_end)
455 abstol = s.set.abs_tol # max error in m/s and m
456 prob = DAEProblem{true}(residual!, yd0, y0, tspan, s, differential_vars=differential_vars)
457 integrator = Sundials.init(prob, solver, abstol=abstol, reltol=s.set.rel_tol)
458 end
459
460 """
461 next_step!(s::AKM, integrator; v_ro = 0.0, v_wind_gnd=s.set.v_wind, wind_dir=0.0, dt=1/s.set.sample_freq)
462
463 Calculates the next simulation step.
464
465 Parameters:
466 - s: an instance of an abstract kite model
467 - integrator: an integrator instance as returned by the function [`init_sim!`](@ref)
468 - v_ro: set value of reel out speed in m/s
469 - `v_wind_gnd`: wind speed at reference height in m/s
470 - wind_dir: wind direction in radians
471 - dt: time step in seconds
472
473 Only the first two parameters are required.
474
475 Returns:
476 The end time of the time step in seconds.
477 """
478 58 (100 %)
116 (200 %) samples spent in next_step!
58 (50 %) (incl.) when called from next_step! line 478
58 (50 %) (incl.) when called from nsteps line 9
58 (100 %) samples spent calling #next_step!#40
function next_step!(s::AKM, integrator; v_ro = 0.0, v_wind_gnd=s.set.v_wind, wind_dir=0.0, dt=1/s.set.sample_freq)
479 KitePodModels.on_timer(s.kcu)
480 KiteModels.set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu))
481 set_v_reel_out!(s, v_ro, integrator.t)
482 set_v_wind_ground!(s, calc_height(s), v_wind_gnd, wind_dir)
483 58 (100 %)
58 (100 %) samples spent calling step!
Sundials.step!(integrator, dt, true)
484 if s.stiffness_factor < 1.0
485 s.stiffness_factor+=0.01
486 if s.stiffness_factor > 1.0
487 s.stiffness_factor = 1.0
488 end
489 end
490 integrator.t
491 end
492
493 """
494 copy_examples()
495
496 Copy the example scripts to the folder "examples"
497 (it will be created if it doesn't exist).
498 """
499 function copy_examples()
500 PATH = "examples"
501 if ! isdir(PATH)
502 mkdir(PATH)
503 end
504 src_path = joinpath(dirname(pathof(@__MODULE__)), "..", PATH)
505 cp(joinpath(src_path, "compare_kps3_kps4.jl"), joinpath(PATH, "compare_kps3_kps4.jl"), force=true)
506 cp(joinpath(src_path, "plot2d.jl"), joinpath(PATH, "plot2d.jl"), force=true)
507 cp(joinpath(src_path, "simulate.jl"), joinpath(PATH, "simulate.jl"), force=true)
508 cp(joinpath(src_path, "reel_out_1p.jl"), joinpath(PATH, "reel_out_1p.jl"), force=true)
509 cp(joinpath(src_path, "reel_out_4p.jl"), joinpath(PATH, "reel_out_4p.jl"), force=true)
510 chmod(joinpath(PATH, "compare_kps3_kps4.jl"), 0o664)
511 chmod(joinpath(PATH, "plot2d.jl"), 0o664)
512 chmod(joinpath(PATH, "simulate.jl"), 0o664)
513 chmod(joinpath(PATH, "reel_out_1p.jl"), 0o664)
514 chmod(joinpath(PATH, "reel_out_4p.jl"), 0o664)
515 end
516
517 """
518 copy_bin()
519
520 Copy the scripts create_sys_image and run_julia to the folder "bin"
521 (it will be created if it doesn't exist).
522 """
523 function copy_bin()
524 PATH = "bin"
525 if ! isdir(PATH)
526 mkdir(PATH)
527 end
528 src_path = joinpath(dirname(pathof(@__MODULE__)), "..", PATH)
529 cp(joinpath(src_path, "create_sys_image2"), joinpath(PATH, "create_sys_image"), force=true)
530 cp(joinpath(src_path, "run_julia2"), joinpath(PATH, "run_julia"), force=true)
531 chmod(joinpath(PATH, "create_sys_image"), 0o774)
532 chmod(joinpath(PATH, "run_julia"), 0o774)
533 PATH = "test"
534 if ! isdir(PATH)
535 mkdir(PATH)
536 end
537 src_path = joinpath(dirname(pathof(@__MODULE__)), "..", PATH)
538 cp(joinpath(src_path, "create_sys_image2.jl"), joinpath(PATH, "create_sys_image.jl"), force=true)
539 cp(joinpath(src_path, "test_for_precompile.jl"), joinpath(PATH, "test_for_precompile.jl"), force=true)
540 cp(joinpath(src_path, "update_packages.jl"), joinpath(PATH, "update_packages.jl"), force=true)
541 chmod(joinpath(PATH, "create_sys_image.jl"), 0o664)
542 chmod(joinpath(PATH, "test_for_precompile.jl"), 0o664)
543 chmod(joinpath(PATH, "update_packages.jl"), 0o664)
544 end
545
546 @setup_workload begin
547 # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
548 # precompile file and potentially make loading faster.
549 # list = [OtherType("hello"), OtherType("world!")]
550 set_data_path()
551 kps4_::KPS4 = KPS4(KCU(se()))
552 kps3_::KPS3 = KPS3(KCU(se()))
553 @compile_workload begin
554 # all calls in this block will be precompiled, regardless of whether
555 # they belong to your package or not (on Julia 1.8 and higher)
556 integrator = KiteModels.init_sim!(kps3_; stiffness_factor=0.035, prn=false)
557 integrator = KiteModels.init_sim!(kps4_; stiffness_factor=0.035, prn=false)
558 nothing
559 end
560 end
561 end