# Reasoning about Exponential Decay and Growth in KeYmaera X

--
Author: Nathan Fulton
Category: Automated Reasoning
Date Published: 01-14-2017
--

**Update: There's an exciting new paper on differential ghosts from 
[André Platzer and Yong Kiam Tan](https://lfcps.org/pub/diffaxiomatic.pdf)**

**Warning: The Interface to differential ghosts provided by KeYmaera X has changed since I wrote this post.
I maintain working examples from this blog post at [the KeYmaera X projects repository](https://github.com/LS-Lab/KeYmaeraX-projects/tree/master/ghosts_simple)
but haven't had time to update this blog post. All of the conceptual explanations remain relevant, but refer to the projects repository for up-to-date tactic scripts.**

<!--
**Application**: Reachability proofs for Ordinary Differential Equations. 
<br/>
**Abstraction**: Differential Invariants; Differential Ghosts
-->

The purpose of this post is to explain how [KeYmaera X](http://keymaerax.org) automatically proves reachability properties about
differential equations of the form `x' = f(x)`; and, in particular, equations whose solutions are exponential functions.
My intended audience are people who are already familiar with the first 5 or 6 lectures of [Foundations of Cyber-Physical Systems](http://symbolaris.com/course/fcps16-schedule.html).
This is about equivalent to the half-day KeYmaera X tutorials[^1].

[^1]: 
    In particular, I assume the reader is familiar with Differential Dynamic Logic and differential invariants/ghosts, 
    but I'll quickly review the basics. There are also a couple of appendices with incomplete explanations of various concepts.

Ultimately, the goal of this note is to explain **how to find a useful differential ghost in order to prove a property about an exponential system such as x'=x**.

## Background and Motivation: Continuous Rechability Problems in KeYmaera X

Many control problems boil down to ensuring that a continuous system can only reach a safe subset of the overall state space.

For example, consider a model of an adaptive cruise control protocol. The engineering problem is to
design a control algorithm that actuates a follower car's acceleration such that the distance between the
lead car and the follower car, denoted *rx*, is always strictly positive. In this example, the safe
states are all states where *rx > 0*.

More generally, continuous reachability problems are expressed in Differential Dynamic Logic (dL) as formulas of the form:

    P -> [{x' = f & H}]Q

where P, Q, and H are a first-order formulas of real arithmetic, and `P -> [{X' = f}]Q` is true iff
*every* flow of the system *x' = f* restricted to the set `H` and starting in the set *P* stays within the set *Q*.

Returning to the adaptive cruise control example, we can phrase our reachability problem as:

    P -> [{rx' = rv, rv' = ra & rv > 0}] rx > 0

where `rx` is relative position and `& rv > 0` is our syntax for restricting the flow of the system to sets in which *rv > 0*. There are many valid choices for `P`; for example:

    rx > 0 & rv = 12 & ra = 0

In dL, we prove properties like this one by reducing the question to one that is expressible in a decidable fragment of first-order real arithmetic.
In this case, we can integrate the system and arrive at a formula that looks something like:

    P -> \forall t (t>=0->\forall s (0<=s&s<=t->ra*s+rv>0)->ra/2*t^2+rv*t+rx>0)

The truth of this formula is then established using a decision procedure for real arithmetic.

### So why not just solve x' = x?

Because model theory.

Unfortunately, not all systems have explicit closed form solutions expressible in terms of decidable fragments of Real Arithmetic.
In particular, consider the equation `x' = x`, which has the solution `e^t`. Therefore, the arithmetic question corresponding to the reachability question 
`P -> [{x' = x}] x > 0` looks something like:

    P -> \forall t. e^t > 0

*In general*, we can't simply appeal to a arithmetic decision procedure to answer questions of this form -- 
the decidability of real arithmetic with exponentials [remains an open problem](https://en.wikipedia.org/wiki/Tarski%27s_exponential_function_problem).

In this post, I'll explain how you can prove reachability properties about equations of the form `x' = f(x)` in [KeYmaera X](https://keymaerax.org) 
while staying with fragments of arithmetic that are known to be decidable.

## Differential Ghosts

I assume as background familiarity with André Platzer's [lecture notes on differential ghosts](http://symbolaris.com/course/fcps16/12-diffghost.pdf),
or a bit of experience with dL (e.g., having attended one of our KeYmaera X tutorials).
The following review of ghosts is non-comprehensive, but might be a nice refresher.

dL's **differential ghost** axiom augments a continuous system with a new equation that isn't mentioned in the rest of the system: 

    [{ode & H}]P <-> ∃y.[{ode, y'=a*y+b & H}]P

The terms `a` and `b` may not mention `y`.
Additionally, `y` must not occur anywhere in `[{ode & H}]P`.
I.e., `y' = a*y+b` must by a linear system and `y` must be a fresh (new) variable.

**Differential auxiliaries** are a related concept that allow the post-condition to be re-stated in terms of the system's new variable:

    P <-> ∃y.G    G |- [{ode, y' = a*y+b & H}]G
    --------------------------------------------- diffAux(y, a, b, G)
            P |- [{ode & H}]P

This [proof rule](https://en.wikipedia.org/wiki/Sequent_calculus) appears a bit confusing at first.
Seeing why this rule should be *sound* is not too subtle.
`P <-> ∃y.G` establishes that `G` implies `P` for *some* choice of `y`, and `G |- [{ode, y' = a*y+b & H}]G`
establishes that `G` remains invariant for *any* choice of `y`.

Seeing why this rule is *helpful* is a bit more subtle. I could waste some ink, but instead, I'll just jump into examples!

## Ghosts for Open Sets

In this section, we'll consider various systems of the form
`x > 0 -> [{x' = f(x)}] x > 0`.

Notice that we can't prove anything about these systems using differential variants, and we can't solve the system without turning one undecidable problem into another (for now, at least).

**Example 1** `x > 0 -> [{x' = -x}] x > 0`

Before reading further, see if you can come up with a first order formula G that makes the premise of diffAux true;
i.e., find a predicate `G(x,y)` that makes `x>0 <-> ∃y.G` true.


                         ------------------------------------- R
                                |- -xy^2 + 2xy(y/2) = 0
    ----------------- R  ------------------------------------- diffInd
    x>0 <-> ∃y.xy^2=1    xy^2=1 |- [{x' = -x, y' = y/2}]xy^2=1
    ----------------------------------------------------------- diffAux(y, 1/2, 0, xy^2 = 1)
            x > 0 |- [{x' = -x}] x > 0

Hopefully you arrived at `xy^2 = 1` by yourself. 
**Notice that getting to a choice of `y' = ay + b` after fixing `G` is pretty mechanical**; we just compute Lie derivatives of our choice of G:

    (xy^2=1)' <-> (xy^2)' = 0         (def'n Lie operator =)
              <-> x'y^2 + x2yy' = 0   (def'n Lie operator *)
              <-> -xy^2 + x2yy' = 0   (because x' = -x)
              <-> y'(2xy) = xy^2
              <-> y' = xy^2 / 2xy
              <-> y' = y/2
              <-> y' = (1/2)y + 0     (just re-stating so it's clear a=1/2 and b=0)

The same proof generalizes nicely to monomials. 

**Example 2** `x > 0 -> [{x' = -x^2}] x > 0`

The choice of `y` changes slightly: 

    (xy^2=1)' <-> ... 
              <-> -x^2y^2 + x2yy' = 0
              <-> y' = x^2y^2 / 2xy
              <-> y' = (x/2)y + 0

Instead of giving the sequent calculus proof, here's the KeYmaera X tactic:

    implyR(1) ; DA4({`x*y^2=1`}, {`y`}, {`x/2`}, {`0`}, 1) ; <(
      QE,
      implyR(1) ; diffInd(1)
    )

<!-- TODO-nrf find time to add a polynomial example. Maybe include as an appendix all of the files/tactics in ~/dev/Notes/invcheck 
Handling polynomials isn't essentially different.

**Example 3** `x > 0 -> [{x' = -x^2 + 2x + 5}] x > 0`
-->


#### Some Proving Advice

At this point, you've already noticed that the choice of `G` and `y' = ay+b` are closely coupled.
In fact, you only have to be creative once -- there's a systematic way of deriving `a` and `b` from a fixed `G`.
This is true in the other direction as well.

I typically start with G, but as we'll see with equilibrium points, it's somethings helpful to move back and forth.

## Ghosts for Equilibrium Points

In some sense, you would expect the proof of `x=0 -> [{x' = -x}]x=0` to be trivial. But the proof of this property in dL
(without extra proof rules like [Khalil Ghorbal's](http://www.lix.polytechnique.fr/Labo/Khalil.Ghorbal/) [DRI](http://repository.cmu.edu/cgi/viewcontent.cgi?article=3734&context=compsci))
is a tough exercise. Before continuing, see if you can find some candidates for `G` and/or `y'`.

Here's a choice for G: `y>0 & x*y=0`. 
Obviously, `∃y.x=0 <-> y>0 & x*y=0`.
However, a simple differential induction argument doesn't suffice to establish the remaining subgoal:

    y>0 & x*y=0 -> [{x' = -x}](y>0 & x*y=0)

We'll need another ghost. 

Notice that we already know how to prove

    y>0 -> [{x' = -x}](y>0)

via a differential ghost argument. So let's split this subgoal into two cases,
which is a proof we can exploit using the `boxAnd` (`[]^`) axiom:

     (use open set approach)                            ...
    -----------------------------           ----------------------------------
    y>0 & x*y=0 |- [{x' = -x, y'=???}]y>0     y>0 & x*y=0 |- [{x' = -x}]x*y=0
    --------------------------------------------------------------------------
            y>0 & x*y=0 |- [{x' = -x, y'=???}](y>0 & x*y=0)


Because we already know how to prove the first case, it makes sense to 
choose a value for `y' = ax + b` that makes the `x*y=0` case prove by a differential invariance argument:

    (x*y=0)' <-> (x*y)'=0
             <-> x'y + y'x = 0
             <-> -xy + y'x = 0
             <-> y' = xy/x = y

The proof for the case with the `y>0` post-condition is the same as the proof for the close set examples.

## Ghosts for Closed/Clopen Sets

The easiest way of dealing with properties of the form `x>=0` is to think of them as a combination of a closed set and an equilibrium point:

    x >= 0 <-> x > 0 | x = 0

The key to this techniques is to cut `x > 0 | x = 0`, case distinguish, and then use the following proof rule to rewrite the post-condition:

     |- Q -> P    G |- [a]Q
    ------------------------ G[]
           G |- [a]P 

E.g. by taking `P as x >= 0` and `Q as x > 0`.

Here's how that proof goes (starting *after* the cut):

              *             cont'd                      *             cont'd
        --------------  --------------------       -------------  -------------------
        |- x>0 -> x>=0  |- x>0 |- [{ode}]x>0       |- x=0 -> x>0  x=0 |- [{ode}] x=0
        ------------------------------------ G[]   ----------------------------------- G[]
                x>0 |- [{ode}] x>=0                     x=0 |- [{ode}] x>=0
    -------------------------------------------------------------------------------- orL
                              x > 0 | x = 0 |- [{ode}] x>=0

For reference, here's the tactic that does that (the open goals are `skip`'d over):

    implyR(1) ; cut({`x>0 | x=0`}) ; <(
      hideL(-1) ; orL(-1) ; <(
        generalizeb({`x>0`}, 1) ; <(skip, QE),
        generalizeb({`x=0`}, 1) ; <(skip, QE)
      ),
      hideR(1) ; QE
    )

## Conclusion

Proving reachability properties about exponential functions is one of the more difficult tasks for
new KeYmaera X users. 
This may seem surprising because properties such as `x=0 -> [{x'=x}]x=0` are so intuitively simple;
however, given [what we know about real arithmetic](https://en.wikipedia.org/wiki/Tarski%27s_exponential_function_problem),
perhaps this shouldn't be so surprising after all.

<hr/>

## Appendices

### Appendix A Further Resources
 * Source code for tactics that automate [equilibrium point](https://github.com/LS-Lab/KeYmaeraX-release/tree/1c05309f318f15b5469ff068adf4949516bffeaf/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/DifferentialTactics.scala#L709)
   and [open set](https://github.com/LS-Lab/KeYmaeraX-release/tree/1c05309f318f15b5469ff068adf4949516bffeaf/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/DifferentialTactics.scala#L587) proofs.
 * Lecture videos and notes from [15-424 Foundations of Cyber-Physical Systems](http://www.cs.cmu.edu/~aplatzer/course/fcps16-schedule.html) at Carnegie Mellon.
   The lectures on differential invariants and differential ghosts are particularly relevant to this blog post.

<!--

### Appendix B Refresher on ODEs

All undergraduate courses on differential equations teach students how to *solve* a system of differential equations;
that is, how to compute a closed-form representation of the functions `x_i(t)` when given a system of differential equations
`x_1' = s_1, ..., x_n' = s_n`.
For those who haven't yet had a course on ODEs, here are a few examples:

**Example 1: Solve the system {x'=v, v'=a} where a is a constant**

    To solve x'(t) = v(t), v'(t) = a 
    N.T. v(t)  = a*t + v(0)                   by integrating v'(t) wrt t            (1)
    .:.  x'(t) = v(t) = a*t + v(0)            by (1)                                (2)
    .:.  x(t)  = a*t^2/2 + v(0)*t + x(0)      by (2) and integration x'(t) wrt t    (3)

    Therefore, x(t)  = a*t^2/2 + v(0)*t + x(0) and v(t)  = a*t + v(0) is a solution
    to the differential equations {x'=v,v'=a}.

**Example 2: Solve the system {x'=y, y'=-x} with x_0 = 0, y_0=1**

    To solve: x'(t) = y(t), t'(t) = x(t) where x(0) = 0 and y(0) = 1.
    Let
      x(t) = sin(t) and 
      y(t) = cos(t).
    Verify:
      x'(t) = sin'(t) = cos(t)  = y(t)
      y'(t) = cos'(t) = -sin(t) = -x(t)
      x(0)  = sin(0)  = 0
      y(0)  = cos(0)  = 1

    Therefore, x = sin and y = cos is a solution to the system {x'=y,y'=-x}
    with initial values x_0=0 and y_0=1.

Note that in this post we only consider linear homogeneous equations; i.e., we'll discuss equations like the ones in example 1, but not equations like the ones in example 2.
The solutions computed in example 2 are not very interesting from a proof theoretic (i.e., pragmatic) perspective because *sin* and *cos* exist only in an undecidable fragment of
real arithemtic; therefore, using the solutions instead of the ODE itself merely reduces one undcidable question into another undeciable question.

-->

### Appendix B: Axiomatizing Differential Dynamics; Lie Derivatives


**Differential invariants** are the key piece of technology that allow us to obtain 
correct-by-construction solutions without resorting to a verified implementation of a fixed-point procedure.
[KeYmaera X](https://keymaerax.org) axiomatizes continuous dynamical systems
in terms of differential invariants and Lie derivatives.


A differential invariant is just a formula that remains true throughout the entire flow of an ODE.
For example, suppose `x=1` initially.
Then:

 * `x=1` is an invariant of `{x'=0}`,
 * `x>0` is an invariant of `{x'=1}`, but
 * `x<1` is *not* an invariant of the equation `{x'=-1}` (because `x<1` is not true in the first instant of the flow).

In [KeYmaera X](https://keymaerax.org) we can ask questions about invariants of differential equations using
a formula of the form `[ODE]inv` where `ODE` is a (system of) differential equation(s) and `inv` is a formula describing the invariant.
Notice that `[ODE]inv` is a formula, so it can be used with other logical connectives.
For example, we can express the english prose "if `x=1` initially then `x>0` is an invariant of `{x'=1}`"
using the formula `x=1 -> [{x'=1}]x>0`. Notice that `x=0 -> [{x'=-1}]x<0` is also well-formed, but is false.

**Evolution domain constraints** allow us to *assume* invariants of a differential equation by restricting the flow of a differential equation to a particular domain.
For example, if we define our dynamical system as the flow of `x'=-1` restricted to the set `x>0`, then `x>0` is, by force of definition, an invariant of the system.
The [KeYmaera X](https://keymaerax.org) syntax for expressing that `x>0` is an invariant of the system `x'=-1` constrainted to `x>0` is 
`x=1 -> [{x'=-1 & x>0}]x>0`.

**Differential cuts** are a way of embedding an invariant into the defintiion of a contiuous dynamical system.
Think of cuts like lemmas -- we first establish that `I` is an invariant of our system, and then get to assume throughout the 
rest of our proof that `I` is a component of the evolution domain constraint:

    Γ |- [{ode & C}]I   Γ |- [{ode & C & I}]P
    ------------------------------------------ diffCut
            Γ |- [{ode & C}]P


**Differential induction** can be used to prove that a formula is an invariant of a differential equation:


    [{x'=t}]P <-> [x':=t;]P'    (diffInd)

where P' is the [Lie Derivative](https://en.wikipedia.org/wiki/Lie_derivative) of `P`.
The definition of `P'` is straight-forward for terms:

    (s+t)' = s' + t'
    (s*t)' = s'*t + t'*s

and so on. The definition of unquantified formulas of real arithmetic is subtle (there are no typos on the following lines):

    (f=g)'  <-> f' = g'
    (f!=g)' <-> f' = g'
    (f>g)'  <-> f' >= g'
    (f<g)'  <-> f' <= g'
    ...

Explaining differential induction is beyond the scope of this blog post, but there is an excellent [video](https://scs.hosted.panopto.com/Panopto/Podcast/Download/c49e45fa-5498-45e4-9a49-05f50ef56751.mp4)
with accompanying [lecture notes](http://symbolaris.com/course/fcps16/10-diffinv.pdf) for interested readers.

<hr/>

Change Log:

 * 2/26/2017: Fixed some typos and updated Bellerophon code samples to new `;` syntax.
 * July 2017: Added warning about tactics not working anymore. I'll update the note soon.

<hr/>