Viabilitree
This library propose a set of algorithms, based on a particular kdtree structure, in order to compute viability kernels and capture basins.
Motivation
Mathematical viability theory offers concepts and methods that are suitable to study the compatibility between a dynamical system described by a set of differential equations and constraints in the state space. The result sets built during the viability analysis can give very useful information regarding management issues in fields where it is easier to discuss constraints than objective functions. Viabilitree is a framework in which the viability sets are represented and approximated with particular kdtrees. The computation of the viability kernel is seen as an active learning problem. We prove the convergence of the algorithm and assess the approximation it produces for known problems with analytical solution. This framework aims at simplifying the declaration of the viability problem and provides useful methods to assist further use of viability sets produced by the computation.
Viability problem
Simple example
Population Growth Model
This example is taken from [4]. The population model is Malthusian. The population viability problem consists in maintaining the size of a population in a given interval $[a;b]
x(t)
y(t)
$, the population growth rate. The dynamics are described by the following equations:
The dynamics are controlled by taking the growth rate evolution in interval $[c,c]
dt
$ tends toward $0
$, the theoretical viability kernel is defined by:
[Figure 1: Viability kernel of the population viability problem][Figure 1]
[//]: # (Figure 1: Viability kernel of the population viability problem)
The above figure shows an approximation of the viability kernel for the population problem with:
 constraint set $
K=[a=0.2,b=3]\times[d=2,e=2]
$,  parameters $
dt=0.1
$,  control set $
U=[0.5;0.5]
$ with discretization step 0.02. The color stands for the value of a control $u
$ which allows the state to stay in the viability kernel. In black the boundary of the true kernel.
By definition of the viability kernel, starting from any point in the viability kernel, there exists always an evolution that stays in the viability kernel. Starting from any point outside the viability kernel, any evolution will exit the constraint set in finite time. This is why the notion of viability kernel is so useful.
The corresponding code is the following:
For the definition of the model: dynamics, perturbations, etc.
import viabilitree.model.Dynamic
case class Population(integrationStep: Double = 0.01, timeStep: Double = 0.1) {
def dynamic(state: Vector[Double], control: Vector[Double]) = {
def xDot(state: Vector[Double], t: Double) = state(1) * state(0)
def yDot(state: Vector[Double], t: Double) = control(0)
val dynamic = Dynamic(xDot, yDot)
val res = dynamic.integrate(state.toArray, integrationStep, timeStep)
res
}
timeStep stands for to
For the definition of the viability problem:
import scala.util.Random
import viabilitree.viability._
import viabilitree.export._
import viabilitree.kdtree.Tree
import viabilitree.viability.kernel._
import java.io.File
object PopulationViability extends App {
// accuracy parameter
val depth = 20
// algorithm parameter
val rng = new Random(42)
// model definition
val population = Population()
// control parameter
def c = 0.5
// constraint set parameters
def a = 0.2
def b = 3.0
def d = 2.0
def e = 2.0
// definition of the viability problem
val vk = KernelComputation(
dynamic = population.dynamic,
depth = depth,
zone = Vector((a, b), (d, e)),
controls = Vector(c to c by 0.02)
)
// computation of the viability kernel corresponding to problem vk
val (ak, steps) = approximate(vk, rng)
// save viability kernel to file (vtk format, to be processed by paraview)
val f = new File(s"population${steps}depth${depth}.vtk")
saveVTK2D(ak, f)
}
a to e are the same parameters as in the mathematical definition.
The viability problem is defined by class KernelComputation with the following parameters:

depth which defines the accuracy of the approximation. There are $
2^{depth}
$ grid points (here, $2^{\frac{depth}{2}}
$ points per axes).  dynamic: the model dynamic

zone: the area to explore and here it is also the constraint set, $
[a,b]\times[d,e]
$ 
controls: the set of admissible controls, it is the same set for each state,$
[c,c]
$
The computation itself is done by the approximate function.
Mathematical Viability Theory ([2], [3])
In Viabilitree we consider a viability problem defined by a controlled dynamical system $S
U
$ (the set of admissible controls depending on the state of the system), and a compact subset $K
$ of the state space (the set of constraints):
$x(t)
$ is the state of the system $S
x(t)\in {\mathbb R}^p
$ a finite dimensional vector space.
$u(t)
$ is the control, with $u(t)\in {\mathbb{R}}^q
$.
The setvalued map $U : {\mathbb R}^p\leadsto {\mathbb{R}}^q
$ gives the set of admissible control for each state $x
\Phi
$ is a function from $\text{Graph}(U)
$ to ${\mathbb R}^p
$.
$K\subset {\mathbb R}^p
$ is a compact subset of ${\mathbb R}^p
x(t)
$ is supposed to stay.
The viability kernel $viab_S(K)
$ is the largest subset of $K
$ (possibly empty) that gathers the states from which it is possible to find a control function $u(t)
$ such that the evolution $x(.)
$ stays in the compact set $K
$.
In Viabilitree we follow the method described in [5], we consider dynamical system ($S_{dt}
$) discretized in time:
We use the learning algorithm $L
$ of kdtree described in [1] on a discretized grid $K_h
$ to compute an approximation $L(K_h)
$ of the viability kernel $viab_{S_{dt}}(K)
$ of the discretized dynamical system ($S_{dt}
K
S_{dt}
$ verify some conditions, [1], [5] and [6] ensure that $L(K_h)
$ converges to $viab_S(K)
$ when $h
$ and $dt
$ tend to 0.
The convergence conditions are:
 $
F
$ is Marchaud and $\mu
$Lipschitz with closed images  $
K
$ is a compact set.  $
L(K^0)=K
$.  $
viab_S(K)
$ is compact, it is pathconnected  $
viab_S(K)
$ erosion with $B(\epsilon)
$ is pathconnected and points of $viab_S(K)
$ are at most distant from the eroded set by $\epsilon \sqrt{p}
$  $
viab_S(K)
$ complementary set is pathconnected as is its erosion with $B(\epsilon)
$, and its points are at most distant from the eroded set by $ \epsilon \sqrt{p}
$. The two last properties ensure that there are no small tentacles.
In practice points of the grid are removed from the current dilated approximation when following the dynamic with each available control they always leave the dilated approximation at the next step. See KernelComputation for implementation details.
More details on both theoretical and practical aspects in 7
Install
Coming soon
Use
The main purpose of Viabilitree is to help users to declare a viability problem, compute the corresponding viability kernel and possibly its capture basin, in order to perform a viability analysis.
Example of viability problems are gathered in package example. Each example is detailled in its own readme file.
 Bilingual :
 approximation of the viability kernel (BilingualViabDomain)
 approximation of the capture basin of the viability kernel (BilingualBasin)
 Consumer
 approximation of the viability kernel (ConsumerViability)
 approximation of a set (the analytical kernel ) (ConsumerKernel)
 Lake
 approximation of the viability kernel (LakeViability)
 use of output files (OutputLake)
 Population
 approximation of the viability kernel (PopulationViability)
 exploration with OpenMOLE (PopulationViability)
 approximation of a set (the analytical kernel ) (PopulationApproximation)
Set operators
 volume
 intersection of trees
 indicator function
 erosion
 dilation
Examples in circle.
Learning algorithm
Apart from viability study it is possible to use Viabilitree as a simple learning algorithm.
See approximation package for more details.
Examples in population and circle
To be completed
References
[1] Rouquier et al, A kdtree algorithm to discover the boundary of a black box hypervolume, Annals of Mathematics and Artificial Intelligence, vol 75, 3, pp "335350, 2015.
[2] Aubin, Viability theory, Birkhäuser, 1991
[3] Aubin, Bayen, A. and SaintPierre, P. Viability Theory: New Directions, Springer, 2011.
[4] Aubin et SaintPierre, "An introduction to viability theory and management of renewable resources", chapter in Advanced Methods for Decision Making and Risk Management, J. Kropp and J. Scheffran eds., Nova Science Publishers, 2007
[5] Deffuant , Chapel, Martin (2007). Approximating viability kernels with support vector machines. IEEE T. Automat. Contr., 52(5), 933–937.
[6] SaintPierre, P. (1994). Approximation of the viability kernel. Applied Mathematics & Optimisation, 29(2), 187–209.
[7] Alvarez, Reuillon, de Aldama. Viabilitree: A kdtree Framework for Viabilitybased Decision. hal01319738.
Remarks
A setvalued map $F
c > 0
$ such that $\forall x, \ \sup _{y \in F(x)} \left\y\right\ \leq c(\left\x\right\+1)
$.
A setvalued map $F
$ is $\mu
\mu>0
$ if for all $x
$ and $y
F(x)\subset F(y)+B(0,\mu xy)
$