From 39f244a0128ccc2f6a63fbc97b3c9bdac19e13f0 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 17 May 2024 12:08:08 -0400 Subject: [PATCH 01/34] Organized documenation --- docs/make.jl | 16 +++++++- docs/src/demos/damped_SHO.md | 1 + docs/src/index.md | 5 +-- docs/src/man.md | 18 +++++++++ docs/src/man/decrease.md | 24 ++++++++++++ docs/src/man/internals.md | 5 +++ docs/src/man/local_lyapunov.md | 5 +++ docs/src/man/minimization.md | 25 ++++++++++++ docs/src/man/pdesystem.md | 14 +++++++ docs/src/man/policy_search.md | 9 +++++ docs/src/man/roa.md | 9 +++++ docs/src/man/structure.md | 15 +++++++ src/NeuralLyapunovPDESystem.jl | 61 +++++++++++++++-------------- src/conditions_specification.jl | 12 +++++- src/numerical_lyapunov_functions.jl | 61 +++++++++++++++-------------- src/structure_specification.jl | 16 +++++--- 16 files changed, 225 insertions(+), 71 deletions(-) create mode 100644 docs/src/demos/damped_SHO.md create mode 100644 docs/src/man.md create mode 100644 docs/src/man/decrease.md create mode 100644 docs/src/man/internals.md create mode 100644 docs/src/man/local_lyapunov.md create mode 100644 docs/src/man/minimization.md create mode 100644 docs/src/man/pdesystem.md create mode 100644 docs/src/man/policy_search.md create mode 100644 docs/src/man/roa.md create mode 100644 docs/src/man/structure.md diff --git a/docs/make.jl b/docs/make.jl index 5779690..4de1343 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -14,7 +14,21 @@ makedocs(; assets = String[] ), pages = [ - "Home" => "index.md" + "Home" => "index.md", + "Manual" => [ + "man.md", + "man/pdesystem.md", + "man/structure.md", + "man/minimization.md", + "man/decrease.md", + "man/roa.md", + "man/policy_search.md", + "man/local_lyapunov.md", + hide("man/internals.md") + ], + "Demonstrations" => [ + "demos/damped_SHO.md" + ] ] ) diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md new file mode 100644 index 0000000..80eb854 --- /dev/null +++ b/docs/src/demos/damped_SHO.md @@ -0,0 +1 @@ +# Damped Simple Harmonic Oscillator \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 65267f6..ece41fa 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -6,9 +6,8 @@ CurrentModule = NeuralLyapunov Documentation for [NeuralLyapunov](https://github.com/SciML/NeuralLyapunov.jl). -```@index +```@contents ``` -```@autodocs -Modules = [NeuralLyapunov] +```@index ``` diff --git a/docs/src/man.md b/docs/src/man.md new file mode 100644 index 0000000..8fb53c8 --- /dev/null +++ b/docs/src/man.md @@ -0,0 +1,18 @@ +# Components of a Neural Lyapunov Problem + +For a candidate Lyapunov function $V(x)$ to certify the stability of an equilibrium point $x_0$ of the dynamical system $\dot{x}(t) = f(x(t))$, it must satisfy two conditions: +1. The function $V$ must be uniquely minimized at $x_0$, and +2. The function $V$ must decrease along system trajectories (i.e., $V(x(t))$ decreases as long as $x(t)$ is a trajectory of the dynamical system). + +A neural Lyapunov function represents the candidate Lyapunov function $V$ using a neural network, sometimes modifying the output of the network slightly so as to enforce one of the above conditions. + +Thus, we specify our neural Lyapunov problems with three components, each answering a different question: +1. How is $V$ structured in terms of the neural network? +2. How is the minimization condition to be enforced? +3. How is the decrease condition to be enforced? + +These three components are represented by the three fields of a `NeuralLyapunovSpecification` object. + +```@docs +NeuralLyapunovSpecification +``` diff --git a/docs/src/man/decrease.md b/docs/src/man/decrease.md new file mode 100644 index 0000000..c6572f6 --- /dev/null +++ b/docs/src/man/decrease.md @@ -0,0 +1,24 @@ +# Lyapunov Decrease Condition + +```@docs +LyapunovDecreaseCondition +``` + +## Pre-defined decrease conditions + +```@docs +AsymptoticDecrease +ExponentialDecrease +DontCheckDecrease +``` + +## Defining your own decrease condition + +```@docs +NeuralLyapunov.AbstractLyapunovDecreaseCondition +``` + +```@docs +NeuralLyapunov.check_decrease +NeuralLyapunov.get_decrease_condition +``` diff --git a/docs/src/man/internals.md b/docs/src/man/internals.md new file mode 100644 index 0000000..3a377ee --- /dev/null +++ b/docs/src/man/internals.md @@ -0,0 +1,5 @@ +# Internals + +```@docs +NeuralLyapunov.phi_to_net +``` diff --git a/docs/src/man/local_lyapunov.md b/docs/src/man/local_lyapunov.md new file mode 100644 index 0000000..c4804f8 --- /dev/null +++ b/docs/src/man/local_lyapunov.md @@ -0,0 +1,5 @@ +# Local Lyapunov analysis + +```@docs +local_lyapunov +``` diff --git a/docs/src/man/minimization.md b/docs/src/man/minimization.md new file mode 100644 index 0000000..a2b6433 --- /dev/null +++ b/docs/src/man/minimization.md @@ -0,0 +1,25 @@ +# Lyapunov Minimization Condition + +```@docs +LyapunovMinimizationCondition +``` + +## Pre-defined minimization conditions + +```@docs +PositiveSemiDefinite +DontCheckNonnegativity +StrictlyPositiveDefinite +``` + +## Defining your own minimization condition + +```@docs +NeuralLyapunov.AbstractLyapunovMinimizationCondition +``` + +```@docs +NeuralLyapunov.check_nonnegativity +NeuralLyapunov.check_minimal_fixed_point +NeuralLyapunov.get_minimization_condition +``` diff --git a/docs/src/man/pdesystem.md b/docs/src/man/pdesystem.md new file mode 100644 index 0000000..9d0de95 --- /dev/null +++ b/docs/src/man/pdesystem.md @@ -0,0 +1,14 @@ +# Solving a Neural Lyapunov Problem + +NeuralLyapunov.jl represents neural Lyapunov problems as systems of partial differential equations, using the `ModelingToolkit.PDESystem` type. +Such a `PDESystem` can then be solved using a physics-informed neural network through [NeuralPDE.jl](https://github.com/SciML/NeuralPDE.jl). + +```@docs +NeuralLyapunovPDESystem +``` + +## Extracting the numerical Lyapunov function + +```@docs +get_numerical_lyapunov_function +``` diff --git a/docs/src/man/policy_search.md b/docs/src/man/policy_search.md new file mode 100644 index 0000000..c6be6cc --- /dev/null +++ b/docs/src/man/policy_search.md @@ -0,0 +1,9 @@ +# Policy Search + +```@docs +add_policy_search +``` + +```@docs +get_policy +``` diff --git a/docs/src/man/roa.md b/docs/src/man/roa.md new file mode 100644 index 0000000..4795a2f --- /dev/null +++ b/docs/src/man/roa.md @@ -0,0 +1,9 @@ +# Training for Region of Attraction Identification + +```@docs +RoAAwareDecreaseCondition +``` + +```@docs +make_RoA_aware +``` diff --git a/docs/src/man/structure.md b/docs/src/man/structure.md new file mode 100644 index 0000000..b2d8c26 --- /dev/null +++ b/docs/src/man/structure.md @@ -0,0 +1,15 @@ +# Structuring a neural Lyapunov function + +```@docs +NeuralLyapunovStructure +``` + +## Pre-defined structures + +```@docs +UnstructuredNeuralLyapunov +NonnegativeNeuralLyapunov +PositiveSemiDefiniteStructure +``` + +## Defining your own neural Lyapunov structure diff --git a/src/NeuralLyapunovPDESystem.jl b/src/NeuralLyapunovPDESystem.jl index 065c01e..58306de 100644 --- a/src/NeuralLyapunovPDESystem.jl +++ b/src/NeuralLyapunovPDESystem.jl @@ -4,35 +4,38 @@ Construct a `ModelingToolkit.PDESystem` representing the specified neural Lyapunov problem. -# Arguments -- `dynamics`: the dynamical system being analyzed, represented as an `ODESystem` or the - function `f` such that `ẋ = f(x[, u], p, t)`; either way, the ODE should not depend - on time and only `t = 0.0` will be used -- `bounds`: an array of domains, defining the training domain by bounding the states (and - derivatives, when applicable) of `dynamics`; only used when `dynamics isa - ODESystem`, otherwise use `lb` and `ub`. -- `lb` and `ub`: the training domain will be ``[lb_1, ub_1]×[lb_2, ub_2]×...``; not used - when `dynamics isa ODESystem`, then use `bounds`. -- `spec::NeuralLyapunovSpecification`: defines the Lyapunov function structure, as well as - the minimization and decrease conditions. -- `fixed_point`: the equilibrium being analyzed; defaults to the origin. -- `p`: the values of the parameters of the dynamical system being analyzed; defaults to - `SciMLBase.NullParameters()`; not used when `dynamics isa ODESystem`, then use the - default parameter values of `dynamics`. -- `state_syms`: an array of the `Symbol` representing each state; not used when `dynamics - isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics isa - ODEFunction`, symbols stored there are used, unless overridden here; if not provided - here and cannot be inferred, `[:state1, :state2, ...]` will be used. -- `parameter_syms`: an array of the `Symbol` representing each parameter; not used when - `dynamics isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics - isa ODEFunction`, symbols stored there are used, unless overridden here; if not - provided here and cannot be inferred, `[:param1, :param2, ...]` will be used. -- `policy_search::Bool`: whether or not to include a loss term enforcing `fixed_point` to - actually be a fixed point; defaults to `false`; only used when `dynamics isa - 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` +# Positional Arguments + - `dynamics`: the dynamical system being analyzed, represented as an `ODESystem` or the + function `f` such that `ẋ = f(x[, u], p, t)`; either way, the ODE should not depend on + time and only `t = 0.0` will be used. (For an example of when `f` would have a `u` + argument, see [`add_policy_search`](@ref).) + - `bounds`: an array of domains, defining the training domain by bounding the states (and + derivatives, when applicable) of `dynamics`; only used when `dynamics isa + ODESystem`, otherwise use `lb` and `ub`. + - `lb` and `ub`: the training domain will be ``[lb_1, ub_1]×[lb_2, ub_2]×...``; not used + when `dynamics isa ODESystem`, then use `bounds`. + - `spec`: a [`NeuralLyapunovSpecification`](@ref) defining the Lyapunov function + structure, as well as the minimization and decrease conditions. + +# Keyword Arguments + - `fixed_point`: the equilibrium being analyzed; defaults to the origin. + - `p`: the values of the parameters of the dynamical system being analyzed; defaults to + `SciMLBase.NullParameters()`; not used when `dynamics isa ODESystem`, then use the + default parameter values of `dynamics`. + - `state_syms`: an array of the `Symbol` representing each state; not used when `dynamics + isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics isa ODEFunction`, + symbols stored there are used, unless overridden here; if not provided here and cannot + be inferred, `[:state1, :state2, ...]` will be used. + - `parameter_syms`: an array of the `Symbol` representing each parameter; not used when + `dynamics isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics isa + ODEFunction`, symbols stored there are used, unless overridden here; if not provided + here and cannot be inferred, `[:param1, :param2, ...]` will be used. + - `policy_search::Bool`: whether or not to include a loss term enforcing `fixed_point` to + actually be a fixed point; defaults to `false`; only used when `dynamics isa 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, diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl index 9de4922..8460436 100644 --- a/src/conditions_specification.jl +++ b/src/conditions_specification.jl @@ -46,9 +46,17 @@ All concrete `AbstractLyapunovDecreaseCondition` subtypes should define the abstract type AbstractLyapunovDecreaseCondition end """ - NeuralLyapunovSpecification + NeuralLyapunovSpecification(structure, minimzation_condition, decrease_condition) -Specifies a neural Lyapunov problem +Specifies a neural Lyapunov problem. + +# Fields + - `structure`: a [`NeuralLyapunovStructure`](@ref) specifying the relationship between the + neural network and the candidate Lyapunov function. + - `minimzation_condition`: an [`AbstractLyapunovMinimizationCondition`](@ref) specifying + how the minimization condition will be enforced. + - `decrease_condition`: an [`AbstractLyapunovDecreaseCondition`](@ref) specifying how the + decrease condition will be enforced. """ struct NeuralLyapunovSpecification structure::NeuralLyapunovStructure diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl index f915d47..fb1e5c1 100644 --- a/src/numerical_lyapunov_functions.jl +++ b/src/numerical_lyapunov_functions.jl @@ -7,28 +7,30 @@ 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`: 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 - `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`. +# Positional Arguments + - `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 + `structure.f_call`. + - `fixed_point`: the equilibrium point being analyzed by the Lyapunov function. + +# Keyword Arguments + - `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, @@ -94,14 +96,13 @@ end Return the network as a function of state alone. # Arguments - -- `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. -- `idx`: the neural network outputs to include in the returned function; defaults to all and - only applicable when `phi isa Vector`. + - `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. + - `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 diff --git a/src/structure_specification.jl b/src/structure_specification.jl index 0b356c6..aac8f8e 100644 --- a/src/structure_specification.jl +++ b/src/structure_specification.jl @@ -1,7 +1,7 @@ """ UnstructuredNeuralLyapunov() -Create a `NeuralLyapunovStructure` where the Lyapunov function is the neural network +Create a [`NeuralLyapunovStructure`](@ref) 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,8 +21,8 @@ end """ NonnegativeNeuralLyapunov(network_dim, δ, pos_def; grad_pos_def, grad) -Create a `NeuralLyapunovStructure` where the Lyapunov function is the L2 norm of the neural -network output plus a constant δ times a function `pos_def`. +Create a [`NeuralLyapunovStructure`](@ref) 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 be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This @@ -40,6 +40,8 @@ The neural network output has dimension `network_dim`. Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For `f(state, input, p, t)`, consider using `add_policy_search`. + +See also: [`DontCheckNonnegativity`](@ref) """ function NonnegativeNeuralLyapunov( network_dim::Integer; @@ -85,9 +87,9 @@ end """ PositiveSemiDefiniteStructure(network_dim; pos_def, non_neg, grad_pos_def, grad_non_neg, grad) -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. +Create a [`NeuralLyapunovStructure`](@ref) 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. The condition that the Lyapunov function must be minimized uniquely at the fixed point can be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This @@ -106,6 +108,8 @@ The neural network output has dimension `network_dim`. Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For `f(state, input, p, t)`, consider using `add_policy_search`. + +Typically used with [`DontCheckNonnegativity`](@ref). """ function PositiveSemiDefiniteStructure( network_dim::Integer; From a67ca6b9aa17a77e5d9f47ae1b535b2d782ecae1 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 17 May 2024 13:56:58 -0400 Subject: [PATCH 02/34] Clarified documentation for `get_numerical_lyapunov_function` --- docs/src/man/pdesystem.md | 5 +++++ src/NeuralLyapunovPDESystem.jl | 6 +++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/docs/src/man/pdesystem.md b/docs/src/man/pdesystem.md index 9d0de95..4270337 100644 --- a/docs/src/man/pdesystem.md +++ b/docs/src/man/pdesystem.md @@ -3,12 +3,17 @@ NeuralLyapunov.jl represents neural Lyapunov problems as systems of partial differential equations, using the `ModelingToolkit.PDESystem` type. Such a `PDESystem` can then be solved using a physics-informed neural network through [NeuralPDE.jl](https://github.com/SciML/NeuralPDE.jl). +Candidate Lyapunov functions will be trained within a box domain subset of the state space. + ```@docs NeuralLyapunovPDESystem ``` ## Extracting the numerical Lyapunov function +We provide the following convenience function for generating the Lyapunov function after the parameters have been found. +If the `PDESystem` was solved using `NeuralPDE.jl` and `Optimization.jl`, then the argument `phi` is a field of the output of `NeuralPDE.discretize` and the argument `θ` is `res.u.depvar` where `res` is the result of the optimization. + ```@docs get_numerical_lyapunov_function ``` diff --git a/src/NeuralLyapunovPDESystem.jl b/src/NeuralLyapunovPDESystem.jl index 58306de..9a29bb6 100644 --- a/src/NeuralLyapunovPDESystem.jl +++ b/src/NeuralLyapunovPDESystem.jl @@ -32,9 +32,9 @@ Construct a `ModelingToolkit.PDESystem` representing the specified neural Lyapun here and cannot be inferred, `[:param1, :param2, ...]` will be used. - `policy_search::Bool`: whether or not to include a loss term enforcing `fixed_point` to actually be a fixed point; defaults to `false`; only used when `dynamics isa 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. + !(dynamics isa ODEFunction)`; when `dynamics isa ODEFunction`, `policy_search` should + not be supplied (as it must be false); when `dynamics isa ODESystem`, value inferred by + the presence of unbound inputs. - `name`: the name of the constructed `PDESystem` """ function NeuralLyapunovPDESystem( From 7353da3271be8636e30613432f503fbed1c19f8f Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 17 May 2024 15:52:22 -0400 Subject: [PATCH 03/34] Updated documentation on Lyapunov function structures --- docs/src/man/structure.md | 15 +++++- src/conditions_specification.jl | 30 ++++++----- src/structure_specification.jl | 90 +++++++++++++++++++++------------ 3 files changed, 86 insertions(+), 49 deletions(-) diff --git a/docs/src/man/structure.md b/docs/src/man/structure.md index b2d8c26..c476f4b 100644 --- a/docs/src/man/structure.md +++ b/docs/src/man/structure.md @@ -6,10 +6,21 @@ NeuralLyapunovStructure ## Pre-defined structures +Three simple neural Lyapunov function structures are provided, but users can always specify a custom structure using the [`NeuralLyapunovStructure`](@ref) struct. + +The simplest structure is to train the neural network directly to be the Lyapunov function, which can be accomplished using an [`UnstructuredNeuralLyapunov`](@ref). + ```@docs UnstructuredNeuralLyapunov +``` + +The condition that the Lyapunov function ``V(x)`` must be minimized uniquely at the fixed point ``x_0`` is often represented as a requirement for ``V(x)`` to be positive away from the fixed point and zero at the fixed point. +Put mathematically, it is sufficient to require ``V(x) > 0 \, \forall x \ne x_0`` and ``V(x_0) = 0``. +We call such functions positive definite (around the fixed point ``x_0``). + +Two structures are provided which partially or fully enforce the minimization condition: [`NonnegativeNeuralLyapunov`](@ref), which structurally enforces ``V(x) \ge 0`` everywhere, and [`PositiveSemiDefiniteStructure`](@ref), which additionally enforces ``V(x_0) = 0``. + +```@docs NonnegativeNeuralLyapunov PositiveSemiDefiniteStructure ``` - -## Defining your own neural Lyapunov structure diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl index 8460436..ec0a954 100644 --- a/src/conditions_specification.jl +++ b/src/conditions_specification.jl @@ -1,21 +1,23 @@ """ - NeuralLyapunovStructure + NeuralLyapunovStructure(V, ∇V, V̇, f_call, network_dim) Specifies the structure of the neural Lyapunov function and its derivative. -Allows the user to define the Lyapunov in terms of the neural network to structurally -enforce Lyapunov conditions. -`V(phi::Function, state, fixed_point)` takes in the neural network, the state, and the fixed -point, and outputs the value of the Lyapunov function at `state`. -`V̇(phi::Function, J_phi::Function, f::Function, state, params, t, fixed_point)` takes in the -neural network, jacobian of the neural network, dynamics, state, parameters and time (for -calling the dynamics, when relevant), and fixed point, and outputs the time derivative of -the Lyapunov function at `state`. -`f_call(dynamics::Function, phi::Function, state, p, t)` takes in the dynamics, the neural -network, the state, the parameters of the dynamics, and time, and outputs the derivative of -the state; this is useful for making closed-loop dynamics which depend on the neural -network, such as in the policy search case -`network_dim` is the dimension of the output of the neural network. +Allows the user to define the Lyapunov in terms of the neural network, potentially +structurally enforcing some Lyapunov conditions. + +# Fields + - `V(phi::Function, state, fixed_point)`: takes in the neural network, the state, and the + fixed point, and outputs the value of the Lyapunov function at `state`. + - `V̇(phi::Function, J_phi::Function, dynamics::Function, state, params, t, fixed_point)`: + takes in the neural network, jacobian of the neural network, dynamics, state, parameters + and time (for calling the dynamics, when relevant), and fixed point, and outputs the + time derivative of the Lyapunov function at `state`. + - `f_call(dynamics::Function, phi::Function, state, p, t)`: takes in the dynamics, the + neural network, the state, the parameters of the dynamics, and time, and outputs the + derivative of the state; this is useful for making closed-loop dynamics which depend on + the neural network, such as in the policy search case. + - `network_dim`: the dimension of the output of the neural network. """ struct NeuralLyapunovStructure V::Function diff --git a/src/structure_specification.jl b/src/structure_specification.jl index aac8f8e..0e0d6f5 100644 --- a/src/structure_specification.jl +++ b/src/structure_specification.jl @@ -4,8 +4,10 @@ Create a [`NeuralLyapunovStructure`](@ref) where the Lyapunov function is the neural network evaluated at the state. This does not structurally enforce any Lyapunov conditions. +Corresponds to ``V(x) = ϕ(x)``, where ``ϕ`` is the neural network. + Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For -`f(state, input, p, t)`, consider using `add_policy_search`. +`f(state, input, p, t)`, consider using [`add_policy_search`](@ref). """ function UnstructuredNeuralLyapunov()::NeuralLyapunovStructure NeuralLyapunovStructure( @@ -24,22 +26,31 @@ end Create a [`NeuralLyapunovStructure`](@ref) 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 -be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This -structure ensures `V(state) ≥ 0`. Further, if `δ > 0` and -`pos_def(fixed_point, fixed_point) = 0`, but `pos_def(state, fixed_point) > 0` when -`state ≠ fixed_point`, this ensures that `V(state) > 0` when `state != fixed_point`. This -does not enforce `V(fixed_point) = 0`, so that condition must included in the neural -Lyapunov loss function. - -`grad_pos_def(state, fixed_point)` should be the gradient of `pos_def` with respect to -`state` at `state`. If `grad_pos_def` is not defined, it is evaluated using `grad`, which -defaults to `ForwardDiff.gradient`. - -The neural network output has dimension `network_dim`. +Corresponds to ``V(x) = \\lVert ϕ(x) \\rVert^2 + δ \\, \\texttt{pos\\_def}(x, x_0)``, where +``ϕ`` is the neural network and `x_0` is the equilibrium point. + +This structure ensures ``V(x) ≥ 0 \\, ∀ x`` when ``δ ≥ 0`` and `pos_def` is always +nonnegative. Further, if ``δ > 0`` and `pos_def` is strictly positive definite around +`fixed_point`, the structure ensures that ``V(x)`` is strictly positive away from +`fixed_point`. In such cases, the minimization condition reduces to ensuring +``V(\\texttt{fixed\\_point}) = 0``, and so [`DontCheckNonnegativity(true)`](@ref) should be +used. + +# Arguments + - `network_dim`: output dimensionality of the neural network. + - `δ`: weight of `pos_def`, as above. + - `pos_def(state, fixed_point)`: a function that is postive (semi-)definite in `state` + around `fixed_point`. + +# Keyword Arguments + - `grad_pos_def(state, fixed_point)`: the gradient of `pos_def` with respect to `state` at + `state`. If `isnothing(grad_pos_def)` (as is the default), the gradient of `pos_def` + will be evaluated using `grad`. + - `grad`: a function for evaluating gradients to be used when `isnothing(grad_pos_def)`; + defaults to, and expects the same arguments as, `ForwardDiff.gradient`. Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For -`f(state, input, p, t)`, consider using `add_policy_search`. +`f(state, input, p, t)`, consider using [`add_policy_search`](@ref). See also: [`DontCheckNonnegativity`](@ref) """ @@ -89,27 +100,40 @@ end Create a [`NeuralLyapunovStructure`](@ref) 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. - -The condition that the Lyapunov function must be minimized uniquely at the fixed point can -be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This -structure ensures `V(state) ≥ 0`. Further, if `pos_def` is `0` only at `fixed_point` (and -positive elsewhere) and if `non_neg` is strictly positive away from `fixed_point` (as is the -case for the default values of `pos_def` and `non_neg`), then this structure ensures -`V(fixed_point) = 0` and `V(state) > 0` when `state ≠ fixed_point`. - -`grad_pos_def(state, fixed_point)` should be the gradient of `pos_def` with respect to -`state` at `state`. Similarly, `grad_non_neg(net, J_net, state, fixed_point)` should be the -gradient of `non_neg(net, state, fixed_point)` with respect to `state` at `state`. If -`grad_pos_def` or `grad_non_neg` is not defined, it is evaluated using `grad`, which -defaults to `ForwardDiff.gradient`. - -The neural network output has dimension `network_dim`. +nonnegative function `non_neg` which does depend the network. + +Corresponds to ``V(x) = \\texttt{pos\\_def}(x, x_0) * \\texttt{non\\_neg}(ϕ, x, x_0)``, where +``ϕ`` is the neural network and `x_0` is the equilibrium point. + +This structure ensures ``V(x) ≥ 0``. Further, if `pos_def` is strictly positive definite +`fixed_point` and `non_neg` is strictly positive (as is the case for the default values of +`pos_def` and `non_neg`), then this structure ensures ``V(x)`` is strictly positive definite +around `fixed_point`. In such cases, the minimization condition is satisfied structurally, +so [`DontCheckNonnegativity(false)`](@ref) should be used. + +# Arguments + - network_dim: output dimensionality of the neural network. + +# Keyword Arguments + - `pos_def(state, fixed_point)`: a function that is postive (semi-)definite in `state` + around `fixed_point`. + - `non_neg(net, state, fixed_point)`: a nonnegative function of the neural network; note + that `net` is the neural network ``ϕ``, and `net(state)` is the value of the neural + network at a point ``ϕ(x)``. + - `grad_pos_def(state, fixed_point)`: the gradient of `pos_def` with respect to `state` at + `state`. If `isnothing(grad_pos_def)` (as is the default), the gradient of `pos_def` + will be evaluated using `grad`. + - `grad_non_neg(net, J_net, state, fixed_point)`: the gradient of `non_neg` with respect + to `state` at `state`; `J_net` is a function outputting the Jacobian of `net` at the + input. If `isnothing(grad_non_neg)` (as is the default), the gradient of `non_neg` will + be evaluated using `grad`. + - `grad`: a function for evaluating gradients to be used when `isnothing(grad_pos_def) || + isnothing(grad_non_neg)`; defaults to, and expects the same arguments as, `ForwardDiff.gradient`. Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For -`f(state, input, p, t)`, consider using `add_policy_search`. +`f(state, input, p, t)`, consider using [`add_policy_search`](@ref). -Typically used with [`DontCheckNonnegativity`](@ref). +See also: [`DontCheckNonnegativity`](@ref) """ function PositiveSemiDefiniteStructure( network_dim::Integer; From e18ec5763bcbdd15ee10d46595134ef488cdccda Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 17 May 2024 18:18:19 -0400 Subject: [PATCH 04/34] Eliminated the grad V field from `NeuralLyapunovStructure` --- src/conditions_specification.jl | 20 +++++++++----------- src/policy_search.jl | 5 +---- src/structure_specification.jl | 8 -------- 3 files changed, 10 insertions(+), 23 deletions(-) diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl index ec0a954..d47c9be 100644 --- a/src/conditions_specification.jl +++ b/src/conditions_specification.jl @@ -1,5 +1,5 @@ """ - NeuralLyapunovStructure(V, ∇V, V̇, f_call, network_dim) + NeuralLyapunovStructure(V, V̇, f_call, network_dim) Specifies the structure of the neural Lyapunov function and its derivative. @@ -7,21 +7,19 @@ Allows the user to define the Lyapunov in terms of the neural network, potential structurally enforcing some Lyapunov conditions. # Fields - - `V(phi::Function, state, fixed_point)`: takes in the neural network, the state, and the - fixed point, and outputs the value of the Lyapunov function at `state`. + - `V(phi::Function, state, fixed_point)`: outputs the value of the Lyapunov function at + `state`. - `V̇(phi::Function, J_phi::Function, dynamics::Function, state, params, t, fixed_point)`: - takes in the neural network, jacobian of the neural network, dynamics, state, parameters - and time (for calling the dynamics, when relevant), and fixed point, and outputs the - time derivative of the Lyapunov function at `state`. - - `f_call(dynamics::Function, phi::Function, state, p, t)`: takes in the dynamics, the - neural network, the state, the parameters of the dynamics, and time, and outputs the - derivative of the state; this is useful for making closed-loop dynamics which depend on - the neural network, such as in the policy search case. + outputs the time derivative of the Lyapunov function at `state`. + - `f_call(dynamics::Function, phi::Function, state, p, t)`: outputs the derivative of the + state; this is useful for making closed-loop dynamics which depend on the neural + network, such as in the policy search case. - `network_dim`: the dimension of the output of the neural network. + +`phi` and `J_phi` above are both functions of `state` alone. """ struct NeuralLyapunovStructure V::Function - ∇V::Function V̇::Function f_call::Function network_dim::Integer diff --git a/src/policy_search.jl b/src/policy_search.jl index eb099b4..b4bb214 100644 --- a/src/policy_search.jl +++ b/src/policy_search.jl @@ -19,7 +19,7 @@ function add_policy_search( new_dims::Integer; control_structure::Function = identity )::NeuralLyapunovStructure - let V = lyapunov_structure.V, ∇V = lyapunov_structure.∇V, V̇ = lyapunov_structure.V̇, + let V = lyapunov_structure.V, V̇ = lyapunov_structure.V̇, V_dim = lyapunov_structure.network_dim, nd = new_dims, u = control_structure NeuralLyapunovStructure( @@ -34,9 +34,6 @@ function add_policy_search( V(st -> net(st)[1:V_dim, :], state, fixed_point) end end, - function (net, J_net, state, fixed_point) - ∇V(st -> net(st)[1:V_dim], st -> J_net(st)[1:V_dim, :], state, fixed_point) - end, function (net, J_net, f, state, params, t, fixed_point) V̇(st -> net(st)[1:V_dim], st -> J_net(st)[1:V_dim, :], (st, p, t) -> f(st, u(net(st)[(V_dim + 1):end]), p, t), state, params, diff --git a/src/structure_specification.jl b/src/structure_specification.jl index 0e0d6f5..042e2c3 100644 --- a/src/structure_specification.jl +++ b/src/structure_specification.jl @@ -12,7 +12,6 @@ Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For function UnstructuredNeuralLyapunov()::NeuralLyapunovStructure NeuralLyapunovStructure( (net, state, fixed_point) -> net(state), - (net, grad_net, state, fixed_point) -> grad_net(state), (net, grad_net, f, state, params, t, fixed_point) -> grad_net(state) ⋅ f(state, params, t), (f, net, state, p, t) -> f(state, p, t), @@ -66,7 +65,6 @@ function NonnegativeNeuralLyapunov( if δ == 0.0 NeuralLyapunovStructure( (net, state, fixed_point) -> net(state) ⋅ net(state), - (net, J_net, state, fixed_point) -> 2 * transpose(net(state)) * J_net(state), (net, J_net, f, state, params, t, fixed_point) -> 2 * dot( net(state), J_net(state), f(state, params, t)), @@ -82,8 +80,6 @@ function NonnegativeNeuralLyapunov( NeuralLyapunovStructure( (net, state, fixed_point) -> net(state) ⋅ net(state) + δ * pos_def(state, fixed_point), - (net, J_net, state, fixed_point) -> 2 * transpose(net(state)) * J_net(state) + - δ * grad_pos_def(state, fixed_point), (net, J_net, f, state, params, t, fixed_point) -> 2 * dot( net(state), J_net(state), @@ -161,10 +157,6 @@ function PositiveSemiDefiniteStructure( NeuralLyapunovStructure( (net, state, fixed_point) -> pos_def(state, fixed_point) * non_neg(net, state, fixed_point), - (net, J_net, state, fixed_point) -> grad_pos_def(state, fixed_point) * - non_neg(net, state, fixed_point) + - pos_def(state, fixed_point) * - grad_non_neg(net, J_net, state, fixed_point), (net, J_net, f, state, params, t, fixed_point) -> (f(state, params, t) ⋅ grad_pos_def( state, fixed_point)) * From 67f0520d875bcf01daea1a0603422378aa03b008 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 17 May 2024 18:20:14 -0400 Subject: [PATCH 05/34] Added documentation for defining custom Lyapunov function structures --- docs/src/man/structure.md | 39 +++++++++++++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 4 deletions(-) diff --git a/docs/src/man/structure.md b/docs/src/man/structure.md index c476f4b..85cccfa 100644 --- a/docs/src/man/structure.md +++ b/docs/src/man/structure.md @@ -1,9 +1,5 @@ # Structuring a neural Lyapunov function -```@docs -NeuralLyapunovStructure -``` - ## Pre-defined structures Three simple neural Lyapunov function structures are provided, but users can always specify a custom structure using the [`NeuralLyapunovStructure`](@ref) struct. @@ -24,3 +20,38 @@ Two structures are provided which partially or fully enforce the minimization co NonnegativeNeuralLyapunov PositiveSemiDefiniteStructure ``` + +## Defining your own neural Lyapunov function structure + +To define a new structure for a neural Lyapunov function, one must specify the form of the Lyapunov candidate ``V`` and its time derivative along a trajectory ``\dot{V}``, as well as how to call the dynamics ``f``. +Additionally, the dimensionality of the output of the neural network must be known in advance. + +```@docs +NeuralLyapunovStructure +``` + +### Calling the dynamics + +Very generally, the dynamical system can be a system of ODEs ``\dot{x} = f(x, u, p, t)``, where ``u`` is a control input, ``p`` contains parameters, and ``f`` depends on the neural network in some way. +To capture this variety, users must supply the function `f_call(dynamics::Function, phi::Function, state, p, t)`. + +The most common example is when `dynamics` takes the same form as an `ODEFunction`. +i.e., ``\dot{x} = \texttt{dynamics}(x, p, t)``. +In that case, `f_call(dynamics, phi, state, p, t) = dynamics(state, p, t)`. + +Suppose instead, the dynamics were supplied as a function of state alone: ``\dot{x} = \texttt{dynamics}(x)``. +Then, `f_call(dynamics, phi, state, p, t) = dynamics(state)`. + +Finally, suppose ``\dot{x} = \texttt{dynamics}(x, u, p, t)`` has a unidimensional control input that is being trained (via [policy search](policy_search.md)) to be the second output of the neural network. +Then `f_call(dynamics, phi, state, p, t) = dynamics(state, phi(state)[2], p, t)`. + +Note that, despite the inclusion of the time variable ``t``, NeuralLyapunov.jl currently only supports time-invariant systems, so only `t = 0.0` is used. + +### Structuring ``V`` and ``\dot{V}`` + +The Lyapunov candidate function ``V`` gets specified as a function `V(phi, state, fixed_point)`, where `phi` is the neural network as a function `phi(state)`. +Note that this form allows ``V(x)`` to depend on the neural network evaluated at points other than just the input ``x``. + +The time derivative ``\dot{V}`` is similarly defined by a function `V̇(phi, J_phi, dynamics, state, params, t, fixed_point)`. +The function `J_phi(state)` gives the Jacobian of the neural network `phi` at `state`. +The function `dynamics` is as above (with parameters `params`). From 256deb9b963f00bebece4a164afcb72fc91c53f5 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Tue, 21 May 2024 12:58:24 -0400 Subject: [PATCH 06/34] Updated `LyapunovMinimizationCondition` documentation --- docs/make.jl | 2 +- docs/src/man/minimization.md | 7 +++++++ src/minimization_conditions.jl | 33 ++++++++++++++++++++------------- 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 4de1343..5fc7743 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -18,9 +18,9 @@ makedocs(; "Manual" => [ "man.md", "man/pdesystem.md", - "man/structure.md", "man/minimization.md", "man/decrease.md", + "man/structure.md", "man/roa.md", "man/policy_search.md", "man/local_lyapunov.md", diff --git a/docs/src/man/minimization.md b/docs/src/man/minimization.md index a2b6433..0a5e031 100644 --- a/docs/src/man/minimization.md +++ b/docs/src/man/minimization.md @@ -1,5 +1,12 @@ # Lyapunov Minimization Condition +The condition that the Lyapunov function ``V(x)`` must be minimized uniquely at the fixed point ``x_0`` is often represented as a requirement for ``V(x)`` to be positive away from the fixed point and zero at the fixed point. +Put mathematically, it is sufficient to require ``V(x) > 0 \, \forall x \ne x_0`` and ``V(x_0) = 0``. +We call such functions positive definite (around the fixed point ``x_0``). +The weaker condition that ``V(x) \ge 0 \, \forall x \ne x_0`` and ``V(x_0) = 0`` is positive *semi-*definiteness. + +Users specify the form of the minimization condition to enforce through training using a [`LyapunovMinimizationCondition`](@ref). + ```@docs LyapunovMinimizationCondition ``` diff --git a/src/minimization_conditions.jl b/src/minimization_conditions.jl index ad67436..c8afd60 100644 --- a/src/minimization_conditions.jl +++ b/src/minimization_conditions.jl @@ -1,28 +1,35 @@ """ - LyapunovMinimizationCondition + LyapunovMinimizationCondition(check_nonnegativity, strength, rectifier, check_fixed_point) -Specifies the form of the Lyapunov conditions to be used. +Specifies the form of the Lyapunov minimization condition to be used. + +# Fields + - `check_nonnegativity::Bool`: whether or not to train for positivity/non-negativity of ``V(x)`` + - `strength::Function`: level of strictness for positivity training, if `check_nonnegativity == true` + - `rectifier::Function`: positive when the input is positive and (approximately) zero when the input is negative + - `check_fixed_point`: whether or not to train for ``V(x_0) = 0``. + +# Training conditions If `check_nonnegativity` is `true`, training will attempt to enforce - `V(state) ≥ strength(state, fixed_point)`. + ``V(x) ≥ \\texttt{strength}(x, x_0)``. The inequality will be approximated by the equation - `rectifier(strength(state, fixed_point) - V(state)) = 0.0`. + ``\\texttt{rectifier}(\\texttt{strength}(x, x_0) - V(x_0)) = 0``. + If `check_fixed_point` is `true`, then training will also attempt to enforce - `V(fixed_point) = 0`. + ``V(x_0) = 0``. # Examples -The condition that the Lyapunov function must be minimized uniquely at the fixed point can -be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This -could be enfored by `V(fixed_point) ≥ ||state - fixed_point||^2`, which would be -represented, with `check_nonnegativity = true`, by - strength(state, fixed_point) = ||state - fixed_point||^2, -paired with `V(fixed_point) = 0`, which can be enforced with `check_fixed_point = true`. +When training for a strictly positive definite ``V``, an example of an appropriate `strength` +is ``\\texttt{strength}(x, x_0) = \\lVert x - x_0 \\rVert^2``. +This form is used in [`StrictlyPositiveDefinite`](@ref). -If `V` were structured such that it is always nonnegative, then `V(fixed_point) = 0` is all +If ``V`` were structured such that it is always nonnegative, then ``V(x_0) = 0`` is all that must be enforced in training for the Lyapunov function to be uniquely minimized at -`fixed_point`. So, in that case, we would use +``x_0``. So, in that case, we would use `check_nonnegativity = false; check_fixed_point = true`. +This can also be accomplished with [`DontCheckNonnegativity(true)`](@ref). In either case, `rectifier = (t) -> max(0.0, t)` exactly represents the inequality, but differentiable approximations of this function may be employed. From d70582c82c9cbd5c2bf51ece264dbf807c8a9c6a Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Tue, 21 May 2024 15:00:36 -0400 Subject: [PATCH 07/34] Added documentation of minimization conditions --- docs/src/man/minimization.md | 11 +++++-- src/conditions_specification.jl | 8 +++-- src/minimization_conditions.jl | 54 +++++++++++++++++++-------------- 3 files changed, 46 insertions(+), 27 deletions(-) diff --git a/docs/src/man/minimization.md b/docs/src/man/minimization.md index 0a5e031..04de57d 100644 --- a/docs/src/man/minimization.md +++ b/docs/src/man/minimization.md @@ -5,7 +5,7 @@ Put mathematically, it is sufficient to require ``V(x) > 0 \, \forall x \ne x_0` We call such functions positive definite (around the fixed point ``x_0``). The weaker condition that ``V(x) \ge 0 \, \forall x \ne x_0`` and ``V(x_0) = 0`` is positive *semi-*definiteness. -Users specify the form of the minimization condition to enforce through training using a [`LyapunovMinimizationCondition`](@ref). +NeuralLyapunov.jl provides the [`LyapunovMinimizationCondition`](@ref) struct for users to specify the form of the minimization condition to enforce through training. ```@docs LyapunovMinimizationCondition @@ -14,17 +14,24 @@ LyapunovMinimizationCondition ## Pre-defined minimization conditions ```@docs +StrictlyPositiveDefinite PositiveSemiDefinite DontCheckNonnegativity -StrictlyPositiveDefinite ``` ## Defining your own minimization condition +If a user wishes to define their own version of the minimization condition in a form other than +``V(x) ≥ \texttt{strength}(x, x_0)``, +they must define their own subtype of [`NeuralLyapunov.AbstractLyapunovMinimizationCondition`](@ref). + ```@docs NeuralLyapunov.AbstractLyapunovMinimizationCondition ``` +When constructing the PDESystem, [`NeuralLyapunovPDESystem`](@ref) uses [`NeuralLyapunov.check_nonnegativity`](@ref) to determine if it should include an equation equating the result of [`NeuralLyapunov.get_minimization_condition`](@ref) to zero. +It additionally uses [`NeuralLyapunov.check_minimal_fixed_point`](@ref) to determine if it should include the equation ``V(x_0) = 0``. + ```@docs NeuralLyapunov.check_nonnegativity NeuralLyapunov.check_minimal_fixed_point diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl index d47c9be..5481212 100644 --- a/src/conditions_specification.jl +++ b/src/conditions_specification.jl @@ -89,8 +89,12 @@ end """ get_minimization_condition(cond::AbstractLyapunovMinimizationCondition) -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. +Return a function of ``V``, ``x``, and ``x_0`` that equals zero when the Lyapunov +minimization condition is met for the Lyapunov candidate function ``V`` at the point ``x``, +and is greater than zero if it's violated. + +Note that the first input, ``V``, is a function, so the minimization condition can depend on +the value of the candidate Lyapunov function at multiple points. """ function get_minimization_condition(cond::AbstractLyapunovMinimizationCondition) error("get_condition not implemented for AbstractLyapunovMinimizationCondition of " * diff --git a/src/minimization_conditions.jl b/src/minimization_conditions.jl index c8afd60..2575ed9 100644 --- a/src/minimization_conditions.jl +++ b/src/minimization_conditions.jl @@ -4,9 +4,13 @@ Specifies the form of the Lyapunov minimization condition to be used. # Fields - - `check_nonnegativity::Bool`: whether or not to train for positivity/non-negativity of ``V(x)`` - - `strength::Function`: level of strictness for positivity training, if `check_nonnegativity == true` - - `rectifier::Function`: positive when the input is positive and (approximately) zero when the input is negative + - `check_nonnegativity::Bool`: whether or not to train for positivity/nonnegativity of + ``V(x)`` + - `strength::Function`: specifies the level of strictness for positivity training; should + be zero when the two inputs are equal and nonnegative otherwise; used when + `check_nonnegativity == true` + - `rectifier::Function`: positive when the input is positive and (approximately) zero when + the input is negative - `check_fixed_point`: whether or not to train for ``V(x_0) = 0``. # Training conditions @@ -15,6 +19,8 @@ If `check_nonnegativity` is `true`, training will attempt to enforce ``V(x) ≥ \\texttt{strength}(x, x_0)``. The inequality will be approximated by the equation ``\\texttt{rectifier}(\\texttt{strength}(x, x_0) - V(x_0)) = 0``. +Note that the approximate equation and inequality are identical when +``\\texttt{rectifier}(t) = \\max(0, t)``. If `check_fixed_point` is `true`, then training will also attempt to enforce ``V(x_0) = 0``. @@ -27,12 +33,13 @@ This form is used in [`StrictlyPositiveDefinite`](@ref). If ``V`` were structured such that it is always nonnegative, then ``V(x_0) = 0`` is all that must be enforced in training for the Lyapunov function to be uniquely minimized at -``x_0``. So, in that case, we would use +``x_0``. In that case, we would use `check_nonnegativity = false; check_fixed_point = true`. This can also be accomplished with [`DontCheckNonnegativity(true)`](@ref). -In either case, `rectifier = (t) -> max(0.0, t)` exactly represents the inequality, but -differentiable approximations of this function may be employed. +In either case, the rectified linear unit `rectifier = (t) -> max(zero(t), t)` exactly +represents the inequality, but differentiable approximations of this function may be +employed. """ struct LyapunovMinimizationCondition <: AbstractLyapunovMinimizationCondition check_nonnegativity::Bool @@ -60,12 +67,13 @@ end """ StrictlyPositiveDefinite(C; check_fixed_point, rectifier) -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`. +Construct a [`LyapunovMinimizationCondition`](@ref) representing + ``V(x) ≥ C \\lVert x - x_0 \\rVert^2``. +If `check_fixed_point == true`, then training will also attempt to enforce + ``V(x_0) = 0``. -The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. +The inequality is represented by + ``\\texttt{rectifier}(C \\lVert x - x_0 \\rVert^2 - V(x)) = 0``. """ function StrictlyPositiveDefinite(; check_fixed_point = true, @@ -83,13 +91,13 @@ end """ PositiveSemiDefinite(check_fixed_point) -Construct a `LyapunovMinimizationCondition` representing - `V(state) ≥ 0`. -If `check_fixed_point` is `true`, then training will also attempt to enforce - `V(fixed_point) = 0`. +Construct a [`LyapunovMinimizationCondition`](@ref) representing + ``V(x) ≥ 0``. +If `check_fixed_point == true`, then training will also attempt to enforce + ``V(x_0) = 0``. -The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. -""" +The inequality is represented by + ``\\texttt{rectifier}( -V(x) ) = 0``.""" function PositiveSemiDefinite(; check_fixed_point = true, rectifier = (t) -> max(zero(t), t) @@ -105,13 +113,13 @@ end """ DontCheckNonnegativity(check_fixed_point) -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. +Construct a [`LyapunovMinimizationCondition`](@ref) which represents not checking for +nonnegativity of the Lyapunov function. This is appropriate in cases where this condition +has been structurally enforced. -It is still possible to check for `V(fixed_point) = 0`, even in this case, for example if -`V` is structured to be positive for `state ≠ fixed_point`, but it is not guaranteed -structurally that `V(fixed_point) = 0`. +It is still possible to check for ``V(x_0) = 0``, even in this case, for example if `V` is +structured to be positive for ``x ≠ x_0``, but it is not guaranteed structurally that +``V(x_0) = 0`` (such as [`NonnegativeNeuralLyapunov`](@ref)). """ function DontCheckNonnegativity(; check_fixed_point = false)::LyapunovMinimizationCondition LyapunovMinimizationCondition( From 1e4eba9d8b419aeeb3f79381d66386d3d9356454 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 24 May 2024 15:48:57 -0400 Subject: [PATCH 08/34] Refactored `LyapunovDecreaseCondition` and changed meaning of `strength` to better match `LyapunovMinimizationCondition`; updated decrease condition documentation --- docs/src/man.md | 2 +- docs/src/man/decrease.md | 12 ++++ src/decrease_conditions.jl | 90 +++++++++++++++++++--------- src/decrease_conditions_RoA_aware.jl | 18 +++--- src/minimization_conditions.jl | 14 +++-- 5 files changed, 94 insertions(+), 42 deletions(-) diff --git a/docs/src/man.md b/docs/src/man.md index 8fb53c8..9fd7a2e 100644 --- a/docs/src/man.md +++ b/docs/src/man.md @@ -1,6 +1,6 @@ # Components of a Neural Lyapunov Problem -For a candidate Lyapunov function $V(x)$ to certify the stability of an equilibrium point $x_0$ of the dynamical system $\dot{x}(t) = f(x(t))$, it must satisfy two conditions: +For a candidate Lyapunov function $V(x)$ to certify the stability of an equilibrium point $x_0$ of the dynamical system $\frac{dx}{dt} = f(x(t))$, it must satisfy two conditions: 1. The function $V$ must be uniquely minimized at $x_0$, and 2. The function $V$ must decrease along system trajectories (i.e., $V(x(t))$ decreases as long as $x(t)$ is a trajectory of the dynamical system). diff --git a/docs/src/man/decrease.md b/docs/src/man/decrease.md index c6572f6..727988d 100644 --- a/docs/src/man/decrease.md +++ b/docs/src/man/decrease.md @@ -1,5 +1,17 @@ # Lyapunov Decrease Condition +To represent the condition that the Lyapunov function ``V(x)`` must decrease along system trajectories, we typically introduce a new function ``\dot{V}(x) = \nabla_{f(x)} V(x) = \nabla_x V(x) \cdot f(x)``. +It is then sufficient to show that ``\dot{V}(x)`` is negative away from the fixed point and zero at the fixed point, since ``\dot{V}`` represents the rate of change of ``V`` along system trajectories. +i.e., if ``x(t)`` is a trajectory defined by ``\frac{dx}{dt} = f(x)``, then ``\dot{V}(x(t)) = \frac{d}{dt} [ V(x(t)) ]``. + +Put mathematically, it is sufficient to require ``\dot{V}(x) < 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0``. +We call such functions negative definite (around the fixed point ``x_0``). +The weaker condition that ``\dot{V}(x) \le 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0`` is negative *semi-*definiteness. + +The condition that ``\dot{V}(x_0) = 0`` is satisfied by the definition of ``\dot{V}`` and the fact that ``x_0`` is a fixed point, so we do not need to train for it unless the dynamics have some dependence on the neural network (e.g., in [policy search](policy_search.md)). + +NeuralLyapunov.jl provides the [`LyapunovDecreaseCondition`](@ref) struct for users to specify the form of the decrease condition to enforce through training. + ```@docs LyapunovDecreaseCondition ``` diff --git a/src/decrease_conditions.jl b/src/decrease_conditions.jl index 684905e..bce2e15 100644 --- a/src/decrease_conditions.jl +++ b/src/decrease_conditions.jl @@ -1,31 +1,62 @@ """ - LyapunovDecreaseCondition(check_decrease, decrease, strength, rectifier) + LyapunovDecreaseCondition(check_decrease, rate_metric, 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)`. +Specifies the form of the Lyapunov conditions to be used. -The inequality will be approximated by the equation - `rectifier(decrease(V, dVdt) - strength(state, fixed_point)) = 0.0`. +# Fields + - `check_decrease::Bool`: whether or not to train for negativity/nonpositivity of + ``V̇(x)``. + - `rate_metric::Function`: should increase with ``V̇(x)``; used when + `check_decrease == true`. + - `strength::Function`: specifies the level of strictness for negativity training; should + be zero when the two inputs are equal and nonnegative otherwise; used when + `check_decrease == true`. + - `rectifier::Function`: positive when the input is positive and (approximately) zero when + the input is negative. -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 -to enforce `dVdt(fixed_point) = 0`, so `check_fixed_point` defaults to `false`. +If `check_decrease == true`, training will enforce: + +``\\texttt{rate\\_metric}(V(x), V̇(x)) ≤ -\\texttt{strength}(x, x_0).`` + +The inequality will be approximated by the equation: + +``\\texttt{rectifier}(\\texttt{rate\\_metric}(V(x), V̇(x)) + \\texttt{strength}(x, x_0)) = 0.`` + +Note that the approximate equation and inequality are identical when +``\\texttt{rectifier}(t) = \\max(0, t)``. + +If the dynamics truly have a fixed point at ``x_0`` and ``V̇(x)`` is truly the rate of +decrease of ``V(x)`` along the dynamics, then ``V̇(x_0)`` will be ``0`` and there is no need +to train for ``V̇(x_0) = 0``. # Examples: Asymptotic decrease can be enforced by requiring - `dVdt ≤ -C |state - fixed_point|^2`, + ``V̇(x) ≤ -C \\lVert x - x_0 \\rVert^2``, +for some positive ``C``, which corresponds to + + rate_metric = (V, dVdt) -> dVdt + strength = (x, x0) -> C * (x - x0) ⋅ (x - x0) + +This can also be accomplished with [`AsymptoticDecrease`](@ref). + +Exponential decrease of rate ``k`` is proven by + ``V̇(x) ≤ - k * V(x)``, which corresponds to - `decrease = (V, dVdt) -> dVdt` - `strength = (x, x0) -> -C * (x - x0) ⋅ (x - x0)` -Exponential decrease of rate `k` is proven by `dVdt ≤ - k * V`, so corresponds to - `decrease = (V, dVdt) -> dVdt + k * V` - `strength = (x, x0) -> 0.0` + rate_metric = (V, dVdt) -> dVdt + k * V + strength = (x, x0) -> 0.0 + +This can also be accomplished with [`ExponentialDecrease`](@ref). + + +In either case, the rectified linear unit `rectifier = (t) -> max(zero(t), t)` exactly +represents the inequality, but differentiable approximations of this function may be +employed. """ struct LyapunovDecreaseCondition <: AbstractLyapunovDecreaseCondition check_decrease::Bool - decrease::Function + rate_metric::Function strength::Function rectifier::Function end @@ -37,7 +68,7 @@ end function get_decrease_condition(cond::LyapunovDecreaseCondition) if cond.check_decrease return (V, dVdt, x, fixed_point) -> cond.rectifier( - cond.decrease(V(x), dVdt(x)) - cond.strength(x, fixed_point) + cond.rate_metric(V(x), dVdt(x)) + cond.strength(x, fixed_point) ) else return nothing @@ -47,12 +78,15 @@ end """ AsymptoticDecrease(; strict, C, rectifier) -Construct a `LyapunovDecreaseCondition` corresponding to asymptotic decrease. +Construct a [`LyapunovDecreaseCondition`](@ref) 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`. +If `strict == false`, the decrease condition is +``\\dot{V}(x) ≤ 0``, +and if `strict == true`, the condition is +``\\dot{V}(x) ≤ - C \\lVert x - x_0 \\rVert^2``. -The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. +The inequality is represented by +``\\texttt{rectifier}(\\dot{V}(x) + C \\lVert x - x_0 \\rVert^2) = 0``. """ function AsymptoticDecrease(; strict::Bool = false, @@ -60,7 +94,7 @@ function AsymptoticDecrease(; rectifier = (t) -> max(zero(t), t) )::LyapunovDecreaseCondition strength = if strict - (x, x0) -> -C * (x - x0) ⋅ (x - x0) + (x, x0) -> C * (x - x0) ⋅ (x - x0) else (x, x0) -> 0.0 end @@ -76,12 +110,14 @@ end """ ExponentialDecrease(k; strict, C, rectifier) -Construct a `LyapunovDecreaseCondition` corresponding to exponential decrease of rate `k`. +Construct a [`LyapunovDecreaseCondition`](@ref) 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`. +If `strict == false`, the condition is ``\\dot{V}(x) ≤ -k * V(x)``, and if `strict == true`, +the condition is ``\\dot{V}(x) ≤ -k * V(x) - C * \\lVert x - x_0 \\rVert^2``. -The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`. +The inequality is represented by +``\\texttt{rectifier}(\\dot{V}(x) + k V(x) + C \\lVert x - x_0 \\rVert^2) = 0``. """ function ExponentialDecrease( k::Real; @@ -90,7 +126,7 @@ function ExponentialDecrease( rectifier = (t) -> max(zero(t), t) )::LyapunovDecreaseCondition strength = if strict - (x, x0) -> -C * (x - x0) ⋅ (x - x0) + (x, x0) -> C * (x - x0) ⋅ (x - x0) else (x, x0) -> 0.0 end @@ -106,7 +142,7 @@ end """ DontCheckDecrease() -Construct a `LyapunovDecreaseCondition` which represents not checking for +Construct a [`LyapunovDecreaseCondition`](@ref) 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. """ diff --git a/src/decrease_conditions_RoA_aware.jl b/src/decrease_conditions_RoA_aware.jl index c8ae130..7120a46 100644 --- a/src/decrease_conditions_RoA_aware.jl +++ b/src/decrease_conditions_RoA_aware.jl @@ -1,17 +1,17 @@ """ - RoAAwareDecreaseCondition(check_decrease, decrease, strength, rectifier, ρ, + RoAAwareDecreaseCondition(check_decrease, rate_metric, strength, rectifier, ρ, out_of_RoA_penalty) Specifies the form of the Lyapunov conditions to be used, training for a region of attraction estimate of `{ x : V(x) ≤ ρ }` If `check_decrease`, training will enforce -`decrease(V(state), dVdt(state)) ≤ strength(state, fixed_point)` whenever `V(state) ≤ ρ`, +`rate_metric(V(state), dVdt(state)) ≤ -strength(state, fixed_point)` whenever `V(state) ≤ ρ`, 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 - `rectifier(decrease(V(state), dVdt(state)) - strength(state, fixed_point)) = 0.0`. + `rectifier(rate_metric(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 @@ -22,11 +22,11 @@ to enforce `dVdt(fixed_point) = 0`, so `check_fixed_point` defaults to `false`. Asymptotic decrease can be enforced by requiring `dVdt ≤ -C |state - fixed_point|^2`, which corresponds to - `decrease = (V, dVdt) -> dVdt` and - `strength = (x, x0) -> -C * (x - x0) ⋅ (x - x0)`. + `rate_metric = (V, dVdt) -> dVdt` and + `strength = (x, x0) -> C * (x - x0) ⋅ (x - x0)`. Exponential decrease of rate `k` is proven by `dVdt ≤ - k * V`, so corresponds to - `decrease = (V, dVdt) -> dVdt + k * V` and + `rate_metric = (V, dVdt) -> dVdt + k * V` and `strength = (x, x0) -> 0.0`. Enforcing either condition only in the region of attraction and not penalizing any points @@ -38,7 +38,7 @@ Note that this penalty could also depend on values of `V`, `dVdt`, and `ρ`. """ struct RoAAwareDecreaseCondition <: AbstractLyapunovDecreaseCondition check_decrease::Bool - decrease::Function + rate_metric::Function strength::Function rectifier::Function ρ::Real @@ -53,7 +53,7 @@ function get_decrease_condition(cond::RoAAwareDecreaseCondition) if cond.check_decrease return function (V, dVdt, x, fixed_point) (V(x) ≤ cond.ρ) * cond.rectifier( - cond.decrease(V(x), dVdt(x)) - cond.strength(x, fixed_point) + cond.rate_metric(V(x), dVdt(x)) + cond.strength(x, fixed_point) ) + (V(x) > cond.ρ) * cond.out_of_RoA_penalty(V(x), dVdt(x), x, fixed_point, cond.ρ) @@ -81,7 +81,7 @@ function make_RoA_aware( )::RoAAwareDecreaseCondition RoAAwareDecreaseCondition( cond.check_decrease, - cond.decrease, + cond.rate_metric, cond.strength, cond.rectifier, ρ, diff --git a/src/minimization_conditions.jl b/src/minimization_conditions.jl index 2575ed9..d7fd104 100644 --- a/src/minimization_conditions.jl +++ b/src/minimization_conditions.jl @@ -15,15 +15,19 @@ Specifies the form of the Lyapunov minimization condition to be used. # Training conditions -If `check_nonnegativity` is `true`, training will attempt to enforce - ``V(x) ≥ \\texttt{strength}(x, x_0)``. -The inequality will be approximated by the equation - ``\\texttt{rectifier}(\\texttt{strength}(x, x_0) - V(x_0)) = 0``. +If `check_nonnegativity` is `true`, training will attempt to enforce: + +``V(x) ≥ \\texttt{strength}(x, x_0).`` + +The inequality will be approximated by the equation: + +``\\texttt{rectifier}(\\texttt{strength}(x, x_0) - V(x_0)) = 0.`` + Note that the approximate equation and inequality are identical when ``\\texttt{rectifier}(t) = \\max(0, t)``. If `check_fixed_point` is `true`, then training will also attempt to enforce - ``V(x_0) = 0``. +``V(x_0) = 0``. # Examples From 4e955184c6b9f938cf0bf4ebdb70348959b3ae86 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 24 May 2024 16:38:57 -0400 Subject: [PATCH 09/34] Updated decrease documentation --- docs/src/man.md | 2 +- docs/src/man/decrease.md | 23 +++++++++++++++++------ docs/src/man/minimization.md | 18 +++++++++++------- src/conditions_specification.jl | 7 +++++-- src/decrease_conditions.jl | 4 ++-- 5 files changed, 36 insertions(+), 18 deletions(-) diff --git a/docs/src/man.md b/docs/src/man.md index 9fd7a2e..29de61b 100644 --- a/docs/src/man.md +++ b/docs/src/man.md @@ -11,7 +11,7 @@ Thus, we specify our neural Lyapunov problems with three components, each answer 2. How is the minimization condition to be enforced? 3. How is the decrease condition to be enforced? -These three components are represented by the three fields of a `NeuralLyapunovSpecification` object. +These three components are represented by the three fields of a [`NeuralLyapunovSpecification`](@ref) object. ```@docs NeuralLyapunovSpecification diff --git a/docs/src/man/decrease.md b/docs/src/man/decrease.md index 727988d..5294425 100644 --- a/docs/src/man/decrease.md +++ b/docs/src/man/decrease.md @@ -1,14 +1,15 @@ # Lyapunov Decrease Condition -To represent the condition that the Lyapunov function ``V(x)`` must decrease along system trajectories, we typically introduce a new function ``\dot{V}(x) = \nabla_{f(x)} V(x) = \nabla_x V(x) \cdot f(x)``. +To represent the condition that the Lyapunov function ``V(x)`` must decrease along system trajectories, we typically introduce a new function ``\dot{V}(x) = \nabla V(x) \cdot f(x)``. It is then sufficient to show that ``\dot{V}(x)`` is negative away from the fixed point and zero at the fixed point, since ``\dot{V}`` represents the rate of change of ``V`` along system trajectories. -i.e., if ``x(t)`` is a trajectory defined by ``\frac{dx}{dt} = f(x)``, then ``\dot{V}(x(t)) = \frac{d}{dt} [ V(x(t)) ]``. +That is to say, if ``x(t)`` is a trajectory defined by ``\frac{dx}{dt} = f(x)``, then ``\dot{V}(x(t)) = \frac{d}{dt} [ V(x(t)) ]``. Put mathematically, it is sufficient to require ``\dot{V}(x) < 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0``. We call such functions negative definite (around the fixed point ``x_0``). The weaker condition that ``\dot{V}(x) \le 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0`` is negative *semi-*definiteness. -The condition that ``\dot{V}(x_0) = 0`` is satisfied by the definition of ``\dot{V}`` and the fact that ``x_0`` is a fixed point, so we do not need to train for it unless the dynamics have some dependence on the neural network (e.g., in [policy search](policy_search.md)). +The condition that ``\dot{V}(x_0) = 0`` is satisfied by the definition of ``\dot{V}`` and the fact that ``x_0`` is a fixed point, so we do not need to train for it. +In cases where the dynamics have some dependence on the neural network (e.g., in [policy search](policy_search.md)), we should rather train directly for ``f(x_0) = 0``, since the minimization condition will also guarantee ``[\nabla V](x_0) = 0``, so ``\dot{V}(x_0) = 0``. NeuralLyapunov.jl provides the [`LyapunovDecreaseCondition`](@ref) struct for users to specify the form of the decrease condition to enforce through training. @@ -26,11 +27,21 @@ DontCheckDecrease ## Defining your own decrease condition +```@meta +CurrentModule = NeuralLyapunov +``` + +If a user wishes to define their own version of the decrease condition in a form other than +``\texttt{rate\_metric}(V(x), \dot{V}(x)) \le - \texttt{strength}(x, x_0)``, +they must define their own subtype of [`AbstractLyapunovDecreaseCondition`](@ref). + ```@docs -NeuralLyapunov.AbstractLyapunovDecreaseCondition +AbstractLyapunovDecreaseCondition ``` +When constructing the PDESystem, [`NeuralLyapunovPDESystem`](@ref) uses [`check_decrease`](@ref) to determine if it should include an equation equating the result of [`get_decrease_condition`](@ref) to zero. + ```@docs -NeuralLyapunov.check_decrease -NeuralLyapunov.get_decrease_condition +check_decrease +get_decrease_condition ``` diff --git a/docs/src/man/minimization.md b/docs/src/man/minimization.md index 04de57d..c7d74b7 100644 --- a/docs/src/man/minimization.md +++ b/docs/src/man/minimization.md @@ -21,19 +21,23 @@ DontCheckNonnegativity ## Defining your own minimization condition +```@meta +CurrentModule = NeuralLyapunov +``` + If a user wishes to define their own version of the minimization condition in a form other than ``V(x) ≥ \texttt{strength}(x, x_0)``, -they must define their own subtype of [`NeuralLyapunov.AbstractLyapunovMinimizationCondition`](@ref). +they must define their own subtype of [`AbstractLyapunovMinimizationCondition`](@ref). ```@docs -NeuralLyapunov.AbstractLyapunovMinimizationCondition +AbstractLyapunovMinimizationCondition ``` -When constructing the PDESystem, [`NeuralLyapunovPDESystem`](@ref) uses [`NeuralLyapunov.check_nonnegativity`](@ref) to determine if it should include an equation equating the result of [`NeuralLyapunov.get_minimization_condition`](@ref) to zero. -It additionally uses [`NeuralLyapunov.check_minimal_fixed_point`](@ref) to determine if it should include the equation ``V(x_0) = 0``. +When constructing the PDESystem, [`NeuralLyapunovPDESystem`](@ref) uses [`check_nonnegativity`](@ref) to determine if it should include an equation equating the result of [`get_minimization_condition`](@ref) to zero. +It additionally uses [`check_minimal_fixed_point`](@ref) to determine if it should include the equation ``V(x_0) = 0``. ```@docs -NeuralLyapunov.check_nonnegativity -NeuralLyapunov.check_minimal_fixed_point -NeuralLyapunov.get_minimization_condition +check_nonnegativity +check_minimal_fixed_point +get_minimization_condition ``` diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl index 5481212..9b62677 100644 --- a/src/conditions_specification.jl +++ b/src/conditions_specification.jl @@ -115,8 +115,11 @@ end """ get_decrease_condition(cond::AbstractLyapunovDecreaseCondition) -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. +Return a function of ``V``, ``V̇``, ``x``, and ``x_0`` that returns zero when the Lyapunov +decrease condition is met and a value greater than zero when it is violated. + +Note that the first two inputs, ``V`` and ``V̇``, are functions, so the decrease condition +can depend on the value of these functions at multiple points. """ 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 bce2e15..83e31c6 100644 --- a/src/decrease_conditions.jl +++ b/src/decrease_conditions.jl @@ -113,8 +113,8 @@ end Construct a [`LyapunovDecreaseCondition`](@ref) corresponding to exponential decrease of rate ``k``. -If `strict == false`, the condition is ``\\dot{V}(x) ≤ -k * V(x)``, and if `strict == true`, -the condition is ``\\dot{V}(x) ≤ -k * V(x) - C * \\lVert x - x_0 \\rVert^2``. +If `strict == false`, the condition is ``\\dot{V}(x) ≤ -k V(x)``, and if `strict == true`, +the condition is ``\\dot{V}(x) ≤ -k V(x) - C \\lVert x - x_0 \\rVert^2``. The inequality is represented by ``\\texttt{rectifier}(\\dot{V}(x) + k V(x) + C \\lVert x - x_0 \\rVert^2) = 0``. From e1c55fc3ae111434afdda20b0dc2a1c397040f06 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 24 May 2024 16:57:31 -0400 Subject: [PATCH 10/34] Updated documentation for `local_lyapunov` --- docs/src/man/local_lyapunov.md | 2 ++ docs/src/man/structure.md | 4 ++-- src/local_lyapunov.jl | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/src/man/local_lyapunov.md b/docs/src/man/local_lyapunov.md index c4804f8..16aa55b 100644 --- a/docs/src/man/local_lyapunov.md +++ b/docs/src/man/local_lyapunov.md @@ -1,5 +1,7 @@ # Local Lyapunov analysis +For comparison with neural Lyapunov methods, we also provide a function for local Lyapunov analysis by linearization. + ```@docs local_lyapunov ``` diff --git a/docs/src/man/structure.md b/docs/src/man/structure.md index 85cccfa..688e9b1 100644 --- a/docs/src/man/structure.md +++ b/docs/src/man/structure.md @@ -1,9 +1,9 @@ # Structuring a neural Lyapunov function -## Pre-defined structures - Three simple neural Lyapunov function structures are provided, but users can always specify a custom structure using the [`NeuralLyapunovStructure`](@ref) struct. +## Pre-defined structures + The simplest structure is to train the neural network directly to be the Lyapunov function, which can be accomplished using an [`UnstructuredNeuralLyapunov`](@ref). ```@docs diff --git a/src/local_lyapunov.jl b/src/local_lyapunov.jl index a65e571..25f9cac 100644 --- a/src/local_lyapunov.jl +++ b/src/local_lyapunov.jl @@ -1,5 +1,5 @@ """ - get_local_lyapunov(dynamics, state_dim, optimizer_factory[, jac]; fixed_point, p) + 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)`. From 421113b2e5af1e0f019cdfb67199289475b205fc Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 24 May 2024 18:11:06 -0400 Subject: [PATCH 11/34] Updated documentation for RoA-aware training --- docs/src/man/roa.md | 10 ++- src/decrease_conditions.jl | 7 ++- src/decrease_conditions_RoA_aware.jl | 92 ++++++++++++++++++++-------- 3 files changed, 77 insertions(+), 32 deletions(-) diff --git a/docs/src/man/roa.md b/docs/src/man/roa.md index 4795a2f..80c2ef0 100644 --- a/docs/src/man/roa.md +++ b/docs/src/man/roa.md @@ -1,9 +1,13 @@ # Training for Region of Attraction Identification -```@docs -RoAAwareDecreaseCondition -``` +Satisfying the [minimization](minimization.md) and [decrease](decrease.md) conditions within the training region (or any region around the fixed point, however small) is sufficient for proving local stability. +In many cases, however, we desire an estimate of the region of attraction, rather than simply a guarantee of local stability. + +Any compact sublevel set wherein the minimization and decrease conditions are satisfied is an inner estimate of the region of attraction. +Therefore, we can restrict training for those conditions to only within a predetermined sublevel set ``\{ x : V(x) \le \rho \}``. +To do so, define a [`LyapunovDecreaseCondition`](@ref) as usual and then pass it through the [`make_RoA_aware`](@ref) function, which returns an analogous [`RoAAwareDecreaseCondition`](@ref). ```@docs make_RoA_aware +RoAAwareDecreaseCondition ``` diff --git a/src/decrease_conditions.jl b/src/decrease_conditions.jl index 83e31c6..6b74e83 100644 --- a/src/decrease_conditions.jl +++ b/src/decrease_conditions.jl @@ -1,7 +1,7 @@ """ LyapunovDecreaseCondition(check_decrease, rate_metric, strength, rectifier) -Specifies the form of the Lyapunov conditions to be used. +Specifies the form of the Lyapunov decrease condition to be used. # Fields - `check_decrease::Bool`: whether or not to train for negativity/nonpositivity of @@ -14,6 +14,8 @@ Specifies the form of the Lyapunov conditions to be used. - `rectifier::Function`: positive when the input is positive and (approximately) zero when the input is negative. +# Training conditions + If `check_decrease == true`, training will enforce: ``\\texttt{rate\\_metric}(V(x), V̇(x)) ≤ -\\texttt{strength}(x, x_0).`` @@ -41,7 +43,7 @@ for some positive ``C``, which corresponds to This can also be accomplished with [`AsymptoticDecrease`](@ref). Exponential decrease of rate ``k`` is proven by - ``V̇(x) ≤ - k * V(x)``, + ``V̇(x) ≤ - k V(x)``, which corresponds to rate_metric = (V, dVdt) -> dVdt + k * V @@ -49,7 +51,6 @@ which corresponds to This can also be accomplished with [`ExponentialDecrease`](@ref). - In either case, the rectified linear unit `rectifier = (t) -> max(zero(t), t)` exactly represents the inequality, but differentiable approximations of this function may be employed. diff --git a/src/decrease_conditions_RoA_aware.jl b/src/decrease_conditions_RoA_aware.jl index 7120a46..321394a 100644 --- a/src/decrease_conditions_RoA_aware.jl +++ b/src/decrease_conditions_RoA_aware.jl @@ -1,40 +1,77 @@ """ - RoAAwareDecreaseCondition(check_decrease, rate_metric, strength, rectifier, ρ, - out_of_RoA_penalty) + RoAAwareDecreaseCondition(check_decrease, rate_metric, strength, rectifier, ρ, out_of_RoA_penalty) -Specifies the form of the Lyapunov conditions to be used, training for a region of -attraction estimate of `{ x : V(x) ≤ ρ }` +Specifies the form of the Lyapunov decrease condition to be used, training for a region of +attraction estimate of ``\\{ x : V(x) ≤ ρ \\}``. -If `check_decrease`, training will enforce -`rate_metric(V(state), dVdt(state)) ≤ -strength(state, fixed_point)` whenever `V(state) ≤ ρ`, -and will instead apply -`|out_of_RoA_penalty(V(state), dVdt(state), state, fixed_point, ρ)|^2` when `V(state) > ρ`. +# Fields + - `check_decrease::Bool`: whether or not to train for negativity/nonpositivity of + ``V̇(x)``. + - `rate_metric::Function`: should increase with ``V̇(x)``; used when + `check_decrease == true`. + - `strength::Function`: specifies the level of strictness for negativity training; should + be zero when the two inputs are equal and nonnegative otherwise; used when + `check_decrease == true`. + - `rectifier::Function`: positive when the input is positive and (approximately) zero when + the input is negative. + - `ρ`: the level of the sublevel set forming the estimate of the region of attraction. + - `out_of_RoA_penalty::Function`: a loss function to be applied penalizing points outside + the sublevel set forming the region of attraction estimate. + +# Training conditions + +If `check_decrease == true`, training will enforce + +``\\texttt{rate\\_metric}(V(x), V̇(x)) ≤ - \\texttt{strength}(x, x_0)`` + +whenever ``V(x) ≤ ρ``, and will instead apply a loss of + +``\\lvert \\texttt{out\\_of\\_RoA\\_penalty}(V(x), V̇(x), x, x_0, ρ) \\rvert^2`` + +when ``V(x) > ρ``. The inequality will be approximated by the equation - `rectifier(rate_metric(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 -to enforce `dVdt(fixed_point) = 0`, so `check_fixed_point` defaults to `false`. +``\\texttt{rectifier}(\\texttt{rate\\_metric}(V(x), V̇(x)) + \\texttt{strength}(x, x_0)) = 0``. + +Note that the approximate equation and inequality are identical when +``\\texttt{rectifier}(t) = \\max(0, t)``. + +If the dynamics truly have a fixed point at ``x_0`` and ``V̇(x)`` is truly the rate of +decrease of ``V(x)`` along the dynamics, then ``V̇(x_0)`` will be ``0`` and there is no need +to train for ``V̇(x_0) = 0``. # Examples: Asymptotic decrease can be enforced by requiring - `dVdt ≤ -C |state - fixed_point|^2`, + ``V̇(x) ≤ -C \\lVert x - x_0 \\rVert^2``, +for some positive ``C``, which corresponds to + + rate_metric = (V, dVdt) -> dVdt + strength = (x, x0) -> C * (x - x0) ⋅ (x - x0) + +Exponential decrease of rate ``k`` is proven by + ``V̇(x) ≤ - k V(x)``, which corresponds to - `rate_metric = (V, dVdt) -> dVdt` and - `strength = (x, x0) -> C * (x - x0) ⋅ (x - x0)`. -Exponential decrease of rate `k` is proven by `dVdt ≤ - k * V`, so corresponds to - `rate_metric = (V, dVdt) -> dVdt + k * V` and - `strength = (x, x0) -> 0.0`. + rate_metric = (V, dVdt) -> dVdt + k * V + strength = (x, x0) -> 0.0 Enforcing either condition only in the region of attraction and not penalizing any points outside that region would correspond to - `out_of_RoA_penalty = (V, dVdt, state, fixed_point, ρ) -> 0.0`, + + out_of_RoA_penalty = (V, dVdt, state, fixed_point, ρ) -> 0.0 + whereas an example of a penalty that decays farther in state space from the fixed point is - `out_of_RoA_penalty = (V, dVdt, state, fixed_point, ρ) -> 1.0 / ((x - x0) ⋅ (x - x0))`. -Note that this penalty could also depend on values of `V`, `dVdt`, and `ρ`. + + out_of_RoA_penalty = (V, dVdt, x, x0, ρ) -> 1.0 / ((x - x0) ⋅ (x - x0)) + +Note that this penalty could also depend on values of ``V`` and ``V̇`` at various points, as +well as ``ρ``. + +In any of these cases, the rectified linear unit `rectifier = (t) -> max(zero(t), t)` +exactly represents the inequality, but differentiable approximations of this function may be +employed. """ struct RoAAwareDecreaseCondition <: AbstractLyapunovDecreaseCondition check_decrease::Bool @@ -67,12 +104,15 @@ end make_RoA_aware(cond; ρ, out_of_RoA_penalty) Add awareness of the region of attraction (RoA) estimation task to the supplied -`LyapunovDecreaseCondition` +[`LyapunovDecreaseCondition`](@ref). + +# Arguments + - `cond::LyapunovDecreaseCondition`: specifies the loss to be applied when ``V(x) ≤ ρ``. + - `ρ`: the target level such that the RoA will be ``\\{ x : V(x) ≤ ρ \\}``. + - `out_of_RoA_penalty::Function`: specifies the loss to be applied when ``V(x) > ρ``. -`ρ` is the target level such that the RoA will be `{ x : V(x) ≤ ρ }`. -`cond` specifies the loss applied when `V(state) ≤ ρ`, and -`|out_of_RoA_penalty(V(state), dVdt(state), state, fixed_point, ρ)|^2` is the loss from -`state` values such that `V(state) > ρ`. +The loss applied to samples ``x`` such that ``V(x) > ρ`` is +``\\lvert \\texttt{out\\_of\\_RoA\\_penalty}(V(x), V̇(x), x, x_0, ρ) \\rvert^2``. """ function make_RoA_aware( cond::LyapunovDecreaseCondition; From 0314e66748a022a3de3013f84ab225b366302b1e Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Mon, 27 May 2024 13:26:37 -0400 Subject: [PATCH 12/34] Added warning about symbolic tracing to documentation --- docs/src/man/decrease.md | 3 ++- docs/src/man/pdesystem.md | 9 ++++++++- src/numerical_lyapunov_functions.jl | 2 +- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/docs/src/man/decrease.md b/docs/src/man/decrease.md index 5294425..6152caa 100644 --- a/docs/src/man/decrease.md +++ b/docs/src/man/decrease.md @@ -1,8 +1,9 @@ # Lyapunov Decrease Condition To represent the condition that the Lyapunov function ``V(x)`` must decrease along system trajectories, we typically introduce a new function ``\dot{V}(x) = \nabla V(x) \cdot f(x)``. -It is then sufficient to show that ``\dot{V}(x)`` is negative away from the fixed point and zero at the fixed point, since ``\dot{V}`` represents the rate of change of ``V`` along system trajectories. +This function represents the rate of change of ``V`` along system trajectories. That is to say, if ``x(t)`` is a trajectory defined by ``\frac{dx}{dt} = f(x)``, then ``\dot{V}(x(t)) = \frac{d}{dt} [ V(x(t)) ]``. +It is then sufficient to show that ``\dot{V}(x)`` is negative away from the fixed point and zero at the fixed point, since a negative derivative means a decreasing function. Put mathematically, it is sufficient to require ``\dot{V}(x) < 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0``. We call such functions negative definite (around the fixed point ``x_0``). diff --git a/docs/src/man/pdesystem.md b/docs/src/man/pdesystem.md index 4270337..8cd353b 100644 --- a/docs/src/man/pdesystem.md +++ b/docs/src/man/pdesystem.md @@ -9,10 +9,17 @@ Candidate Lyapunov functions will be trained within a box domain subset of the s NeuralLyapunovPDESystem ``` +!!! warning + + When using [`NeuralLyapunovPDESystem`](@ref), the Lyapuonv function structure, minimization and decrease conditions, and dynamics will all be symbolically traced to generate the resulting `PDESystem` equations. + In some cases tracing results in more efficient code, but in others it can result in inefficiencies or even errors. + + If the generated `PDESystem` is then used with NeuralPDE.jl, that library's parser will convert the equations into Julia functions representing the loss, which presents another opportunity for unexpected errors. + ## Extracting the numerical Lyapunov function We provide the following convenience function for generating the Lyapunov function after the parameters have been found. -If the `PDESystem` was solved using `NeuralPDE.jl` and `Optimization.jl`, then the argument `phi` is a field of the output of `NeuralPDE.discretize` and the argument `θ` is `res.u.depvar` where `res` is the result of the optimization. +If the `PDESystem` was solved using NeuralPDE.jl and Optimization.jl, then the argument `phi` is a field of the output of `NeuralPDE.discretize` and the argument `θ` is `res.u.depvar` where `res` is the result of the optimization. ```@docs get_numerical_lyapunov_function diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl index fb1e5c1..47c2402 100644 --- a/src/numerical_lyapunov_functions.jl +++ b/src/numerical_lyapunov_functions.jl @@ -22,7 +22,7 @@ These functions can operate on a state vector or columnwise on a matrix of state # Keyword Arguments - `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``; + ``V̇(x)`` is calculated using `deriv` as ``\\frac{∂}{∂t} 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`. From 1d6fa4474764610ad3bd818b6a38db917c42ab4d Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Mon, 27 May 2024 16:07:49 -0400 Subject: [PATCH 13/34] Added policy search documentation --- docs/src/man/policy_search.md | 12 +++++--- src/numerical_lyapunov_functions.jl | 22 +++++++------- src/policy_search.jl | 46 +++++++++++++++++++---------- 3 files changed, 49 insertions(+), 31 deletions(-) diff --git a/docs/src/man/policy_search.md b/docs/src/man/policy_search.md index c6be6cc..aaa60ce 100644 --- a/docs/src/man/policy_search.md +++ b/docs/src/man/policy_search.md @@ -1,9 +1,13 @@ -# Policy Search +# Policy Search and Network-dependent Dynamics -```@docs -add_policy_search -``` +At times, we wish to model a component of the dynamics with a neural network. +A common example is the policy search case, when the closed-loop dynamics include a neural network controller. +In such cases, we consider the dynamics to take the form of ``\frac{dx}{dt} = f(x, u, p, t)``, where ``u`` is the control input/the contribution to the dynamics from the neural network. +We provide the [`add_policy_search`](@ref) function to transform a [`NeuralLyapunovStructure`](@ref) to include training the neural network to represent not just the Lyapunov function, but also the relevant part of the dynamics. + +Similar to [`get_numerical_lyapunov_function`](@ref), we provide the [`get_policy`](@ref) convenience function to construct ``u(x)`` that can be combined with the open-loop dynamics ``f(x, u, p, t)`` to create closed loop dynamics ``f_{cl}(x, p, t) = f(x, u(x), p, t)``. ```@docs +add_policy_search get_policy ``` diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl index 47c2402..7ff37fe 100644 --- a/src/numerical_lyapunov_functions.jl +++ b/src/numerical_lyapunov_functions.jl @@ -10,20 +10,20 @@ These functions can operate on a state vector or columnwise on a matrix of state # Positional Arguments - `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 + - `θ`: 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. + - `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. # Keyword Arguments - `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{∂}{∂t} V(x + t f(x))`` at ``t = 0``; - defaults to `false`, as it is more efficient in many cases. + - `use_V̇_structure`: when `true`, ``V̇(x)`` is calculated using `structure.V̇`; when ` + false`, ``V̇(x)`` is calculated using `deriv` as ``\\frac{∂}{∂t} 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) @@ -98,11 +98,11 @@ Return the network as a function of state alone. # Arguments - `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 + - `θ`: 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`. + - `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 diff --git a/src/policy_search.jl b/src/policy_search.jl index b4bb214..38a9d9c 100644 --- a/src/policy_search.jl +++ b/src/policy_search.jl @@ -1,18 +1,23 @@ """ - add_policy_search(lyapunov_structure, new_dims, control_structure) + add_policy_search(lyapunov_structure, new_dims; control_structure) -Add dependence on the neural network to the dynamics in a `NeuralLyapunovStructure`. +Add dependence on the neural network to the dynamics in a [`NeuralLyapunovStructure`](@ref). -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. +# Arguments + - `lyapunov_structure::NeuralLyapunovStructure`: provides structure for ``V, V̇``; should + assume dynamics take a form of `f(x, p, t)`. + - `new_dims::Integer`: number of outputs of the neural network to pass into the dynamics + through `control_structure`. + - `control_structure::Function`: transforms the final `new_dims` outputs of the neural net + before passing them into the dynamics; defaults to passing in the neural network outputs + unchanged. -`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 -contribution from the neural network. Therefore, this structure cannot be used with a -`NeuralLyapunovPDESystem` method that requires an `ODEFunction`, `ODESystem`, or -`ODEProblem`. +The returned `NeuralLyapunovStructure` expects dynamics of the form `f(x, u, p, t)`, where +`u` captures the dependence of dynamics on the neural network (e.g., through a control +input). When evaluating the dynamics, it uses `u = control_structure(phi_end(x))` where +`phi_end` is a function that returns the final `new_dims` outputs of the neural network. +The other `lyapunov_structure.network_dim` outputs are used for calculating ``V`` and ``V̇``, +as specified originally by `lyapunov_structure`. """ function add_policy_search( lyapunov_structure::NeuralLyapunovStructure, @@ -46,16 +51,25 @@ function add_policy_search( end """ - get_policy(phi, θ, network_func, dim; control_structure) + get_policy(phi, θ, network_dim, control_dim; control_structure) -Generate a Julia function representing the control policy as a function of the state +Generate a Julia function representing the control policy/unmodeled portion of the dynamics +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 `θ`. The network should have `network_dim` -outputs, the last `control_dim` of which will be passed into `control_structure` to create -the policy output. +# 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. + - `network_dim`: total number of neural network outputs. + - `control_dim`: number of neural network outputs used in the control policy. + - `control_structure`: transforms the final `control_dim` outputs of the neural net before + passing them into the dynamics; defaults to passing in the neural network outputs + unchanged. """ function get_policy( phi, From ae4da03a5011eea9f6d2a6a38594188d4d2a20ee Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Fri, 31 May 2024 18:37:40 -0400 Subject: [PATCH 14/34] Added damped SHO demo to docs --- docs/Project.toml | 7 + docs/src/demos/damped_SHO.md | 269 ++++++++++++++++++++++++++++++++++- test/damped_sho.jl | 44 +++--- 3 files changed, 297 insertions(+), 23 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 6118f9b..34a5ac0 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,10 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Lux = "b2108857-7c20-44ae-9111-449ecde12c47" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NeuralLyapunov = "61252e1c-87f0-426a-93a7-3cdb1ed5a867" +NeuralPDE = "315f7962-48a3-4962-8226-d0f33b1235f0" +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" +OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md index 80eb854..0bb1a07 100644 --- a/docs/src/demos/damped_SHO.md +++ b/docs/src/demos/damped_SHO.md @@ -1 +1,268 @@ -# Damped Simple Harmonic Oscillator \ No newline at end of file +# Damped Simple Harmonic Oscillator + +Let's train a neural network to prove the exponential stability of the damped simple harmonic oscillator (SHO). + +The damped SHO is represented by the system of differential equations +```math +\begin{align*} + \frac{dx}{dt} &= v \\ + \frac{dv}{dt} &= -2 \zeta \omega_0 v - \omega_0^2 x +\end{align*} +``` +where ``x`` is the position, ``v`` is the velocity, ``t`` is time, and ``\zeta, \omega_0`` are parameters. + +We'll consider just the box domain ``x \in [-5, 5], v \in [-2, 2]``. + +## Copy-Pastable Code + +```julia +using LinearAlgebra +using NeuralPDE, Lux, ModelingToolkit +using Optimization, OptimizationOptimisers, OptimizationOptimJL +using NeuralLyapunov +using Random + +Random.seed!(200) + +######################### Define dynamics and domain ########################## + +"Simple Harmonic Oscillator Dynamics" +function f(state, p, t) + ζ, ω_0 = p + pos = state[1] + vel = state[2] + vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos) +end +lb = [-5.0, -2.0]; +ub = [ 5.0, 2.0]; +p = [0.5, 1.0]; +fixed_point = zeros(2); +dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) + +####################### Specify neural Lyapunov problem ####################### + +# Define neural network discretization +dim_state = length(lb) +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) + ) for _ in 1:dim_output] + +# Define neural network discretization +strategy = GridTraining(0.05) +discretization = PhysicsInformedNN(chain, strategy) + +# Define neural Lyapunov structure +structure = NonnegativeNeuralLyapunov( + dim_output; + δ = 1e-6 +) +minimization_condition = DontCheckNonnegativity(check_fixed_point = true) + +# Define Lyapunov decrease condition +# 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( + structure, + minimization_condition, + decrease_condition +) + +############################# Construct PDESystem ############################# + +@named pde_system = NeuralLyapunovPDESystem( + dynamics, + lb, + ub, + spec; + p = p +) + +######################## Construct OptimizationProblem ######################## + +prob = discretize(pde_system, discretization) +sym_prob = symbolic_discretize(pde_system, discretization) + +########################## Solve OptimizationProblem ########################## + +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500) +prob = Optimization.remake(prob, u0 = res.u) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 500) + +###################### Get numerical numerical functions ###################### +net = discretization.phi +θ = res.u.depvar + +V, V̇ = get_numerical_lyapunov_function( + net, + θ, + structure, + f, + fixed_point; + p = p +) +``` + +## Detailed description + +In this example, we set the dynamics up as an `ODEFunction` and use a `SciMLBase.SymbolCache` to tell the ultimate `PDESystem` what to call our state and parameter variables. + +```@example SHO +using NeuralPDE # for ODEFunction and SciMLBase.SymbolCache + +"Simple Harmonic Oscillator Dynamics" +function f(state, p, t) + ζ, ω_0 = p + pos = state[1] + vel = state[2] + vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos) +end +lb = [-5.0, -2.0]; +ub = [ 5.0, 2.0]; +p = [0.5, 1.0]; +fixed_point = zeros(2); +dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) +``` + +Setting up the neural network using Lux and NeuralPDE training strategy is no different from any other physics-informed neural network problem. +For more on that aspect, see the [NeuralPDE documentation](https://docs.sciml.ai/NeuralPDE/stable/). + +```@example SHO +using Lux + +# Define neural network discretization +dim_state = length(lb) +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) + ) for _ in 1:dim_output] +``` + +```@example SHO +# Define training strategy +strategy = GridTraining(0.05) +discretization = PhysicsInformedNN(chain, strategy) +``` + +We now define our Lyapunov candidate structure along with the form of the Lyapunov conditions we'll be using. + +For this example, let's use a Lyapunov candidate +```math +V(x) = \lVert \phi(x) \rVert^2 + \delta \log \left( 1 + \lVert x \rVert^2 \right), +``` +which structurally enforces nonnegativity, but doesn't guarantee ``V([0, 0]) = 0``. +We therefore don't need a term in the loss function enforcing ``V(x) > 0 \, \forall x \ne 0``, but we do need something enforcing ``V([0, 0]) = 0``. +So, we use [`DontCheckNonnegativity(check_fixed_point = true)`](@ref). + +To train for exponential decrease we use [`ExponentialDecrease`](@ref), but we must specify the rate of exponential decrease, which we know in this case to be ``\zeta \omega_0``. + +```@example SHO +using NeuralLyapunov + +# Define neural Lyapunov structure +structure = NonnegativeNeuralLyapunov( + dim_output; + δ = 1e-6 +) +minimization_condition = DontCheckNonnegativity(check_fixed_point = true) + +# Define Lyapunov decrease condition +# 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( + structure, + minimization_condition, + decrease_condition +) + +# Construct PDESystem +@named pde_system = NeuralLyapunovPDESystem( + dynamics, + lb, + ub, + spec; + p = p +); +``` + +Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem. + +```@example SHO +prob = discretize(pde_system, discretization) + +using Optimization, OptimizationOptimisers, OptimizationOptimJL + +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500) +prob = Optimization.remake(prob, u0 = res.u) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 500) + +net = discretization.phi +θ = res.u.depvar +``` + +We can use the result of the optimization problem to build the Lyapunov candidate as a Julia function. + +```@example SHO +V, V̇ = get_numerical_lyapunov_function( + net, + θ, + structure, + f, + fixed_point; + p = p +) +``` + +Now let's see how we did. +We'll evaluate both ``V`` and ``\dot{V}`` on a finer grid than we trained on: + +```@example SHO +xs, ys = [lb[i]:0.02:ub[i] for i in eachindex(lb)] +states = Iterators.map(collect, Iterators.product(xs, ys)) +V_samples = vec(V(hcat(states...))) +V̇_samples = vec(V̇(hcat(states...))) + +# Print statistics +println("V(0.,0.) = ", V(fixed_point)) +println("V̇(0.,0.) = ", V̇(fixed_point)) +println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]") +println( + "V̇ ∋ [", + minimum(V̇_samples), + ", ", + max(V̇(fixed_point), maximum(V̇_samples)), + "]", +) +``` + +At least at these validation samples, ``\dot{V}`` is at least negative semi-definite and ``V`` is almost positive definite and is minimized (though perhaps not uniquely) at the origin. + +```@example SHO +using Plots + +p1 = plot(xs, ys, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "v"); +p1 = scatter!([0], [0], label = "Equilibrium"); +p2 = plot( + xs, + ys, + V̇_samples, + linetype = :contourf, + title = "V̇", + xlabel = "x", + ylabel = "v", +); +p2 = scatter!([0], [0], label = "Equilibrium"); +plot(p1, p2) +``` + +Each sublevel set of ``V`` completely contained in the plot above has been verified as a subset of the region of attraction. \ No newline at end of file diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 429d810..1f0231f 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -1,6 +1,6 @@ using LinearAlgebra using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt +using Optimization, OptimizationOptimisers, OptimizationOptimJL using NeuralLyapunov using Random using Test @@ -18,8 +18,8 @@ function f(state, p, t) vel = state[2] vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos) end -lb = [-2 * pi, -2.0]; -ub = [2 * pi, 2.0]; +lb = [-5.0, -2.0]; +ub = [5.0, 2.0]; p = [0.5, 1.0] dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) @@ -35,7 +35,7 @@ chain = [Lux.Chain( Dense(dim_hidden, 1) ) for _ in 1:dim_output] -# Define neural network discretization +# Define training strategy strategy = GridTraining(0.05) discretization = PhysicsInformedNN(chain, strategy) @@ -76,10 +76,10 @@ sym_prob = symbolic_discretize(pde_system, discretization) res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, BFGS(); maxiters = 500) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 500) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = get_numerical_lyapunov_function( +V, V̇ = get_numerical_lyapunov_function( discretization.phi, res.u.depvar, structure, @@ -92,49 +92,49 @@ V_func, V̇_func = get_numerical_lyapunov_function( ################################## Simulate ################################### xs, ys = [lb[i]:0.02:ub[i] for i in eachindex(lb)] states = Iterators.map(collect, Iterators.product(xs, ys)) -V_predict = vec(V_func(hcat(states...))) -dVdt_predict = vec(V̇_func(hcat(states...))) +V_samples = vec(V(hcat(states...))) +V̇_samples = vec(V̇(hcat(states...))) #################################### Tests #################################### # Network structure should enforce nonegativeness of V -@test min(V_func([0.0, 0.0]), minimum(V_predict)) ≥ 0.0 +@test min(V([0.0, 0.0]), minimum(V_samples)) ≥ 0.0 # Trained for V's minimum to be at the fixed point -@test V_func([0.0, 0.0])≈minimum(V_predict) atol=1e-4 -@test V_func([0.0, 0.0]) < 1e-4 +@test V([0.0, 0.0])≈minimum(V_samples) atol=1e-4 +@test V([0.0, 0.0]) < 1e-4 # Dynamics should result in a fixed point at the origin -@test V̇_func([0.0, 0.0]) == 0.0 +@test V̇([0.0, 0.0]) == 0.0 # V̇ should be negative almost everywhere -@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 1e-5 +@test sum(V̇_samples .> 0) / length(V̇_samples) < 1e-5 #= # Get RoA Estimate -data = reshape(V_predict, (length(xs), length(ys))); +data = reshape(V_samples, (length(xs), length(ys))); data = vcat(data[1, :], data[end, :], data[:, 1], data[:, end]); ρ = minimum(data) # Print statistics -println("V(0.,0.) = ", V_func([0.0, 0.0])) -println("V ∋ [", min(V_func([0.0, 0.0]), minimum(V_predict)), ", ", maximum(V_predict), "]") +println("V(0.,0.) = ", V([0.0, 0.0])) +println("V ∋ [", min(V([0.0, 0.0]), minimum(V_samples)), ", ", maximum(V_samples), "]") println( "V̇ ∋ [", - minimum(dVdt_predict), + minimum(V̇_samples), ", ", - max(V̇_func([0.0, 0.0]), maximum(dVdt_predict)), + max(V̇([0.0, 0.0]), maximum(V̇_samples)), "]", ) # Plot results -p1 = plot(xs, ys, V_predict, linetype = :contourf, title = "V", xlabel = "x", ylabel = "ẋ"); +p1 = plot(xs, ys, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "ẋ"); p1 = scatter!([0], [0], label = "Equilibrium"); p2 = plot( xs, ys, - dVdt_predict, + V̇_samples, linetype = :contourf, title = "dV/dt", xlabel = "x", @@ -144,7 +144,7 @@ p2 = scatter!([0], [0], label = "Equilibrium"); p3 = plot( xs, ys, - V_predict .< ρ, + V_samples .< ρ, linetype = :contourf, title = "Estimated RoA", xlabel = "x", @@ -154,7 +154,7 @@ p3 = plot( p4 = plot( xs, ys, - dVdt_predict .< 0, + V̇_samples .< 0, linetype = :contourf, title = "dV/dt < 0", xlabel = "x", From e036d63127c62eefb043e9633793a84b607f3c98 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 5 Jun 2024 17:26:24 -0400 Subject: [PATCH 15/34] Changed `[0.0, 0.0]` to `fixed_point` in SHO test --- test/damped_sho.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 1f0231f..37d6c95 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -20,7 +20,8 @@ function f(state, p, t) end lb = [-5.0, -2.0]; ub = [5.0, 2.0]; -p = [0.5, 1.0] +p = [0.5, 1.0]; +fixed_point = [0.0, 0.0]; dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) ####################### Specify neural Lyapunov problem ####################### @@ -84,7 +85,7 @@ V, V̇ = get_numerical_lyapunov_function( res.u.depvar, structure, f, - zeros(2); + fixed_point; p = p, use_V̇_structure = true ) @@ -98,14 +99,14 @@ V̇_samples = vec(V̇(hcat(states...))) #################################### Tests #################################### # Network structure should enforce nonegativeness of V -@test min(V([0.0, 0.0]), minimum(V_samples)) ≥ 0.0 +@test min(V(fixed_point), minimum(V_samples)) ≥ 0.0 # Trained for V's minimum to be at the fixed point -@test V([0.0, 0.0])≈minimum(V_samples) atol=1e-4 -@test V([0.0, 0.0]) < 1e-4 +@test V(fixed_point)≈minimum(V_samples) atol=1e-4 +@test V(fixed_point) < 1e-4 # Dynamics should result in a fixed point at the origin -@test V̇([0.0, 0.0]) == 0.0 +@test V̇(fixed_point) == 0.0 # V̇ should be negative almost everywhere @test sum(V̇_samples .> 0) / length(V̇_samples) < 1e-5 @@ -117,13 +118,13 @@ data = vcat(data[1, :], data[end, :], data[:, 1], data[:, end]); ρ = minimum(data) # Print statistics -println("V(0.,0.) = ", V([0.0, 0.0])) -println("V ∋ [", min(V([0.0, 0.0]), minimum(V_samples)), ", ", maximum(V_samples), "]") +println("V(0.,0.) = ", V(fixed_point)) +println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]") println( "V̇ ∋ [", minimum(V̇_samples), ", ", - max(V̇([0.0, 0.0]), maximum(V̇_samples)), + max(V̇(fixed_point), maximum(V̇_samples)), "]", ) From ac6460a7f09c0ca7d18130a326ec28fa5793bd74 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 12:30:43 -0400 Subject: [PATCH 16/34] Sped up damped SHO test --- test/damped_sho.jl | 67 ++++++++++++++++++---------------------------- 1 file changed, 26 insertions(+), 41 deletions(-) diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 37d6c95..c57ed95 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -28,8 +28,8 @@ dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) # Define neural network discretization dim_state = length(lb) -dim_hidden = 20 -dim_output = 5 +dim_hidden = 10 +dim_output = 3 chain = [Lux.Chain( Dense(dim_state, dim_hidden, tanh), Dense(dim_hidden, dim_hidden, tanh), @@ -37,7 +37,7 @@ chain = [Lux.Chain( ) for _ in 1:dim_output] # Define training strategy -strategy = GridTraining(0.05) +strategy = QuasiRandomTraining(1000) discretization = PhysicsInformedNN(chain, strategy) # Define neural Lyapunov structure @@ -91,35 +91,40 @@ V, V̇ = get_numerical_lyapunov_function( ) ################################## Simulate ################################### -xs, ys = [lb[i]:0.02:ub[i] for i in eachindex(lb)] -states = Iterators.map(collect, Iterators.product(xs, ys)) +Δx = (ub[1] - lb[1]) / 100 +Δv = (ub[2] - lb[2]) / 100 +xs = lb[1]:Δx:ub[1] +vs = lb[2]:Δv:ub[2] +states = Iterators.map(collect, Iterators.product(xs, vs)) V_samples = vec(V(hcat(states...))) V̇_samples = vec(V̇(hcat(states...))) #################################### Tests #################################### # Network structure should enforce nonegativeness of V -@test min(V(fixed_point), minimum(V_samples)) ≥ 0.0 - -# Trained for V's minimum to be at the fixed point -@test V(fixed_point)≈minimum(V_samples) atol=1e-4 -@test V(fixed_point) < 1e-4 +V_min, i_min = findmin(V_samples) +state_min = collect(states)[i_min] +V_min, state_min = if V(fixed_point) ≤ V_min + V(fixed_point), fixed_point + else + V_min, state_min + end +@test V_min ≥ 0.0 + +# Trained for V's minimum to be near the fixed point +@test all(abs.(state_min .- fixed_point) .≤ [Δx, Δv]) # Dynamics should result in a fixed point at the origin @test V̇(fixed_point) == 0.0 # V̇ should be negative almost everywhere -@test sum(V̇_samples .> 0) / length(V̇_samples) < 1e-5 +@test sum(V̇_samples .> 0) / length(V̇_samples) < 1e-4 #= -# Get RoA Estimate -data = reshape(V_samples, (length(xs), length(ys))); -data = vcat(data[1, :], data[end, :], data[:, 1], data[:, end]); -ρ = minimum(data) - # Print statistics println("V(0.,0.) = ", V(fixed_point)) -println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]") +println("V ∋ [", V_min, ", ", maximum(V_samples), "]") +println("Minimial sample of V is at ", state_min) println( "V̇ ∋ [", minimum(V̇_samples), @@ -129,12 +134,13 @@ println( ) # Plot results +using Plots -p1 = plot(xs, ys, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "ẋ"); +p1 = plot(xs, vs, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "ẋ"); p1 = scatter!([0], [0], label = "Equilibrium"); p2 = plot( xs, - ys, + vs, V̇_samples, linetype = :contourf, title = "dV/dt", @@ -142,26 +148,5 @@ p2 = plot( ylabel = "ẋ", ); p2 = scatter!([0], [0], label = "Equilibrium"); -p3 = plot( - xs, - ys, - V_samples .< ρ, - linetype = :contourf, - title = "Estimated RoA", - xlabel = "x", - ylabel = "ẋ", - colorbar = false, -); -p4 = plot( - xs, - ys, - V̇_samples .< 0, - linetype = :contourf, - title = "dV/dt < 0", - xlabel = "x", - ylabel = "ẋ", - colorbar = false, -); -p4 = scatter!([0], [0], label = "Equilibrium"); -plot(p1, p2, p3, p4) +plot(p1, p2) =# From 574f53b3f396378d675d43d298cc74e4331386d5 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 12:50:08 -0400 Subject: [PATCH 17/34] Sped up dampd SHO documentation demo --- docs/src/demos/damped_SHO.md | 45 ++++++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md index 0bb1a07..c3430c4 100644 --- a/docs/src/demos/damped_SHO.md +++ b/docs/src/demos/damped_SHO.md @@ -36,23 +36,23 @@ end lb = [-5.0, -2.0]; ub = [ 5.0, 2.0]; p = [0.5, 1.0]; -fixed_point = zeros(2); +fixed_point = [0.0, 0.0]; dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) ####################### Specify neural Lyapunov problem ####################### # Define neural network discretization dim_state = length(lb) -dim_hidden = 20 -dim_output = 5 +dim_hidden = 10 +dim_output = 3 chain = [Lux.Chain( Dense(dim_state, dim_hidden, tanh), Dense(dim_hidden, dim_hidden, tanh), Dense(dim_hidden, 1) ) for _ in 1:dim_output] -# Define neural network discretization -strategy = GridTraining(0.05) +# Define training strategy +strategy = QuasiRandomTraining(1000) discretization = PhysicsInformedNN(chain, strategy) # Define neural Lyapunov structure @@ -125,7 +125,7 @@ end lb = [-5.0, -2.0]; ub = [ 5.0, 2.0]; p = [0.5, 1.0]; -fixed_point = zeros(2); +fixed_point = [0.0, 0.0]; dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0])) ``` @@ -137,8 +137,8 @@ using Lux # Define neural network discretization dim_state = length(lb) -dim_hidden = 20 -dim_output = 5 +dim_hidden = 10 +dim_output = 3 chain = [Lux.Chain( Dense(dim_state, dim_hidden, tanh), Dense(dim_hidden, dim_hidden, tanh), @@ -148,7 +148,7 @@ chain = [Lux.Chain( ```@example SHO # Define training strategy -strategy = GridTraining(0.05) +strategy = QuasiRandomTraining(1000) discretization = PhysicsInformedNN(chain, strategy) ``` @@ -224,18 +224,29 @@ V, V̇ = get_numerical_lyapunov_function( ``` Now let's see how we did. -We'll evaluate both ``V`` and ``\dot{V}`` on a finer grid than we trained on: +We'll evaluate both ``V`` and ``\dot{V}`` on a ``101 \times 101`` grid: ```@example SHO -xs, ys = [lb[i]:0.02:ub[i] for i in eachindex(lb)] -states = Iterators.map(collect, Iterators.product(xs, ys)) +Δx = (ub[1] - lb[1]) / 100 +Δv = (ub[2] - lb[2]) / 100 +xs = lb[1]:Δx:ub[1] +vs = lb[2]:Δv:ub[2] +states = Iterators.map(collect, Iterators.product(xs, vs)) V_samples = vec(V(hcat(states...))) V̇_samples = vec(V̇(hcat(states...))) # Print statistics +V_min, i_min = findmin(V_samples) +state_min = collect(states)[i_min] +V_min, state_min = if V(fixed_point) ≤ V_min + V(fixed_point), fixed_point + else + V_min, state_min + end + println("V(0.,0.) = ", V(fixed_point)) -println("V̇(0.,0.) = ", V̇(fixed_point)) -println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]") +println("V ∋ [", V_min, ", ", maximum(V_samples), "]") +println("Minimial sample of V is at ", state_min) println( "V̇ ∋ [", minimum(V̇_samples), @@ -245,16 +256,16 @@ println( ) ``` -At least at these validation samples, ``\dot{V}`` is at least negative semi-definite and ``V`` is almost positive definite and is minimized (though perhaps not uniquely) at the origin. +At least at these validation samples, the conditions that ``\dot{V}`` be negative semi-definite and ``V`` be minimized at the origin are nearly sastisfied. ```@example SHO using Plots -p1 = plot(xs, ys, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "v"); +p1 = plot(xs, vs, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "v"); p1 = scatter!([0], [0], label = "Equilibrium"); p2 = plot( xs, - ys, + vs, V̇_samples, linetype = :contourf, title = "V̇", From 20b725bb220e92fd1ddbb60b8a0a55dcdae19fdd Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 15:15:27 -0400 Subject: [PATCH 18/34] Fixed minimization condition documentation not showing when arguments are keyword arguments --- src/minimization_conditions.jl | 12 +++++++----- test/roa_estimation.jl | 4 ++-- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/minimization_conditions.jl b/src/minimization_conditions.jl index d7fd104..9f4355e 100644 --- a/src/minimization_conditions.jl +++ b/src/minimization_conditions.jl @@ -69,7 +69,7 @@ function get_minimization_condition(cond::LyapunovMinimizationCondition) end """ - StrictlyPositiveDefinite(C; check_fixed_point, rectifier) + StrictlyPositiveDefinite(; C, check_fixed_point, rectifier) Construct a [`LyapunovMinimizationCondition`](@ref) representing ``V(x) ≥ C \\lVert x - x_0 \\rVert^2``. @@ -93,7 +93,7 @@ function StrictlyPositiveDefinite(; end """ - PositiveSemiDefinite(check_fixed_point) + PositiveSemiDefinite(; check_fixed_point) Construct a [`LyapunovMinimizationCondition`](@ref) representing ``V(x) ≥ 0``. @@ -115,7 +115,7 @@ function PositiveSemiDefinite(; end """ - DontCheckNonnegativity(check_fixed_point) + DontCheckNonnegativity(; check_fixed_point) Construct a [`LyapunovMinimizationCondition`](@ref) which represents not checking for nonnegativity of the Lyapunov function. This is appropriate in cases where this condition @@ -123,9 +123,11 @@ has been structurally enforced. It is still possible to check for ``V(x_0) = 0``, even in this case, for example if `V` is structured to be positive for ``x ≠ x_0``, but it is not guaranteed structurally that -``V(x_0) = 0`` (such as [`NonnegativeNeuralLyapunov`](@ref)). +``V(x_0) = 0`` (such as [`NonnegativeNeuralLyapunov`](@ref)). `check_fixed_point` defaults +to `true`, since in cases where ``V(x_0) = 0`` is enforced structurally, the equation will +reduce to `0.0 ~ 0.0`, which gets removed in most cases. """ -function DontCheckNonnegativity(; check_fixed_point = false)::LyapunovMinimizationCondition +function DontCheckNonnegativity(; check_fixed_point = true)::LyapunovMinimizationCondition LyapunovMinimizationCondition( false, (state, fixed_point) -> 0.0, diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index c3d2418..e1d317a 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -1,6 +1,6 @@ using LinearAlgebra using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt +using Optimization, OptimizationOptimisers, OptimizationOptimJL using NeuralLyapunov using Random using Test @@ -63,7 +63,7 @@ sym_prob = symbolic_discretize(pde_system, discretization) res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, BFGS(); maxiters = 300) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### V_func, V̇_func = get_numerical_lyapunov_function( From c6779342c5cb4d5356ba7ce3fd99c524fafaaa7d Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 16:03:24 -0400 Subject: [PATCH 19/34] Add default behavior to structure documentation --- docs/src/demos/damped_SHO.md | 2 +- src/structure_specification.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md index c3430c4..7f3ac41 100644 --- a/docs/src/demos/damped_SHO.md +++ b/docs/src/demos/damped_SHO.md @@ -192,7 +192,7 @@ spec = NeuralLyapunovSpecification( ub, spec; p = p -); +) ``` Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem. diff --git a/src/structure_specification.jl b/src/structure_specification.jl index 042e2c3..fabe589 100644 --- a/src/structure_specification.jl +++ b/src/structure_specification.jl @@ -39,7 +39,7 @@ used. - `network_dim`: output dimensionality of the neural network. - `δ`: weight of `pos_def`, as above. - `pos_def(state, fixed_point)`: a function that is postive (semi-)definite in `state` - around `fixed_point`. + around `fixed_point`; defaults to ``log(1 + \\lVert x - x_0 \\rVert^2)``. # Keyword Arguments - `grad_pos_def(state, fixed_point)`: the gradient of `pos_def` with respect to `state` at @@ -112,10 +112,10 @@ so [`DontCheckNonnegativity(false)`](@ref) should be used. # Keyword Arguments - `pos_def(state, fixed_point)`: a function that is postive (semi-)definite in `state` - around `fixed_point`. + around `fixed_point`; defaults to ``log(1 + \\lVert x - x_0 \\rVert^2)``. - `non_neg(net, state, fixed_point)`: a nonnegative function of the neural network; note that `net` is the neural network ``ϕ``, and `net(state)` is the value of the neural - network at a point ``ϕ(x)``. + network at a point ``ϕ(x)``; defaults to ``1 + \\lVert ϕ(x) \\rVert^2``. - `grad_pos_def(state, fixed_point)`: the gradient of `pos_def` with respect to `state` at `state`. If `isnothing(grad_pos_def)` (as is the default), the gradient of `pos_def` will be evaluated using `grad`. From 1b1b2a2a167b0a713809033f9650bac55e213ae5 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 16:31:33 -0400 Subject: [PATCH 20/34] Added RoA estimation demo to documentation --- docs/make.jl | 3 +- docs/src/demos/roa_estimation.md | 245 +++++++++++++++++++++++++++++++ test/roa_estimation.jl | 39 ++--- 3 files changed, 267 insertions(+), 20 deletions(-) create mode 100644 docs/src/demos/roa_estimation.md diff --git a/docs/make.jl b/docs/make.jl index 5fc7743..ca3812d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,7 +27,8 @@ makedocs(; hide("man/internals.md") ], "Demonstrations" => [ - "demos/damped_SHO.md" + "demos/damped_SHO.md", + "demos/roa_estimation.md" ] ] ) diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md new file mode 100644 index 0000000..fd547ed --- /dev/null +++ b/docs/src/demos/roa_estimation.md @@ -0,0 +1,245 @@ +# Estimating the region of attraction + +In this demonstration, we add awareness of the region of attraction (RoA) estimation task to our training. + +We'll be examining the simple one-dimensional differential equation: +```math +\frac{dx}{dt} = - x + x^3. +``` +This system has a fixed point at ``x = 0`` which has a RoA of ``x \in (-1, 1)``, which we will attempt to identify. + +We'll train in the larger domain ``x \in [-2, 2]``. + +## Copy-Pastable Code + +```julia +using LinearAlgebra +using NeuralPDE, Lux, ModelingToolkit +using Optimization, OptimizationOptimisers, OptimizationOptimJL +using NeuralLyapunov +using Random + +Random.seed!(250) + +######################### Define dynamics and domain ########################## + +f(x, p, t) = -x .+ x .^ 3 +lb = [-2.0]; +ub = [ 2.0]; +fixed_point = [0.0]; + +####################### Specify neural Lyapunov problem ####################### + +# Define neural network discretization +dim_state = length(lb) +dim_hidden = 5 +dim_output = 2 +chain = [Lux.Chain( + Dense(dim_state, dim_hidden, tanh), + Dense(dim_hidden, dim_hidden, tanh), + Dense(dim_hidden, 1, use_bias = false) + ) for _ in 1:dim_output] + +# Define training strategy +strategy = GridTraining(0.1) +discretization = PhysicsInformedNN(chain, strategy) + +# Define neural Lyapunov structure +structure = PositiveSemiDefiniteStructure(dim_output) +minimization_condition = DontCheckNonnegativity() + +# Define Lyapunov decrease condition +decrease_condition = make_RoA_aware(AsymptoticDecrease(strict = true)) + +# Construct neural Lyapunov specification +spec = NeuralLyapunovSpecification( + structure, + minimization_condition, + decrease_condition +) + +############################# Construct PDESystem ############################# + +@named pde_system = NeuralLyapunovPDESystem( + f, + lb, + ub, + spec +) + +######################## Construct OptimizationProblem ######################## + +prob = discretize(pde_system, discretization) +sym_prob = symbolic_discretize(pde_system, discretization) + +########################## Solve OptimizationProblem ########################## + +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) +prob = Optimization.remake(prob, u0 = res.u) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) + +###################### Get numerical numerical functions ###################### +net = discretization.phi +θ = res.u.depvar + +V, V̇ = get_numerical_lyapunov_function( + net, + θ, + structure, + f, + fixed_point +) + +################################## Simulate ################################### +states = lb[]:0.001:ub[] +V_samples = vec(V(states')) +V̇_samples = vec(V̇(states')) + +# Calculated RoA estimate +ρ = decrease_condition.ρ +RoA_states = states[vec(V(transpose(states))) .≤ ρ] +RoA = (first(RoA_states), last(RoA_states)) +``` + +## Detailed description + +In this example, we set up the dynamics as a Julia function and don't bother specifying the symbols for the variables (so ``x`` will be called the default `state1` in the PDESystem). + +```@example RoA +f(x, p, t) = -x .+ x .^ 3 +lb = [-2.0]; +ub = [ 2.0]; +fixed_point = [0.0]; +nothing # hide +``` + +Setting up the neural network using Lux and NeuralPDE training strategy is no different from any other physics-informed neural network problem. +For more on that aspect, see the [NeuralPDE documentation](https://docs.sciml.ai/NeuralPDE/stable/). +Since we're only considering one dimension, training on a grid isn't so bad in this case. + +```@example RoA +using Lux + +# Define neural network discretization +dim_state = length(lb) +dim_hidden = 5 +dim_output = 2 +chain = [Lux.Chain( + Dense(dim_state, dim_hidden, tanh), + Dense(dim_hidden, dim_hidden, tanh), + Dense(dim_hidden, 1, use_bias = false) + ) for _ in 1:dim_output] +``` + +```@example RoA +using NeuralPDE + +# Define training strategy +strategy = GridTraining(0.1) +discretization = PhysicsInformedNN(chain, strategy) +``` + +We now define our Lyapunov candidate structure along with the form of the Lyapunov conditions we'll be using. + +For this example, let's use the default Lyapunov candidate from [`PositiveSemiDefiniteStructure`](@ref): +```math +V(x) = \left( 1 + \lVert \phi(x) \rVert^2 \right) \log \left( 1 + \lVert x \rVert^2 \right), +``` +which structurally enforces positive definiteness. +We therefore use [`DontCheckNonnegativity()`](@ref). + +We only require asymptotic decrease in this example, but we use `make_RoA_aware` to only penalize positive values of ``\dot{V}(x)`` when ``V(x) \le 1``. + +```@example RoA +using NeuralLyapunov + +# Define neural Lyapunov structure +structure = PositiveSemiDefiniteStructure(dim_output) +minimization_condition = DontCheckNonnegativity() + +# Define Lyapunov decrease condition +decrease_condition = make_RoA_aware(AsymptoticDecrease(strict = true)) + +# Construct neural Lyapunov specification +spec = NeuralLyapunovSpecification( + structure, + minimization_condition, + decrease_condition +) + +# Construct PDESystem +@named pde_system = NeuralLyapunovPDESystem( + f, + lb, + ub, + spec +) +``` + +Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem. + +```@example RoA +prob = discretize(pde_system, discretization) + +using Optimization, OptimizationOptimisers, OptimizationOptimJL + +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) +prob = Optimization.remake(prob, u0 = res.u) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) + +net = discretization.phi +θ = res.u.depvar +``` + +We can use the result of the optimization problem to build the Lyapunov candidate as a Julia function, then sample on a finer grid than we trained on to find the estimated region of attraction. + +```@example RoA +V, V̇ = get_numerical_lyapunov_function( + net, + θ, + structure, + f, + fixed_point +) + +# Sample +states = lb[]:0.001:ub[] +V_samples = vec(V(states')) +V̇_samples = vec(V̇(states')) + +# Calculated RoA estimate +ρ = decrease_condition.ρ +RoA_states = states[vec(V(transpose(states))) .≤ ρ] +RoA = (first(RoA_states), last(RoA_states)) + +# Print statistics +println("V(0.,0.) = ", V(fixed_point)) +println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]") +println( + "V̇ ∋ [", + minimum(V̇_samples), + ", ", + max(V̇(fixed_point), maximum(V̇_samples)), + "]", +) +println("True region of attraction: (-1, 1)") +println("Estimated region of attraction: ", RoA) +``` + +The estimated region of attraction is within the true region of attraction. + +```@example RoA +using Plots + +p_V = plot(states, V_samples, label = "V", xlabel = "x", linewidth=2); +p_V = hline!([ρ], label = "V = ρ", legend = :top); +p_V = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/); +p_V = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green); + +p_V̇ = plot(states, V̇_samples, label = "dV/dt", xlabel = "x", linewidth=2); +p_V̇ = hline!([0.0], label = "dV/dt = 0", legend = :top); +p_V̇ = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/); +p_V̇ = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green); + +plt = plot(p_V, p_V̇) +``` \ No newline at end of file diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index e1d317a..805b338 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -12,8 +12,9 @@ println("Region of Attraction Estimation") ######################### Define dynamics and domain ########################## f(x, p, t) = -x .+ x .^ 3 -lb = [-2]; -ub = [2]; +lb = [-2.0]; +ub = [2.0]; +fixed_point = [0.0]; ####################### Specify neural Lyapunov problem ####################### @@ -27,7 +28,7 @@ chain = [Lux.Chain( Dense(dim_hidden, 1, use_bias = false) ) for _ in 1:dim_output] -# Define neural network discretization +# Define training strategy strategy = GridTraining(0.1) discretization = PhysicsInformedNN(chain, strategy) @@ -66,48 +67,48 @@ prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### -V_func, V̇_func = get_numerical_lyapunov_function( +V, V̇ = get_numerical_lyapunov_function( discretization.phi, res.u.depvar, structure, f, - zeros(length(lb)) + fixed_point ) ################################## Simulate ################################### -states = first(lb):0.001:first(ub) -V_predict = vec(V_func(states')) -dVdt_predict = vec(V̇_func(states')) +states = lb[]:0.001:ub[] +V_samples = vec(V(states')) +V̇_samples = vec(V̇(states')) # Calculated RoA estimate ρ = decrease_condition.ρ -RoA_states = states[vec(V_func(transpose(states))) .≤ ρ] +RoA_states = states[vec(V(transpose(states))) .≤ ρ] RoA = (first(RoA_states), last(RoA_states)) #################################### Tests #################################### # Network structure should enforce positive definiteness of V -@test min(V_func([0.0]), minimum(V_predict)) ≥ 0.0 -@test V_func([0.0]) == 0.0 +@test min(V(fixed_point), minimum(V_samples)) ≥ 0.0 +@test V(fixed_point) == 0.0 # Dynamics should result in a fixed point at the origin -@test V̇_func([0.0]) == 0.0 +@test V̇(fixed_point) == 0.0 # V̇ should be negative everywhere in the region of attraction except the fixed point -@test all(V̇_func(transpose(RoA_states[RoA_states .!= 0.0])) .< 0) +@test all(V̇(transpose(RoA_states[ RoA_states .!= fixed_point[] ])) .< 0) # The estimated region of attraction should be a subset of the real region of attraction @test first(RoA) ≥ -1.0 && last(RoA) ≤ 1.0 #= # Print statistics -println("V(0.,0.) = ", V_func([0.0])) -println("V ∋ [", min(V_func([0.0]), minimum(V_predict)), ", ", maximum(V_predict), "]") +println("V(0.,0.) = ", V(fixed_point)) +println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]") println( "V̇ ∋ [", - minimum(dVdt_predict), + minimum(V̇_samples), ", ", - max(V̇_func([0.0]), maximum(dVdt_predict)), + max(V̇(fixed_point), maximum(V̇_samples)), "]", ) println("True region of attraction: (-1, 1)") @@ -116,12 +117,12 @@ println("Estimated region of attraction: ", RoA) # Plot results using Plots -p_V = plot(states, V_predict, label = "V", xlabel = "x", linewidth=2); +p_V = plot(states, V_samples, label = "V", xlabel = "x", linewidth=2); p_V = hline!([ρ], label = "V = ρ", legend = :top); p_V = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/); p_V = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green); -p_V̇ = plot(states, dVdt_predict, label = "dV/dt", xlabel = "x", linewidth=2); +p_V̇ = plot(states, V̇_samples, label = "dV/dt", xlabel = "x", linewidth=2); p_V̇ = hline!([0.0], label = "dV/dt = 0", legend = :top); p_V̇ = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/); p_V̇ = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green); From 8f40a912a48462cd192b8131f6c2a27b64532840 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 16:32:47 -0400 Subject: [PATCH 21/34] Added link to `make_RoA_aware` docs in RoA estimation demo --- docs/src/demos/roa_estimation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md index fd547ed..e5d2347 100644 --- a/docs/src/demos/roa_estimation.md +++ b/docs/src/demos/roa_estimation.md @@ -148,7 +148,7 @@ V(x) = \left( 1 + \lVert \phi(x) \rVert^2 \right) \log \left( 1 + \lVert x \rVer which structurally enforces positive definiteness. We therefore use [`DontCheckNonnegativity()`](@ref). -We only require asymptotic decrease in this example, but we use `make_RoA_aware` to only penalize positive values of ``\dot{V}(x)`` when ``V(x) \le 1``. +We only require asymptotic decrease in this example, but we use [`make_RoA_aware`](@ref) to only penalize positive values of ``\dot{V}(x)`` when ``V(x) \le 1``. ```@example RoA using NeuralLyapunov From 44fed18e911ea478afc42a6c2be820497101c382 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 16:33:23 -0400 Subject: [PATCH 22/34] Fixed typo --- docs/src/demos/roa_estimation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md index e5d2347..4b67074 100644 --- a/docs/src/demos/roa_estimation.md +++ b/docs/src/demos/roa_estimation.md @@ -207,7 +207,7 @@ states = lb[]:0.001:ub[] V_samples = vec(V(states')) V̇_samples = vec(V̇(states')) -# Calculated RoA estimate +# Calculate RoA estimate ρ = decrease_condition.ρ RoA_states = states[vec(V(transpose(states))) .≤ ρ] RoA = (first(RoA_states), last(RoA_states)) From af9edf71279333b4f2cbcabb2cf80cdff8515076 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 17:40:50 -0400 Subject: [PATCH 23/34] Cleaned up test import statements --- test/damped_pendulum.jl | 5 ++--- test/damped_sho.jl | 4 +--- test/inverted_pendulum.jl | 5 ++--- test/inverted_pendulum_ODESystem.jl | 5 ++--- test/roa_estimation.jl | 4 +--- 5 files changed, 8 insertions(+), 15 deletions(-) diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl index df60974..b4c59f3 100644 --- a/test/damped_pendulum.jl +++ b/test/damped_pendulum.jl @@ -1,6 +1,5 @@ -using LinearAlgebra using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt +using Optimization, OptimizationOptimisers, OptimizationOptimJL using NeuralLyapunov using Random using Test @@ -99,7 +98,7 @@ prob = discretize(pde_system, discretization) res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, BFGS(); maxiters = 300) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### diff --git a/test/damped_sho.jl b/test/damped_sho.jl index c57ed95..2e3a3a8 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -1,7 +1,5 @@ -using LinearAlgebra -using NeuralPDE, Lux, ModelingToolkit +using NeuralPDE, Lux, NeuralLyapunov using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov using Random using Test diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl index ef9e274..06ba4f4 100644 --- a/test/inverted_pendulum.jl +++ b/test/inverted_pendulum.jl @@ -1,6 +1,5 @@ -using LinearAlgebra using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt +using Optimization, OptimizationOptimisers, OptimizationOptimJL using NeuralLyapunov using Random using Test @@ -94,7 +93,7 @@ prob = discretize(pde_system, discretization) res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, BFGS(); maxiters = 300) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index 03f02bb..c893fa5 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -1,6 +1,5 @@ -using LinearAlgebra using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt +using Optimization, OptimizationOptimisers, OptimizationOptimJL using NeuralLyapunov using Random using Test @@ -102,7 +101,7 @@ prob = discretize(pde_system, discretization) res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, BFGS(); maxiters = 300) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index 805b338..24d2211 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -1,7 +1,5 @@ -using LinearAlgebra -using NeuralPDE, Lux, ModelingToolkit +using NeuralPDE, Lux, NeuralLyapunov using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov using Random using Test From 6660d34afad7ffba764d19a42739dce05d7d5a80 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 17:46:04 -0400 Subject: [PATCH 24/34] Added setup to damped SHO and RoA estimation demos; cleaned up imports --- docs/src/demos/damped_SHO.md | 12 ++++++++---- docs/src/demos/roa_estimation.md | 12 ++++++++---- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md index 7f3ac41..2456914 100644 --- a/docs/src/demos/damped_SHO.md +++ b/docs/src/demos/damped_SHO.md @@ -1,4 +1,4 @@ -# Damped Simple Harmonic Oscillator +# Damped simple harmonic oscillator Let's train a neural network to prove the exponential stability of the damped simple harmonic oscillator (SHO). @@ -16,10 +16,8 @@ We'll consider just the box domain ``x \in [-5, 5], v \in [-2, 2]``. ## Copy-Pastable Code ```julia -using LinearAlgebra -using NeuralPDE, Lux, ModelingToolkit +using NeuralPDE, Lux, NeuralLyapunov using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov using Random Random.seed!(200) @@ -112,6 +110,12 @@ V, V̇ = get_numerical_lyapunov_function( In this example, we set the dynamics up as an `ODEFunction` and use a `SciMLBase.SymbolCache` to tell the ultimate `PDESystem` what to call our state and parameter variables. +```@setup SHO +using Random + +Random.seed!(200) +``` + ```@example SHO using NeuralPDE # for ODEFunction and SciMLBase.SymbolCache diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md index 4b67074..87dee7e 100644 --- a/docs/src/demos/roa_estimation.md +++ b/docs/src/demos/roa_estimation.md @@ -13,13 +13,11 @@ We'll train in the larger domain ``x \in [-2, 2]``. ## Copy-Pastable Code ```julia -using LinearAlgebra -using NeuralPDE, Lux, ModelingToolkit +using NeuralPDE, Lux, NeuralLyapunov using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov using Random -Random.seed!(250) +Random.seed!(200) ######################### Define dynamics and domain ########################## @@ -105,6 +103,12 @@ RoA = (first(RoA_states), last(RoA_states)) In this example, we set up the dynamics as a Julia function and don't bother specifying the symbols for the variables (so ``x`` will be called the default `state1` in the PDESystem). +```@setup RoA +using Random + +Random.seed!(200) +``` + ```@example RoA f(x, p, t) = -x .+ x .^ 3 lb = [-2.0]; From 895ce789f2ca18928ada2516b5c308a916c6f585 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 18:12:41 -0400 Subject: [PATCH 25/34] Removed symbolic discretization from demos --- docs/src/demos/damped_SHO.md | 1 - docs/src/demos/roa_estimation.md | 1 - 2 files changed, 2 deletions(-) diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md index 2456914..c568c7b 100644 --- a/docs/src/demos/damped_SHO.md +++ b/docs/src/demos/damped_SHO.md @@ -84,7 +84,6 @@ spec = NeuralLyapunovSpecification( ######################## Construct OptimizationProblem ######################## prob = discretize(pde_system, discretization) -sym_prob = symbolic_discretize(pde_system, discretization) ########################## Solve OptimizationProblem ########################## diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md index 87dee7e..5b8a82b 100644 --- a/docs/src/demos/roa_estimation.md +++ b/docs/src/demos/roa_estimation.md @@ -68,7 +68,6 @@ spec = NeuralLyapunovSpecification( ######################## Construct OptimizationProblem ######################## prob = discretize(pde_system, discretization) -sym_prob = symbolic_discretize(pde_system, discretization) ########################## Solve OptimizationProblem ########################## From 07ea72761c80904ce3a9928dc975f1ab4dc66d4f Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Thu, 6 Jun 2024 18:32:58 -0400 Subject: [PATCH 26/34] Replaced `GridTraining` with `QuasiRandomTraining` in inverted pendulum ODESystem test case --- test/inverted_pendulum_ODESystem.jl | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index c893fa5..5d46867 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -1,6 +1,5 @@ -using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov +using NeuralPDE, Lux, ModelingToolkit, NeuralLyapunov +import Optimization, OptimizationOptimisers, OptimizationOptimJL using Random using Test @@ -43,9 +42,9 @@ p = [defaults[param] for param in p_order] # Define neural network discretization # We use an input layer that is periodic with period 2π with respect to θ dim_state = length(bounds) -dim_hidden = 15 -dim_phi = 2 -dim_u = 1 +dim_hidden = 20 +dim_phi = 3 +dim_u = 2 dim_output = dim_phi + dim_u chain = [Lux.Chain( PeriodicEmbedding([1], [2π]), @@ -55,7 +54,7 @@ chain = [Lux.Chain( ) for _ in 1:dim_output] # Define neural network discretization -strategy = GridTraining(0.1) +strategy = QuasiRandomTraining(12_630)# GridTraining(0.1) discretization = PhysicsInformedNN(chain, strategy) # Define neural Lyapunov structure @@ -105,16 +104,19 @@ res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### +net = discretization.phi +_θ = res.u.depvar + V_func, V̇_func = get_numerical_lyapunov_function( - discretization.phi, - res.u.depvar, + net, + _θ, structure, open_loop_pendulum_dynamics, upright_equilibrium; p = p ) -u = get_policy(discretization.phi, res.u.depvar, dim_output, dim_u) +u = get_policy(net, _θ, dim_output, dim_u) ################################## Simulate ################################### @@ -139,17 +141,17 @@ x0 = (ub .- lb) .* rand(2, 100) .+ lb # Training should result in a fixed point at the upright equilibrium @test all(isapprox.( open_loop_pendulum_dynamics(upright_equilibrium, u(upright_equilibrium), p, 0.0), - 0.0; atol = 1e-8)) + 0.0; atol = 1e-5)) @test V̇_func(upright_equilibrium) == 0.0 # V̇ should be negative almost everywhere -@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 1.04e-3 +@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 3e-3 ################################## Simulate ################################### using DifferentialEquations -state_order = map(st -> istree(st) ? operation(st) : st, state_order) +state_order = map(st -> SymbolicUtils.iscall(st) ? operation(st) : st, state_order) state_syms = Symbol.(state_order) closed_loop_dynamics = ODEFunction( From 3116eee0b8a0af84ec8471a87cb8ba891ce39c0b Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Mon, 10 Jun 2024 14:17:00 -0400 Subject: [PATCH 27/34] Sped up inverted pendulum ODESystem test --- test/inverted_pendulum_ODESystem.jl | 52 ++++++++++++++--------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index 5d46867..9b09860 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -44,7 +44,7 @@ p = [defaults[param] for param in p_order] dim_state = length(bounds) dim_hidden = 20 dim_phi = 3 -dim_u = 2 +dim_u = 1 dim_output = dim_phi + dim_u chain = [Lux.Chain( PeriodicEmbedding([1], [2π]), @@ -54,7 +54,7 @@ chain = [Lux.Chain( ) for _ in 1:dim_output] # Define neural network discretization -strategy = QuasiRandomTraining(12_630)# GridTraining(0.1) +strategy = QuasiRandomTraining(1_250)# GridTraining(0.1) discretization = PhysicsInformedNN(chain, strategy) # Define neural Lyapunov structure @@ -98,7 +98,7 @@ prob = discretize(pde_system, discretization) ########################## Solve OptimizationProblem ########################## -res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 400) prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) @@ -125,14 +125,14 @@ ub = [2π, 10.0]; xs = (-2π):0.1:(2π) ys = lb[2]:0.1:ub[2] states = Iterators.map(collect, Iterators.product(xs, ys)) -V_predict = vec(V_func(hcat(states...))) -dVdt_predict = vec(V̇_func(hcat(states...))) +V_samples = vec(V_func(hcat(states...))) +V̇_samples = vec(V̇_func(hcat(states...))) #################################### Tests #################################### # Network structure should enforce positive definiteness @test V_func(upright_equilibrium) == 0.0 -@test min(V_func(upright_equilibrium), minimum(V_predict)) ≥ 0.0 +@test min(V_func(upright_equilibrium), minimum(V_samples)) ≥ 0.0 # Network structure should enforce periodicity in θ x0 = (ub .- lb) .* rand(2, 100) .+ lb @@ -141,11 +141,11 @@ x0 = (ub .- lb) .* rand(2, 100) .+ lb # Training should result in a fixed point at the upright equilibrium @test all(isapprox.( open_loop_pendulum_dynamics(upright_equilibrium, u(upright_equilibrium), p, 0.0), - 0.0; atol = 1e-5)) + 0.0; atol = 1e-3)) @test V̇_func(upright_equilibrium) == 0.0 # V̇ should be negative almost everywhere -@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 3e-3 +@test sum(V̇_samples .> 0) / length(V_samples) < 0.01 ################################## Simulate ################################### @@ -159,24 +159,24 @@ closed_loop_dynamics = ODEFunction( sys = SciMLBase.SymbolCache(state_syms, Symbol.(p_order)) ) -# Starting still at bottom +# Starting still at bottom ... downward_equilibrium = zeros(2) -ode_prob = ODEProblem(closed_loop_dynamics, downward_equilibrium, [0.0, 70.0], p) +ode_prob = ODEProblem(closed_loop_dynamics, downward_equilibrium, [0.0, 120.0], p) sol = solve(ode_prob, Tsit5()) # plot(sol) -# Should make it to the top +# ...the system should make it to the top θ_end, ω_end = sol.u[end] x_end, y_end = sin(θ_end), -cos(θ_end) @test all(isapprox.([x_end, y_end, ω_end], [0.0, 1.0, 0.0]; atol = 1e-3)) -# Starting at a random point +# Starting at a random point ... x0 = lb .+ rand(2) .* (ub .- lb) -ode_prob = ODEProblem(closed_loop_dynamics, x0, [0.0, 100.0], p) +ode_prob = ODEProblem(closed_loop_dynamics, x0, [0.0, 150.0], p) sol = solve(ode_prob, Tsit5()) # plot(sol) -# Should make it to the top +# ...the system should make it to the top θ_end, ω_end = sol.u[end] x_end, y_end = sin(θ_end), -cos(θ_end) @test all(isapprox.([x_end, y_end, ω_end], [0.0, 1.0, 0.0]; atol = 1e-3)) @@ -191,16 +191,16 @@ println( println( "V ∋ [", min(V_func(upright_equilibrium), - minimum(V_predict)), + minimum(V_samples)), ", ", - maximum(V_predict), + maximum(V_samples), "]" ) println( "V̇ ∋ [", - minimum(dVdt_predict), + minimum(V̇_samples), ", ", - max(V̇_func(upright_equilibrium), maximum(dVdt_predict)), + max(V̇_func(upright_equilibrium), maximum(V̇_samples)), "]" ) @@ -210,7 +210,7 @@ using Plots p1 = plot( xs / pi, ys, - V_predict, + V_samples, linetype = :contourf, title = "V", @@ -219,13 +219,13 @@ p1 = plot( c = :bone_1 ); p1 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0], - label = "Downward Equilibria", color = :green, markershape = :+); + label = "Downward Equilibria", color = :red, markershape = :x); p1 = scatter!( - [-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :red, markershape = :x); + [-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green, markershape = :+); p2 = plot( xs / pi, ys, - dVdt_predict, + V̇_samples, linetype = :contourf, title = "dV/dt", xlabel = "θ/π", @@ -233,13 +233,13 @@ p2 = plot( c = :binary ); p2 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0], - label = "Downward Equilibria", color = :green, markershape = :+); -p2 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria", - color = :red, markershape = :x, legend = false); + label = "Downward Equilibria", color = :red, markershape = :x); +p2 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green, + markershape = :+, legend = false); p3 = plot( xs / pi, ys, - dVdt_predict .< 0, + V̇_samples .< 0, linetype = :contourf, title = "dV/dt < 0", xlabel = "θ/π", From a1135d2d915bee994906ef7aeb905917fd37c250 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Mon, 10 Jun 2024 16:51:27 -0400 Subject: [PATCH 28/34] Added policy search demo to docs --- docs/Project.toml | 1 + docs/make.jl | 3 +- docs/src/demos/policy_search.md | 426 ++++++++++++++++++++++++++++ docs/src/demos/roa_estimation.md | 6 +- test/inverted_pendulum_ODESystem.jl | 10 +- 5 files changed, 437 insertions(+), 9 deletions(-) create mode 100644 docs/src/demos/policy_search.md diff --git a/docs/Project.toml b/docs/Project.toml index 34a5ac0..57f966a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Lux = "b2108857-7c20-44ae-9111-449ecde12c47" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" diff --git a/docs/make.jl b/docs/make.jl index ca3812d..cf1e395 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,7 +28,8 @@ makedocs(; ], "Demonstrations" => [ "demos/damped_SHO.md", - "demos/roa_estimation.md" + "demos/roa_estimation.md", + "demos/policy_search.md" ] ] ) diff --git a/docs/src/demos/policy_search.md b/docs/src/demos/policy_search.md new file mode 100644 index 0000000..ef2aea3 --- /dev/null +++ b/docs/src/demos/policy_search.md @@ -0,0 +1,426 @@ +# Policy search on the driven inverted pendulum + +In this demonstration, we'll search for a neural network policy to stabilize the upright equilibrium of the inverted pendulum. + +The governing differential equation for the driven pendulum is: +```math +\frac{d^2 \theta}{dt^2} + 2 \zeta \omega_0 \frac{d \theta}{dt} + \omega_0^2 \sin(\theta) = \tau, +``` +where ``\theta`` is the counterclockwise angle from the downward equilibrium, ``\zeta`` and ``\omega_0`` are system parameters, and ``\tau`` is our control input—the torque. + +We'll jointly train a neural controller ``\tau = u \left( \theta, \frac{d\theta}{dt} \right)`` and neural Lyapunov function ``V`` which will certify the stability of the closed-loop system. + +## Copy-Pastable Code + +```julia +using NeuralPDE, Lux, ModelingToolkit, NeuralLyapunov +import Optimization, OptimizationOptimisers, OptimizationOptimJL +using Random + +Random.seed!(200) + +######################### Define dynamics and domain ########################## + +@parameters ζ ω_0 +defaults = Dict([ζ => 0.5, ω_0 => 1.0]) + +@variables t θ(t) τ(t) [input = true] +Dt = Differential(t) +DDt = Dt^2 + +eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ τ] + +@named driven_pendulum = ODESystem( + eqs, + t, + [θ, τ], + [ζ, ω_0]; + defaults = defaults +) + +bounds = [ + θ ∈ (0, 2π), + Dt(θ) ∈ (-10.0, 10.0) +] + +upright_equilibrium = [π, 0.0] + +####################### Specify neural Lyapunov problem ####################### + +# Define neural network discretization +# We use an input layer that is periodic with period 2π with respect to θ +dim_state = length(bounds) +dim_hidden = 20 +dim_phi = 3 +dim_u = 1 +dim_output = dim_phi + dim_u +chain = [Lux.Chain( + PeriodicEmbedding([1], [2π]), + Dense(3, dim_hidden, tanh), + Dense(dim_hidden, dim_hidden, tanh), + Dense(dim_hidden, 1) + ) for _ in 1:dim_output] + +# Define neural network discretization +strategy = QuasiRandomTraining(1_250) +discretization = PhysicsInformedNN(chain, strategy) + +# Define neural Lyapunov structure +structure = PositiveSemiDefiniteStructure( + dim_phi; + pos_def = function (state, fixed_point) + θ, ω = state + θ_eq, ω_eq = fixed_point + log(1.0 + (sin(θ) - sin(θ_eq))^2 + (cos(θ) - cos(θ_eq))^2 + (ω - ω_eq)^2) + end +) +structure = add_policy_search( + structure, + dim_u +) +minimization_condition = DontCheckNonnegativity(check_fixed_point = false) + +# Define Lyapunov decrease condition +decrease_condition = AsymptoticDecrease(strict = true) + +# Construct neural Lyapunov specification +spec = NeuralLyapunovSpecification( + structure, + minimization_condition, + decrease_condition +) + +############################# Construct PDESystem ############################# + +@named pde_system = NeuralLyapunovPDESystem( + driven_pendulum, + bounds, + spec; + fixed_point = upright_equilibrium +) + +######################## Construct OptimizationProblem ######################## + +prob = discretize(pde_system, discretization) + +########################## Solve OptimizationProblem ########################## + +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 400) +prob = Optimization.remake(prob, u0 = res.u) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) + +###################### Get numerical numerical functions ###################### + +net = discretization.phi +_θ = res.u.depvar + +(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function( + driven_pendulum; simplify = true) +p = [defaults[param] for param in p_order] + +V_func, V̇_func = get_numerical_lyapunov_function( + net, + _θ, + structure, + open_loop_pendulum_dynamics, + upright_equilibrium; + p = p +) + +u = get_policy(net, _θ, dim_output, dim_u) +``` + +## Detailed description + +```@setup policy_search +using Random + +Random.seed!(200) +``` + +In this example, we'll set up the dynamics using a ModelingToolkit `ODESystem`. +Since the torque ``\tau`` is our control input, we use the `[input = true]` flag for it. + +Since the angle ``\theta`` is periodic with period ``2\pi``, our box domain will be one period in ``\theta`` and an interval in ``\frac{d\theta}{dt}``. + +```@example policy_search +using ModelingToolkit + +@parameters ζ ω_0 +defaults = Dict([ζ => 0.5, ω_0 => 1.0]) + +@variables t θ(t) τ(t) [input = true] +Dt = Differential(t) +DDt = Dt^2 + +eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ τ] + +bounds = [ + θ ∈ (0, 2π), + Dt(θ) ∈ (-10.0, 10.0) +] + +upright_equilibrium = [π, 0.0] + +@named driven_pendulum = ODESystem( + eqs, + t, + [θ, τ], + [ζ, ω_0]; + defaults = defaults +) +``` + +We'll use an architecture that's ``2\pi``-periodic in ``\theta`` so that we can train on just one period of ``\theta`` and don't need to add any periodic boundary conditions. +To achieve that, we use `Lux.PeriodicEmbedding([1], [2pi])`, enforces `2pi`-periodicity in input number `1`. +Additionally, we include output dimensions for both the neural Lyapunov function and the neural controller. + +Other than that, setting up the neural network using Lux and NeuralPDE training strategy is no different from any other physics-informed neural network problem. +For more on that aspect, see the [NeuralPDE documentation](https://docs.sciml.ai/NeuralPDE/stable/). + +```@example policy_search +using Lux + +# Define neural network discretization +# We use an input layer that is periodic with period 2π with respect to θ +dim_state = length(bounds) +dim_hidden = 20 +dim_phi = 3 +dim_u = 1 +dim_output = dim_phi + dim_u +chain = [Lux.Chain( + PeriodicEmbedding([1], [2π]), + Dense(3, dim_hidden, tanh), + Dense(dim_hidden, dim_hidden, tanh), + Dense(dim_hidden, 1) + ) for _ in 1:dim_output] +``` + +```@example policy_search +using NeuralPDE + +# Define neural network discretization +strategy = QuasiRandomTraining(1250) +discretization = PhysicsInformedNN(chain, strategy) +``` + +We now define our Lyapunov candidate structure along with the form of the Lyapunov conditions we'll be using. + +The default Lyapunov candidate from [`PositiveSemiDefiniteStructure`](@ref) is: +```math +V(x) = \left( 1 + \lVert \phi(x) \rVert^2 \right) \log \left( 1 + \lVert x \rVert^2 \right), +``` +which structurally enforces positive definiteness. +We'll modify the second factor to be ``2\pi``-periodic in ``\theta``: + +```@example policy_search +using NeuralLyapunov + +# Define neural Lyapunov structure +structure = PositiveSemiDefiniteStructure( + dim_phi; + pos_def = function (state, fixed_point) + θ, ω = state + θ_eq, ω_eq = fixed_point + log(1.0 + (sin(θ) - sin(θ_eq))^2 + (cos(θ) - cos(θ_eq))^2 + (ω - ω_eq)^2) + end +) +``` + +In addition to representing the neural Lyapunov function, our neural network must also represent the controller. +For this, we use the [`add_policy_search`](@ref) function, which tells NeuralLyapunov to expect dynamics with a control input and to treat the last `dim_u` dimensions of the neural network as the output of our controller. + +```@example policy_search +structure = add_policy_search( + structure, + dim_u +) +``` + +Since our Lyapunov candidate structurally enforces positive definiteness, we use [`DontCheckNonnegativity`](@ref). + +```@example policy_search +minimization_condition = DontCheckNonnegativity(check_fixed_point = false) + +# Define Lyapunov decrease condition +decrease_condition = AsymptoticDecrease(strict = true) + +# Construct neural Lyapunov specification +spec = NeuralLyapunovSpecification( + structure, + minimization_condition, + decrease_condition +) + +# Construct PDESystem +@named pde_system = NeuralLyapunovPDESystem( + driven_pendulum, + bounds, + spec; + fixed_point = upright_equilibrium +) +``` + +Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem. + +```@example policy_search +prob = discretize(pde_system, discretization) + +import Optimization, OptimizationOptimisers, OptimizationOptimJL + +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 400) +prob = Optimization.remake(prob, u0 = res.u) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) + +net = discretization.phi +_θ = res.u.depvar +``` + +We can use the result of the optimization problem to build the Lyapunov candidate as a Julia function, as well as extract our controller, using the [`get_policy`](@ref) function. + +```@example policy_search +(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function(driven_pendulum; simplify = true) +p = [defaults[param] for param in p_order] + +V_func, V̇_func = get_numerical_lyapunov_function( + net, + _θ, + structure, + open_loop_pendulum_dynamics, + upright_equilibrium; + p = p +) + +u = get_policy(net, _θ, dim_output, dim_u) +``` + +Now, let's evaluate our controller. +First, we'll get the usual summary statistics on the Lyapunov function and plot ``V``, ``\dot{V}``, and the violations of the decrease condition. + +```@example policy_search +lb = [0.0, -10.0]; +ub = [2π, 10.0]; +xs = (-2π):0.1:(2π) +ys = lb[2]:0.1:ub[2] +states = Iterators.map(collect, Iterators.product(xs, ys)) +V_samples = vec(V_func(hcat(states...))) +V̇_samples = vec(V̇_func(hcat(states...))) + +# Print statistics +println("V(π, 0) = ", V_func(upright_equilibrium)) +println( + "f([π, 0], u([π, 0])) = ", + open_loop_pendulum_dynamics(upright_equilibrium, u(upright_equilibrium), p, 0.0) +) +println( + "V ∋ [", + min(V_func(upright_equilibrium), + minimum(V_samples)), + ", ", + maximum(V_samples), + "]" +) +println( + "V̇ ∋ [", + minimum(V̇_samples), + ", ", + max(V̇_func(upright_equilibrium), maximum(V̇_samples)), + "]" +) +``` + +```@example policy_search +using Plots + +p1 = plot( + xs / pi, + ys, + V_samples, + linetype = + :contourf, + title = "V", + xlabel = "θ/π", + ylabel = "ω", + c = :bone_1 +); +p1 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0], + label = "Downward Equilibria", color = :red, markershape = :x); +p1 = scatter!( + [-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green, markershape = :+); +p2 = plot( + xs / pi, + ys, + V̇_samples, + linetype = :contourf, + title = "dV/dt", + xlabel = "θ/π", + ylabel = "ω", + c = :binary +); +p2 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0], + label = "Downward Equilibria", color = :red, markershape = :x); +p2 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green, + markershape = :+, legend = false); +p3 = plot( + xs / pi, + ys, + V̇_samples .< 0, + linetype = :contourf, + title = "dV/dt < 0", + xlabel = "θ/π", + ylabel = "ω", + colorbar = false, + linewidth = 0 +); +p3 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0], + label = "Downward Equilibria", color = :green, markershape = :+); +p3 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria", + color = :red, markershape = :x, legend = false); +plot(p1, p2, p3) +``` + +Now, let's simulate the closed-loop dynamics to verify that the controller can get our system to the upward equilibrium. + +First, we'll start at the downward equilibrium: + +```@example policy_search +state_order = map(st -> SymbolicUtils.iscall(st) ? operation(st) : st, state_order) +state_syms = Symbol.(state_order) + +closed_loop_dynamics = ODEFunction( + (x, p, t) -> open_loop_pendulum_dynamics(x, u(x), p, t); + sys = SciMLBase.SymbolCache(state_syms, Symbol.(p_order)) +) + +using DifferentialEquations + +# Starting still at bottom ... +downward_equilibrium = zeros(2) +ode_prob = ODEProblem(closed_loop_dynamics, downward_equilibrium, [0.0, 120.0], p) +sol = solve(ode_prob, Tsit5()) +plot(sol) +``` + +```@example policy_search +# ...the system should make it to the top +θ_end, ω_end = sol.u[end] +x_end, y_end = sin(θ_end), -cos(θ_end) +[x_end, y_end, ω_end] # Should be approximately [0.0, 1.0, 0.0] +``` + +Then, we'll start at a random state: + +```@example policy_search +# Starting at a random point ... +x0 = lb .+ rand(2) .* (ub .- lb) +ode_prob = ODEProblem(closed_loop_dynamics, x0, [0.0, 150.0], p) +sol = solve(ode_prob, Tsit5()) +plot(sol) +``` + +```@example policy_search +# ...the system should make it to the top +θ_end, ω_end = sol.u[end] +x_end, y_end = sin(θ_end), -cos(θ_end) +[x_end, y_end, ω_end] # Should be approximately [0.0, 1.0, 0.0] +``` diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md index 5b8a82b..4d8b35d 100644 --- a/docs/src/demos/roa_estimation.md +++ b/docs/src/demos/roa_estimation.md @@ -100,14 +100,14 @@ RoA = (first(RoA_states), last(RoA_states)) ## Detailed description -In this example, we set up the dynamics as a Julia function and don't bother specifying the symbols for the variables (so ``x`` will be called the default `state1` in the PDESystem). - ```@setup RoA using Random Random.seed!(200) ``` +In this example, we set up the dynamics as a Julia function and don't bother specifying the symbols for the variables (so ``x`` will be called the default `state1` in the PDESystem). + ```@example RoA f(x, p, t) = -x .+ x .^ 3 lb = [-2.0]; @@ -184,7 +184,7 @@ Now, we solve the PDESystem using NeuralPDE the same way we would any PINN probl ```@example RoA prob = discretize(pde_system, discretization) -using Optimization, OptimizationOptimisers, OptimizationOptimJL +import Optimization, OptimizationOptimisers, OptimizationOptimJL res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300) prob = Optimization.remake(prob, u0 = res.u) diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index 9b09860..dd4967d 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -31,11 +31,7 @@ bounds = [ Dt(θ) ∈ (-10.0, 10.0) ] -(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function( - driven_pendulum; simplify = true) - upright_equilibrium = [π, 0.0] -p = [defaults[param] for param in p_order] ####################### Specify neural Lyapunov problem ####################### @@ -54,7 +50,7 @@ chain = [Lux.Chain( ) for _ in 1:dim_output] # Define neural network discretization -strategy = QuasiRandomTraining(1_250)# GridTraining(0.1) +strategy = QuasiRandomTraining(1_250) discretization = PhysicsInformedNN(chain, strategy) # Define neural Lyapunov structure @@ -107,6 +103,10 @@ res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) net = discretization.phi _θ = res.u.depvar +(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function( + driven_pendulum; simplify = true) +p = [defaults[param] for param in p_order] + V_func, V̇_func = get_numerical_lyapunov_function( net, _θ, From c34be64bdc784e67dbbe6fdcdb7e75c58865dbf4 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Tue, 11 Jun 2024 10:50:18 -0400 Subject: [PATCH 29/34] Sped up damped SHO test by decreasing training iterations --- test/damped_sho.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 2e3a3a8..8a2f881 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -73,9 +73,9 @@ sym_prob = symbolic_discretize(pde_system, discretization) ########################## Solve OptimizationProblem ########################## -res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500) +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 450) prob = Optimization.remake(prob, u0 = res.u) -res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 500) +res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) ###################### Get numerical numerical functions ###################### V, V̇ = get_numerical_lyapunov_function( From 7be0f1b98970f7ceb5611ebfb5bbda422b616f59 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Tue, 11 Jun 2024 11:32:40 -0400 Subject: [PATCH 30/34] Increasing iterations in policy search test 2, since it fails on GitHub (but not locally) --- test/inverted_pendulum_ODESystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl index dd4967d..de161d3 100644 --- a/test/inverted_pendulum_ODESystem.jl +++ b/test/inverted_pendulum_ODESystem.jl @@ -94,7 +94,7 @@ prob = discretize(pde_system, discretization) ########################## Solve OptimizationProblem ########################## -res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 400) +res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 450) prob = Optimization.remake(prob, u0 = res.u) res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300) From 2cdcc666710ce5199ffc0adce310d27e48574879 Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 12 Jun 2024 13:34:58 -0400 Subject: [PATCH 31/34] Wrote documentation homepage --- docs/make.jl | 34 +++++++++++++++++--------------- docs/src/demos/damped_SHO.md | 2 +- docs/src/demos/policy_search.md | 2 +- docs/src/demos/roa_estimation.md | 2 +- docs/src/index.md | 18 ++++++++++++++--- docs/src/man/policy_search.md | 2 +- docs/src/man/structure.md | 2 +- 7 files changed, 38 insertions(+), 24 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index cf1e395..5c82102 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -4,6 +4,22 @@ using Documenter DocMeta.setdocmeta!( NeuralLyapunov, :DocTestSetup, :(using NeuralLyapunov); recursive = true) +MANUAL_PAGES = [ + "man.md", + "man/pdesystem.md", + "man/minimization.md", + "man/decrease.md", + "man/structure.md", + "man/roa.md", + "man/policy_search.md", + "man/local_lyapunov.md" +] +DEMONSTRATION_PAGES = [ + "demos/damped_SHO.md", + "demos/roa_estimation.md", + "demos/policy_search.md" +] + makedocs(; modules = [NeuralLyapunov], authors = "Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> and contributors", @@ -15,22 +31,8 @@ makedocs(; ), pages = [ "Home" => "index.md", - "Manual" => [ - "man.md", - "man/pdesystem.md", - "man/minimization.md", - "man/decrease.md", - "man/structure.md", - "man/roa.md", - "man/policy_search.md", - "man/local_lyapunov.md", - hide("man/internals.md") - ], - "Demonstrations" => [ - "demos/damped_SHO.md", - "demos/roa_estimation.md", - "demos/policy_search.md" - ] + "Manual" => vcat(MANUAL_PAGES, hide("man/internals.md")), + "Demonstrations" => DEMONSTRATION_PAGES ] ) diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md index c568c7b..46f7104 100644 --- a/docs/src/demos/damped_SHO.md +++ b/docs/src/demos/damped_SHO.md @@ -1,4 +1,4 @@ -# Damped simple harmonic oscillator +# Damped Simple Harmonic Oscillator Let's train a neural network to prove the exponential stability of the damped simple harmonic oscillator (SHO). diff --git a/docs/src/demos/policy_search.md b/docs/src/demos/policy_search.md index ef2aea3..569a6c9 100644 --- a/docs/src/demos/policy_search.md +++ b/docs/src/demos/policy_search.md @@ -1,4 +1,4 @@ -# Policy search on the driven inverted pendulum +# Policy Search on the Driven Inverted Pendulum In this demonstration, we'll search for a neural network policy to stabilize the upright equilibrium of the inverted pendulum. diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md index 4d8b35d..2bbed30 100644 --- a/docs/src/demos/roa_estimation.md +++ b/docs/src/demos/roa_estimation.md @@ -1,4 +1,4 @@ -# Estimating the region of attraction +# Estimating the Region of Attraction In this demonstration, we add awareness of the region of attraction (RoA) estimation task to our training. diff --git a/docs/src/index.md b/docs/src/index.md index ece41fa..9393401 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,12 +2,24 @@ CurrentModule = NeuralLyapunov ``` -# NeuralLyapunov +# NeuralLyapunov.jl -Documentation for [NeuralLyapunov](https://github.com/SciML/NeuralLyapunov.jl). +[NeuralLyapunov.jl](https://github.com/SciML/NeuralLyapunov.jl) is a library for searching for neural Lyapunov functions in Julia. + +This package provides an API for setting up the search for a neural Lyapunov function. Such a search can be formulated as a partial differential inequality, and this library generates a [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) PDESystem to be solved using [NeuralPDE.jl](https://github.com/SciML/NeuralPDE.jl). Since the Lyapunov conditions can be formulated in several different ways and a neural Lyapunov function can be set up in many different forms, this library presents an extensible interface for users to choose how they wish to set up the search, with useful pre-built options for common setups. + +## Getting Started + +If this is your first time using the library, start by familiarizing yourself with the [components of a neural Lyapunov problem](man.md) in NeuralLyapunov.jl. +Then, you can dive in with any of the following demonstrations (the [damped simple harmonic oscillator](demos/damped_SHO.md) is recommended to begin): ```@contents +Pages = Main.DEMONSTRATION_PAGES +Depth = 1 ``` -```@index +When you begin to write your own neural Lyapunov code, especially if you hope to define your own neural Lyapunov formulation, you may find any of the following manual pages useful: + +```@contents +Pages = Main.MANUAL_PAGES ``` diff --git a/docs/src/man/policy_search.md b/docs/src/man/policy_search.md index aaa60ce..2f8214b 100644 --- a/docs/src/man/policy_search.md +++ b/docs/src/man/policy_search.md @@ -1,4 +1,4 @@ -# Policy Search and Network-dependent Dynamics +# Policy Search and Network-Sependent Dynamics At times, we wish to model a component of the dynamics with a neural network. A common example is the policy search case, when the closed-loop dynamics include a neural network controller. diff --git a/docs/src/man/structure.md b/docs/src/man/structure.md index 688e9b1..fc071eb 100644 --- a/docs/src/man/structure.md +++ b/docs/src/man/structure.md @@ -1,4 +1,4 @@ -# Structuring a neural Lyapunov function +# Structuring a Neural Lyapunov function Three simple neural Lyapunov function structures are provided, but users can always specify a custom structure using the [`NeuralLyapunovStructure`](@ref) struct. From 8aacd7e1c8b829135ca41e3e98305b8cd068ccbd Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 12 Jun 2024 14:47:52 -0400 Subject: [PATCH 32/34] Formatting changes --- src/local_lyapunov.jl | 15 ++++++--------- test/damped_sho.jl | 8 ++++---- test/local_lyapunov.jl | 8 ++++---- test/roa_estimation.jl | 2 +- 4 files changed, 15 insertions(+), 18 deletions(-) diff --git a/src/local_lyapunov.jl b/src/local_lyapunov.jl index 25f9cac..61f48ae 100644 --- a/src/local_lyapunov.jl +++ b/src/local_lyapunov.jl @@ -19,9 +19,8 @@ If `fixed_point` is not specified, it defaults to the origin, i.e., `zeros(state Parameters `p` for the dynamics should be supplied when the dynamics depend on them. """ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, - jac::AbstractMatrix{T}; fixed_point = zeros(T, state_dim), - p = SciMLBase.NullParameters()) where T <: Number - + jac::AbstractMatrix{T}; fixed_point = zeros(T, state_dim), + p = SciMLBase.NullParameters()) where {T <: Number} model = JuMP.Model(optimizer_factory) JuMP.set_silent(model) @@ -45,8 +44,8 @@ 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) @@ -56,8 +55,7 @@ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, end function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, jac::Function; - fixed_point = zeros(state_dim), p = SciMLBase.NullParameters()) - + fixed_point = zeros(state_dim), p = SciMLBase.NullParameters()) A::AbstractMatrix = jac(fixed_point, p) return local_lyapunov( dynamics, @@ -70,8 +68,7 @@ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, jac::F end function local_lyapunov(dynamics::Function, state_dim, optimizer_factory; - fixed_point = zeros(state_dim), p = SciMLBase.NullParameters()) - + fixed_point = zeros(state_dim), p = SciMLBase.NullParameters()) A::AbstractMatrix = ForwardDiff.jacobian(x -> dynamics(x, p, 0.0), fixed_point) return local_lyapunov( dynamics, diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 8a2f881..62e50ba 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -103,10 +103,10 @@ V̇_samples = vec(V̇(hcat(states...))) V_min, i_min = findmin(V_samples) state_min = collect(states)[i_min] V_min, state_min = if V(fixed_point) ≤ V_min - V(fixed_point), fixed_point - else - V_min, state_min - end + V(fixed_point), fixed_point +else + V_min, state_min +end @test V_min ≥ 0.0 # Trained for V's minimum to be near the fixed point diff --git a/test/local_lyapunov.jl b/test/local_lyapunov.jl index 84e70f4..7a2c0aa 100644 --- a/test/local_lyapunov.jl +++ b/test/local_lyapunov.jl @@ -5,7 +5,7 @@ using NeuralLyapunov, CSDP, ForwardDiff, Test, LinearAlgebra f1(x, p, t) = -x # Set up Jacobian -J1(x, p) = (-1.0*I)(3) +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) @@ -25,7 +25,7 @@ function f2(state::AbstractVector, p, t) end function f2(states::AbstractMatrix, p, t) ζ, ω_0 = p - pos, vel = states[1,:], states[2,:] + pos, vel = states[1, :], states[2, :] vcat(transpose(vel), transpose(-2ζ * ω_0 * vel - ω_0^2 * pos)) end p2 = [3.2, 5.1] @@ -37,9 +37,9 @@ dV2dt = (state) -> ForwardDiff.derivative((δt) -> V2(state + δt * f2(state, p2 # Test V̇2 = d/dt V2 xs2 = -1:0.03:1 states2 = hcat(Iterators.map(collect, Iterators.product(xs2, xs2))...) -errs2 = abs.( V̇2(states2) .- dV2dt(states2)) +errs2 = abs.(V̇2(states2) .- dV2dt(states2)) @test all(errs2 .< 1e-10) # Test V̇2 ≺ 0 -@test all( V̇2(states2) .< 0) +@test all(V̇2(states2) .< 0) @test V̇2(zeros(2)) == 0.0 diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index 24d2211..8e27dd8 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -93,7 +93,7 @@ RoA = (first(RoA_states), last(RoA_states)) @test V̇(fixed_point) == 0.0 # V̇ should be negative everywhere in the region of attraction except the fixed point -@test all(V̇(transpose(RoA_states[ RoA_states .!= fixed_point[] ])) .< 0) +@test all(V̇(transpose(RoA_states[RoA_states .!= fixed_point[]])) .< 0) # The estimated region of attraction should be a subset of the real region of attraction @test first(RoA) ≥ -1.0 && last(RoA) ≤ 1.0 From db9e358919114e42984fe8926fbe30b2017421ba Mon Sep 17 00:00:00 2001 From: Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> Date: Wed, 12 Jun 2024 15:14:42 -0400 Subject: [PATCH 33/34] Homogenized and pruned test imports --- test/damped_pendulum.jl | 5 ++--- test/damped_sho.jl | 2 +- test/inverted_pendulum.jl | 5 ++--- test/roa_estimation.jl | 2 +- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl index b4c59f3..a61472d 100644 --- a/test/damped_pendulum.jl +++ b/test/damped_pendulum.jl @@ -1,6 +1,5 @@ -using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov +using NeuralPDE, Lux, ModelingToolkit, NeuralLyapunov +import Optimization, OptimizationOptimisers, OptimizationOptimJL using Random using Test diff --git a/test/damped_sho.jl b/test/damped_sho.jl index 62e50ba..83fb33c 100644 --- a/test/damped_sho.jl +++ b/test/damped_sho.jl @@ -1,5 +1,5 @@ using NeuralPDE, Lux, NeuralLyapunov -using Optimization, OptimizationOptimisers, OptimizationOptimJL +import Optimization, OptimizationOptimisers, OptimizationOptimJL using Random using Test diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl index 06ba4f4..419b050 100644 --- a/test/inverted_pendulum.jl +++ b/test/inverted_pendulum.jl @@ -1,6 +1,5 @@ -using NeuralPDE, Lux, ModelingToolkit -using Optimization, OptimizationOptimisers, OptimizationOptimJL -using NeuralLyapunov +using NeuralPDE, Lux, NeuralLyapunov +import Optimization, OptimizationOptimisers, OptimizationOptimJL using Random using Test diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl index 8e27dd8..e8ab570 100644 --- a/test/roa_estimation.jl +++ b/test/roa_estimation.jl @@ -1,5 +1,5 @@ using NeuralPDE, Lux, NeuralLyapunov -using Optimization, OptimizationOptimisers, OptimizationOptimJL +import Optimization, OptimizationOptimisers, OptimizationOptimJL using Random using Test From fc41720aa5fc0e9cc0e317aad739937cc32524de Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 13 Jun 2024 05:31:22 -0400 Subject: [PATCH 34/34] Create Documentation.yml Signed-off-by: Christopher Rackauckas --- .github/workflows/Documentation.yml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 .github/workflows/Documentation.yml diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml new file mode 100644 index 0000000..030c5a4 --- /dev/null +++ b/.github/workflows/Documentation.yml @@ -0,0 +1,18 @@ +name: "Documentation" + +on: + push: + branches: + - master + tags: '*' + pull_request: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ github.ref_name != github.event.repository.default_branch || github.ref != 'refs/tags/v*' }} + +jobs: + build-and-deploy-docs: + name: "Documentation" + uses: "SciML/.github/.github/workflows/documentation.yml@v1" + secrets: "inherit"