<jupyter_start><jupyter_text>Introduction to Models and Reachability in Kaa We first import the necessary modules and classes from Kaa. To start defining models, we at least need the sympy and numpy packages alongside the Model and Bundle classes.<jupyter_code>import sympy as sp
import numpy as np

from kaa.bundle import Bundle
from kaa.model import Model<jupyter_output><empty_output><jupyter_text>Basic Parallelotope  A model file in Kaa consists of first defining the variables and the dynamics. Kaa is designed for the analysis of **polynomial discrete dynamics**. The following model definition illustates how to define the unitbox shifting diagonally up to the right.<jupyter_code>class Basic(Model):

    def __init__(self):

        x,y = sp.Symbol('x'), sp.Symbol('y')

        dx = x + 1
        dy = y + 1

        dyns  = [dx, dy]
        vars = [x, y]

        L = np.empty([2,2])
        T = np.empty(2)

        L[0] = [1, 0]
        L[1] = [0, 1]

        T[0][0] = 0
        T[0][1] = 1

        offu = np.empty(2)
        offl = np.empty(2)

        offu[0] = 1
        offu[1] = 1
    
        offl[0] = 1
        offl[1] = 1
  
        b = Bundle(T, L, offu, offl, vars)

        super().__init__(b, dyns, vars)<jupyter_output><empty_output><jupyter_text>The numpy matrix $L$ represents the directions matrix of the model. Each row of the matrix reprsents a normal vector to the plane used to define the positive halfspace i.e the $i^{th}$ row of $L$ defines the halfspaces $L_i{\bf x} \leq \text{offu}_i$ and $-L_i{\bf x} \leq \text{offl}_i$ where $\text{offu}$ is the vector of offsets for the upper facets (corresponding to the $L_i$) and $\text{offl}$ is the vector of offsets for the lower facets (corresponding to the $-L_i$). The matrix $T$ represents the template matrix defining the parallelotope bundle of the model. Each row $T_i$ represents a vector of row-indices in $L$ defining parallelotope $i$ in the bundle. This defines $n$ hyperplanes and their parallel counterparts from the directions matrix. In other words, element $T_{ij}$ indicates that the postive and negative halfspaces defined by normal vector $L_i, -L_i$ with offsets in $\text{offu}_i$ and $\text{offl}_i$ respectively is contained in parallelotope $j$. The intersection of the positive and negative halfspaces will be the parallelotope in question. Note that the number of columns of $T$ must exactly be the dimension of the system.  Thus, observe that definitions of $L$, $T$, $\text{offu}$. and $\text{offl}$ designate the initial set as the **unit box** in $\mathbb{R}^2$. To initialize the reachability computation, we must first create an initial bundle denoted by the variable $b$ in this case. The bundle is a data structure which contains the directions matrix $L$ along with the upper and lower offsets $\text{offu}$, $\text{offl}$ and the template matrix $T$. It additionally encapsulates methods describing the parallelotope bundle in terms of linear constraints (half-spaces).In this case, there is only one row in the template matrix. Examples below will leverage multiple parallelotopes to generate tighter over-approximations. We now construct a Reach object and a FlowPipePlotter object to calculate the reachable set and plot the projection to the $x$-axis for 150 steps.First, we import the necessary modules and classes as usual:<jupyter_code>from kaa.reach import ReachSet
from kaa.flowpipe import FlowPipePlotter<jupyter_output><empty_output><jupyter_text>kaa.reach.ReachSet is the object responsible for computing the reachale set for a supplied number of time steps. The constructor takes in a Model object defining the dynamics and the initial bundle. Calling **ReachSet.computeReachSet(num)** will calculate the reachable set for $num$ time steps. To plot the projection on the $x$-axis, we first create a **FlowPipePlotter** object supplied with the **Flowpipe** returned by **ReachSet.computeReachSet(num)** and feed it into **FlowPipePlotter.plot2DProj(var)** where $var$ is the index of the desired variable. The relevant code is posted below:<jupyter_code>basic_model = Basic()
basic_reach = ReachSet(basic_model)
basic_flowpipe = basic_reach.computeReachSet(150)
FlowPipePlotter(basic_flowpipe).plot2DProj(0)<jupyter_output><empty_output><jupyter_text>We can also plot the phase diagram for this model as follows. This should just be the unit box shifting along the positive diagonal for 150 steps.  Phase Plot for Basic (x,y)<jupyter_code>basic_model = basic.Basic()
basic_reach = ReachSet(basic_model)
FlowPipePlotter(basic_reach.computeReachSet(150)).plot2DPhase(0,1)<jupyter_output><empty_output><jupyter_text>SIR Epidemic Model We now investigate a slightly more complicated model, the SIR epidemic model. Following the format from the example above, we define the following model:<jupyter_code>class SIR(Model):

  def __init__(self):

      s, i, r = sp.Symbol('s'), sp.Symbol('i'), sp.Symbol('r')

      ds = s - (0.34*s*i)*0.05;
      di = i + (0.34*s*i - 0.05*i)*0.05;
      dr = r + 0.05*i*0.05;

      dyns = [ds, di, dr]
      vars = [s, i, r] 
      sys_dim = len(vars)

      num_direct = 3
      num_temps = 1

      L = np.zeros([num_direct, sys_dim])
      T = np.zeros([num_temps, sys_dim])

      # Directions matrix
      L[0][0] = 1  #[1 0 0 ]^T
      L[1][1] = 1  #[0 1 0 ]^T
      L[2][2] = 1  #[0 0 1 ]^T

      # Template matrix
      T[0][0] = 0 
      T[0][1] = 1 
      T[0][2] = 2

      offu = np.zeros(num_direct)
      offl = np.zeros(num_direct)

      offu[0] = 0.8
      offl[0] = -0.79

      offu[1] = 0.2
      offl[1] = -0.19

      offu[2] = 0.001
      offl[2] = -0.00099

      b = Bundle(T, L, offu, offl, vars)
      super().__init__(b, dyns, vars)<jupyter_output><empty_output><jupyter_text>The initial set is taken from Dreossi et. al. Note that we only use one parallelotope here as $T$ only has one row. Furthermore, this parallelotope uses the axis-aligned unit vectors as its directions. Based on the definitions of $\text{offu}$ and $\text{offl}$, observe that initial set will be:$$ s \in [0.79, 0.8] \quad i \in [0.19, 0.2]  \quad r = 0 $$We can now calculate the flowpipe for 150 steps, and plot the resulting projections and phase plots.   SIR Model projection to r<jupyter_code>sir_model = SIR()
sir_reach = ReachSet(sir_model)
sir_flow = sir_reach.computeReachSet(150)
FlowPipePlotter(sir_flow).plot2DProj(2)<jupyter_output><empty_output><jupyter_text>SIR Model projection to i<jupyter_code>FlowPipePlotter(sir_flow).plot2DProj(1)<jupyter_output><empty_output><jupyter_text>Phase Plot for SIR (S,I)<jupyter_code>FlowPipePlotter(flowpipe).plot2DPhase(0,1)<jupyter_output><empty_output><jupyter_text>Phase Plot for SIR (I,R)<jupyter_code>FlowPipePlotter(flowpipe).plot2DPhase(1,2)<jupyter_output><empty_output><jupyter_text>Rossler Model We finally introduce the Rossler model with parameters set by Dreossi et. al as follows:<jupyter_code>class Rossler(Model):

    def __init__(self):

        x, y, z = sp.Symbol('x'), sp.Symbol('y'), sp.Symbol('z')
        vars = [x, y, z]

        dim_sys = len(vars)

        dx = x + (-y-z)*0.025
        dy = y + (x + 0.1*y)*0.025
        dz = z + (0.1 + z*(x-14))*0.02

        dyns = [dx, dy ,dz]

        num_direct = 5
        num_temps = 3

        L = np.zeros([num_direct, dim_sys])
        T = np.zeros([num_temps, dim_sys])

        L[0][0] = 1
        L[1][1] = 1
        L[2][2] = 1

        L[3][0] = 1
        L[3][1] = 0.5

        L[4][0] = 0.5
        L[4][2] = 0.5

        T[0][0] = 0; T[0][1] = 1; T[0][2] = 2;
        T[1][0] = 1; T[1][1] = 2; T[1][2] = 3;
        T[2][0] = 2; T[2][1] = 3; T[2][2] = 4;

        offu = np.zeros(num_direct)
        offl = np.zeros(num_direct)

        offu[0] = 0.1; offl[0] = -0.09;
        offu[1] = 5; offl[1] = -4.99;
        offu[2] = 0.1; offl[2] = -0.09;
        offu[3] = 10; offl[3] = 0;
        offu[4] = 10; offl[4] = 0;

        b = Bundle(T, L, offu, offl, vars)

        super().__init__(b, dyns, vars)<jupyter_output><empty_output><jupyter_text>Rossler Model projection to x<jupyter_code>rossler_model = rossler.Rossler()
rossler_reach = ReachSet(rossler_model)
rossler_flow = rossler_reach.computeReachSet(150)
FlowPipePlotter(rossler_flow).plot2DProj(0)<jupyter_output><empty_output><jupyter_text>Rossler Model projection to z<jupyter_code>FlowPipePlotter(rossler_flow).plot2DProj(1)<jupyter_output><empty_output><jupyter_text>Kaa also supports phase plotting. Run the examples below to see the phase plots for the Rossler model.  Phase Plot for Rossler (x,y)<jupyter_code>FlowPipePlotter(rossler_flow).plot2DPhase(0,1)<jupyter_output><empty_output><jupyter_text>Phase Plot for Rossler (y,z)<jupyter_code>FlowPipePlotter(rossler_flow).plot2DPhase(0,2)<jupyter_output><empty_output><jupyter_text>Quadcopter modelThe Quadcopter model is a system of 17 dimensions modeling the dynamics of a quadcopter. We define its dynamics below:<jupyter_code>class Quadcopter(Model):

    def __init__(self):

           dim_sys = 17
           pn, pe, h, u, v, w, q0v, q1v = sp.Symbol("pn"), sp.Symbol("pe"), sp.Symbol("h"), sp.Symbol("u"), sp.Symbol("v"), sp.Symbol("w"), sp.Symbol("q0v"), sp.Symbol("q1v")
           q2v, q3v, p, q, r, hI, uI, vI, psiI = sp.Symbol("q2v"), sp.Symbol("q3v"), sp.Symbol("p"),sp.Symbol("q"),sp.Symbol("r"), sp.Symbol("hI"), sp.Symbol("uI"), sp.Symbol("vI"), sp.Symbol("psiI")

           vars = [pn, pe, h, u, v, w, q0v, q1v, q2v, q3v, p, q, r, hI, uI, vI, psiI]

           M = 0.0015;
           mr = 0.001;
           R = 0.020;
           l = 0.045;
           g = 9.81;
           m = M + 4*mr;
           Jx = (2*M*R**2)/5 + (2*l**2)*mr;
           Jy = (2*M*R**2)/5 + (2*l**2)*mr;
           Jz = (2*M*R**2)/5 + (4*l**2)*mr;

           ur = 0;
           vr = 0;
           psir = 0;
           hr = 1;

           phi = 2*q1v;
           theta = 2*q2v;
           psi = 2*q3v;

           delta = 0.01;

           F = 0.0361*hI + 0.0694*h + 0.0603*w;
           tauphi = -0.0003*vI - 0.0005*v - 0.0018*phi - 0.0004*p;
           tautheta = 0.0003*uI + 0.0005*u - 0.0018*theta - 0.0004*q;
           taupsi = -0.0003*psiI - 0.0006*psi - 0.0003*r;

           dpn = pn + (u*(2*q0v**2 + 2*q1v**2 - 1) - v*(2*q0v*q3v - 2*q1v*q2v ) + w*(2*q0v*q2v + 2*q1v*q3v ))*delta;
           dpe = pe + (v*(2*q0v**2 + 2*(q2v**2) - 1) + u*(2*q0v*q3v + 2*q1v*q2v ) - w*(2*q0v*q1v - 2*q2v*q3v ))*delta;
           dh = h + (w*(2*(q0v**2) + 2*(q3v**2) - 1) - u*(2*q0v*q2v - 2*q1v*q3v ) + v*(2*q0v*q1v + 2*q2v*q3v ))*delta;

           du = u + (r*v - q*w - g*(2*q0v*q2v - 2*q1v*q3v ))*delta;
           dv = v + (p*w - r*u + g*(2*q0v*q1v + 2*q2v*q3v ))*delta;
           dw = w + (q*u - p*v - F/m + g*(2*(q0v**2) + 2*(q3v**2) - 1 ))*delta;

           dq0v = q0v +(-(q1v/2)*p - (q2v/2)*q - (q3v/2)*r)*delta;
           dq1v = q1v + ((q0v/2)*p - (q3v/2)*q + (q2v/2)*r)*delta;
           dq2v = q2v + ((q3v/2)*p + (q0v/2)*q - (q1v/2)*r)*delta;
           dq3v = q3v + ((q1v/2)*q - (q2v/2)*p + (q0v/2)*r)*delta;

           dp = p + ((1/Jx)*tauphi + ((Jy - Jz)/Jx)*q*r)*delta;
           dq = q + ((1/Jy)*tautheta - ((Jx - Jz)/Jy)*p*r)*delta;
           dr = r + ((1/Jz)*taupsi + ((Jx - Jy)/Jz)*p*q)*delta;

           dhI = hI + (h - hr)*delta;
           duI = uI +(u - ur)*delta;
           dvI = vI + (v - vr)*delta;
           dpsiI = psiI + (psi - psir)*delta;

           dyns = [dpn,dpe,dh,du,dv,dw,dq0v,dq1v,dq2v,dq3v,dp,dq,dr,dhI,duI,dvI,dpsiI]

           num_dirs = 18
           num_temp = 2

           L = np.zeros([num_dirs,dim_sys])

           for i in range(dim_sys):
               L[i][i] = 1

           L[17][2] = 0.5; L[17][5] = 0.5; L[17][6] = 0.5; L[17][15] = 0.25;

           T = np.zeros([num_temp, dim_sys]);
           for i in range(dim_sys):
               T[0][i] = i
               T[1][i] = i

           T[1][5] = 17

           offu = np.zeros(num_dirs);
           offl = np.zeros(num_dirs);

           offu[2] = 0.21; offl[2] = -0.20;
           offu[6] = 1; offl[6] = -1;
           offu[17] = 100; offl[17] = 100;

           b = Bundle(T, L, offu, offl, vars)
           super().__init__(b, dyns, vars)<jupyter_output><empty_output><jupyter_text>We now run the model. Note that this is an example where multiple parallelotopes increase the quality of the over-approximation.As of writing this, the model takes a while to compute.<jupyter_code>model = Quadcopter()
mod_reach = ReachSet(model)
mod_flow = mod_reach.computeReachSet(300)

FlowPipePlotter(mod_flow).plot2DProj(2,5,13)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14465408805031446, 'fraction_numerical': 0.04569699864660457, 'mean_word_length': 2.960277427490542, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '286668', 'license': 'MIT', 'n_tokens_mistral': 5115, 'n_tokens_neox': 4727, 'n_words': 1295, 'path': 'kaa-intro.ipynb', 'repo_name': 'Tarheel-Formal-Methods/kaa', 'size': '12561'}