diff --git a/examples/training.jl b/examples/training.jl
index fe5da594085eb21cb72e3d414cfd307aaf801e0f..b108b3f6e5769c8536bf0f4f4b53042b7606d25f 100644
--- a/examples/training.jl
+++ b/examples/training.jl
@@ -51,6 +51,9 @@ import ProgressMeter
using Printf
using Plots
+# %%
+product = Iterators.product
+
# %%
using Logging
@@ -392,6 +395,9 @@ plot!(ks[end:end], as[end:end], line = nothing, marker = :dot, label = "")
plot(plt_a, plt_b, plt_c, plt_d)
+# %% [markdown]
+# **TODO**: Can the R-dependence be obtained directly?
+
# %%
ks = range(a+1, 100, length=200)
plot()
@@ -508,12 +514,10 @@ plot(kps, xaxis = "Steps", yaxis = "k_p")
plot(kds, xaxis = "Steps", yaxis = "k_d")
# %%
-function train(kp, kd)
+function train(kp, kd; n_iters = 250)
λs = []
kps = []
kds = []
- n_iters = 250
- # pm = Progress(n_iters)
@time for iter in 1:n_iters
(λ, ∂kp, ∂kd) = conv_grad(ζ, ω, kp, kd)
step = 1e-1 * 1/sqrt(∂kp^2 + ∂kd^2) * 1/sqrt(iter)
@@ -522,15 +526,14 @@ function train(kp, kd)
push!(λs, λ)
push!(kps, kp)
push!(kds, kd)
- # next!(pm, showvalues=[(:λ, λ), (:kp, kp), (:kd, kd)])
end
(λs, kps, kds)
end
# %%
ress = Dict([])
-@time for kp in [1, 0, -1], kd in [1, 0, -1]
- (λs, kps, kds) = train(kp, kd)
+for (kp, kd) in product([1, 0, -1], [1, 0, -1])
+ @time (λs, kps, kds) = train(kp, kd, n_iters = 50)
ress[(kp,kd)] = (kps, kds)
plt_a = plot(λs, label = "", yaxis = "conv_rate", xaxis = "steps")
plt_b = plot(kps, label = "", yaxis = "kp", xaxis = "steps")
@@ -570,8 +573,7 @@ ys = []
us = []
vs = []
-@showprogress for kp in kps
- for kd in kds
+@showprogress for (kp, kd) in product(kps, kds)
(λ, ∂kp, ∂kd) = conv_grad(ζ, ω, kp, kd)
push!(xs, kp)
push!(ys, kd)
@@ -877,7 +879,7 @@ function real_conv_rate(ζ, ω, kp, kd)
(2ζ*ω+kd)/2
end
end
-
+
# %%
@time contourf(kps, kds, (kp, kd) -> real_conv_rate(.1, 1, kp, kd))
@@ -1212,4 +1214,239 @@ plot!(kps[end:end], kds[end:end], line=nothing, marker=:dot, label = "")
# %%
plot(τs)
+# %% [markdown]
+# # Gain Bound of PD-controlled Second Order System, with shift
+
+# %%
+function _gain_grad(ζ_val, ω_val, kp_val, kd_val; margin = 1e-3, R = 0.0)
+
+ sys = DynamicSystem(opt = Hypatia.Optimizer)
+ kp = @IQC.parameter sys kp kp_val
+ kd = @IQC.parameter sys kd kd_val
+ ζ = @IQC.parameter sys ζ ζ_val
+ ω = @IQC.parameter sys ω ω_val
+
+ @signals sys begin
+ x, 1
+ dx, 1
+ u, 1
+ w, 1
+ end
+
+ @equalities sys begin
+ 0 = u+kp*x+kd*dx
+ end
+
+ @tfs sys begin
+ x = (1/(s^2+2*ζ*ω*s+ω^2))*(u+w)
+ dx = s*x
+ end
+
+ τ = gen_τ(sys, IQC.Scalar(), minimize=true, nonneg=true)
+ @target sys begin
+ x^2 + R*u^2 <= τ*w^2
+ end
+
+
+ λ = IQC.convergence_rate(sys, margin=margin, ϵ = 1e-2)
+ λ = clamp(λ, -Inf, 0.0)
+ (success, P, feas_prob) = solve!(sys, λ=λ, margin=margin, req_PSD_P=true, silent=true)
+ IQC.set_param!(sys, τ, IQC.value(sys, τ))
+ @assert success
+
+ grads = IQC.gradient(sys, τ, λ=λ, margin=margin)
+ (IQC.value(sys, τ), grads[kp], grads[kd])
+end
+
+
+function gain_grad(ζ_val, ω_val, kp_val, kd_val; margin = 1e-3, R=0.0, n_tries = 10)
+ n_failed = 0
+ done = false
+ ret_val = nothing
+ while n_failed < n_tries && !done
+ try
+ ret_val = _gain_grad(ζ_val, ω_val, kp_val, kd_val, margin=margin, R=0.0)
+ done = true
+ catch e
+ n_failed += 1
+ end
+ end
+ ret_val
+end
+
+# %%
+ζ = .1
+ω = 1
+kp = 0.5
+kd = 0.5
+
+for _ in 1:20
+ try
+ @show _gain_grad(ζ, ω, kp, kd, margin=1e-3)
+ catch e
+ @show e
+ end
+end
+
+# %%
+ζ = .1
+ω = 1
+kp = 0.5
+kd = 0.5
+
+
+for _ in 1:20
+ try
+ @show _gain_grad(ζ, ω, kp, kd, margin=1e-6)
+ catch e
+ @show e
+ end
+end
+
+# %%
+ζ = .1
+ω = 1
+kp = 0.0
+kd = 0.0
+
+
+for _ in 1:20
+ try
+ @show _gain_grad(ζ, ω, kp, kd, margin=1e-6)
+ catch e
+ @show e
+ end
+end
+
+# %% [markdown]
+# Now, the gain bound varies a lot between runs. This is definitely a problem.
+
+# %%
+ζ = -.1
+ω = 1
+kp = -.5
+kd = -.1
+∂kp = 0
+∂kd = 0
+τ = 0
+
+τs = []
+kps = []
+kds = []
+n_iters = 10^2
+pm = Progress(n_iters)
+for iter in 1:n_iters
+ (τ, ∂kp, ∂kd) = gain_grad(ζ, ω, kp, kd, margin=1e-6)
+ step = 1e-1/sqrt(iter)
+ kp -= step*∂kp
+ kd -= step*∂kd
+ push!(τs, τ)
+ push!(kps, kp)
+ push!(kds, kd)
+ next!(pm, showvalues=[(:τ, τ), (:kp, kp), (:kd, kd)])
+end
+
+# %%
+plot(τs, yaxis = ("Gain Bound", :log10), xaxis = "Gradient Steps")
+
+# %%
+plot(kps, kds, xaxis = "kp", yaxis = "kd", label = "")
+plot!(kps[end:end], kds[end:end], line=nothing, marker=:dot, label = "")
+
+# %%
+ζ = .1
+ω = 1
+kp = .5
+kd = .5
+∂kp = 0
+∂kd = 0
+τ = 0
+
+τs = []
+kps = []
+kds = []
+n_iters = 10^2
+pm = Progress(n_iters)
+for iter in 1:n_iters
+ (τ, ∂kp, ∂kd) = gain_grad(ζ, ω, kp, kd, margin=1e-6, R = 1.0)
+ step = 1e-1/sqrt(iter)
+ kp -= step*sign(∂kp)
+ kd -= step*sign(∂kd)
+ push!(τs, τ)
+ push!(kps, kp)
+ push!(kds, kd)
+ next!(pm, showvalues=[(:τ, τ), (:kp, kp), (:kd, kd)])
+end
+
+# %%
+plot(τs, yaxis = ("Gain Bound", :log10), xaxis = "Gradient Steps")
+
+# %%
+plot(plot(kps, xaxis = "Steps", yaxis = "kp", label=""),plot(kds, xaxis = "Steps", yaxis = "kp", label = ""))
+
+# %%
+plot(kps, kds, xaxis = "kp", yaxis = "kd", label = "")
+plot!(kps[end:end], kds[end:end], line=nothing, marker=:dot, label = "")
+
+# %% [markdown]
+# For sufficiently large $R$, I expect the optimal $kp$ to be finite, but I seem to get large $kp$. Matlab should be able to tell me what the answer is, and then I should try to interpret the result.
+
+# %%
+
+# %%
+ζ = .1
+ω = 1
+kp = 1
+kd = 1
+∂kp = 0
+∂kd = 0
+τ = 0
+
+τs = []
+kps = []
+kds = []
+n_iters = 10^3
+pm = Progress(n_iters)
+for iter in 1:n_iters
+ (τ, ∂kp, ∂kd) = conv_grad(ζ, ω, kp, kd, margin=1e-6)
+ step = 1e-1 * 1/sqrt(∂kp^2 + ∂kd^2) * 1/sqrt(iter)
+ kp += step*∂kp
+ kd += step*∂kd
+ push!(τs, τ)
+ push!(kps, kp)
+ push!(kds, kd)
+ next!(pm, showvalues=[(:τ, τ), (:kp, kp), (:kd, kd)])
+end
+
+# %%
+plot(kps, kds, xaxis = "kp", yaxis = "kd", label = "")
+plot!(kps[end:end], kds[end:end], line=nothing, marker=:dot, label = "")
+
+# %%
+plot(τs)
+
+# %%
+τs = []
+
+for iter in 1:n_iters
+ (τ, ∂kp, ∂kd) = gain_grad(ζ, ω, kp, kd, R = 1.0)
+ step = 1e-1/sqrt(iter)
+ kp -= step*∂kp
+ kd -= step*∂kd
+ push!(τs, τ)
+ push!(kps, kp)
+ push!(kds, kd)
+ next!(pm, showvalues=[(:τ, τ), (:kp, kp), (:kd, kd)])
+end
+
+# %%
+plot(kps, kds, xaxis = "kp", yaxis = "kd", label = "")
+plot!(kps[end:end], kds[end:end], line=nothing, marker=:dot, label = "")
+
+# %%
+plot(τs)
+
+# %% [markdown]
+# # Shift, and then do gain bound gradients
+
# %%