Merge pull request #2 from Mohit-coder-droid/main

Adding Julia and PyTorch Scripts

Author: jbramburger

Autoencoder Neural Networks/PyTorch Code/Gissinger_conj.ipynb

@@ -0,0 +1,138 @@
+# -------------------------------------------------------------------------
+# Identifying Lyapunov Functions from Data with EDMD
+#
+# The goal of this script is to use EDMD and semidefinite programming to
+# identify a Lyapunov function for the planar discrete-time system
+#      x1 --> 0.3*x1 ,
+#      x2 --> -x1 + 0.5*x2 + (7/18)*x1^2,
+# using only data gathered from the system.
+#
+# To run this notebook one requires paths to the freely available software
+# package MOSEK, which can be downloaded at:
+#           https://www.mosek.com/downloads/
+#
+# This script accompanies Section 4.3 of Data-Driven Methods for
+# Dynamic Systems.
+#
+# This script was adapted from the MATLAB version by Jason J. Bramburger
+# 
+# Author: Mohit Sahu
+# -------------------------------------------------------------------------
+
+
+using DynamicPolynomials, SumOfSquares, MosekTools, LinearAlgebra
+include("utils.jl")
+
+# Method parameters
+# Dpₘₐₓ = max degree of Dp dictionary of obserables
+# Dqₘₐₓ = max degree of Dq dictionary of obserables
+# ϵ = hyperparameter specific to Lyapunov function for sharp bounds
+# cleanVal = remove coefficients smaller than this value from Lyap function
+
+Dpₘₐₓ = 4;
+Dqₘₐₓ = Dpₘₐₓ + 2;
+ϵ = 1;
+cleanVal = 1e-4;
+
+# Generate Synthetic Data
+# Number of data points
+N = 100;
+
+# Random points in [-2,2]x[-2,2]
+xdat = 2 .*rand(N,2) .- 1;
+ydat = zeros(N,2);
+
+# Images of random points under map dynamics
+ydat[:,1] = 0.3*xdat[:,1];
+ydat[:,2] = -xdat[:,1] + 0.5*xdat[:,2] + (7/18)*xdat[:,1].^2;
+
+# Create Q matrix
+pow = monpowers(2, Dqₘₐₓ)
+ell = size(pow)[1]  # number of nontrivial monomials in Q
+Q = zeros(ell, N)
+for i = 1:ell
+    zx = xdat.^permutedims(pow[i]);
+    Q[i,:] = prod(zx,dims=2);
+end
+
+# Create P matrix
+pow = monpowers(2, Dpₘₐₓ)
+ell = size(pow)[1]  # number of nontrivial monomials in P
+P = zeros(ell, N)
+for i = 1:ell
+    zy = ydat.^permutedims(pow[i]);
+    P[i,:] = prod(zy,dims=2);
+end
+
+# Create Koopman and Lie derivative matrix approximations
+# Koopman approximations
+K = P * pinv(Q)
+
+# Symbolic variables
+@polyvar x[1:2] # 2d state variables
+
+z = monomials(x, 0:Dpₘₐₓ)  # monomials that make up span(Dp)
+w = monomials(x, 0:Dqₘₐₓ)  # monomials that make up span(Dq)
+
+model = SOSModel(Mosek.Optimizer)
+
+@variable(model, c[1:length(z)])  # coeffs to build Lyapunov function in span(Dp)
+
+# As julia doesn't handle L1 norm, so t here is a dummy variable to solve this problem
+# See Tutorials/Linear programs/Tips and Tricks of Jump.jl documentation to know more about this
+@variable(model, t) 
+
+# Lie approximation
+# ---> Lyapunov function = c.'*z in span(Dp)
+L = c'*(K*w) - c'*z;
+
+# Identify Lyapunov function with SOS programming
+# Inequality constraints posed relaxed to SOS conditions
+@constraints(model, begin
+    c'*z - ϵ*dot(x,x) in SOSCone()
+    -L - ϵ*dot(x, x) in SOSCone()
+    [t; c] in MOI.NormOneCone(1+length(c))  # Minimizing L1 Norm of c
+end)
+
+# Objective function
+@objective(model, Min, t)
+
+# Solve SOS problem
+optimize!(model)
+
+# Remove coefficients smaller than cleanVal
+c = [abs(c_)>cleanVal ? c_ : zero(c_) for c_ in value.(c)]
+
+# Check identified Lyapunov function is a real Lyapunov function 
+# This time we maximize ϵ and check if it is positive
+model = SOSModel(Mosek.Optimizer)
+@variable(model, ϵ)  # ϵ is now a variable to be maximized
+
+# Discovered Lyapunov function 
+v = c'*z
+
+# Reverse the order of x₁ and x₂ to match the theory, as julia create polynomials in lexical order, so it's necessary to do to get to right answer
+v = subs(v, x[1]=>x[2], x[2]=>x[1])
+
+# True RHS of map
+f = [0.3*x[1]; -x[1] + 0.5*x[2] + (7/18)*x[1]^2]
+
+# True Lie Derivative
+Lie_exact = subs(v, x=>f)
+
+# Find Lyapunov function
+@constraints(model, begin
+    (v - ϵ * dot(x,x)) in SOSCone()
+    (-(Lie_exact - v) - ϵ*dot(x,x)) in SOSCone()
+end)
+
+@objective(model, Max, ϵ)
+optimize!(model)
+solution_summary(model)
+println("Maximized value of ϵ is $(value(ϵ))")
+
+if value(ϵ) > 0
+    println("We have discovered a true Lyapunov function from the data! \n") 
+else
+    println("Something went wrong - we cannot verify that this is a Lyapunov function. \n")
+end

Autoencoder Neural Networks/PyTorch Code/Global_Linearization.ipynb

@@ -0,0 +1,86 @@
+# -------------------------------------------------------------------------
+# Identifying Lyapunov Functions with Semidefinite Programming
+#
+# The goal of this script is to use semidefinite programming to identify a
+# Lyapunov function for the planar Moore-Greitzer system
+#      x1' = -x2 - 1.5x1^2 - 0.5x1^3,
+#      x2' = 3x1 - x2.
+# The script expands the unknown Lyapunov function as a polynomial up to a
+# given degree and tunes the coefficients in order to satisfy both of the
+# constraints:
+#      V - x1^2 - x2^2 is SOS
+#      -LV - x1^2 - x1^2 is SOS
+# where LV is the Lie derivative associated to the Moore-Greitzer system.
+#
+# To run this notebook one requires paths to the freely available software
+# package MOSEK, which can be downloaded at:
+#           https://www.mosek.com/downloads/
+#
+# This script accompanies Section 4.2 of Data-Driven Methods for
+# Dynamic Systems.
+#
+# This script was adapted from the MATLAB version by Jason J. Bramburger
+# 
+# Author: Mohit Sahu
+# -------------------------------------------------------------------------
+
+
+using DynamicPolynomials, SumOfSquares
+using MosekTools
+using CSDP
+
+@polyvar x[1:2]
+
+f = [-x[2] - 1.5*x[1]^2 - 0.5*x[1]^3;
+      3*x[1] - x[2]]
+
+# solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
+model = Model(Mosek.Optimizer)
+V_basis = monomials(x, 0:6)  # monomials for polynomial V
+@variable(model, c[1:length(V_basis)])
+
+# As julia doesn't handle L1 norm, so t here is a dummy variable to solve this problem
+# See Tutorials/Linear programs/Tips and Tricks of Jump.jl documentation to know more about this
+@variable(model, t)
+
+V = c'*V_basis
+
+using LinearAlgebra
+Dv = differentiate(V, x)
+Lie = dot(Dv, f)
+
+# Add SOS constraints
+@constraints(model, begin
+    V .- sum(x.^2) in SOSCone()
+    -Lie - sum(x.^2) in SOSCone()
+    [t; c] in MOI.NormOneCone(1+length(c))   # L1 Norm
+end)
+
+@objective(model, Min, t)
+optimize!(model)
+
+# Check if a solution was found
+if termination_status(model) == MOI.OPTIMAL
+    println("Lyapunov function identified!")
+    V_sol = value(V)
+    threshold = 1e-1
+    filtered_coeffs = [abs(c) > threshold ? c : zero(c) for c in coefficients(V_sol)]
+    println("Filtered V(x1, x2) = ", sum(filtered_coeffs .* monomials(V_sol)))
+else
+    println("Sorry, unable to identify a Lyapunov function")
+end
+
+# Plot the Lyapunov function
+using Plots
+plotlyjs()
+N = 1000;
+x1_range = LinRange(-5, 5, N)
+
+x = x1_range' .* ones(length(x1_range))
+y = ones(length(x1_range))' .* x1_range
+V_discovered = value.(c)'*V_basis;
+
+plot(x, y, V_discovered.(x,y), st=:surface)
+xlabel!("x1")
+ylabel!("x2")
+zlabel!("V(x1, x2)")

Autoencoder Neural Networks/PyTorch Code/NormalForm_pd.ipynb

@@ -0,0 +1,93 @@
+# -------------------------------------------------------------------------
+# Plotting upper and lower bounds on heteroclinic existence
+#
+# This script plots the upper and lower bounds over a range of the
+# hyperparameter lambda for the existence of heteroclinic orbit. Data is
+# loaded in from the .mat files heteroclinic_lambda_ + 'upper' or 'lower' for
+# upper or lower bounds + 'm2' or 'm3' for m = 2 or 3, respectively. Each
+# file containts the lambda values, the value of m for reference, and
+# arrays c2, c3, c4, c5 representing the values of upper or lower bounds
+# on mu for degree 2,3,4, or 5 auxiliary functions, respectively.
+#
+# This script accompanies Section 4.2 of Data-Driven Methods for
+# Dynamic Systems.
+#
+# This script was adapted from the MATLAB version by Jason J. Bramburger
+# 
+# Author: Mohit Sahu
+# -------------------------------------------------------------------------
+
+
+using MAT, Plots, Parameters
+
+# Upper Bound Data
+# m = 2 data
+file_path = "DataDrivenDynSyst/data_driven_poly_optimiazation/"
+vars = matread("$(file_path)heteroclinic_lambda_upper_m2.mat")
+@unpack c2, c3, c4, c5, lambda, m = vars
+vars2 = matread("$(file_path)heteroclinic_lambda_upper_m2_part2.mat")
+@unpack c3_2, c4_2, c5_2, lambda_2, m = vars2
+
+# m = 3 data
+# vars = matread("$(file_path)heteroclinic_lambda_upper_m3.mat")
+# @unpack c2, c3, c4, c5, lambda, m = vars
+# vars2 = matread("$(file_path)heteroclinic_lambda_upper_m3_part2.mat")
+# @unpack c2_2,c3_2, c4_2, c5_2, lambda_2, m = vars2
+
+# Unpack data and aggregate
+λ = [lambda lambda_2]
+c3 = [c3 c3_2]
+c4 = [c4 c4_2]
+c5 = [c5[:,1:181] c5_2]
+
+# Move failures off the figure
+if m == 2
+    c2 = [c2 2*ones(1,length(lambda_2))];
+    c2[1:46] .= 2;
+    c2[174:end] .= 2;
+    c3[260:end] .= 2;
+else
+    c2 = [c2 c2_2];
+end
+
+# Plot Upper bounds
+figure1 = plot()
+
+plot!(λ[1:3:end], c2[1:3:end], c=RGB(1, 69/255, 79/255), lw=4, label="c2")
+plot!(λ[1:3:end], c3[1:3:end], ls=:dash, c=RGB(36/255, 122/255, 254/255), lw=4, label="c3")
+plot!(λ[1:3:end], c4[1:3:end], ls=:dot, c=RGB(0, 168/255, 0), lw=4, label="c4")
+plot!(λ[1:3:end], c5[1:3:end], ls=:dashdot, c=RGB(255/255, 66/255, 161/255), lw=4, label="c5")
+
+display(figure1)
+
+# Load Lower bound data
+if m == 2
+    vars = matread("$(file_path)heteroclinic_lambda_lower_m2.mat")
+    @unpack c3, c4, c5, lambda, m = vars
+elseif m == 3
+    vars = matread("$(file_path)heteroclinic_lambda_lower_m3.mat")
+    @unpack c3, c4, c5, lambda, m = vars
+end
+
+# Plot Lower Bounds
+figure2 = plot()
+if m == 2
+    plot!(lambda[7:37], c3[7:37], ls=:dash, c=RGB(36/255, 122/255, 254/255), lw=4, label="c3")
+    plot!(lambda[7:37], c4[7:37], ls=:dashdot, c=RGB(0, 168/255, 0), lw=4, label="c4")
+    plot!(lambda[7:2:39], c5[7:2:39], ls=:dot, c=RGB(255/255, 66/255, 161/255), lw=4, label="c5")
+    xlims!(0, lambda[37])
+    ylims!(0.5, 1)
+    yticks!(0.5:0.1:1)
+elseif m == 3
+    plot!(lambda[7:end], c3[7:end], ls=:dash, c=RGB(36/255, 122/255, 254/255), lw=4, label="c3")
+    plot!(lambda[7:end], c4[7:end], ls=:dashdot, c=RGB(0, 168/255, 0), lw=4, label="c4")
+    plot!(lambda[7:2:end], c5[7:2:end], ls=:dot, c=RGB(255/255, 66/255, 161/255), lw=4, label="c5")
+    xlims!(0, lambda[37])
+    ylims!(0.30, 0.7)
+    yticks!(0.3:0.1:0.765)
+end
+
+xlabel!("λ", fontsize=24, fontweight=:bold)
+ylabel!("μ", fontsize=24, fontweight=:bold)
+
+display(figure2)

Autoencoder Neural Networks/PyTorch Code/NormalForm_sn.ipynb

@@ -0,0 +1,70 @@
+using DynamicPolynomials, SumOfSquares, MosekTools
+include("utils.jl")
+
+# ODE mu parameter
+m = 2
+
+# Bounding method parameters
+degV = 20  # Degree of auxiliary function
+λ = 10^3
+
+# Bisection Method for Finding Upper Bound
+function volume(μ, m, d, λ)
+    @polyvar u v
+
+    ϵ = 1e-4;
+
+    model = SOSModel(Mosek.Optimizer)
+    set_silent(model)
+
+    #Auxiliary function
+    @variable(model, V, createPoly([u,v], d))
+
+    # S procedure
+    d2 = d + m 
+    @variables(model, begin
+        s1, createPoly([u,v], d2)
+        s2 ,createPoly([u,v], d2)
+        s3 , createPoly(u, d2)
+    end)
+
+    # Derivatives
+    dVdu = differentiate(V, u)
+    dVdv = differentiate(V, v)
+
+    # Replacements
+    Vv0 = subs(V, v=>0)
+
+    @constraints(model, begin
+        subs(V,u=>0, v=>-ϵ) == 0
+        subs(V,u=>1,v=>0) == 0
+        λ*(dVdu*v - μ*dVdv*v - dVdv*(1-u)*u^m) -V -u*(1-u)*s1 + v*s2 in SOSCone()
+        -Vv0 -ϵ*(1-u) -u*(1-u)*s3 in SOSCone()
+        s1 in SOSCone()
+        s2 in SOSCone()
+        s3 in SOSCone()
+    end)
+
+    optimize!(model)
+
+    # Return whether we are able to solve the problem or not
+    return primal_status(model)
+end
+
+#Bisection method
+μ_left = 0
+μ_right = 2
+μ_mid = 0
+
+while abs((μ_right - μ_left)) >= 1e-5
+    μ_mid = 0.5 * (μ_left + μ_right)
+    flag = volume(μ_mid, m, degV, λ)
+
+    if flag==FEASIBLE_POINT
+        μ_left = μ_mid
+    else
+        μ_right = μ_mid
+    end
+end
+
+println("The lower bound on the minimum speed for m = $m is $μ_mid found using degree $degV polynomials.\n")