From 72924e1aef25b823e1cb76de6e60cafdaee513a3 Mon Sep 17 00:00:00 2001
From: Johan Gronqvist <johan.gronqvist@control.lth.se>
Date: Sat, 10 Jun 2023 12:03:50 +0000
Subject: [PATCH] Test more validation metrics for gradients

---
 TestConstantP.jl | 119 +++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 115 insertions(+), 4 deletions(-)

diff --git a/TestConstantP.jl b/TestConstantP.jl
index 6f0bbad..2f38256 100644
--- a/TestConstantP.jl
+++ b/TestConstantP.jl
@@ -75,6 +75,7 @@ function bisect(; Δ=1e-8, kwargs...)
             λu = λc
             λc *= 2
         end
+        λℓ = λc
     end
     
     while λu-λℓ > Δ
@@ -88,9 +89,18 @@ function bisect(; Δ=1e-8, kwargs...)
     λℓ
 end
 
+# %%
+ϵ = 1e-4
+
 # %%
 bisect(ζ=.1)
 
+# %%
+bisect(ζ=.1+ϵ)
+
+# %%
+bisect(ζ=.2-ϵ)
+
 # %%
 bisect(ζ=.2)
 
@@ -104,16 +114,16 @@ bisect(ζ=.2)
 valid, P = validate(λ=λ, ζ=.1, ret_P=true)
 
 # %%
-validate(λ=λ, ζ=.1, P=P)
+@show validate(λ=λ, ζ=.1, P=P) validate(λ=λ, ζ=.1+ϵ, P=P)
 
 # %%
-validate(λ=λ, ζ=.1, P=P)
+bisect(ζ=.1, P=P)
 
 # %%
-validate(λ=λ, ζ=.2, P=P)
+bisect(ζ=.1+ϵ, P=P)
 
 # %%
-bisect(ζ=.1, P=P)
+bisect(ζ=.2-ϵ, P=P)
 
 # %%
 bisect(ζ=.2, P=P)
@@ -129,9 +139,110 @@ P
 # %%
 valid, P = validate(λ=λ, ζ=.2, ret_P=true)
 
+# %%
+@show validate(λ=λ, ζ=.2-ϵ, P=P) validate(λ=λ, ζ=.2, P=P)
+
+# %%
+bisect(ζ=.1, P=P)
+
+# %%
+bisect(ζ=.1+ϵ, P=P)
+
+# %%
+bisect(ζ=.2-ϵ, P=P)
+
+# %%
+bisect(ζ=.2, P=P)
+
 # %%
 P
 
 # %%
 
 # %%
+function validate_all(;λ, ζs, ω=1.0, ϵ=1e-8, P = nothing, ret_P = false)
+    model = JuMP.Model(MosekTools.Optimizer)
+    set_silent(model)
+
+    
+    M = [ 1+ϵ 0; 0 ϵ ]
+    A(ζ) = [ λ 1; -ω^2 λ-2*ζ*ω ]
+
+    if isnothing(P)
+        lhss = []
+        @variable model P_var[1:2, 1:2] PSD
+        for ζ in ζs
+            lhs = A(ζ)'*P_var+P_var*A(ζ)+M
+            @constraint model lhs<=0 PSDCone()
+            push!(lhss, lhs)
+        end
+            
+        optimize!(model)
+
+        # @show termination_status(model) value.(lhs) eigvals(value.(lhs))
+        valid = maximum([ maximum(eigvals(value.(lhs))) for lhs in lhss ]) < 0
+        if ret_P
+            (valid, value.(P_var))
+        else
+            valid
+        end
+    else
+        lhss = []
+        for ζ in ζs
+            lhs = A(ζ)'*P+P*A(ζ)+M
+            push!(lhss, lhs)
+        end
+        valid = maximum([ maximum(eigvals(lhs)) for lhs in lhss]) < 0
+        valid
+            
+    end
+
+end
+
+# %%
+function bisect_all(; Δ=1e-8, kwargs...)
+    if validate_all(λ=0.0; kwargs...)
+        λℓ = 0.0
+        λc = 1.0
+        while validate_all(λ = λc; kwargs...)
+            λℓ = λc
+            λc *= 2
+        end
+        λu = λc
+    else
+        λu = 0.0
+        λc = -1.0
+        while !validate_all(λ=λc; kwargs...)
+            λu = λc
+            λc *= 2
+        end
+        λℓ = λc
+    end
+    
+    while λu-λℓ > Δ
+        λc = (λu+λℓ)/2
+        if validate_all(λ=λc; kwargs...)
+            λℓ = λc
+        else
+            λu = λc
+        end
+    end
+    λℓ
+end
+
+# %%
+λ = bisect(ζ = .1)
+
+# %%
+valid, P = validate(λ=λ, ζ=.1, ret_P=true)
+
+# %%
+validate_all(λ=λ, ζs = [.1, .1+ϵ])
+
+# %%
+validate_all(λ=λ, ζs = [.2-ϵ, .2])
+
+# %%
+bisect_all(ζs=[.1, .2])
+
+# %%
-- 
GitLab