From 516898362890e4313edeb3bb27005a5e0a0a4b74 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Tue, 26 Mar 2024 18:42:02 -0400 Subject: [PATCH 01/10] Added `name` keyword to and removed function output of `NeuralLyapunovPDESystem`; made dVdt calculation more efficient in `NumericalNeuralLyapunovFunctions` --- src/NeuralLyapunovPDESystem.jl | 116 ++++++++++++++++------------ src/policy_search.jl | 19 +++-- test/damped_pendulum.jl | 6 +- test/damped_sho.jl | 6 +- test/inverted_pendulum.jl | 7 +- test/inverted_pendulum_ODESystem.jl | 7 +- test/roa_estimation.jl | 5 +- 7 files changed, 94 insertions(+), 72 deletions(-) diff --git a/src/NeuralLyapunovPDESystem.jl b/src/NeuralLyapunovPDESystem.jl index 81a059a..2d35ed3 100644 --- a/src/NeuralLyapunovPDESystem.jl +++ b/src/NeuralLyapunovPDESystem.jl @@ -37,6 +37,7 @@ inputs. Function && !(dynamics isa ODEFunction)`; when `dynamics isa ODEFunction`, `policy_search` must be `false`, so should not be supplied; when `dynamics isa ODESystem`, value inferred by the presence of unbound inputs. +- `name`: the name of the constructed `PDESystem` """ function NeuralLyapunovPDESystem( dynamics::Function, @@ -47,8 +48,9 @@ function NeuralLyapunovPDESystem( p = SciMLBase.NullParameters(), state_syms = [], parameter_syms = [], - policy_search::Bool = false -)::Tuple{PDESystem, Function} + policy_search::Bool = false, + name +)::PDESystem ########################## Define state symbols ########################### state_dim = length(lb) @@ -93,7 +95,8 @@ function NeuralLyapunovPDESystem( state, params, defaults, - policy_search + policy_search, + name ) end @@ -103,8 +106,9 @@ function NeuralLyapunovPDESystem( ub, spec::NeuralLyapunovSpecification; fixed_point = zeros(length(lb)), - p = SciMLBase.NullParameters() -)::Tuple{PDESystem, Function} + p = SciMLBase.NullParameters(), + name +)::PDESystem if dynamics.mass_matrix !== I throw(ErrorException("DAEs are not supported at this time. Please supply dynamics" * " without a mass matrix.")) @@ -135,7 +139,8 @@ function NeuralLyapunovPDESystem( p = p, state_syms = s_syms, parameter_syms = p_syms, - policy_search = false + policy_search = false, + name ) end @@ -143,8 +148,9 @@ function NeuralLyapunovPDESystem( dynamics::ODESystem, bounds, spec::NeuralLyapunovSpecification; - fixed_point = zeros(length(bounds)) -)::Tuple{PDESystem, Function} + fixed_point = zeros(length(bounds)), + name +)::PDESystem ######################### Check for policy search ######################### f, x, p, policy_search = if isempty(ModelingToolkit.unbound_inputs(dynamics)) (ODEFunction(dynamics), states(dynamics), parameters(dynamics), false) @@ -180,7 +186,8 @@ function NeuralLyapunovPDESystem( state, p, ModelingToolkit.get_defaults(dynamics), - policy_search + policy_search, + name ) end @@ -192,8 +199,9 @@ function _NeuralLyapunovPDESystem( state, params, defaults, - policy_search::Bool -)::Tuple{PDESystem, Function} + policy_search::Bool, + name +)::PDESystem ########################## Unpack specifications ########################## structure = spec.structure minimzation_condition = spec.minimzation_condition @@ -266,67 +274,66 @@ function _NeuralLyapunovPDESystem( end ########################### Construct PDESystem ########################### - @named lyapunov_pde_system = PDESystem( + lyapunov_pde_system = PDESystem( eqs, bcs, domains, state, φ(state), params; - defaults = defaults + defaults = defaults, + name = name ) - ################### Return PDESystem and neural network ################### - # φ_func is the numerical form of neural network output - function φ_func(phi, θ, x) - reduce( - vcat, - Array(phi[i](x, θ.depvar[net_syms[i]])) for i in 1:output_dim - ) - end - - return lyapunov_pde_system, φ_func + return lyapunov_pde_system end """ - NumericalNeuralLyapunovFunctions(phi, θ, network_func, structure, dynamics, fixed_point; - jac, J_net) + NumericalNeuralLyapunovFunctions(phi, θ, structure, dynamics, fixed_point; jac, J_net) Returns the Lyapunov function, its time derivative, and its gradient: `V(state)`, `V̇(state)`, and `∇V(state)` These functions can operate on a state vector or columnwise on a matrix of state vectors. -`phi` is the neural network with parameters `θ`. `network_func(phi, θ, state)` is an output -of `NeuralLyapunovPDESystem`, which evaluates the neural network represented by `phi` with -parameters `θ` at `state`. +`phi` is the neural network with parameters `θ`. The Lyapunov function structure is specified in structure, which is a `NeuralLyapunovStructure`. The Jacobian of the network is either specified via -`J_net(_phi, _θ, state)` or calculated using `jac`, which defaults to -`ForwardDiff.jacobian`. +`J_net(phi, θ, state)` or calculated using `jac`, which defaults to `ForwardDiff.jacobian`. """ function NumericalNeuralLyapunovFunctions( phi, θ, - network_func::Function, structure::NeuralLyapunovStructure, dynamics::Function, fixed_point; p = SciMLBase.NullParameters(), jac = ForwardDiff.jacobian, - J_net = (_phi, _θ, x) -> jac((y) -> network_func(_phi, _θ, y), x) + J_net = Nothing )::Tuple{Function, Function, Function} - # Make Network function - _net_func = (x) -> network_func(phi, θ, x) - _J_net = (x) -> J_net(phi, θ, x) + # network_func is the numerical form of neural network output + output_dim = structure.network_dim + function network_func(x) + reduce( + vcat, + Array(phi[i](x, θ.depvar[Symbol(:φ, i)])) for i in 1:output_dim + ) + end + + # Make Jacobian of networkfunction + _J_net = if isnothing(J_net) + (x) -> jac(network_func, x) + else + (x) -> J_net(phi, θ, x) + end # Numerical form of Lyapunov function - V_func(state::AbstractVector) = structure.V(_net_func, state, fixed_point) + V_func(state::AbstractVector) = structure.V(network_func, state, fixed_point) V_func(state::AbstractMatrix) = mapslices(V_func, state, dims = [1]) # Numerical gradient of Lyapunov function ∇V_func(state::AbstractVector) = structure.∇V( - _net_func, + network_func, _J_net, state, fixed_point @@ -336,7 +343,7 @@ function NumericalNeuralLyapunovFunctions( # Numerical time derivative of Lyapunov function function V̇_func(state::AbstractVector) structure.V̇( - _net_func, + network_func, _J_net, dynamics, state, @@ -351,35 +358,45 @@ function NumericalNeuralLyapunovFunctions( end """ - NumericalNeuralLyapunovFunctions(phi, θ, network_func, V_structure, dynamics, - fixed_point, grad) + NumericalNeuralLyapunovFunctions(phi, θ, network_dim, V_structure, dynamics, + fixed_point; p, grad, deriv) Returns the Lyapunov function, its time derivative, and its gradient: `V(state)`, `V̇(state)`, and `∇V(state)`. These functions can operate on a state vector or columnwise on a matrix of state vectors. -`phi` is the neural network with parameters `θ`. `network_func` is an output of -`NeuralLyapunovPDESystem`. +`phi` is the neural network with parameters `θ`. `network_dim` is the output dimension of +the neural network. The Lyapunov function structure is defined by `V_structure(_network_func, state, fixed_point)` Its gradient is calculated using `grad`, which defaults to `ForwardDiff.gradient`. + +The Lyapunov function's time derivative is calculated using `deriv`, which defaults to +`ForwardDiff.derivative`, along with `dynamics`, which should be the closed-loop dynamics +`ẋ = dynamics(x, p, t)`; only `t = 0.0` is used. Parameters `p` can be supplied. """ function NumericalNeuralLyapunovFunctions( phi, θ, - network_func::Function, + network_dim, V_structure::Function, dynamics::Function, fixed_point; p = SciMLBase.NullParameters(), - grad = ForwardDiff.gradient + grad = ForwardDiff.gradient, + deriv = ForwardDiff.derivative )::Tuple{Function, Function, Function} - # Make network function - _net_func = (x) -> network_func(phi, θ, x) + # network_func is the numerical form of neural network output + function network_func(x) + reduce( + vcat, + Array(phi[i](x, θ.depvar[Symbol(:φ, i)])) for i in 1:network_dim + ) + end # Numerical form of Lyapunov function - V_func(state::AbstractVector) = V_structure(_net_func, state, fixed_point) + V_func(state::AbstractVector) = V_structure(network_func, state, fixed_point) V_func(state::AbstractMatrix) = mapslices(V_func, state, dims = [1]) # Numerical gradient of Lyapunov function @@ -387,7 +404,10 @@ function NumericalNeuralLyapunovFunctions( ∇V_func(state::AbstractMatrix) = mapslices(∇V_func, state, dims = [1]) # Numerical time derivative of Lyapunov function - V̇_func(state::AbstractVector) = dynamics(state, p, 0.0) ⋅ ∇V_func(state) + V̇_func(state::AbstractVector) = deriv( + (δt) -> V_func(state + δt * dynamics(state, p, 0.0)), + 0.0 + ); V̇_func(state::AbstractMatrix) = mapslices(V̇_func, state, dims = [1]) return V_func, V̇_func, ∇V_func diff --git a/src/policy_search.jl b/src/policy_search.jl index 80b055e..93a4e2d 100644 --- a/src/policy_search.jl +++ b/src/policy_search.jl @@ -56,20 +56,25 @@ Returns the control policy as a function of the state The returned function can operate on a state vector or columnwise on a matrix of state vectors. -`phi` is the neural network with parameters `θ`. `network_func` is an output of -`NeuralLyapunovPDESystem`. -The control policy is `control_structure` composed with the last `dim` outputs of the -neural network, as set up by `add_policy_search`. +`phi` is the neural network with parameters `θ`. The network should have `network_dim` +outputs, the last `control_dim` of which will be passed into `control_structure` to create +the policy output. """ function get_policy( phi, θ, - network_func::Function, - dim::Integer; + network_dim::Integer, + control_dim::Integer; control_structure::Function = identity ) function policy(state::AbstractVector) - control_structure(network_func(phi, θ, state)[(end - dim + 1):end]) + control_structure( + reduce( + vcat, + Array(phi[i](state, θ.depvar[Symbol(:φ, i)])) + for i in (network_dim - control_dim + 1):network_dim + ) + ) end policy(state::AbstractMatrix) = mapslices(policy, state, dims = [1]) diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl index fcda6e9..48cd23d 100644 --- a/test/damped_pendulum.jl +++ b/test/damped_pendulum.jl @@ -82,7 +82,7 @@ spec = NeuralLyapunovSpecification( ############################# Construct PDESystem ############################# -pde_system, network_func = NeuralLyapunovPDESystem( +@named pde_system = NeuralLyapunovPDESystem( ODEFunction(dynamics), lb, ub, @@ -103,10 +103,10 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func, ∇V_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = NumericalNeuralLyapunovFunctions( discretization.phi, res.u, - network_func, + dim_output, structure.V, ODEFunction(dynamics), zeros(length(bounds)); diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 5f06e36..57e4def 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -62,7 +62,7 @@ spec = NeuralLyapunovSpecification( ############################# Construct PDESystem ############################# -pde_system, network_func = NeuralLyapunovPDESystem( +@named pde_system = NeuralLyapunovPDESystem( dynamics, lb, ub, @@ -82,10 +82,10 @@ prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func, ∇V_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = NumericalNeuralLyapunovFunctions( discretization.phi, res.u, - network_func, + dim_output, structure.V, ODEFunction(dynamics), zeros(2); diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl index 54dfbf0..b28ac6d 100644 --- a/test/inverted_pendulum.jl +++ b/test/inverted_pendulum.jl @@ -77,7 +77,7 @@ spec = NeuralLyapunovSpecification( ############################# Construct PDESystem ############################# -pde_system, network_func = NeuralLyapunovPDESystem( +@named pde_system = NeuralLyapunovPDESystem( open_loop_pendulum_dynamics, lb, ub, @@ -102,17 +102,16 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func, ∇V_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = NumericalNeuralLyapunovFunctions( discretization.phi, res.u, - network_func, structure, open_loop_pendulum_dynamics, upright_equilibrium; p = p ) -u = get_policy(discretization.phi, res.u, network_func, dim_u) +u = get_policy(discretization.phi, res.u, dim_output, dim_u) ################################## Simulate ################################### diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index 5738c82..ce70a76 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -90,7 +90,7 @@ spec = NeuralLyapunovSpecification( ############################# Construct PDESystem ############################# -pde_system, network_func = NeuralLyapunovPDESystem( +@named pde_system = NeuralLyapunovPDESystem( driven_pendulum, bounds, spec; @@ -110,17 +110,16 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func, ∇V_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = NumericalNeuralLyapunovFunctions( discretization.phi, res.u, - network_func, structure, open_loop_pendulum_dynamics, upright_equilibrium; p = p ) -u = get_policy(discretization.phi, res.u, network_func, dim_u) +u = get_policy(discretization.phi, res.u, dim_output, dim_u) ################################## Simulate ################################### diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index 8c341d4..cd57d85 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -47,7 +47,7 @@ spec = NeuralLyapunovSpecification( ############################# Construct PDESystem ############################# -pde_system, network_func = NeuralLyapunovPDESystem( +@named pde_system = NeuralLyapunovPDESystem( f, lb, ub, @@ -66,10 +66,9 @@ prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func, ∇V_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = NumericalNeuralLyapunovFunctions( discretization.phi, res.u, - network_func, structure, f, zeros(length(lb)) From ca58f8f49da6a8bf12b32eee3c24df7fcea2218f Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 27 Mar 2024 10:56:49 -0400 Subject: [PATCH 02/10] Update `local_Lyapunov` to use `f(x, p, t)` instead of `f(x)` --- src/local_Lyapunov.jl | 36 +++++++++++++----------------------- 1 file changed, 13 insertions(+), 23 deletions(-) diff --git a/src/local_Lyapunov.jl b/src/local_Lyapunov.jl index 45418f0..7cbd9e4 100644 --- a/src/local_Lyapunov.jl +++ b/src/local_Lyapunov.jl @@ -1,20 +1,21 @@ """ get_local_Lyapunov(dynamics, state_dim; fixed_point, dynamics_jac) -Uses semidefinite programming to derive a quadratic Lyapunov function for the -linearization of dynamics around fixed_point. -Returns (V, dV/dt, ∇V). +Use semidefinite programming to derive a quadratic Lyapunov function for the +linearization of `dynamics` around `fixed_point`. +Return `(V, dV/dt, ∇V)`. -If dynamics_jac is nothing, the Jacobian of the dynamics is calculated using -ForwardDiff. Other allowable forms are a function which takes in the state and -outputs the jacobian of dynamics or an AbstractMatrix representing the Jacobian -at fixed_point. If fixed_point is not specified, it defaults to the origin. +If `isnothing(dynamics_jac)`, the Jacobian of the `dynamics(x, p, t)` with respect to `x` is +calculated using `ForwardDiff`. Other allowable forms are a function which takes in the +state and outputs the jacobian of `dynamics` or an `AbstractMatrix` representing the +Jacobian at `fixed_point`. If `fixed_point` is not specified, it defaults to the origin. """ function local_Lyapunov(dynamics::Function, state_dim, optimizer_factory; - fixed_point = zeros(state_dim), dynamics_jac = nothing) + fixed_point = zeros(state_dim), dynamics_jac = nothing, + p = SciMLBase.NullParameters()) # Linearize the dynamics A = if isnothing(dynamics_jac) - ForwardDiff.jacobian(dynamics, fixed_point) + ForwardDiff.jacobian(x -> dynamics(x, p, 0.0), fixed_point) elseif dynamics_jac isa AbstractMatrix dynamics_jac elseif dynamics_jac isa Function @@ -38,23 +39,12 @@ function local_Lyapunov(dynamics::Function, state_dim, optimizer_factory; V(states::AbstractMatrix) = mapslices(V, states, dims = [1]) # Numerical gradient of Lyapunov function - ∇V(state::AbstractVector) = 2 * Psol * (state - fixed_point) + ∇V(state::AbstractVector) = 2 * ( Psol * (state - fixed_point) ) ∇V(states::AbstractMatrix) = mapslices(∇V, states, dims = [1]) # Numerical time derivative of Lyapunov function - V̇(state::AbstractVector) = dynamics(state) ⋅ ∇V(state) - function V̇(states::AbstractMatrix) - reshape( - map( - x -> x[1] ⋅ x[2], - zip( - eachslice(dynamics(states), dims = 2), - eachslice(∇V(states), dims = 2) - ) - ), - (1, :) - ) - end + V̇(state::AbstractVector) = 2 * dot(dynamics(state, p, 0.0), Psol, state - fixed_point) + V̇(states::AbstractMatrix) = mapslices(V̇, states, dims = [1]) return V, V̇, ∇V end From ae0069b3195430eac1c29acc0c301a7084a332f4 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 27 Mar 2024 12:38:09 -0400 Subject: [PATCH 03/10] Change documentation from indicative to imperative mood and change "relu" to "rectifier" --- src/conditions_specification.jl | 22 +++++++++--------- src/decrease_conditions.jl | 34 ++++++++++++++-------------- src/decrease_conditions_RoA_aware.jl | 12 +++++----- src/minimization_conditions.jl | 32 +++++++++++++------------- src/policy_search.jl | 12 +++++----- src/structure_specification.jl | 6 ++--- 6 files changed, 59 insertions(+), 59 deletions(-) diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl index 010e7b8..9de4922 100644 --- a/src/conditions_specification.jl +++ b/src/conditions_specification.jl @@ -57,10 +57,10 @@ struct NeuralLyapunovSpecification end """ -check_nonnegativity(cond::AbstractLyapunovMinimizationCondition) + check_nonnegativity(cond::AbstractLyapunovMinimizationCondition) -`true` if `cond` specifies training to meet the Lyapunov minimization condition, `false` if -`cond` specifies no training to meet this condition. +Return `true` if `cond` specifies training to meet the Lyapunov minimization condition, and +`false` if `cond` specifies no training to meet this condition. """ function check_nonnegativity(cond::AbstractLyapunovMinimizationCondition)::Bool error("check_nonnegativity not implemented for AbstractLyapunovMinimizationCondition " * @@ -70,8 +70,8 @@ end """ check_minimal_fixed_point(cond::AbstractLyapunovMinimizationCondition) -`true` if `cond` specifies training for the Lyapunov function to equal zero at the -fixed point, `false` if `cond` specifies no training to meet this condition. +Return `true` if `cond` specifies training for the Lyapunov function to equal zero at the +fixed point, and `false` if `cond` specifies no training to meet this condition. """ function check_minimal_fixed_point(cond::AbstractLyapunovMinimizationCondition)::Bool error("check_minimal_fixed_point not implemented for " * @@ -81,8 +81,8 @@ end """ get_minimization_condition(cond::AbstractLyapunovMinimizationCondition) -Returns a function of `V`, `state`, and `fixed_point` that equals zero when the -Lyapunov minimization condition is met and greater than zero when it's violated. +Return a function of `V`, `state`, and `fixed_point` that equals zero when the Lyapunov +minimization condition is met and is greater than zero when it's violated. """ function get_minimization_condition(cond::AbstractLyapunovMinimizationCondition) error("get_condition not implemented for AbstractLyapunovMinimizationCondition of " * @@ -92,8 +92,8 @@ end """ check_decrease(cond::AbstractLyapunovDecreaseCondition) -`true` if `cond` specifies training to meet the Lyapunov decrease condition, `false` -if `cond` specifies no training to meet this condition. +Return `true` if `cond` specifies training to meet the Lyapunov decrease condition, and +`false` if `cond` specifies no training to meet this condition. """ function check_decrease(cond::AbstractLyapunovDecreaseCondition)::Bool error("check_decrease not implemented for AbstractLyapunovDecreaseCondition of type " * @@ -103,8 +103,8 @@ end """ get_decrease_condition(cond::AbstractLyapunovDecreaseCondition) -Returns a function of `V`, `dVdt`, `state`, and `fixed_point` that is equal to zero -when the Lyapunov decrease condition is met and greater than zero when it is violated. +Return a function of `V`, `dVdt`, `state`, and `fixed_point` that is equal to zero +when the Lyapunov decrease condition is met and is greater than zero when it is violated. """ function get_decrease_condition(cond::AbstractLyapunovDecreaseCondition) error("get_condition not implemented for AbstractLyapunovDecreaseCondition of type " * diff --git a/src/decrease_conditions.jl b/src/decrease_conditions.jl index f173092..684905e 100644 --- a/src/decrease_conditions.jl +++ b/src/decrease_conditions.jl @@ -1,11 +1,11 @@ """ - LyapunovDecreaseCondition(check_decrease, decrease, strength, relu) + LyapunovDecreaseCondition(check_decrease, decrease, strength, rectifier) Specifies the form of the Lyapunov conditions to be used; if `check_decrease`, training will enforce `decrease(V, dVdt) ≤ strength(state, fixed_point)`. The inequality will be approximated by the equation - `relu(decrease(V, dVdt) - strength(state, fixed_point)) = 0.0`. + `rectifier(decrease(V, dVdt) - strength(state, fixed_point)) = 0.0`. If the dynamics truly have a fixed point at `fixed_point` and `dVdt` has been defined properly in terms of the dynamics, then `dVdt(fixed_point)` will be `0` and there is no need @@ -27,7 +27,7 @@ struct LyapunovDecreaseCondition <: AbstractLyapunovDecreaseCondition check_decrease::Bool decrease::Function strength::Function - relu::Function + rectifier::Function end function check_decrease(cond::LyapunovDecreaseCondition)::Bool @@ -36,7 +36,7 @@ end function get_decrease_condition(cond::LyapunovDecreaseCondition) if cond.check_decrease - return (V, dVdt, x, fixed_point) -> cond.relu( + return (V, dVdt, x, fixed_point) -> cond.rectifier( cond.decrease(V(x), dVdt(x)) - cond.strength(x, fixed_point) ) else @@ -45,19 +45,19 @@ function get_decrease_condition(cond::LyapunovDecreaseCondition) end """ - AsymptoticDecrease(; strict, C, relu) + AsymptoticDecrease(; strict, C, rectifier) -Constructs a `LyapunovDecreaseCondition` corresponding to asymptotic decrease. +Construct a `LyapunovDecreaseCondition` corresponding to asymptotic decrease. If `strict` is `false`, the condition is `dV/dt ≤ 0`, and if `strict` is `true`, the condition is `dV/dt ≤ - C | state - fixed_point |^2`. -The inequality is represented by `a ≥ b` <==> `relu(b-a) = 0.0`. +The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. """ function AsymptoticDecrease(; strict::Bool = false, C::Real = 1e-6, - relu = (t) -> max(0.0, t) + rectifier = (t) -> max(zero(t), t) )::LyapunovDecreaseCondition strength = if strict (x, x0) -> -C * (x - x0) ⋅ (x - x0) @@ -69,25 +69,25 @@ function AsymptoticDecrease(; true, (V, dVdt) -> dVdt, strength, - relu + rectifier ) end """ - ExponentialDecrease(k; strict, C, relu) + ExponentialDecrease(k; strict, C, rectifier) -Constructs a `LyapunovDecreaseCondition` corresponding to exponential decrease of rate `k`. +Construct a `LyapunovDecreaseCondition` corresponding to exponential decrease of rate `k`. If `strict` is `false`, the condition is `dV/dt ≤ -k * V`, and if `strict` is `true`, the condition is `dV/dt ≤ -k * V - C * ||state - fixed_point||^2`. -The inequality is represented by `a ≥ b` <==> `relu(b-a) = 0.0`. +The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. """ function ExponentialDecrease( k::Real; strict::Bool = false, C::Real = 1e-6, - relu = (t) -> max(0.0, t) + rectifier = (t) -> max(zero(t), t) )::LyapunovDecreaseCondition strength = if strict (x, x0) -> -C * (x - x0) ⋅ (x - x0) @@ -99,22 +99,22 @@ function ExponentialDecrease( true, (V, dVdt) -> dVdt + k * V, strength, - relu + rectifier ) end """ DontCheckDecrease() -Constructs a `LyapunovDecreaseCondition` which represents not checking for +Construct a `LyapunovDecreaseCondition` which represents not checking for decrease of the Lyapunov function along system trajectories. This is appropriate in cases when the Lyapunov decrease condition has been structurally enforced. """ function DontCheckDecrease()::LyapunovDecreaseCondition return LyapunovDecreaseCondition( false, - (V, dVdt) -> 0.0, + (V, dVdt) -> zero(V), (state, fixed_point) -> 0.0, - (t) -> 0.0 + (t) -> zero(t) ) end diff --git a/src/decrease_conditions_RoA_aware.jl b/src/decrease_conditions_RoA_aware.jl index fba7e1f..c8ae130 100644 --- a/src/decrease_conditions_RoA_aware.jl +++ b/src/decrease_conditions_RoA_aware.jl @@ -1,5 +1,5 @@ """ - RoAAwareDecreaseCondition(check_decrease, decrease, strength, relu, ρ, + RoAAwareDecreaseCondition(check_decrease, decrease, strength, rectifier, ρ, out_of_RoA_penalty) Specifies the form of the Lyapunov conditions to be used, training for a region of @@ -11,7 +11,7 @@ and will instead apply `|out_of_RoA_penalty(V(state), dVdt(state), state, fixed_point, ρ)|^2` when `V(state) > ρ`. The inequality will be approximated by the equation - `relu(decrease(V, dVdt) - strength(state, fixed_point)) = 0.0`. + `rectifier(decrease(V(state), dVdt(state)) - strength(state, fixed_point)) = 0.0`. If the dynamics truly have a fixed point at `fixed_point` and `dVdt` has been defined properly in terms of the dynamics, then `dVdt(fixed_point)` will be `0` and there is no need @@ -40,7 +40,7 @@ struct RoAAwareDecreaseCondition <: AbstractLyapunovDecreaseCondition check_decrease::Bool decrease::Function strength::Function - relu::Function + rectifier::Function ρ::Real out_of_RoA_penalty::Function end @@ -52,7 +52,7 @@ end function get_decrease_condition(cond::RoAAwareDecreaseCondition) if cond.check_decrease return function (V, dVdt, x, fixed_point) - (V(x) ≤ cond.ρ) * cond.relu( + (V(x) ≤ cond.ρ) * cond.rectifier( cond.decrease(V(x), dVdt(x)) - cond.strength(x, fixed_point) ) + (V(x) > cond.ρ) * cond.out_of_RoA_penalty(V(x), dVdt(x), x, fixed_point, @@ -66,7 +66,7 @@ end """ make_RoA_aware(cond; ρ, out_of_RoA_penalty) -Adds awareness of the region of attraction (RoA) estimation task to the supplied +Add awareness of the region of attraction (RoA) estimation task to the supplied `LyapunovDecreaseCondition` `ρ` is the target level such that the RoA will be `{ x : V(x) ≤ ρ }`. @@ -83,7 +83,7 @@ function make_RoA_aware( cond.check_decrease, cond.decrease, cond.strength, - cond.relu, + cond.rectifier, ρ, out_of_RoA_penalty ) diff --git a/src/minimization_conditions.jl b/src/minimization_conditions.jl index e11af3e..ad67436 100644 --- a/src/minimization_conditions.jl +++ b/src/minimization_conditions.jl @@ -6,7 +6,7 @@ Specifies the form of the Lyapunov conditions to be used. If `check_nonnegativity` is `true`, training will attempt to enforce `V(state) ≥ strength(state, fixed_point)`. The inequality will be approximated by the equation - `relu(strength(state, fixed_point) - V(state)) = 0.0`. + `rectifier(strength(state, fixed_point) - V(state)) = 0.0`. If `check_fixed_point` is `true`, then training will also attempt to enforce `V(fixed_point) = 0`. @@ -24,13 +24,13 @@ that must be enforced in training for the Lyapunov function to be uniquely minim `fixed_point`. So, in that case, we would use `check_nonnegativity = false; check_fixed_point = true`. -In either case, `relu = (t) -> max(0.0, t)` exactly represents the inequality, but -approximations of this function are allowed. +In either case, `rectifier = (t) -> max(0.0, t)` exactly represents the inequality, but +differentiable approximations of this function may be employed. """ struct LyapunovMinimizationCondition <: AbstractLyapunovMinimizationCondition check_nonnegativity::Bool strength::Function - relu::Function + rectifier::Function check_fixed_point::Bool end @@ -44,31 +44,31 @@ end function get_minimization_condition(cond::LyapunovMinimizationCondition) if cond.check_nonnegativity - return (V, x, fixed_point) -> cond.relu(cond.strength(x, fixed_point) - V(x)) + return (V, x, fixed_point) -> cond.rectifier(cond.strength(x, fixed_point) - V(x)) else return nothing end end """ - StrictlyPositiveDefinite(C; check_fixed_point, relu) + StrictlyPositiveDefinite(C; check_fixed_point, rectifier) -Constructs a `LyapunovMinimizationCondition` representing +Construct a `LyapunovMinimizationCondition` representing `V(state) ≥ C * ||state - fixed_point||^2`. If `check_fixed_point` is `true`, then training will also attempt to enforce `V(fixed_point) = 0`. -The inequality is represented by `a ≥ b` <==> `relu(b-a) = 0.0`. +The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. """ function StrictlyPositiveDefinite(; check_fixed_point = true, C::Real = 1e-6, - relu = (t) -> max(0.0, t) + rectifier = (t) -> max(zero(t), t) )::LyapunovMinimizationCondition LyapunovMinimizationCondition( true, (state, fixed_point) -> C * (state - fixed_point) ⋅ (state - fixed_point), - relu, + rectifier, check_fixed_point ) end @@ -76,21 +76,21 @@ end """ PositiveSemiDefinite(check_fixed_point) -Constructs a `LyapunovMinimizationCondition` representing +Construct a `LyapunovMinimizationCondition` representing `V(state) ≥ 0`. If `check_fixed_point` is `true`, then training will also attempt to enforce `V(fixed_point) = 0`. -The inequality is represented by `a ≥ b` <==> `relu(b-a) = 0.0`. +The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. """ function PositiveSemiDefinite(; check_fixed_point = true, - relu = (t) -> max(0.0, t) + rectifier = (t) -> max(zero(t), t) )::LyapunovMinimizationCondition LyapunovMinimizationCondition( true, (state, fixed_point) -> 0.0, - relu, + rectifier, check_fixed_point ) end @@ -98,7 +98,7 @@ end """ DontCheckNonnegativity(check_fixed_point) -Constructs a `LyapunovMinimizationCondition` which represents not checking for nonnegativity +Construct a `LyapunovMinimizationCondition` which represents not checking for nonnegativity of the Lyapunov function. This is appropriate in cases where this condition has been structurally enforced. @@ -110,7 +110,7 @@ function DontCheckNonnegativity(; check_fixed_point = false)::LyapunovMinimizati LyapunovMinimizationCondition( false, (state, fixed_point) -> 0.0, - (t) -> 0.0, + (t) -> zero(t), check_fixed_point ) end diff --git a/src/policy_search.jl b/src/policy_search.jl index 93a4e2d..2d519f3 100644 --- a/src/policy_search.jl +++ b/src/policy_search.jl @@ -1,12 +1,12 @@ """ add_policy_search(lyapunov_structure, new_dims, control_structure) -Adds dependence on the neural network to the dynamics in a `NeuralLyapunovStructure` +Add dependence on the neural network to the dynamics in a `NeuralLyapunovStructure`. -Adds `new_dims` outputs to the neural network and feeds them through `control_structure` to -calculatethe contribution of the neural network to the dynamics. -The existing `lyapunov_structure.network_dim` dimensions are used as in `lyapunov_structure` -to calculate the Lyapunov function. +Add `new_dims` outputs to the neural network and feeds them through `control_structure` to +calculate the contribution of the neural network to the dynamics. +Use the existing `lyapunov_structure.network_dim` dimensions as in `lyapunov_structure` to +calculate the Lyapunov function. `lyapunov_structure` should assume in its `V̇` that the dynamics take a form `f(x, p, t)`. The returned `NeuralLyapunovStructure` will assume instead `f(x, u, p, t)`, where `u` is the @@ -51,7 +51,7 @@ end """ get_policy(phi, θ, network_func, dim; control_structure) -Returns the control policy as a function of the state +Generate a Julia function representing the control policy as a function of the state The returned function can operate on a state vector or columnwise on a matrix of state vectors. diff --git a/src/structure_specification.jl b/src/structure_specification.jl index b1f6af5..0b356c6 100644 --- a/src/structure_specification.jl +++ b/src/structure_specification.jl @@ -1,7 +1,7 @@ """ UnstructuredNeuralLyapunov() -Creates a `NeuralLyapunovStructure` where the Lyapunov function is the neural network +Create a `NeuralLyapunovStructure` where the Lyapunov function is the neural network evaluated at the state. This does not structurally enforce any Lyapunov conditions. Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For @@ -21,7 +21,7 @@ end """ NonnegativeNeuralLyapunov(network_dim, δ, pos_def; grad_pos_def, grad) -Creates a `NeuralLyapunovStructure` where the Lyapunov function is the L2 norm of the neural +Create a `NeuralLyapunovStructure` where the Lyapunov function is the L2 norm of the neural network output plus a constant δ times a function `pos_def`. The condition that the Lyapunov function must be minimized uniquely at the fixed point can @@ -85,7 +85,7 @@ end """ PositiveSemiDefiniteStructure(network_dim; pos_def, non_neg, grad_pos_def, grad_non_neg, grad) -Creates a `NeuralLyapunovStructure` where the Lyapunov function is the product of a positive +Create a `NeuralLyapunovStructure` where the Lyapunov function is the product of a positive (semi-)definite function `pos_def` which does not depend on the network and a nonnegative function non_neg which does depend the network. From 97c24e38591a19a4d578778874c04636d01ad55e Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 28 Mar 2024 12:22:54 -0400 Subject: [PATCH 04/10] Renamed `NumericalNeuralLyapunovFunctions`, moved it to its own file, and condensed to one method, with performance improvements --- src/NeuralLyapunov.jl | 33 +++-- src/NeuralLyapunovPDESystem.jl | 138 +------------------ src/{local_Lyapunov.jl => local_lyapunov.jl} | 4 +- src/numerical_lyapunov_functions.jl | 93 +++++++++++++ test/damped_pendulum.jl | 5 +- test/damped_sho.jl | 12 +- test/inverted_pendulum.jl | 2 +- test/inverted_pendulum_ODESystem.jl | 2 +- test/roa_estimation.jl | 2 +- 9 files changed, 133 insertions(+), 158 deletions(-) rename src/{local_Lyapunov.jl => local_lyapunov.jl} (94%) create mode 100644 src/numerical_lyapunov_functions.jl diff --git a/src/NeuralLyapunov.jl b/src/NeuralLyapunov.jl index b6e41b7..9952245 100644 --- a/src/NeuralLyapunov.jl +++ b/src/NeuralLyapunov.jl @@ -12,16 +12,31 @@ include("minimization_conditions.jl") include("decrease_conditions.jl") include("decrease_conditions_RoA_aware.jl") include("NeuralLyapunovPDESystem.jl") -include("local_Lyapunov.jl") +include("numerical_lyapunov_functions.jl") +include("local_lyapunov.jl") include("policy_search.jl") -export NeuralLyapunovPDESystem, NumericalNeuralLyapunovFunctions -export local_Lyapunov -export NeuralLyapunovSpecification, NeuralLyapunovStructure, UnstructuredNeuralLyapunov, - NonnegativeNeuralLyapunov, PositiveSemiDefiniteStructure, - LyapunovMinimizationCondition, StrictlyPositiveDefinite, PositiveSemiDefinite, - DontCheckNonnegativity, LyapunovDecreaseCondition, AsymptoticDecrease, - ExponentialDecrease, DontCheckDecrease, RoAAwareDecreaseCondition, make_RoA_aware, - add_policy_search, get_policy +# Lyapunov function structures +export NeuralLyapunovStructure, UnstructuredNeuralLyapunov, NonnegativeNeuralLyapunov, + PositiveSemiDefiniteStructure, get_numerical_lyapunov_function + +# Minimization conditions +export LyapunovMinimizationCondition, StrictlyPositiveDefinite, PositiveSemiDefinite, + DontCheckNonnegativity + +# Decrease conditions +export LyapunovDecreaseCondition, AsymptoticDecrease, ExponentialDecrease, DontCheckDecrease + +# Setting up the PDESystem for NeuralPDE +export NeuralLyapunovSpecification, NeuralLyapunovPDESystem + +# Region of attraction handling +export RoAAwareDecreaseCondition, make_RoA_aware + +# Policy search +export add_policy_search, get_policy + +# Local Lyapunov analysis +export local_lyapunov end diff --git a/src/NeuralLyapunovPDESystem.jl b/src/NeuralLyapunovPDESystem.jl index 2d35ed3..90aa0bc 100644 --- a/src/NeuralLyapunovPDESystem.jl +++ b/src/NeuralLyapunovPDESystem.jl @@ -2,12 +2,7 @@ NeuralLyapunovPDESystem(dynamics::ODESystem, bounds, spec; ) NeuralLyapunovPDESystem(dynamics::Function, lb, ub, spec; ) -Construct and return a `PDESystem` representing the specified neural Lyapunov problem, along -with a function representing the neural network. - -The returned neural network function takes three inputs: the neural network structure `phi`, -the trained network parameters, and a matrix of inputs, then operates columnwise on the -inputs. +Construct a `ModelingToolkit.PDESystem` representing the specified neural Lyapunov problem. # Arguments - `dynamics`: the dynamical system being analyzed, represented as an `ODESystem` or the @@ -178,7 +173,7 @@ function NeuralLyapunovPDESystem( end ########################### Construct PDESystem ########################### - _NeuralLyapunovPDESystem( + return _NeuralLyapunovPDESystem( f, domains, spec, @@ -274,7 +269,7 @@ function _NeuralLyapunovPDESystem( end ########################### Construct PDESystem ########################### - lyapunov_pde_system = PDESystem( + return PDESystem( eqs, bcs, domains, @@ -284,131 +279,4 @@ function _NeuralLyapunovPDESystem( defaults = defaults, name = name ) - - return lyapunov_pde_system -end - -""" - NumericalNeuralLyapunovFunctions(phi, θ, structure, dynamics, fixed_point; jac, J_net) - -Returns the Lyapunov function, its time derivative, and its gradient: `V(state)`, -`V̇(state)`, and `∇V(state)` - -These functions can operate on a state vector or columnwise on a matrix of state vectors. -`phi` is the neural network with parameters `θ`. - -The Lyapunov function structure is specified in structure, which is a -`NeuralLyapunovStructure`. The Jacobian of the network is either specified via -`J_net(phi, θ, state)` or calculated using `jac`, which defaults to `ForwardDiff.jacobian`. -""" -function NumericalNeuralLyapunovFunctions( - phi, - θ, - structure::NeuralLyapunovStructure, - dynamics::Function, - fixed_point; - p = SciMLBase.NullParameters(), - jac = ForwardDiff.jacobian, - J_net = Nothing -)::Tuple{Function, Function, Function} - # network_func is the numerical form of neural network output - output_dim = structure.network_dim - function network_func(x) - reduce( - vcat, - Array(phi[i](x, θ.depvar[Symbol(:φ, i)])) for i in 1:output_dim - ) - end - - # Make Jacobian of networkfunction - _J_net = if isnothing(J_net) - (x) -> jac(network_func, x) - else - (x) -> J_net(phi, θ, x) - end - - # Numerical form of Lyapunov function - V_func(state::AbstractVector) = structure.V(network_func, state, fixed_point) - V_func(state::AbstractMatrix) = mapslices(V_func, state, dims = [1]) - - # Numerical gradient of Lyapunov function - ∇V_func(state::AbstractVector) = structure.∇V( - network_func, - _J_net, - state, - fixed_point - ) - ∇V_func(state::AbstractMatrix) = mapslices(∇V_func, state, dims = [1]) - - # Numerical time derivative of Lyapunov function - function V̇_func(state::AbstractVector) - structure.V̇( - network_func, - _J_net, - dynamics, - state, - p, - 0.0, - fixed_point - ) - end - V̇_func(state::AbstractMatrix) = mapslices(V̇_func, state, dims = [1]) - - return V_func, V̇_func, ∇V_func -end - -""" - NumericalNeuralLyapunovFunctions(phi, θ, network_dim, V_structure, dynamics, - fixed_point; p, grad, deriv) - -Returns the Lyapunov function, its time derivative, and its gradient: `V(state)`, -`V̇(state)`, and `∇V(state)`. - -These functions can operate on a state vector or columnwise on a matrix of state vectors. -`phi` is the neural network with parameters `θ`. `network_dim` is the output dimension of -the neural network. - -The Lyapunov function structure is defined by - `V_structure(_network_func, state, fixed_point)` -Its gradient is calculated using `grad`, which defaults to `ForwardDiff.gradient`. - -The Lyapunov function's time derivative is calculated using `deriv`, which defaults to -`ForwardDiff.derivative`, along with `dynamics`, which should be the closed-loop dynamics -`ẋ = dynamics(x, p, t)`; only `t = 0.0` is used. Parameters `p` can be supplied. -""" -function NumericalNeuralLyapunovFunctions( - phi, - θ, - network_dim, - V_structure::Function, - dynamics::Function, - fixed_point; - p = SciMLBase.NullParameters(), - grad = ForwardDiff.gradient, - deriv = ForwardDiff.derivative -)::Tuple{Function, Function, Function} - # network_func is the numerical form of neural network output - function network_func(x) - reduce( - vcat, - Array(phi[i](x, θ.depvar[Symbol(:φ, i)])) for i in 1:network_dim - ) - end - - # Numerical form of Lyapunov function - V_func(state::AbstractVector) = V_structure(network_func, state, fixed_point) - V_func(state::AbstractMatrix) = mapslices(V_func, state, dims = [1]) - - # Numerical gradient of Lyapunov function - ∇V_func(state::AbstractVector) = grad(V_func, state) - ∇V_func(state::AbstractMatrix) = mapslices(∇V_func, state, dims = [1]) - - # Numerical time derivative of Lyapunov function - V̇_func(state::AbstractVector) = deriv( - (δt) -> V_func(state + δt * dynamics(state, p, 0.0)), - 0.0 - ); - V̇_func(state::AbstractMatrix) = mapslices(V̇_func, state, dims = [1]) - - return V_func, V̇_func, ∇V_func end diff --git a/src/local_Lyapunov.jl b/src/local_lyapunov.jl similarity index 94% rename from src/local_Lyapunov.jl rename to src/local_lyapunov.jl index 7cbd9e4..6d4f0ea 100644 --- a/src/local_Lyapunov.jl +++ b/src/local_lyapunov.jl @@ -1,5 +1,5 @@ """ - get_local_Lyapunov(dynamics, state_dim; fixed_point, dynamics_jac) + get_local_lyapunov(dynamics, state_dim; fixed_point, dynamics_jac) Use semidefinite programming to derive a quadratic Lyapunov function for the linearization of `dynamics` around `fixed_point`. @@ -10,7 +10,7 @@ calculated using `ForwardDiff`. Other allowable forms are a function which takes state and outputs the jacobian of `dynamics` or an `AbstractMatrix` representing the Jacobian at `fixed_point`. If `fixed_point` is not specified, it defaults to the origin. """ -function local_Lyapunov(dynamics::Function, state_dim, optimizer_factory; +function local_lyapunov(dynamics::Function, state_dim, optimizer_factory; fixed_point = zeros(state_dim), dynamics_jac = nothing, p = SciMLBase.NullParameters()) # Linearize the dynamics diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl new file mode 100644 index 0000000..d449b2b --- /dev/null +++ b/src/numerical_lyapunov_functions.jl @@ -0,0 +1,93 @@ +""" + get_numerical_lyapunov_function(phi, θ, structure, dynamics, fixed_point; + ) + +Combine Lyapunov function structure, dynamics, and neural network weights to generate Julia +functions representing the Lyapunov function and its time derivative: ``V(x), V̇(x)``. + +These functions can operate on a state vector or columnwise on a matrix of state vectors. + +# Arguments +- `phi`, `θ`: `phi` is the neural network with parameters `θ`. +- `structure`: a [`NeuralLyapunovStructure`](@ref) representing the structure of the neural + Lyapunov function. +- `dynamics`: the system dynamics, as a function to be used in conjunction with + `structure.f_call`. +- `fixed_point`: the equilibrium point being analyzed by the Lyapunov function. +- `p`: parameters to be passed into `dynamics`; defaults to `SciMLBase.NullParameters()`. +- `use_V̇_structure`: when `true`, ``V̇(x)`` is calculated using `structure.V̇`; when `false`, + ``V̇(x)`` is calculated using `deriv` as ``\\frac{d}{dt} V(x + t f(x))`` at + ``t = 0``; defaults to `false`, as it is more efficient in many cases. +- `deriv`: a function for calculating derivatives; defaults to (and expects same arguments + as) `ForwardDiff.derivative`; only used when `use_V̇_structure` is `false`. +- `jac`: a function for calculating Jacobians; defaults to (and expects same arguments as) + `ForwardDiff.jacobian`; only used when `use_V̇_structure` is `true`. +- `J_net`: the Jacobian of the neural network, specified as a function + `J_net(phi, θ, state)`; if `isnothing(J_net)` (as is the default), `J_net` will be + calculated using `jac`; only used when `use_V̇_structure` is `true`. +""" +function get_numerical_lyapunov_function( + phi, + θ, + structure::NeuralLyapunovStructure, + dynamics::Function, + fixed_point; + p = SciMLBase.NullParameters(), + use_V̇_structure = false, + deriv = ForwardDiff.derivative, + jac = ForwardDiff.jacobian, + J_net = nothing +)::Tuple{Function, Function} + # network_func is the numerical form of neural network output + output_dim = structure.network_dim + network_func = let φ = phi, _θ = θ, dim = output_dim + function (x) + reduce( + vcat, + Array(φ[i](x, _θ.depvar[Symbol(:φ, i)])) for i in 1:dim + ) + end + end + + # V is the numerical form of Lyapunov function + V = let V_structure = structure.V, net = network_func, x0 = fixed_point + V(state::AbstractVector) = V_structure(net, state, x0) + V(state::AbstractMatrix) = mapslices(V, state, dims = [1]) + end + + if use_V̇_structure + # Make Jacobian of network_func + network_jacobian = if isnothing(J_net) + let net = network_func, J = jac + (x) -> J(net, x) + end + else + let _J_net = J_net, φ = phi, _θ = θ + (x) -> _J_net(φ, _θ, x) + end + end + + let _V = V, V̇_structure = structure.V̇, net = network_func, x0 = fixed_point, + f = dynamics, params = p, _J_net = network_jacobian + + # Numerical time derivative of Lyapunov function + V̇(state::AbstractVector) = V̇_structure(net, _J_net, f, state, params, 0.0, x0) + V̇(state::AbstractMatrix) = mapslices(V̇, state, dims = [1]) + + return _V, V̇ + end + else + let f_call = structure.f_call, _V = V, net = network_func, f = dynamics, params = p, + _deriv = deriv + + # Numerical time derivative of Lyapunov function + V̇(state::AbstractVector) = _deriv( + (δt) -> _V(state + δt * f_call(f, net, state, params, 0.0)), + 0.0 + ) + V̇(state::AbstractMatrix) = mapslices(V̇, state, dims = [1]) + + return _V, V̇ + end + end +end diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl index 48cd23d..ef5aa35 100644 --- a/test/damped_pendulum.jl +++ b/test/damped_pendulum.jl @@ -103,11 +103,10 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, res.u, - dim_output, - structure.V, + structure, ODEFunction(dynamics), zeros(length(bounds)); p = p diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 57e4def..2b1d397 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -50,7 +50,7 @@ minimization_condition = DontCheckNonnegativity(check_fixed_point = true) κ = 20.0 decrease_condition = AsymptoticDecrease( strict = true, - relu = (t) -> log(1.0 + exp(κ * t)) / κ + rectifier = (t) -> log(1.0 + exp(κ * t)) / κ ) # Construct neural Lyapunov specification @@ -82,14 +82,14 @@ prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, res.u, - dim_output, - structure.V, - ODEFunction(dynamics), + structure, + f, zeros(2); - p = p + p = p, + use_V̇_structure = true ) ################################## Simulate ################################### diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl index b28ac6d..2122103 100644 --- a/test/inverted_pendulum.jl +++ b/test/inverted_pendulum.jl @@ -102,7 +102,7 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, res.u, structure, diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index ce70a76..ab03bf6 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -110,7 +110,7 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, res.u, structure, diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index cd57d85..347d58c 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -66,7 +66,7 @@ prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = NumericalNeuralLyapunovFunctions( +V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, res.u, structure, From d104c83dacef60fc282dea29f2b44f5e62c15f71 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 28 Mar 2024 20:04:01 -0400 Subject: [PATCH 05/10] Cleaned up and added test for `local_lyapunov` --- Project.toml | 3 +- src/local_lyapunov.jl | 90 +++++++++++++++++++++++++++++------------- test/local_lyapunov.jl | 40 +++++++++++++++++++ test/runtests.jl | 3 ++ 4 files changed, 107 insertions(+), 29 deletions(-) create mode 100644 test/local_lyapunov.jl diff --git a/Project.toml b/Project.toml index e950ba6..2ac2161 100644 --- a/Project.toml +++ b/Project.toml @@ -34,6 +34,7 @@ OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +CSDP = "0a46da34-8e4b-519e-b418-48813639ff34" [targets] -test = ["SafeTestsets", "Test", "Lux", "Optimization", "OptimizationOptimJL", "OptimizationOptimisers", "NLopt", "Random", "NeuralPDE", "DifferentialEquations"] +test = ["SafeTestsets", "Test", "Lux", "Optimization", "OptimizationOptimJL", "OptimizationOptimisers", "NLopt", "Random", "NeuralPDE", "DifferentialEquations", "CSDP"] diff --git a/src/local_lyapunov.jl b/src/local_lyapunov.jl index 6d4f0ea..a65e571 100644 --- a/src/local_lyapunov.jl +++ b/src/local_lyapunov.jl @@ -1,36 +1,42 @@ """ - get_local_lyapunov(dynamics, state_dim; fixed_point, dynamics_jac) + get_local_lyapunov(dynamics, state_dim, optimizer_factory[, jac]; fixed_point, p) -Use semidefinite programming to derive a quadratic Lyapunov function for the -linearization of `dynamics` around `fixed_point`. -Return `(V, dV/dt, ∇V)`. +Use semidefinite programming to derive a quadratic Lyapunov function for the linearization +of `dynamics` around `fixed_point`. Return `(V, dV/dt)`. -If `isnothing(dynamics_jac)`, the Jacobian of the `dynamics(x, p, t)` with respect to `x` is -calculated using `ForwardDiff`. Other allowable forms are a function which takes in the -state and outputs the jacobian of `dynamics` or an `AbstractMatrix` representing the -Jacobian at `fixed_point`. If `fixed_point` is not specified, it defaults to the origin. +To solve the semidefinite program, `JuMP.Model` requires an `optimizer_factory` capable of +semidefinite programming (SDP). See the +[JuMP documentation](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) for +examples. + +If `jac` is not supplied, the Jacobian of the `dynamics(x, p, t)` with respect to `x` is +calculated using `ForwardDiff`. Otherwise, `jac` is expected to be either a function or an +`AbstractMatrix`. If `jac isa Function`, it should take in the state and parameters and +output the Jacobian of `dynamics` with respect to the state `x`. If `jac isa +AbstractMatrix`, it should be the value of the Jacobian at `fixed_point`. + +If `fixed_point` is not specified, it defaults to the origin, i.e., `zeros(state_dim)`. +Parameters `p` for the dynamics should be supplied when the dynamics depend on them. """ -function local_lyapunov(dynamics::Function, state_dim, optimizer_factory; - fixed_point = zeros(state_dim), dynamics_jac = nothing, - p = SciMLBase.NullParameters()) - # Linearize the dynamics - A = if isnothing(dynamics_jac) - ForwardDiff.jacobian(x -> dynamics(x, p, 0.0), fixed_point) - elseif dynamics_jac isa AbstractMatrix - dynamics_jac - elseif dynamics_jac isa Function - dynamics_jac(fixed_point) - else - throw(ErrorException("Unable to calculate Jacobian from dynamics_jac.")) - end - - # Use quadratic semidefinite programming to calculate a Lyapunov function - # for the linearized system +function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, + jac::AbstractMatrix{T}; fixed_point = zeros(T, state_dim), + p = SciMLBase.NullParameters()) where T <: Number + model = JuMP.Model(optimizer_factory) JuMP.set_silent(model) + + # If jac is a matrix A, then the linearization is ẋ = A (x - x0), where + # x is the state and x0 is the fixed point. + # A quadratic Lyapunov function is V(x) = (x - x0) ⋅ (P * (x - x0)) for some positive + # definite matrix P, where V̇(x) = -(x - x0) ⋅ (Q * (x - x0)) for some positive definite + # matrix Q JuMP.@variable(model, P[1:state_dim, 1:state_dim], PSD) JuMP.@variable(model, Q[1:state_dim, 1:state_dim], PSD) - JuMP.@constraint(model, P * A + transpose(A) * P.==-Q) + + # V̇(x) = (x - x0) ⋅ ( (P A + A^T P) * (x - x0)), so we require P A + A^T P = -Q + JuMP.@constraint(model, P * jac + transpose(jac) * P.==-Q) + + # Solve the semidefinite program and get the value of P JuMP.optimize!(model) Psol = JuMP.value.(P) @@ -39,12 +45,40 @@ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory; V(states::AbstractMatrix) = mapslices(V, states, dims = [1]) # Numerical gradient of Lyapunov function - ∇V(state::AbstractVector) = 2 * ( Psol * (state - fixed_point) ) - ∇V(states::AbstractMatrix) = mapslices(∇V, states, dims = [1]) +# ∇V(state::AbstractVector) = 2 * ( Psol * (state - fixed_point) ) +# ∇V(states::AbstractMatrix) = mapslices(∇V, states, dims = [1]) # Numerical time derivative of Lyapunov function V̇(state::AbstractVector) = 2 * dot(dynamics(state, p, 0.0), Psol, state - fixed_point) V̇(states::AbstractMatrix) = mapslices(V̇, states, dims = [1]) - return V, V̇, ∇V + return V, V̇ +end + +function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, jac::Function; + fixed_point = zeros(state_dim), p = SciMLBase.NullParameters()) + + A::AbstractMatrix = jac(fixed_point, p) + return local_lyapunov( + dynamics, + state_dim, + optimizer_factory, + A; + fixed_point = fixed_point, + p = p + ) +end + +function local_lyapunov(dynamics::Function, state_dim, optimizer_factory; + fixed_point = zeros(state_dim), p = SciMLBase.NullParameters()) + + A::AbstractMatrix = ForwardDiff.jacobian(x -> dynamics(x, p, 0.0), fixed_point) + return local_lyapunov( + dynamics, + state_dim, + optimizer_factory, + A; + fixed_point = fixed_point, + p = p + ) end diff --git a/test/local_lyapunov.jl b/test/local_lyapunov.jl new file mode 100644 index 0000000..644cd81 --- /dev/null +++ b/test/local_lyapunov.jl @@ -0,0 +1,40 @@ +using NeuralLyapunov, CSDP, ForwardDiff, Test, LinearAlgebra + +################################### ẋ = -x ################################### +# Set up dynamics with a known Lyapunov function: V(x) = C x ⋅ x +f1(x, p, t) = -x + +# Set up Jacobian +J1(x, p) = (-1.0*I)(3) + +# Find Lyapunov function and rescale for C = 1 +_V1, _V̇1 = local_lyapunov(f1, 3, CSDP.Optimizer, J1) +V1(x) = _V1(x) ./ _V1([0, 0, 1]) + +# Test equality +xs = -1:0.1:1 +states = Iterators.map(collect, Iterators.product(xs, xs, xs)) +errs1 = @. abs(V1(states) - states ⋅ states) +@test all(errs1 .< 1e-10) + +########################## Simple harmonic oscillator ######################### +function f2(state, p, t) + ζ, ω_0 = p + pos, vel = state + vcat(vel, -2ζ * vel - ω_0^2 * pos) +end +p2 = [3.2, 5.1] + +# Find Lyapunov function and derivative +V2, V̇2 = local_lyapunov(f2, 2, CSDP.Optimizer; p = p2) +dV2dt = (state) -> ForwardDiff.derivative((δt) -> V2(state + δt * f2(state, p2, 0.0)), 0.0) + +# Test V̇2 = d/dt V2 +xs = -1:0.03:1 +states = Iterators.map(collect, Iterators.product(xs, xs)) +errs2 = @. abs(dV2dt(states) - V̇2(states)) +@test all(errs2 .< 1e-10) + +# Test V̇2 ≺ 0 +@test all( @. V̇2(states) < 0) +@test V̇2(zeros(2)) == 0.0 diff --git a/test/runtests.jl b/test/runtests.jl index ef1120d..07ce606 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,4 +16,7 @@ using SafeTestsets @time @safetestset "Policy search - inverted pendulum 2" begin include("inverted_pendulum_ODESystem.jl") end + @time @safetestset "Local Lyapunov function search" begin + include("local_lyapunov.jl") + end end From 4f0fef1cff48f31298648b3f806dec3ac810aeb1 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 29 Mar 2024 15:29:40 -0400 Subject: [PATCH 06/10] Tweak `local_lyapunov` tests to improve code coverage --- test/local_lyapunov.jl | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/test/local_lyapunov.jl b/test/local_lyapunov.jl index 644cd81..91bc5d4 100644 --- a/test/local_lyapunov.jl +++ b/test/local_lyapunov.jl @@ -12,17 +12,22 @@ _V1, _V̇1 = local_lyapunov(f1, 3, CSDP.Optimizer, J1) V1(x) = _V1(x) ./ _V1([0, 0, 1]) # Test equality -xs = -1:0.1:1 -states = Iterators.map(collect, Iterators.product(xs, xs, xs)) -errs1 = @. abs(V1(states) - states ⋅ states) +xs1 = -1:0.1:1 +states1 = Iterators.map(collect, Iterators.product(xs1, xs1, xs1)) +errs1 = @. abs(V1(states1) - states1 ⋅ states1) @test all(errs1 .< 1e-10) ########################## Simple harmonic oscillator ######################### -function f2(state, p, t) +function f2(state::AbstractVector, p, t) ζ, ω_0 = p pos, vel = state vcat(vel, -2ζ * vel - ω_0^2 * pos) end +function f2(states::AbstractMatrix, p, t) + ζ, ω_0 = p + pos, vel = states[1,:], states[2,:] + vcat(transpose(vel), transpose(-2ζ * vel - ω_0^2 * pos)) +end p2 = [3.2, 5.1] # Find Lyapunov function and derivative @@ -30,11 +35,11 @@ V2, V̇2 = local_lyapunov(f2, 2, CSDP.Optimizer; p = p2) dV2dt = (state) -> ForwardDiff.derivative((δt) -> V2(state + δt * f2(state, p2, 0.0)), 0.0) # Test V̇2 = d/dt V2 -xs = -1:0.03:1 -states = Iterators.map(collect, Iterators.product(xs, xs)) -errs2 = @. abs(dV2dt(states) - V̇2(states)) +xs2 = -1:0.03:1 +states2 = hcat(Iterators.map(collect, Iterators.product(xs2, xs2))...) +errs2 = abs.( V̇2(states2) .- dV2dt(states2)) @test all(errs2 .< 1e-10) # Test V̇2 ≺ 0 -@test all( @. V̇2(states) < 0) +@test all( V̇2(states2) .< 0) @test V̇2(zeros(2)) == 0.0 From a9dd80aea34464879ffb4971c9984ad8550cb936 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 29 Mar 2024 18:32:01 -0400 Subject: [PATCH 07/10] Updated SHO/pendulum parameter definitions to match convention and modified damped SHO and damped pendulum tests for greater code coverage --- test/damped_pendulum.jl | 8 ++++++-- test/damped_sho.jl | 27 ++++++++++++--------------- test/inverted_pendulum.jl | 2 +- test/inverted_pendulum_ODESystem.jl | 2 +- test/local_lyapunov.jl | 4 ++-- 5 files changed, 22 insertions(+), 21 deletions(-) diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl index ef5aa35..20c0df1 100644 --- a/test/damped_pendulum.jl +++ b/test/damped_pendulum.jl @@ -18,7 +18,7 @@ defaults = Dict([ζ => 0.5, ω_0 => 1.0]) Dt = Differential(t) DDt = Dt^2 -eqs = [DDt(θ) + 2ζ * Dt(θ) + ω_0^2 * sin(θ) ~ 0.0] +eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ 0.0] @named dynamics = ODESystem( eqs, @@ -71,7 +71,11 @@ structure = PositiveSemiDefiniteStructure( minimization_condition = DontCheckNonnegativity(check_fixed_point = false) # Define Lyapunov decrease condition -decrease_condition = AsymptoticDecrease(strict = true) +κ = 20.0 +decrease_condition = AsymptoticDecrease( + strict = true, + rectifier = (t) -> log(1.0 + exp(κ * t)) / κ +) # Construct neural Lyapunov specification spec = NeuralLyapunovSpecification( diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 2b1d397..e2e5044 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -16,10 +16,10 @@ function f(state, p, t) ζ, ω_0 = p pos = state[1] vel = state[2] - vcat(vel, -2ζ * vel - ω_0^2 * pos) + vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos) end -lb = [-2 * pi, -10.0]; -ub = [2 * pi, 10.0]; +lb = [-2 * pi, -2.0]; +ub = [2 * pi, 2.0]; p = [0.5, 1.0] dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) @@ -27,16 +27,16 @@ dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) # Define neural network discretization dim_state = length(lb) -dim_hidden = 15 -dim_output = 2 +dim_hidden = 20 +dim_output = 5 chain = [Lux.Chain( Dense(dim_state, dim_hidden, tanh), Dense(dim_hidden, dim_hidden, tanh), - Dense(dim_hidden, 1, use_bias = false) + Dense(dim_hidden, 1) ) for _ in 1:dim_output] # Define neural network discretization -strategy = GridTraining(0.1) +strategy = GridTraining(0.05) discretization = PhysicsInformedNN(chain, strategy) # Define neural Lyapunov structure @@ -47,11 +47,8 @@ structure = NonnegativeNeuralLyapunov( minimization_condition = DontCheckNonnegativity(check_fixed_point = true) # Define Lyapunov decrease condition -κ = 20.0 -decrease_condition = AsymptoticDecrease( - strict = true, - rectifier = (t) -> log(1.0 + exp(κ * t)) / κ -) +# Damped SHO has exponential decrease at a rate of k = ζ * ω_0, so we train to certify that +decrease_condition = ExponentialDecrease(prod(p)) # Construct neural Lyapunov specification spec = NeuralLyapunovSpecification( @@ -77,9 +74,9 @@ sym_prob = symbolic_discretize(pde_system, discretization) ########################## Solve OptimizationProblem ########################## -res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, BFGS(); maxiters = 300) +res = Optimization.solve(prob, BFGS(); maxiters = 500) ###################### Get numerical numerical functions ###################### V_func, V̇_func = get_numerical_lyapunov_function( @@ -111,7 +108,7 @@ dVdt_predict = vec(V̇_func(hcat(states...))) @test V̇_func([0.0, 0.0]) == 0.0 # V̇ should be negative almost everywhere -@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 1e-4 +@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 1e-5 #= # Get RoA Estimate diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl index 2122103..7d46784 100644 --- a/test/inverted_pendulum.jl +++ b/test/inverted_pendulum.jl @@ -16,7 +16,7 @@ function open_loop_pendulum_dynamics(x, u, p, t) ζ, ω_0 = p τ = u[] return [ω - -2ζ * ω - ω_0^2 * sin(θ) + τ] + -2ζ * ω_0 * ω - ω_0^2 * sin(θ) + τ] end lb = [0.0, -10.0]; diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index ab03bf6..48ccab7 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -18,7 +18,7 @@ defaults = Dict([ζ => 0.5, ω_0 => 1.0]) Dt = Differential(t) DDt = Dt^2 -eqs = [DDt(θ) + 2ζ * Dt(θ) + ω_0^2 * sin(θ) ~ τ] +eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ τ] @named driven_pendulum = ODESystem( eqs, diff --git a/test/local_lyapunov.jl b/test/local_lyapunov.jl index 91bc5d4..84e70f4 100644 --- a/test/local_lyapunov.jl +++ b/test/local_lyapunov.jl @@ -21,12 +21,12 @@ errs1 = @. abs(V1(states1) - states1 ⋅ states1) function f2(state::AbstractVector, p, t) ζ, ω_0 = p pos, vel = state - vcat(vel, -2ζ * vel - ω_0^2 * pos) + vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos) end function f2(states::AbstractMatrix, p, t) ζ, ω_0 = p pos, vel = states[1,:], states[2,:] - vcat(transpose(vel), transpose(-2ζ * vel - ω_0^2 * pos)) + vcat(transpose(vel), transpose(-2ζ * ω_0 * vel - ω_0^2 * pos)) end p2 = [3.2, 5.1] From 2cd47204dd5895ba68641a2c123059f6bab0cfe6 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 29 Mar 2024 19:45:42 -0400 Subject: [PATCH 08/10] ModelingToolkit v9.8 renamed `states` accessor to `unknowns` --- src/NeuralLyapunovPDESystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NeuralLyapunovPDESystem.jl b/src/NeuralLyapunovPDESystem.jl index 90aa0bc..065c01e 100644 --- a/src/NeuralLyapunovPDESystem.jl +++ b/src/NeuralLyapunovPDESystem.jl @@ -110,7 +110,7 @@ function NeuralLyapunovPDESystem( end s_syms, p_syms = if dynamics.sys isa ODESystem - s_syms = Symbol.(operation.(states(dynamics.sys))) + s_syms = Symbol.(operation.(unknowns(dynamics.sys))) p_syms = Symbol.(parameters(dynamics.sys)) (s_syms, p_syms) elseif dynamics.sys isa SciMLBase.SymbolCache @@ -148,7 +148,7 @@ function NeuralLyapunovPDESystem( )::PDESystem ######################### Check for policy search ######################### f, x, p, policy_search = if isempty(ModelingToolkit.unbound_inputs(dynamics)) - (ODEFunction(dynamics), states(dynamics), parameters(dynamics), false) + (ODEFunction(dynamics), unknowns(dynamics), parameters(dynamics), false) else (f, _), x, p = ModelingToolkit.generate_control_function(dynamics; simplify = true) (f, x, p, true) From d25f53441ff2d701cbabbc844f186cc618448af3 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 3 Apr 2024 13:06:54 -0400 Subject: [PATCH 09/10] Created `phi_to_net` and changed parameters to `res.u.depvar` instead of `res.u`. --- src/numerical_lyapunov_functions.jl | 48 +++++++++++++++++++++++------ src/policy_search.jl | 17 +++++----- test/damped_pendulum.jl | 2 +- test/damped_sho.jl | 2 +- test/inverted_pendulum.jl | 4 +-- test/inverted_pendulum_ODESystem.jl | 4 +-- test/roa_estimation.jl | 2 +- 7 files changed, 52 insertions(+), 27 deletions(-) diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl index d449b2b..0285e39 100644 --- a/src/numerical_lyapunov_functions.jl +++ b/src/numerical_lyapunov_functions.jl @@ -8,7 +8,11 @@ functions representing the Lyapunov function and its time derivative: ``V(x), V These functions can operate on a state vector or columnwise on a matrix of state vectors. # Arguments -- `phi`, `θ`: `phi` is the neural network with parameters `θ`. +- `phi`: the neural network, represented as `phi(x, θ)` if the neural network has a single + output, or a `Vector` of the same with one entry per neural network output. +- `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the first + neural network output (even if there is only one), `θ[:φ2]` the parameters of the + second (if there are multiple), and so on. - `structure`: a [`NeuralLyapunovStructure`](@ref) representing the structure of the neural Lyapunov function. - `dynamics`: the system dynamics, as a function to be used in conjunction with @@ -39,15 +43,7 @@ function get_numerical_lyapunov_function( J_net = nothing )::Tuple{Function, Function} # network_func is the numerical form of neural network output - output_dim = structure.network_dim - network_func = let φ = phi, _θ = θ, dim = output_dim - function (x) - reduce( - vcat, - Array(φ[i](x, _θ.depvar[Symbol(:φ, i)])) for i in 1:dim - ) - end - end + network_func = phi_to_net(phi, θ) # V is the numerical form of Lyapunov function V = let V_structure = structure.V, net = network_func, x0 = fixed_point @@ -91,3 +87,35 @@ function get_numerical_lyapunov_function( end end end + +""" + phi_to_net(phi, θ[; idx]) + +Return the network as a function of state alone. + +# Arguments + +- `phi`: the neural network, represented as `phi(state, θ)` if the neural network has a + single output, or a `Vector` of the same with one entry per neural network output. +- `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the first + neural network output (even if there is only one), `θ[:φ2]` the parameters of the + second (if there are multiple), and so on. +- `idx`: the neural network outputs to include in the returned function; defaults to all and + only applicable when `phi isa Vector`. +""" +function phi_to_net(phi, θ) + let _θ = θ, φ = phi + return (state) -> φ(state, _θ[:φ1]) + end +end + +function phi_to_net(phi::Vector, θ; idx = eachindex(phi)) + let _θ = θ, φ = phi, _idx = idx + return function (x) + reduce( + vcat, + Array(φ[i](x, _θ[Symbol(:φ, i)])) for i in _idx + ) + end + end +end diff --git a/src/policy_search.jl b/src/policy_search.jl index 2d519f3..eb099b4 100644 --- a/src/policy_search.jl +++ b/src/policy_search.jl @@ -67,17 +67,14 @@ function get_policy( control_dim::Integer; control_structure::Function = identity ) - function policy(state::AbstractVector) - control_structure( - reduce( - vcat, - Array(phi[i](state, θ.depvar[Symbol(:φ, i)])) - for i in (network_dim - control_dim + 1):network_dim - ) - ) - end + network_func = phi_to_net(phi, θ; idx = (network_dim - control_dim + 1):network_dim) - policy(state::AbstractMatrix) = mapslices(policy, state, dims = [1]) + policy(state::AbstractVector) = control_structure(network_func(state)) + policy(states::AbstractMatrix) = mapslices( + control_structure, + network_func(states), + dims = [1] + ) return policy end diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl index 20c0df1..b5ef2dc 100644 --- a/test/damped_pendulum.jl +++ b/test/damped_pendulum.jl @@ -109,7 +109,7 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, - res.u, + res.u.depvar, structure, ODEFunction(dynamics), zeros(length(bounds)); diff --git a/test/damped_sho.jl b/test/damped_sho.jl index e2e5044..429d810 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -81,7 +81,7 @@ res = Optimization.solve(prob, BFGS(); maxiters = 500) ###################### Get numerical numerical functions ###################### V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, - res.u, + res.u.depvar, structure, f, zeros(2); diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl index 7d46784..6433995 100644 --- a/test/inverted_pendulum.jl +++ b/test/inverted_pendulum.jl @@ -104,14 +104,14 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, - res.u, + res.u.depvar, structure, open_loop_pendulum_dynamics, upright_equilibrium; p = p ) -u = get_policy(discretization.phi, res.u, dim_output, dim_u) +u = get_policy(discretization.phi, res.u.depvar, dim_output, dim_u) ################################## Simulate ################################### diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index 48ccab7..caff5e4 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -112,14 +112,14 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, - res.u, + res.u.depvar, structure, open_loop_pendulum_dynamics, upright_equilibrium; p = p ) -u = get_policy(discretization.phi, res.u, dim_output, dim_u) +u = get_policy(discretization.phi, res.u.depvar, dim_output, dim_u) ################################## Simulate ################################### diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index 347d58c..c3d2418 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -68,7 +68,7 @@ res = Optimization.solve(prob, BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### V_func, V̇_func = get_numerical_lyapunov_function( discretization.phi, - res.u, + res.u.depvar, structure, f, zeros(length(lb)) From 7e1a4c6bb920c15f3dfe51b726f1143196e165d2 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 3 Apr 2024 15:23:31 -0400 Subject: [PATCH 10/10] Renamed "state" to "states" when it's a matrix --- src/numerical_lyapunov_functions.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl index 0285e39..f915d47 100644 --- a/src/numerical_lyapunov_functions.jl +++ b/src/numerical_lyapunov_functions.jl @@ -48,7 +48,7 @@ function get_numerical_lyapunov_function( # V is the numerical form of Lyapunov function V = let V_structure = structure.V, net = network_func, x0 = fixed_point V(state::AbstractVector) = V_structure(net, state, x0) - V(state::AbstractMatrix) = mapslices(V, state, dims = [1]) + V(states::AbstractMatrix) = mapslices(V, states, dims = [1]) end if use_V̇_structure @@ -68,7 +68,7 @@ function get_numerical_lyapunov_function( # Numerical time derivative of Lyapunov function V̇(state::AbstractVector) = V̇_structure(net, _J_net, f, state, params, 0.0, x0) - V̇(state::AbstractMatrix) = mapslices(V̇, state, dims = [1]) + V̇(states::AbstractMatrix) = mapslices(V̇, states, dims = [1]) return _V, V̇ end @@ -81,7 +81,7 @@ function get_numerical_lyapunov_function( (δt) -> _V(state + δt * f_call(f, net, state, params, 0.0)), 0.0 ) - V̇(state::AbstractMatrix) = mapslices(V̇, state, dims = [1]) + V̇(states::AbstractMatrix) = mapslices(V̇, states, dims = [1]) return _V, V̇ end @@ -95,8 +95,8 @@ Return the network as a function of state alone. # Arguments -- `phi`: the neural network, represented as `phi(state, θ)` if the neural network has a - single output, or a `Vector` of the same with one entry per neural network output. +- `phi`: the neural network, represented as `phi(x, θ)` if the neural network has a single + output, or a `Vector` of the same with one entry per neural network output. - `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the first neural network output (even if there is only one), `θ[:φ2]` the parameters of the second (if there are multiple), and so on. @@ -105,7 +105,7 @@ Return the network as a function of state alone. """ function phi_to_net(phi, θ) let _θ = θ, φ = phi - return (state) -> φ(state, _θ[:φ1]) + return (x) -> φ(x, _θ[:φ1]) end end