Use different formula for damage calculation

Project.toml

@@ -8,7 +8,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
-NOMAD = "02130f1c-4665-5b79-af82-ff1385104aa0"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -27,7 +26,6 @@ KiteUtils = "0.10.7"
 LinearAlgebra = "1.10, 1.11"
 NLsolve = "4.5.1"
 NOMAD = "2.4.1"
-PRIMA = "0.2.2"
 Parameters = "0.12.3"
 Pkg = "1.10, 1.11"
 StaticArrays = "1.9.13"
@@ -41,10 +39,10 @@ julia = "1.10, 1.11"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+NOMAD = "02130f1c-4665-5b79-af82-ff1385104aa0"
 ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ControlPlots", "Documenter", "Aqua", "PRIMA"]
+test = ["Test", "ControlPlots", "Documenter", "Aqua", "NOMAD"]

data/wc_settings.yaml

@@ -1,38 +1,38 @@
 wc_settings:
-    dt: 0.02           # time step of the winch controller
+    dt: 0.02            # time step of the winch controller
     test: false
-    fac: 1.0           # factor for I and P of lower force controller
-    max_iter: 100      # max iterations limit for the PID solvers
-    iter: 0            # actual max iterations of the PID solvers
-    t_startup: 0.25    # startup time for soft start  
-    t_blend: 0.1       # blending time of the mixers in seconds
-    v_sat_error: 1.0   # limit of reel-out speed error, used by the input sat block of the speed controller
-    v_sat: 8.0         # limit of reel-out speed, used by the output sat block of the speed controller
-    v_ri_max: 8.0      # maximal reel-in speed [m/s]
-    p_speed: 0.125     # P value of the speed controller
-    i_speed:  6.0      # I value of the speed controller
-    kb_speed: 4.0      # back calculation constant for the anti-windup loop of the speed controller
-    kt_speed: 5.0      # tracking constant of the speed controller
-    vf_max: 2.75       # reel-out velocity where the set force should reach it's maximum
-    pf_low: 0.00014    # P constant of the lower force controller 0.013,  0.00014 also works
-    if_low: 0.01125    # I constant of the lower force controller 0.0225, 0.01125 also works
+    fac: 1.0            # factor for I and P of lower force controller
+    max_iter: 100       # max iterations limit for the PID solvers
+    iter: 0             # actual max iterations of the PID solvers
+    t_startup: 0.25     # startup time for soft start  
+    t_blend: 0.1        # blending time of the mixers in seconds
+    v_sat_error: 1.0    # limit of reel-out speed error, used by the input sat block of the speed controller
+    v_sat: 8.0          # limit of reel-out speed, used by the output sat block of the speed controller
+    v_ri_max: 8.0       # maximal reel-in speed [m/s]
+    p_speed: 0.125      # P value of the speed controller
+    i_speed:  6.0       # I value of the speed controller
+    kb_speed: 4.0       # back calculation constant for the anti-windup loop of the speed controller
+    kt_speed: 5.0       # tracking constant of the speed controller
+    vf_max: 2.75        # reel-out velocity where the set force should reach it's maximum
+    pf_low: 0.00014     # P constant of the lower force controller 0.013,  0.00014 also works
+    if_low: 0.01125     # I constant of the lower force controller 0.0225, 0.01125 also works
 
-    df_low: 0.0        # D constant of lower force controller 0.000017
-    nf_low:  0         # filter constant n of lower force controller
-    kbf_low: 8.0       # back calculation constant for the anti-windup loop of the lower force controller
-    ktf_low: 8.0       # tracking constant of the lower force controller
-    f_low: 350         # lower force limit [N]
-    f_reelin: 700      # set force for reel-in phase [N]
-    f_high: 3800       # upper force limit [N]
-    pf_high: 0.0002304 # P constant of upper force controller
-    if_high: 0.012     # I constant of upper force controller
-    df_high: 2.4e-5    # D constant of upper force controller
-    nf_high: 15.0      # filter constant n of upper force controller
-    kbf_high: 1.0      # back calculation constant for the anti-windup loop of the upper force controller
-    ktf_high: 10.0     # tracking constant of the upper force controller
-    winch_iter: 10     # iterations of the winch model
-    max_acc: 8.0       # maximal acceleration of the winch (derivative of the set value of the reel-out speed)
-    damage_factor: 0.4 # damage at max acceleration for jerk_factor=0
-    jerk_factor: 0.9   # factor for the jerk (derivative of the acceleration), 0..1 
-    kv: 0.06           # proportional factor of the square root law, see function calc_vro
+    df_low: 0.0         # D constant of lower force controller 0.000017
+    nf_low:  0          # filter constant n of lower force controller
+    kbf_low: 8.0        # back calculation constant for the anti-windup loop of the lower force controller
+    ktf_low: 8.0        # tracking constant of the lower force controller
+    f_low: 350          # lower force limit [N]
+    f_reelin: 700       # set force for reel-in phase [N]
+    f_high: 3800        # upper force limit [N]
+    pf_high: 0.0002304  # P constant of upper force controller
+    if_high: 0.012      # I constant of upper force controller
+    df_high: 2.4e-5     # D constant of upper force controller
+    nf_high: 15.0       # filter constant n of upper force controller
+    kbf_high: 1.0       # back calculation constant for the anti-windup loop of the upper force controller
+    ktf_high: 10.0      # tracking constant of the upper force controller
+    winch_iter: 10      # iterations of the winch model
+    max_acc: 8.0        # maximal acceleration of the winch (derivative of the set value of the reel-out speed)
+    damage_factor: 0.05 # damage at max acceleration for jerk_factor=0
+    jerk_factor: 0.90   # factor for the jerk (derivative of the acceleration), 0..1 
+    kv: 0.06            # proportional factor of the square root law, see function calc_vro
 

docs/make.jl

@@ -22,6 +22,7 @@ makedocs(;
         "Winchcontroller" => "winchcontroller.md",
         "Winchcontroller Settings" => "settings.md",
         "Performance Indicators" => "performance_indicators.md",
+        "Autotuning" => "autotuning.md",
         "Tests" => "tests.md",
     ],
 )

docs/src/autotuning.md

@@ -0,0 +1,53 @@
+```@meta
+CurrentModule = WinchControllers
+```
+
+# Autotuning
+
+## Introduction
+The winch controller has 34 parameters (see [Winchcontroller Settings](@ref)). To determine and optimize all of them manually is time consuming and error prone, in particular if the physical system is changing. Therefore some form of automated tuning is desirable.
+
+## Methodology
+When the [Performance Indicators](@ref) are defined, an optimizer can be used to determine the controller parameters. The controller is using three PID controllers and a mixer. The following parameters are tuned automatically:
+```julia
+    wcs.i_speed  # the speed controller gain
+    wcs.p_speed  # the speed controller proportional gain
+    wcs.t_blend  # the blending time for switching between controllers
+    wcs.pf_low   # the lower force controller gain
+    wcs.if_low   # the lower force controller proportional gain
+    wcs.pf_high  # the upper force controller gain
+    wcs.if_high  # the upper force controller integral gain
+    wcs.df_high  # the differential gain of the upper force controller
+```
+The global, blackbox optimizer package [NOMAD.jl](https://github.com/bbopt/NOMAD.jl) is used for the optimization process, together with a test case that mimics extreme wind conditions.
+With a very simple system model it allows to find an optimal solution within less than one minute.
+
+## Example
+To run the auto-tuning script, launch Julia and execute:
+```julia
+include("scripts/autotune.jl")
+```
+It will use the initial conditions provided in the file "data/wc_settings.yaml" and save the result in "data/wc_settings_tuned.yaml".
+
+All values are limited to $\hat x \le 2.5x~\land~\hat x \le x/2.5$, where $x$ is the original value and $\hat x$ the optimized value. Therefore make sure that none of the original values is zero.
+
+Example output:
+```julia
+Blackbox evaluations:         897
+Total model evaluations:      116337
+Cache hits:                   107
+Total number of evaluations:  1004
+Autotuning results: [14.999895, 0.125, 0.1090168, 2.24464e-5, 0.04164095, 3.83195e-5, 0.0271663, 4.0008673e-6]
+Running simulation with tuned parameters...
+
+Performance of force controllers: 92.86 %
+Performance of speed controller:  93.82 %
+Damage with jerk:                 1.19 %
+Combined performance γ: 92.06 %
+
+[ Info: Tuned settings saved to data/wc_settings_tuned.yaml
+```
+You can plot the result with the command:
+```julia
+plot(lg)
+```

docs/src/performance_indicators.md

@@ -22,19 +22,21 @@ $\bar{v} = \frac{1}{n} \sum_{i=1}^{n} \lvert v_\mathrm{set} \rvert$
 The reason to use the mean square error for the velocity is that in this mode of operation the error only reduces the power output. Short peak errors are usually not a problem.
 
 ## Damage
-While the performance shall be maximized, the damage shall be minimized. Here, a very simple
+While the performance shall be maximized, the damage shall be minimized. Here, a simple
 damage model is presented,
 
-$\epsilon = \frac{\Delta}{a_\mathrm{max}}~{(\max |a|)^2}$,
+$\epsilon = \Delta~\left((1-\zeta)\left(\frac{\max |a|}{a_\mathrm{max}}\right)^2~+~\zeta \left(\frac{\mathrm{rms}(j)}{a_\mathrm{max}^2}\right)^4\right)$,
 
-where $a$ is the actual acceleration, $a_{max}$ the specified, maximal acceleration and $\Delta$ the damage that occurs at $a=a_\mathrm{max}$. The default value used is $\Delta=0.2$, but in the end this value needs to be determined based on the specification of the Winch. Control-codesign can be used to determine this parameter.
+where $a$ is the actual acceleration, $a_{max}$ the specified, maximal acceleration, $\mathrm{rms}(j)$ the root-mean-square of the jerk, the derivative of the acceleration, $\Delta$ the damage that occurs at $a=a_\mathrm{max}$ and $\zeta$ the jerk factor. The default values used are $\Delta=0.05$ and $\zeta = 0.9$, but in the end this value needs to be determined based on the specification of the Winch. Control-codesign can be used to determine this parameter.
+
+Using the forth potence of the jerk has proven to be effective in suppressing any oscillations and reducing overshoot of the controller while keeping it fast and stable.
 
 ## Combined performance
 The combined performance indicator $\gamma$ in the range of 0..1 is defined as
 
-$\gamma = 1 - \frac{1}{2}(F_\mathrm{err} + v_\mathrm{err})~-~\epsilon~$,
+$\gamma = 1 - \frac{1}{2.5}(1.5F_\mathrm{err} + v_\mathrm{err})~-~\epsilon~$,
 
-the average of the error of the speed and the force controllers minus the damage caused by the acceleration of the winch.
+the average of the error of the speed and the force controllers minus the damage caused by the acceleration of the winch. The force error has a slightly higher weight, because it is more critical for the controller performance.
 
 ## TODO
 - Quantify the robustness. This could be done by linearizing the system and checking gain and phase margin (for example), or by varying the model parameters (e.g. inertia of the drum) and checking if the combined performance stays above a required minimum. Possible robustness requirements:  

scripts/autotune.jl

@@ -1,11 +1,12 @@
 # this script tunes the controller parameters (well, eight of them)
 using Pkg
-if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies))
+if ! ("NOMAD" ∈ keys(Pkg.project().dependencies))
     using TestEnv; TestEnv.activate()
     using Test
 end
-using WinchControllers, KiteUtils, PRIMA, NOMAD, ControlPlots
+using WinchControllers, KiteUtils, NOMAD, ControlPlots
 
+LF = 2.5 # limit factor
 TUNED::Bool = false
 load_settings("system.yaml")
 
@@ -116,8 +117,8 @@ function autotune(max_iter=1000)
                       1,         # number of outputs of the blackbox
                       ["OBJ"],   # type of outputs of the blackbox
                       eval_fct;
-                      lower_bound = 0.5 .* x0,
-                      upper_bound = 2.0 .* x0)
+                      lower_bound = 1/LF .* x0,
+                      upper_bound = LF   .* x0)
     pb.options.max_bb_eval = max_iter 
     result = solve(pb, x0)
     x = result.x_best_feas
@@ -140,7 +141,7 @@ function copy_settings()
     load_settings("system_tuned.yaml")
 end
 function change_value(lines, varname, value::Union{Integer, Float64})
-    KiteUtils.change_value(lines, varname, repr(round(value, digits=5)))
+    KiteUtils.change_value(lines, varname, repr(round(value, digits=6)))
 end
 function update_settings(wcs::WCSettings)
     lines = KiteUtils.readfile("data/wc_settings_tuned.yaml") 

src/logging.jl

@@ -166,8 +166,8 @@
 
 function damage(logger::WCLogger)
     jerk_factor= logger.jerk_factor
-    logger.damage_factor*((jerk_factor*(rms(logger.jerk) / logger.max_acc^2)^2 + 
-                             (1-jerk_factor)*(maximum(norm.(logger.acc)))) / logger.max_acc)^2
+    logger.damage_factor*(jerk_factor*(rms(logger.jerk) / logger.max_acc^2)^4 + 
+                         (1-jerk_factor)*((maximum(norm.(logger.acc)))) / logger.max_acc)^2
 end
 
 """
@@ -183,5 +183,5 @@
 - The gamma value associated with the log of the used test case.
 """
 function gamma(logger::WCLogger)
-    1 - 0.5 * (f_err(logger) + v_err(logger)) - damage(logger)
+    1 - 1/2.5 * (1.5*f_err(logger) + v_err(logger)) - damage(logger)
 end

test/test_winchcontroller.jl

@@ -86,7 +86,7 @@ if @isdefined __TEST__
     @test V_WIND[Int64(2/0.02)+1] ≈ 9.0 rtol=1e-4
     @test V_WIND[Int64(4/0.02)+1] ≈ 0.0 atol=1e-4
     @test V_WIND[Int64(6/0.02)+1] ≈ 9.0 rtol=1e-4
-    @test mean(FORCE) ≈ 1783.2979123136008 rtol=1e-4
+    @test mean(FORCE) ≈ 1782.9912136687296 rtol=1e-4
     @test maximum(FORCE) < 4295
     @test minimum(FORCE[10:end]) > 10.0
     @test STATE[Int64(0.5/dt)] == 0