From a265544218e1fb686726d4add1c6700d333c414f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 8 Feb 2023 13:44:56 +0100 Subject: [PATCH 1/4] Add example of eigenfunction of Koopman operator --- .../koopman_eigenfunctions.jl | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl diff --git a/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl b/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl new file mode 100644 index 000000000..acc978f77 --- /dev/null +++ b/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl @@ -0,0 +1,89 @@ +# # Eigenfunctions of the Koopman operator + +#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Systems and Control/koopman_eigenfunctions.ipynb) +#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Systems and Control/koopman_eigenfunctions.ipynb) +# **Adapted from**: Example 2.6 of [MAI20] +# +# [MAI20] Mauroy, Alexandre, Aivar Sootla, and Igor Mezić. +# *Koopman framework for global stability analysis.* +# The Koopman Operator in Systems and Control: Concepts, Methodologies, and Applications (2020): 35-58. + +using Test #src +using DynamicPolynomials +@polyvar x[1:2] + +a = 1 +I = 0.05 +ε = 0.08 +γ = 1 +F0 = [-x[2] - x[1] * (x[1] - 1) * (x[1] - a) + I, ε * (x[1] - γ * x[2])] + +# We move equilibrium `(0.0256, 0.0256)` to the origin + +x1 = x2 = 0.0256 +F = [f(x => [x[1] + x1, x[2] + x2]) for f in F0] + +# We compute the Jacobian at the equilibrium + +J = [j(x => zeros(2)) for j in differentiate(F, x)] + +# We see that its eigenvalues indeed have negative real part: +using LinearAlgebra +E = eigen(J) + +# We set `w` as its dominant eigenvector: + +λ = E.values[end] +w = E.vectors[:, end] + +using SumOfSquares +r = 0.3 +X = @set x[1]^2 + x[2]^2 ≤ r^2 + +import CSDP +model = SOSModel(CSDP.Optimizer) +@variable(model, γ) +@objective(model, Min, γ) +@variable(model, ϕN, Poly(monomials(x, 2:10))) +ϕ = w ⋅ x + ϕN +∇ϕ = differentiate(ϕ, x) +@constraint(model, -γ ≤ F ⋅ ∇ϕ - λ * ϕ, domain = X) +@constraint(model, F ⋅ ∇ϕ - λ * ϕ ≤ γ, domain = X) +optimize!(model) +solution_summary(model) + +ϕ_opt = value(ϕ) + +using Plots +x1s = x2s = range(-0.3, stop = 0.3, length = 40) +ϕs = ϕ_opt.(x1s', x2s) +contour(x1s, x2s, ϕs, levels=[-0.2, -0.1, 0, 0.1, 0.2, 0.5, 0.7]) +θ = range(0, stop = 2π, length = 100) +plot!(r * cos.(θ), r * sin.(θ), label = "") +scatter!([0], [0], label = "") + +# We can compute the Laplace average as follows: + +using DifferentialEquations +function S(t, x1, x2, solver = DifferentialEquations.Tsit5()) + tspan = (0.0, t) + prob = DifferentialEquations.ODEProblem((v, p, t) -> [f(x => v) for f in F], [x1, x2], tspan) + traj = DifferentialEquations.solve(prob, solver, reltol=1e-4, abstol=1e-4) + return traj[end] +end + +using QuadGK +function laplace_average(f, λ, x1, x2, T = 10, args...) + v, _ = quadgk(0, T, rtol=1e-3) do t + s = S(t, x1, x2, args...) + return f(S(t, x1, x2, args...)) * exp(-λ * t) + end + return v / T +end + +lap(x1, x2) = laplace_average(v -> ϕ_opt(x => v), λ, x1, x2, 10) +laplace = lap.(x1s', x2s) + +# The error is given by: + +norm(laplace - ϕs) From 814b80005ca2a87aa28c3d4cf7460d517e0a51ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 16 Feb 2023 11:52:16 +0100 Subject: [PATCH 2/4] Add QuadGK --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index 68359e058..56add1b7a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,6 +18,7 @@ MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" From d09fd9634dc31d920662233e7fa1543e87c42405 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 20 Feb 2023 09:02:05 +0100 Subject: [PATCH 3/4] Split code --- .../tutorials/Systems and Control/koopman_eigenfunctions.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl b/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl index acc978f77..4abb1057d 100644 --- a/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl +++ b/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl @@ -52,8 +52,12 @@ model = SOSModel(CSDP.Optimizer) optimize!(model) solution_summary(model) +# The optimal value of `ϕ` is obtained as follows: + ϕ_opt = value(ϕ) +# Its plot is given below: + using Plots x1s = x2s = range(-0.3, stop = 0.3, length = 40) ϕs = ϕ_opt.(x1s', x2s) From efa3fd0c8f2143a6ae7ed62dbbe0c8d793356b70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 22 Feb 2023 14:30:04 +0100 Subject: [PATCH 4/4] Create function --- .../koopman_eigenfunctions.jl | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl b/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl index 4abb1057d..a2bc20db3 100644 --- a/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl +++ b/docs/src/tutorials/Systems and Control/koopman_eigenfunctions.jl @@ -40,16 +40,30 @@ using SumOfSquares r = 0.3 X = @set x[1]^2 + x[2]^2 ≤ r^2 +# We define the the program for the FitzHugh-Nagumo problem below: +# `N` is the degree of `ϕN` as defined in [MAI20, p. 50] and `M` is the +# degree of the multipliers for the constraints of `X`. +# As `maxdegree` corresponds to the degree of the multiplier multiplied +# by `r^2 - x[1]^2 - x[2]^2`, we set it to `M + 2`. + +function fitzhugh_nagumo(solver, N, M) + model = SOSModel(solver) + @variable(model, γ) + @objective(model, Min, γ) + @variable(model, ϕN, Poly(monomials(x, 2:N))) + ϕ = w ⋅ x + ϕN + ∇ϕ = differentiate(ϕ, x) + @constraint(model, -γ ≤ F ⋅ ∇ϕ - λ * ϕ, domain = X, maxdegree = M + 2) + @constraint(model, F ⋅ ∇ϕ - λ * ϕ ≤ γ, domain = X, maxdegree = M + 2) + optimize!(model) + return ϕ, model +end + +# In [MAI20, p. 50], we read that the result is obtained with `N = 10` and +# `M = 20`. + import CSDP -model = SOSModel(CSDP.Optimizer) -@variable(model, γ) -@objective(model, Min, γ) -@variable(model, ϕN, Poly(monomials(x, 2:10))) -ϕ = w ⋅ x + ϕN -∇ϕ = differentiate(ϕ, x) -@constraint(model, -γ ≤ F ⋅ ∇ϕ - λ * ϕ, domain = X) -@constraint(model, F ⋅ ∇ϕ - λ * ϕ ≤ γ, domain = X) -optimize!(model) +ϕ, model = fitzhugh_nagumo(CSDP.Optimizer, 10, 20) solution_summary(model) # The optimal value of `ϕ` is obtained as follows: