Skip to content

Latest commit

 

History

History
1705 lines (1193 loc) · 50.7 KB

smt-solving.md

File metadata and controls

1705 lines (1193 loc) · 50.7 KB
title author date output
SAT Modulo Theories for Fun and Profit
Matt Peddie, Kittyhawk <[email protected]>
10 July 2018
slidy_presentation
highlight
pygments

Overview

  • What is

    • SAT
    • SMT
    • SBV

  • Proofs

    • Simple proofs
    • Stability proof for a feedback controller

  • Code generation

  • Constraint solving

  • Optimization

SAT

Boolean SATisfiability

For some formula involving logical variables and logical connectives, is there an assignment of boolean values to the variables that makes the formula true?

SAT

Boolean SATisfiability

For some formula involving logical variables and logical connectives, is there an assignment of boolean values to the variables that makes the formula true?

Example:

$\Huge \displaystyle{A \land \left(B \lor \left(\neg C \right)\right)}$

SAT

Boolean SATisfiability

For some formula involving logical variables and logical connectives, is there an assignment of boolean values to the variables that makes the formula true?

Example:

$\Huge \displaystyle{A \land \left(B \lor \left(\neg C \right)\right)}$

  • Variables: A, B, C . . .
  • Connectives are always $\Huge \land$ (AND), $\Huge \lor$ (OR) and $\Huge \neg$ (NOT).

SMT

SAT Modulo Theories

Take a SAT formula and replace some or all of the variables with predicates over a more complicated theory. Example:

$\Huge \left(2*x + y &lt; 0) \land (x \geq 0)$

SMT

SAT Modulo Theories

Take a SAT formula and replace some or all of the variables with predicates over a more complicated theory. Example:

$\Huge \left(2*x + y &lt; 0) \land (x \geq 0)$

"More complicated theories" include

  • algebraic real numbers
  • floating-point arithmetic
  • integer arithmetic
  • arrays
  • bitvectors
  • "uninterpreted functions"
  • strings (Z3)

SMTLIB

http://smtlib.cs.uiowa.edu/Logics/logics.png

SMTLIB

http://smtlib.cs.uiowa.edu/Logics/logics.png

Examples:

  • QF_BV is "quantifier-free formulae over the theory of fixed-size bit-vectors."

SMTLIB

http://smtlib.cs.uiowa.edu/Logics/logics.png

Examples:

  • QF_BV is "quantifier-free formulae over the theory of fixed-size bit-vectors."

  • QF_NIA is "quantifier-free formulae over the theory of nonlinear integer arithmetic."

SMTLIB

http://smtlib.cs.uiowa.edu/Logics/logics.png

Examples:

  • QF_BV is "quantifier-free formulae over the theory of fixed-size bit-vectors."

  • QF_NIA is "quantifier-free formulae over the theory of nonlinear integer arithmetic."

  • QF_FPBV is a combination of floating point (FP) and array and bit-vector (BV) theories

These specify the underlying theories from which the predicates in the boolean formula may be drawn.

SBV

SMT-Based Verification -- a Haskell DSL for working with SMT solvers.

Write your programs in Haskell, run them in SMTLIB.

Additional features:

  • Constraint solving
  • Overflow/underflow checking
  • C code generation
  • Optimization
  • Parallel solvers (race-mode or check-mode)

SBV uses Z3 by default, but it supports several others as well.

Basic SBV interface

Simple DSL:

data SBV a

a ranges over types we reason about:

  • Double
  • Bool
  • Int32
  • Word8

etc. These must be instances of a class SymWord.

type SDouble = SBV Double
type SBool   = SBV Bool
...

Basic SBV interface

Some new interfaces:

class Boolean b where
  true :: b
  bnot :: b -> b
  (&&&) :: b -> b -> b
  ...

class EqSymbolic a where
  (.==) :: a -> a -> SBool
  (./=) :: a -> a -> SBool
  ...

class Mergeable a where
  select :: (SymWord b, Num b) => [a] -> a -> SBV b -> a
  ...

ite :: Mergeable a => SBool -> a -> a -> a

Basic SBV interface

Symbolic reasoning type:

data Symbolic a

instance Functor Symbolic
instance Applicative Symbolic
instance Monad Symbolic

Introducing symbolic variables:

sDouble :: String -> Symbolic (SBV Double)
sBools :: [String] -> Symbolic [SBV Bool]

Constraints:

constrain :: SBool -> Symbolic ()

Proof of associative property of addition

additionAssoc x y z = (x + y) + z .== x + (y + z)

Proof of associative property of addition

additionAssoc x y z = (x + y) + z .== x + (y + z)

prove $ do
  [x, y, z] <- sDoubles (pure <$> "xyz")
  pure $ additionAssoc x y z

Proof of associative property of addition

additionAssoc x y z = (x + y) + z .== x + (y + z)

prove $ do
  [x, y, z] <- sDoubles (pure <$> "xyz")
  pure $ additionAssoc x y z

Falsifiable. Counter-example:
  x =               Infinity :: Double
  y =              -Infinity :: Double
  z = -1.78607059478183e-310 :: Double

Proof of associative property of addition (take 2)

prove $ do
  numbers@[x, y, z] <- sDoubles (pure <$> "xyz")
  constrain $ foldr (\n s -> fpIsPoint n &&& s) true numbers
  pure $ additionAssoc x y z

Proof of associative property of addition (take 2)

prove $ do
  numbers@[x, y, z] <- sDoubles (pure <$> "xyz")
  constrain $ foldr (\n s -> fpIsPoint n &&& s) true numbers
  pure $ additionAssoc x y z

Falsifiable. Counter-example:
  x = -4.4479747933244543e-308 :: Double
  y =   3.785766995733679e-270 :: Double
  z =  -3.785766995733679e-270 :: Double

Are IEEE754 numbers associative under addition?

> (x + y)
3.785766995733679e-270
> (x + y) + z
0.0
> y + z
0.0
> x + (y + z)
-4.4479747933244543e-308

How about integers?

prove $ do
  [x, y, z] <- sIntegers (pure <$> "xyz")
  pure $ additionAssoc x y z

How about integers?

prove $ do
  [x, y, z] <- sIntegers (pure <$> "xyz")
  pure $ additionAssoc x y z

Q.E.D.

Other proofs

Please check out the SBV documentation for additional cool things to prove, including:

  • correctness of mergeSort
  • correctness of a bit-hacky assembly algorithm for integer multiplication
    • this proof includes a formalization of the archaic architecture for which the clever algorithm was first written
  • correctness of a parallel algorithm for scanl1

A much more interesting proof

Proof that a feedback controller stabilizes a physical system, just by stating the theorem and applying an SMT solver.

A much more interesting proof

Proof that a feedback controller stabilizes a physical system, just by stating the theorem and applying an SMT solver.

Assumptions:

  • our model of the physical system captures all the relevant behavior
  • we know the relevant model parameters
  • we can directly measure the system's state

A much more interesting proof

Proof that a feedback controller stabilizes a physical system, just by stating the theorem and applying an SMT solver.

Assumptions:

  • our model of the physical system captures all the relevant behavior
  • we know the relevant model parameters
  • we can directly measure the system's state

Strategy:

  • Examine the dynamics of the system
    • Translate dynamics into Haskell
  • Change the system's dynamics with feedback control
    • Translate the controller into Haskell
  • Introduce Lyapunov's direct method for proving stability
    • Translate the relevant theorem into Haskell
  • Hook these parts up to SBV

A much more interesting proof

Proof that a feedback controller stabilizes a physical system, just by stating the theorem and applying an SMT solver.

Assumptions:

  • our model of the physical system captures all the relevant behavior
  • we know the relevant model parameters
  • we can directly measure the system's state

Strategy:

  • Examine the dynamics of the system
    • Translate dynamics into Haskell
  • Change the system's dynamics with feedback control
    • Translate the controller into Haskell
  • Introduce Lyapunov's direct method for proving stability
    • Translate the relevant theorem into Haskell
  • Hook these parts up to SBV
  • ???
  • Pro(o)fit!

The inverted pendulum

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

  • $\Huge \theta$ - angle

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

  • $\Huge \theta$ - angle
  • $\Huge \omega$ - velocity

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

  • $\Huge \theta$ - angle
  • $\Huge \omega$ - velocity
  • $\Huge \alpha$ - acceleration

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

  • $\Huge \theta$ - angle
  • $\Huge \omega$ - velocity
  • $\Huge \alpha$ - acceleration
  • $\Huge \tau$ - torque

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

  • $\Huge \theta$ - angle
  • $\Huge \omega$ - velocity
  • $\Huge \alpha$ - acceleration
  • $\Huge \tau$ - torque
  • $\Huge b$ - damping coefficient
  • $\Huge m$ - mass at the tip
  • $\Huge l$ - length of arm
  • $\Huge g$ - gravitational constant

The inverted pendulum

Here are the dynamics of the pendulum, written as a differential equation:

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

  • $\Huge \theta$ - angle
  • $\Huge \omega$ - velocity
  • $\Huge \alpha$ - acceleration
  • $\Huge \tau$ - torque
  • $\Huge b$ - damping coefficient
  • $\Huge m$ - mass at the tip
  • $\Huge l$ - length of arm
  • $\Huge g$ - gravitational constant

$\Huge \omega = \frac{d}{dt} \theta$

$\Huge \alpha = \frac{d}{dt} \omega$

The inverted pendulum (data types)

data Pendulum a = Pendulum
  { pendulumLength          :: a
  , pendulumDampingConstant :: a
  , pendulumMass            :: a
  , pendulumGravity         :: a
  } deriving (Eq, Show, Functor, Foldable, Traversable)

data State a = State
  { stateθ :: a
  , stateω :: a
  } deriving (Eq, Show, Functor, Foldable, Traversable)

newtype Controller a = Controller
  { controllerDamping :: a
  } deriving (Eq, Show, Functor, Foldable, Traversable)

The inverted pendulum (implementing the dynamics)

pendulumInertia :: Fractional a => Pendulum a -> a
pendulumInertia (Pendulum l _ m _) = m * l * l

The inverted pendulum (implementing the dynamics)

pendulumInertia :: Fractional a => Pendulum a -> a
pendulumInertia (Pendulum l _ m _) = m * l * l

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

The inverted pendulum (implementing the dynamics)

pendulumInertia :: Fractional a => Pendulum a -> a
pendulumInertia (Pendulum l _ m _) = m * l * l

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

$\Huge \frac{d}{dt}\begin{bmatrix}\theta \ \omega\end{bmatrix} = \begin{bmatrix} \omega \ \frac{b}{ml^2}\omega + \frac{g}{l} \sin \theta + \frac{1}{ml^2}\tau\end{bmatrix}$

The inverted pendulum (implementing the dynamics)

pendulumInertia :: Fractional a => Pendulum a -> a
pendulumInertia (Pendulum l _ m _) = m * l * l

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

$\Huge \frac{d}{dt}\begin{bmatrix}\theta \ \omega\end{bmatrix} = \begin{bmatrix} \omega \ \frac{b}{ml^2}\omega + \frac{g}{l} \sin \theta + \frac{1}{ml^2}\tau\end{bmatrix}$

pendulum :: Fractional a
         => Pendulum a  -- ^ System specification
         -> a           -- ^ Input torque
         -> State a     -- ^ Current state vector
         -> State a     -- ^ Time-derivative of state vector
pendulum sys@(Pendulum l b _ g) τ (State θ ω) =
  State ω $
  (g * taylorSin θ) / l + (b * ω) / inertia + τ / inertia
  where
    inertia = pendulumInertia sys

The inverted pendulum (simulation results -- angle)

$\Huge \left(\theta, \omega\right) \in \begin{bmatrix}(10^{-3}, 10^{-3}), & (-0.5, 0.1), & (0.3, 0.3)\end{bmatrix}$

The inverted pendulum (simulation results -- velocity)

$\Huge \left(\theta, \omega\right) \in \begin{bmatrix}(10^{-3}, 10^{-3}), & (-0.5, 0.1), & (0.3, 0.3)\end{bmatrix}$

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$ as a function of the state

$\Huge \vec{x} = \begin{bmatrix}\theta & \omega \end{bmatrix}^{T}$

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$ as a function of the state

$\Huge \vec{x} = \begin{bmatrix}\theta & \omega \end{bmatrix}^{T} \sim$ State a.

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$ as a function of the state

$\Huge \vec{x} = \begin{bmatrix}\theta & \omega \end{bmatrix}^{T} \sim ~$ State a.

Goals:

  • The pendulum should point straight up ($\Huge \theta = 0$).

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$ as a function of the state

$\Huge \vec{x} = \begin{bmatrix}\theta & \omega \end{bmatrix}^{T} \sim ~$ State a.

Goals:

  • The pendulum should point straight up ($\Huge \theta = 0$).
  • The pendulum shouldn't be moving ($\Huge \omega = 0$).

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$ as a function of the state

$\Huge \vec{x} = \begin{bmatrix}\theta & \omega \end{bmatrix}^{T} \sim ~$ State a.

Goals:

  • The pendulum should point straight up ($\Huge \theta = 0$).
  • The pendulum shouldn't be moving ($\Huge \omega = 0$).
  • No matter where you start it out or push it, it should eventually return to this state.

Controlling the pendulum

$\Huge ml^2\alpha = b\omega + mgl \sin \theta + \tau$

Feedback control: specify input torque $\Huge \tau$ as a function of the state

$\Huge \vec{x} = \begin{bmatrix}\theta & \omega \end{bmatrix}^{T} \sim ~$ State a.

Goals:

  • The pendulum should point straight up ($\Huge \theta = 0$).
  • The pendulum shouldn't be moving ($\Huge \omega = 0$).
  • No matter where you start it out or push it, it should eventually return to this state.
  • Less wobbly would be nice.

Controlling the pendulum (the control law)

Proposed feedback law:

$\Huge \tau = -2mgl \sin \theta - k_{d} \omega$

There are two parts:

  • $\Huge -2mgl \sin \theta$
  • $\Huge -k_{d} \omega$

Controlling the pendulum (the control law)

Proposed feedback law:

$\Huge \tau = -2mgl \sin \theta - k_{d} \omega$

There are two parts:

  • $\Huge -2mgl \sin \theta$

  • $\Huge -k_{d} \omega$

In Haskell:

runController :: Fractional a => Controller a -> Pendulum a -> State a -> a
runController (Controller kd) (Pendulum l _ m g) (State θ ω) =
  (-2) * m * g * l * taylorSin θ - kd * ω

Controlling the pendulum (the control law)

Proposed feedback law:

$\Huge \tau = -2mgl \sin \theta - k_{d} \omega$

There are two parts:

  • $\Huge -2mgl \sin \theta$

  • $\Huge -k_{d} \omega$

In Haskell:

runController :: Fractional a => Controller a -> Pendulum a -> State a -> a
runController (Controller kd) (Pendulum l _ m g) (State θ ω) =
  (-2) * m * g * l * taylorSin θ - kd * ω

closedLoop ::
  Fractional a => Controller a -> Pendulum a -> State a -> State a
closedLoop ctrl system initialState =
  pendulum system torque initialState
  where
    torque = runController ctrl system initialState

Controlling the pendulum (simulation results -- angle)

$\Huge \left(\theta, \omega\right) \in \begin{bmatrix}(10^{-3}, 10^{-3}), & (-0.5, 0.1), & (0.3, 0.3)\end{bmatrix}$

Controlling the pendulum (simulation results -- velocity)

$\Huge \left(\theta, \omega\right) \in \begin{bmatrix}(10^{-3}, 10^{-3}), & (-0.5, 0.1), & (0.3, 0.3)\end{bmatrix}$

Lyapunov's direct method

One of several theorems by Lyapunov:

$\Huge \forall V \in \mathbb{R}^{n} \to \mathbb{R}$

v :: Fractional a => State a -> a

$\Huge \vec{x} \sim$ State a

Lyapunov's direct method

One of several theorems by Lyapunov:

$\Huge \forall V \in \mathbb{R}^{n} \to \mathbb{R}$

v :: Fractional a => State a -> a

$\Huge \vec{x} \sim$ State a

If

  • $\Huge V({0}) = 0$
  • $\Huge V\left(\vec{x}\right) &gt; 0, \vec{x} \neq \vec{0}$
  • $\Huge \frac{d}{dt}V\left(\vec{x}\right) \leq 0$

Then the system is stable at $\Huge \vec{x} = 0$

Lyapunov's direct method

One of several theorems by Lyapunov:

$\Huge \forall V \in \mathbb{R}^{n} \to \mathbb{R}$

v :: Fractional a => State a -> a

$\Huge \vec{x} \sim$ State a

If

  • $\Huge V({0}) = 0$
  • $\Huge V\left(\vec{x}\right) &gt; 0, \vec{x} \neq \vec{0}$
  • $\Huge \frac{d}{dt}V\left(\vec{x}\right) \leq 0$

Then the system is stable at $\Huge \vec{x} = 0$

How can this prove stability if it doesn't mention the system dynamics $\Huge \left(\alpha = \texttt{\ldots}\right)$?

Lyapunov's direct method

One of several theorems by Lyapunov:

$\Huge \forall V \in \mathbb{R}^{n} \to \mathbb{R}$

v :: Fractional a => State a -> a

$\Huge \vec{x} \sim$ State a

If

  • $\Huge V({0}) = 0$
  • $\Huge V\left(\vec{x}\right) &gt; 0, \vec{x} \neq \vec{0}$
  • $\Huge \frac{d}{dt}V\left(\vec{x}\right) \leq 0$

Then the system is stable at $\Huge \vec{x} = 0$

How can this prove stability if it doesn't mention the system dynamics $\Huge \left(\alpha = \texttt{\ldots}\right)$?

$\Huge \frac{d}{dt}\begin{bmatrix}\theta \ \omega\end{bmatrix} = \begin{bmatrix} \omega \ \frac{b}{ml^2}\omega + \frac{g}{l} \sin \theta + \frac{1}{ml^2}\tau\end{bmatrix}$

Implementation of the theorem (Lyapunov function)

Kinetic energy: $\Huge \frac{1}{2} ml^2\omega^2$

kineticEnergy :: Fractional a => Pendulum a -> State a -> a
kineticEnergy system (State _ ω) =
  0.5 * pendulumInertia system * ω * ω

Implementation of the theorem (Lyapunov function)

Kinetic energy: $\Huge \frac{1}{2} ml^2\omega^2$

kineticEnergy :: Fractional a => Pendulum a -> State a -> a
kineticEnergy system (State _ ω) =
  0.5 * pendulumInertia system * ω * ω

Potential energy: $\Huge mgl \left(\cos \theta - 1\right)$

potentialEnergy :: Fractional a => Pendulum a -> State a -> a
potentialEnergy (Pendulum l _ m g) (State θ _) =
  m * g * l * (taylorCos θ - 1)

spans $\Huge mgl \cdot [-2, 0]$

Implementation of the theorem (Lyapunov function)

Kinetic energy: $\Huge \frac{1}{2}ml^2\omega^2$

kineticEnergy :: Fractional a => Pendulum a -> State a -> a
kineticEnergy system (State _ ω) =
  0.5 * pendulumInertia system * ω * ω

Potential energy: $\Huge mgl \left(\cos \theta - 1\right)$

potentialEnergy :: Fractional a => Pendulum a -> State a -> a
potentialEnergy (Pendulum l _ m g) (State θ _) =
  m * g * l * (taylorCos θ - 1)

$\Huge V(\vec{x}) = T - U$ (so that $\Huge V &gt; 0$)

v :: Fractional a => Pendulum a -> State a -> a
v system st =
  kineticEnergy system st - potentialEnergy system st

Implementation of the theorem (Lyapunov function derivative)

$\Huge \begin{matrix}\frac{d}{dt}V(\vec{x}) &= \frac{d}{dt}T - \frac{d}{dt}U \ ~ &= ml^2 \cdot \omega \cdot \alpha + mgl \sin\theta\end{matrix}$

dkedt :: Fractional a =>
  Controller a -> Pendulum a -> State a -> a
dkedt ctrl system state@(State _ ω) =
  pendulumInertia system * ω * stateω (closedLoop ctrl system state)

dpedt :: Fractional a =>
  Pendulum a -> State a -> a
dpedt (Pendulum l _ m g) (State θ ω) =
  m * g * l * (- taylorSin θ) * ω

dvdt :: Fractional a =>
  Controller a -> Pendulum a -> State a -> a
dvdt ctrl system st =
  dkedt ctrl system st - dpedt system st

Implementation of the theorem (proof statement)

lyapunov'sTheorem gen f dfdt = do
  st <- traverse gen stateLabels
  constrainPi st




  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

Implementation of the theorem (proof statement)

lyapunov'sTheorem gen f dfdt = do
  st <- traverse gen stateLabels
  constrainPi st
  eq <- lyapunovEquilibrium st



  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

$\Huge V(\vec{0}) = 0$

    lyapunovEquilibrium _ = pure $
      f (State 0 0) .== 0.0

Implementation of the theorem (proof statement)

lyapunov'sTheorem gen f dfdt = do
  st <- traverse gen stateLabels
  constrainPi st
  eq <- lyapunovEquilibrium st
  nn <- lyapunovNonNegative st


  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

$\Huge V(\vec{0}) = 0$

    lyapunovEquilibrium _ = pure $
      f (State 0 0) .== 0.0

$\Huge V(\vec{x}) &gt; 0, \vec{x} \neq \vec{0}$

    lyapunovNonNegative st = do
      constrain $ st ./= State 0 0
      pure $ f st .> 0.0

Implementation of the theorem (proof statement)

lyapunov'sTheorem gen f dfdt = do
  st <- traverse gen stateLabels
  constrainPi st
  eq <- lyapunovEquilibrium st
  nn <- lyapunovNonNegative st
  gn <- lyapunovGradNegative st

  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

$\Huge V(\vec{0}) = 0$

    lyapunovEquilibrium _ = pure $
      f (State 0 0) .== 0.0

$\Huge V(\vec{x}) &gt; 0, \vec{x} \neq \vec{0}$

    lyapunovNonNegative st = do
      constrain $ st ./= State 0 0
      pure $ f st .> 0.0

$\Huge \dot{V}(\vec{x}) &lt;= 0$

    lyapunovGradNegative st = pure $
      dfdt st .<= 0.0 &&& dfdt (State 0 0) .<= 0.0

Implementation of the theorem (proof statement)

lyapunov'sTheorem gen f dfdt = do
  st <- traverse gen stateLabels
  constrainPi st
  eq <- lyapunovEquilibrium st
  nn <- lyapunovNonNegative st
  gn <- lyapunovGradNegative st
  pure $ eq &&& nn &&& gn
  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

$\Huge V(\vec{0}) = 0$

    lyapunovEquilibrium _ = pure $
      f (State 0 0) .== 0.0

$\Huge V(\vec{x}) &gt; 0, \vec{x} \neq \vec{0}$

    lyapunovNonNegative st = do
      constrain $ st ./= State 0 0
      pure $ f st .> 0.0

$\Huge \dot{V}(\vec{x}) &lt;= 0$

    lyapunovGradNegative st = pure $
      dfdt st .<= 0.0 &&& dfdt (State 0 0) .<= 0.0

Proving Lyapunov's theorem for the stabilized pendulum

nominalController = Controller
  { controllerDamping = 0.3
  }

nominalSystem = Pendulum
  { pendulumLength = 0.5
  , pendulumDampingConstant = -0.03
  , pendulumMass = 0.1
  , pendulumGravity = 9.81
  }

Proving Lyapunov's theorem for the stabilized pendulum

nominalController = Controller
  { controllerDamping = 0.3
  }

nominalSystem = Pendulum
  { pendulumLength = 0.5
  , pendulumDampingConstant = -0.03
  , pendulumMass = 0.1
  , pendulumGravity = 9.81
  }

proveStability =
  prove $ lyapunov'sTheorem sReal v' dvdt'
  where
    v' = v nominalSystem
    dvdt' = dvdt nominalController nominalSystem

Proving Lyapunov's theorem for the stabilized pendulum

nominalController = Controller
  { controllerDamping = 0.3
  }

nominalSystem = Pendulum
  { pendulumLength = 0.5
  , pendulumDampingConstant = -0.03
  , pendulumMass = 0.1
  , pendulumGravity = 9.81
  }

proveStability =
  prove $ lyapunov'sTheorem sReal v' dvdt'
  where
    v' = v nominalSystem
    dvdt' = dvdt nominalController nominalSystem

Pendulum> proveStability

Proving Lyapunov's theorem for the stabilized pendulum

nominalController = Controller
  { controllerDamping = 0.3
  }

nominalSystem = Pendulum
  { pendulumLength = 0.5
  , pendulumDampingConstant = -0.03
  , pendulumMass = 0.1
  , pendulumGravity = 9.81
  }

proveStability =
  prove $ lyapunov'sTheorem sReal v' dvdt'
  where
    v' = v nominalSystem
    dvdt' = dvdt nominalController nominalSystem

Pendulum> proveStability
Q.E.D.

Implementation concerns

  • taylorSin and taylorCos
  • m * l * l
  • π = 3.14159...

Further implementation considerations (floating-point implementation)

Next theorem: if the state vector we input to the control law f

  • contains no NaN or Infinity values
  • has been correctly wrapped to the interval $\Huge \left(-\pi, \pi\right]$

then the output will also contain no NaN or Infinity.

Further implementation considerations (floating-point implementation)

Next theorem: if the state vector we input to the control law f

  • contains no NaN or Infinity values
  • has been correctly wrapped to the interval $\Huge \left(-\pi, \pi\right]$

then the output will also contain no NaN or Infinity.

nanFree f = do
  st <- State <$> sFloat "θ" <*> sFloat "ω"
  constrainPi st
  constrainFP st
  pure . fpIsPoint $ f st
  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

    constrainFP = constrain . allIsPoint

Further implementation considerations (floating-point implementation)

Next theorem: if the state vector we input to the control law f

  • contains no NaN or Infinity values
  • has been correctly wrapped to the interval $\Huge \left(-\pi, \pi\right]$

then the output will also contain no NaN or Infinity.

nanFree f = do
  st <- State <$> sFloat "θ" <*> sFloat "ω"
  constrainPi st
  constrainFP st
  pure . fpIsPoint $ f st
  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

    constrainFP = constrain . allIsPoint

Pendulum> proveNanSafety

Further implementation considerations (floating-point implementation)

Next theorem: if the state vector we input to the control law f

  • contains no NaN or Infinity values
  • has been correctly wrapped to the interval $\Huge \left(-\pi, \pi\right]$

then the output will also contain no NaN or Infinity.

nanFree f = do
  st <- State <$> sFloat "θ" <*> sFloat "ω"
  constrainPi st
  constrainFP st
  pure . fpIsPoint $ f st
  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

    constrainFP = constrain . allIsPoint

Pendulum> proveNanSafety
Q.E.D.

Real implementation considerations

emitController gen = compileToC Nothing "runController" $ do
  system <- Pendulum
    <$> gen "length"
    <*> gen "damping"
    <*> gen "mass"
    <*> gen "gravity"
  controller <- Controller <$> gen "kd"
  state <- State <$> gen "theta" <*> gen "omega"
  cgReturn $ runController controller system state

genCCode :: IO ()
genCCode = emitController cgGen
  where
    cgGen :: String -> SBVCodeGen SDouble
    cgGen = cgInput

Generated C code

SDouble runController(const SDouble length,
                      const SDouble damping, const SDouble mass, const SDouble gravity,
                      const SDouble kd, const SDouble theta, const SDouble omega)
{
  const SDouble s0 = length;
  const SDouble s2 = mass;
  const SDouble s3 = gravity;
  const SDouble s4 = kd;
  const SDouble s5 = theta;
  const SDouble s6 = omega;
  const SDouble s8 = s2 * 2.0;
  const SDouble s9 = s3 * s8;
  const SDouble s10 = s0 * s9;
  const SDouble s12 = s5 * s5;
  const SDouble s13 = s5 * s12;
  const SDouble s15 = s13 / 6.0;
  const SDouble s16 = (- s15);
  const SDouble s17 = 0.0 + s16;
  const SDouble s18 = s5 * s15;
  const SDouble s19 = s5 * s18;
  const SDouble s21 = s19 / 20.0;
  const SDouble s22 = s17 + s21;
  const SDouble s23 = s5 * s21;
  const SDouble s24 = s5 * s23;
  const SDouble s26 = s24 / 42.0;
  const SDouble s27 = (- s26);
  const SDouble s28 = s22 + s27;
  const SDouble s29 = s5 * s26;
  const SDouble s30 = s5 * s29;
  const SDouble s32 = s30 / 72.0;
  const SDouble s33 = s28 + s32;
  const SDouble s34 = s5 * s32;
  const SDouble s35 = s5 * s34;
  const SDouble s37 = s35 / 110.0;
  const SDouble s38 = (- s37);
  const SDouble s39 = s33 + s38;
  const SDouble s40 = s5 * s37;
  const SDouble s41 = s5 * s40;
  const SDouble s43 = s41 / 156.0;
  const SDouble s44 = s39 + s43;
  const SDouble s45 = s5 * s43;
  const SDouble s46 = s5 * s45;
  const SDouble s48 = s46 / 210.0;
  const SDouble s49 = (- s48);
  const SDouble s50 = s44 + s49;
  const SDouble s51 = s5 * s48;
  const SDouble s52 = s5 * s51;
  const SDouble s54 = s52 / 272.0;
  const SDouble s55 = s50 + s54;
  const SDouble s56 = s5 * s54;
  const SDouble s57 = s5 * s56;
  const SDouble s59 = s57 / 342.0;
  const SDouble s60 = (- s59);
  const SDouble s61 = s55 + s60;
  const SDouble s62 = s5 * s59;
  const SDouble s63 = s5 * s62;
  const SDouble s65 = s63 / 420.0;
  const SDouble s66 = s61 + s65;
  const SDouble s67 = s5 + s66;
  const SDouble s68 = s10 * s67;
  const SDouble s69 = (- s68);
  const SDouble s70 = (- s6);
  const SDouble s71 = s4 * s70;
  const SDouble s72 = s69 + s71;

  return s72;
}

Separate code-generation for Taylor series functions

emitTaylor f funName gen =
  compileToC Nothing funName $
  gen "x" >>= cgReturn . f

Separate code-generation for Taylor series functions

emitTaylor f funName gen =
  compileToC Nothing funName $
  gen "x" >>= cgReturn . f

...
emitTaylor taylorSin "taylorSin" cgGen
...

Separate code-generation for Taylor series functions

emitTaylor f funName gen =
  compileToC Nothing funName $
  gen "x" >>= cgReturn . f

...
emitTaylor taylorSin "taylorSin" cgGen
...

taylorSin' :: (Fractional a, SymWord a) => SBV a -> SBV a
taylorSin' = cgUninterpret "taylorSin" mempty taylorSin

Separate code-generation for Taylor series functions (sine)

SDouble taylorSin(const SDouble x)
{
  const SDouble s0 = x;
  const SDouble s2 = s0 * s0;
  const SDouble s3 = s0 * s2;
  const SDouble s5 = s3 / 6.0;
  const SDouble s6 = (- s5);
  const SDouble s7 = 0.0 + s6;
  const SDouble s8 = s0 * s5;
  const SDouble s9 = s0 * s8;
  const SDouble s11 = s9 / 20.0;
  const SDouble s12 = s7 + s11;
  const SDouble s13 = s0 * s11;
  const SDouble s14 = s0 * s13;
  const SDouble s16 = s14 / 42.0;
  const SDouble s17 = (- s16);
  const SDouble s18 = s12 + s17;
  const SDouble s19 = s0 * s16;
  const SDouble s20 = s0 * s19;
  const SDouble s22 = s20 / 72.0;
  const SDouble s23 = s18 + s22;
  const SDouble s24 = s0 * s22;
  const SDouble s25 = s0 * s24;
  const SDouble s27 = s25 / 110.0;
  const SDouble s28 = (- s27);
  const SDouble s29 = s23 + s28;
  const SDouble s30 = s0 * s27;
  const SDouble s31 = s0 * s30;
  const SDouble s33 = s31 / 156.0;
  const SDouble s34 = s29 + s33;
  const SDouble s35 = s0 * s33;
  const SDouble s36 = s0 * s35;
  const SDouble s38 = s36 / 210.0;
  const SDouble s39 = (- s38);
  const SDouble s40 = s34 + s39;
  const SDouble s41 = s0 * s38;
  const SDouble s42 = s0 * s41;
  const SDouble s44 = s42 / 272.0;
  const SDouble s45 = s40 + s44;
  const SDouble s46 = s0 * s44;
  const SDouble s47 = s0 * s46;
  const SDouble s49 = s47 / 342.0;
  const SDouble s50 = (- s49);
  const SDouble s51 = s45 + s50;
  const SDouble s52 = s0 * s49;
  const SDouble s53 = s0 * s52;
  const SDouble s55 = s53 / 420.0;
  const SDouble s56 = s51 + s55;
  const SDouble s57 = s0 + s56;

  return s57;
}

Separate code-generation for Taylor series functions (controller)

SDouble runController(const SDouble length,
                      const SDouble damping, const SDouble mass, const SDouble gravity,
                      const SDouble kd, const SDouble theta, const SDouble omega)
{
  const SDouble s0 = length;
  const SDouble s2 = mass;
  const SDouble s3 = gravity;
  const SDouble s4 = kd;
  const SDouble s5 = theta;
  const SDouble s6 = omega;
  const SDouble s8 = s2 * 2.0;
  const SDouble s9 = s3 * s8;
  const SDouble s10 = s0 * s9;
  const SDouble s11 = /* Uninterpreted function */ taylorSin(s5);
  const SDouble s12 = s10 * s11;
  const SDouble s13 = (- s12);
  const SDouble s14 = (- s6);
  const SDouble s15 = s4 * s14;
  const SDouble s16 = s13 + s15;

  return s16;
}

Constraint solving

$\Huge \forall$

$\Huge \exists$

Constraint solving

$\Huge \forall x. p(x) \to \exists x. p(x)$

SEND + MORE = MONEY

sendMoreMoney :: IO AllSatResult
sendMoreMoney = allSat $ do
  ds@[s,e,n,d,m,o,r,y] <- sIntegers $ pure <$> ["sendmory"]
  let
    isDigit x = x .>= 0 &&& x .<= 9
    val xs    = sum $ zipWith (*)
                  (reverse xs)
                  (iterate (*10) 1)
    send      = val [s,e,n,d]
    more      = val [m,o,r,e]
    money     = val [m,o,n,e,y]
  constrain $ bAll isDigit ds
  constrain $ distinct ds
  constrain $ s ./= 0 &&& m ./= 0
  solve [send + more .== money]

SEND + MORE = MONEY

sendMoreMoney :: IO AllSatResult
sendMoreMoney = allSat $ do
  ds@[s,e,n,d,m,o,r,y] <- sIntegers $ pure <$> ["sendmory"]
  let
    isDigit x = x .>= 0 &&& x .<= 9
    val xs    = sum $ zipWith (*)
                  (reverse xs)
                  (iterate (*10) 1)
    send      = val [s,e,n,d]
    more      = val [m,o,r,e]
    money     = val [m,o,n,e,y]
  constrain $ bAll isDigit ds
  constrain $ distinct ds
  constrain $ s ./= 0 &&& m ./= 0
  solve [send + more .== money]

 SendMoreMoney> sendMoreMoney
 Solution #1:
   s = 9 :: Integer
   e = 5 :: Integer
   n = 6 :: Integer
   d = 7 :: Integer
   m = 1 :: Integer
   o = 0 :: Integer
   r = 8 :: Integer
   y = 2 :: Integer

More constraint solving problems

Please check out the SBV documentation for more constraint problem examples, including:

  • n-queens
  • sudoku
  • "Cheryl's birthday problem"
  • knapsack-type problems

Optimization

A simple integer linear programming (ILP) problem:

optimize Lexicographic $ do
  x <- sInteger "x"
  y <- sInteger "y"
  constrain $ x .> 0
  constrain $ x .< 6
  constrain $ y .> 2
  constrain $ y .< 12
  minimize "goal" $ x + 2 * y

Optimization

A simple integer linear programming (ILP) problem:

optimize Lexicographic $ do
  x <- sInteger "x"
  y <- sInteger "y"
  constrain $ x .> 0
  constrain $ x .< 6
  constrain $ y .> 2
  constrain $ y .< 12
  minimize "goal" $ x + 2 * y

Optimal model:
  x    = 1 :: Integer
  y    = 3 :: Integer
  goal = 7 :: Integer

Resources

See Wikipedia for an overview of SAT and SMT

SMTLIB homepage: http://smtlib.cs.uiowa.edu/

SBV homepage: https://leventerkok.github.io/sbv/

SBV hackage docs: https://hackage.haskell.org/package/sbv

Z3 homepage: https://github.com/Z3Prover/z3/wiki

These slides: https://peddie.github.io/smt-solving/

Slide repository, code, etc.: https://github.com/peddie/smt-solving

BACKUPS

Taylor Series functions

10th-order Taylor series approximation to the basic trig functions:

taylorCos :: Fractional a => a -> a
taylorCos x = 1 + sum (take 10 series)
  where
    inc num old =
      let new = old * x * x / (num * (num + 1))
      in new : inc (num + 2) new
    signs = cycle [negate, id]
    series = zipWith ($) signs (inc 1 1)

taylorSin :: Fractional a => a -> a
taylorSin x = x + sum (take 10 series)
  where
    inc num old =
      let new = old * x * x / (num * (num + 1))
      in new : inc (num + 2) new
    signs = cycle [negate, id]
    series = zipWith ($) signs (inc 2 x)

Regexp Crosswords

(Here is another example that ships with SBV.) How many people are familiar with regexp crosswords?

Regular expression crossword

https://regexcrossword.com/challenges/intermediate/puzzles/1

Regexp Crossword solution (part 1)

import qualified Data.SBV.String as S
import qualified Data.SBV.RegExp as R

-- | Solve a given crossword, returning the corresponding rows
solveCrossword :: [R.RegExp] -> [R.RegExp] -> IO [String]
solveCrossword rowRegExps colRegExps = runSMT $ do
        let numRows = genericLength rowRegExps
            numCols = genericLength colRegExps

        -- constrain rows
        let mkRow rowRegExp = do row <- free_
                                 constrain $ row `R.match` rowRegExp
                                 constrain $ S.length row .== literal numCols
                                 return row

        rows <- mapM mkRow rowRegExps

        -- constrain colums
        let mkCol colRegExp = do col <- free_
                                 constrain $ col `R.match` colRegExp
                                 constrain $ S.length col .== literal numRows
                                 return col

        cols <- mapM mkCol colRegExps

Regexp crossword solution (part 2)

. . .

        -- constrain each "cell" as they rows/columns intersect:
        let rowss =           [[r .!! literal i | i <- [0..numCols-1]] | r <- rows]
        let colss = transpose [[c .!! literal i | i <- [0..numRows-1]] | c <- cols]

        constrain $ bAnd $ zipWith (.==) (concat rowss) (concat colss)

        -- Now query to extract the solution
        query $ do cs <- checkSat
                   case cs of
                     Unk   -> error "Solver returned unknown!"
                     Unsat -> error "There are no solutions to this puzzle!"
                     Sat   -> mapM getValue rows

Compiler plugin

The SBV library is also available as a GHC plugin, which runs at compile time and can accept program annotations telling it what parts of the program to prove. This lets you integrate theorem-proving into the compilation phase of your system.