Examples

This page contains some examples of PSY-TaLiRo tests interfacing with different systems. These examples are also available in the example directory of the project repository. Each example will execute a test and create a graph to visualize the outputs.

ODE

This example defines a model using a basic first-order ODE which represents an aircraft in flight.

import logging
import math

import plotly.graph_objects as go

import staliro
import staliro.models as models
import staliro.optimizers as optimizers
import staliro.specifications as specifications


@models.ode()
def nonlinear_model(inputs: models.Ode.Inputs) -> dict[str, float]:
    x1 = inputs.state["x1"]
    x2 = inputs.state["x2"]

    return {
        "x1": x1 - x2 + 0.1 * inputs.time,  # x1_dot
        "x2": x2 * math.cos(2 * math.pi * x1) + 0.1 * inputs.time,  # x2_dot
    }


phi = r"always !(a >= -1.6 and a <= -1.4  and b >= -1.1 and b <= -0.9)"
specification = specifications.rtamt.parse_dense(phi, {"a": 0, "b": 1})
optimizer = optimizers.UniformRandom[float]()
options = staliro.TestOptions(
    runs=1,
    iterations=100,
    tspan=(0, 2),
    static_inputs={
        "x1": (-1, 1),
        "x2": (-1, 1),
    },
)

if __name__ == "__main__":
    logging.basicConfig(level=logging.DEBUG)

    runs = staliro.test(nonlinear_model, specification, optimizer, options)
    run = runs[0]
    min_eval = min(run.evaluations, key=lambda e: e.cost)

    figure = go.Figure()
    figure.add_trace(
        go.Scatter(
            name="Falsification area",
            x=[-1.6, -1.4, -1.4, -1.6],
            y=[-1.1, -1.1, -0.9, -0.9],
            fill="toself",
            fillcolor="red",
            line_color="red",
            mode="lines+markers",
        )
    )

    figure.add_trace(
        go.Scatter(
            name="Initial condition region",
            x=[-1, 1, 1, -1],
            y=[-1, -1, 1, 1],
            fill="toself",
            fillcolor="green",
            line_color="green",
            mode="lines+markers",
        )
    )

    figure.add_trace(
        go.Scatter(
            name="Samples",
            x=[evaluation.sample.static["x1"] for evaluation in run.evaluations],
            y=[evaluation.sample.static["x2"] for evaluation in run.evaluations],
            mode="markers",
            marker=go.scatter.Marker(color="lightblue", symbol="circle"),
        )
    )

    figure.add_trace(
        go.Scatter(
            name="Best evaluation trajectory",
            x=[state[0] for state in min_eval.extra.trace.states],
            y=[state[1] for state in min_eval.extra.trace.states],
            mode="lines+markers",
            line=go.scatter.Line(color="blue", shape="spline"),
        )
    )

    figure.write_image("nonlinear.jpeg")

Blackbox

This example defines a model of the F16 Ground Collision Avoidance System (GCAS) which has been used to save pilots who have become unresponsive from crashing into the ground. The implementation of the GCAS system used in this example can be found here.

import logging
import math
from typing import Final

import numpy as np
import plotly.graph_objects as go
from aerobench.examples.gcas.gcas_autopilot import GcasAutopilot
from aerobench.run_f16_sim import run_f16_sim

from staliro import TestOptions, Trace, staliro
from staliro.models import Blackbox, blackbox
from staliro.optimizers import DualAnnealing
from staliro.specifications import rtamt

TSPAN: Final[tuple[float, float]] = (0, 15)


@blackbox()
def f16_model(inputs: Blackbox.Inputs) -> Trace[list[float]]:
    power = 9
    alpha = np.deg2rad(2.1215)
    beta = 0
    alt = 2330
    vel = 540
    phi = inputs.static["phi"]
    theta = inputs.static["theta"]
    psi = inputs.static["psi"]

    initial_state = [vel, alpha, beta, phi, theta, psi, 0, 0, 0, 0, 0, 0, alt, power]
    step = 1.0 / 30.0
    autopilot = GcasAutopilot(init_mode="roll", stdout=False)
    result = run_f16_sim(initial_state, TSPAN[1], autopilot, step, extended_states=True)
    states = np.vstack(
        (
            np.array([0 if x == "standby" else 1 for x in result["modes"]]),
            result["states"][:, 4],  # roll
            result["states"][:, 5],  # pitch
            result["states"][:, 6],  # yaw
            result["states"][:, 12],  # altitude
        )
    )

    return Trace(times=result["times"], states=np.transpose(states).tolist())


spec = rtamt.parse_dense("always (alt > 0)", {"alt": 4})
optimizer = DualAnnealing()
options = TestOptions(
    runs=1,
    iterations=10,
    tspan=TSPAN,
    static_inputs={
        "phi": math.pi / 4 + np.array([-math.pi / 20, math.pi / 30]),
        "theta": -math.pi / 2 * 0.8 + np.array([0, math.pi / 20]),
        "psi": -math.pi / 4 + np.array([-math.pi / 8, math.pi / 8]),
    },
)

if __name__ == "__main__":
    logging.basicConfig(level=logging.DEBUG)

    runs = staliro(f16_model, spec, optimizer, options)
    run = runs[0]
    min_cost_eval = min(run.evaluations, key=lambda e: e.cost)
    min_cost_trace = min_cost_eval.extra.trace

    figure = go.Figure()
    figure.update_layout(xaxis_title="time (s)", yaxis_title="alt (m)")
    figure.add_hline(y=0, line_color="red")
    figure.add_trace(
        go.Scatter(
            x=list(min_cost_trace.times),
            y=[state[4] for state in min_cost_trace.states],
            mode="lines",
            line_color="green",
            name="altitude",
        )
    )

    figure.write_image("f16.jpeg")

MATLAB/Simulink

This example demonstrates executing a simulink model of an automatic transmission. In order to interface with MATLAB you will need to ensure you have the MATLAB python interface installed according to the MathWorks instructions.

import logging
import pathlib

import matlab
import matlab.engine
import numpy as np
import plotly.graph_objects as go
import plotly.subplots as sp

from staliro import Sample, SignalInput, TestOptions, staliro
from staliro.models import Model, Result
from staliro.optimizers import DualAnnealing
from staliro.specifications import rtamt


class AutotransModel(Model[list[float], None]):
    MODEL_NAME = "Autotrans_shift"

    def __init__(self) -> None:
        script_path = pathlib.Path(__name__)

        self.sampling_step = 0.2
        self.engine = matlab.engine.start_matlab()
        self.engine.addpath(str(script_path.parent))

        model_opts = self.engine.simget(AutotransModel.MODEL_NAME)
        self.model_opts = self.engine.simset(model_opts, "SaveFormat", "Array")

    def simulate(self, sample: Sample) -> Result[list[float], None]:
        assert sample.signals.tspan is not None

        tstart, tend = sample.signals.tspan
        duration = tend - tstart
        sim_t = matlab.double([0, tend])
        n_times = duration // self.sampling_step
        signal_times = np.linspace(tstart, tend, num=int(n_times))
        signal_values = np.array(
            [[signal.at_time(t) for t in signal_times] for signal in sample.signals]
        )

        model_input = matlab.double(np.row_stack((signal_times, signal_values)).T.tolist())
        timestamps, _, data = self.engine.sim(
            self.MODEL_NAME, sim_t, self.model_opts, model_input, nargout=3
        )

        times: list[float] = np.array(timestamps).flatten().tolist()
        states: list[list[float]] = list(data)

        return Result(times=times, states=states, extra=None)


model = AutotransModel()

phi = "always[0,30] (rpm >= 3000) -> (always[0,4] speed >= 35)"
specification = rtamt.parse_discrete(phi, {"rpm": 0, "speed": 1})

optimizer = DualAnnealing()

options = TestOptions(
    runs=1,
    iterations=100,
    tspan=(0, 30),
    signals={
        "throttle": SignalInput(control_points=[(0, 100)] * 7),
        "brake": SignalInput(control_points=[(0, 350)] * 3),
    },
)

if __name__ == "__main__":
    logging.basicConfig(level=logging.DEBUG)

    runs = staliro(model, specification, optimizer, options)
    run = runs[0]
    min_eval = min(run.evaluations, key=lambda e: e.cost)

    times = list(min_eval.extra.trace.times)
    rpm = [state[0] for state in min_eval.extra.trace.states]
    speed = [state[1] for state in min_eval.extra.trace.states]

    figure = sp.make_subplots(rows=2, cols=1, shared_xaxes=True, x_title="Time (s)")
    figure.add_trace(go.Scatter(x=times, y=rpm), row=1, col=1)
    figure.add_trace(go.Scatter(x=times, y=speed), row=2, col=1)
    figure.update_yaxes(title_text="RPM", row=1, col=1)
    figure.update_yaxes(title_text="Speed", row=2, col=1)
    figure.write_image("autotrans.jpeg")

Bi-Level

This example demonstrates a bi-level test where the outer optimization attempt searches over the initial altitude of the aircraft and the inner optimization attempt searches over the initial attitude of the aircraft.

import math

import aerobench.examples.gcas.gcas_autopilot as autopilot
import aerobench.run_f16_sim as f16_sim
import numpy as np

import staliro
import staliro.models as models
import staliro.optimizers as optimizers
import staliro.specifications as specifications

TSPAN = (0, 15)


@staliro.costfunc()
def outer(sample: staliro.Sample) -> float:
    @models.blackbox(step_size=0.1)
    def inner(inputs: models.Blackbox.Inputs) -> models.Trace[list[float]]:
        power = 9
        alpha = np.deg2rad(2.1215)
        beta = 0
        alt = sample.static["alt"]
        vel = 540
        phi = inputs.static["phi"]
        theta = inputs.static["theta"]
        psi = inputs.static["psi"]

        initial_state = [vel, alpha, beta, phi, theta, psi, 0, 0, 0, 0, 0, 0, alt, power]
        step = 1.0 / 30.0
        system = autopilot.GcasAutopilot(init_mode="roll", stdout=False)
        result = f16_sim.run_f16_sim(initial_state, TSPAN[1], system, step, extended_states=True)
        states = np.vstack(
            (
                np.array([0 if x == "standby" else 1 for x in result["modes"]]),
                result["states"][:, 4],  # roll
                result["states"][:, 5],  # pitch
                result["states"][:, 6],  # yaw
                result["states"][:, 12],  # altitude
            )
        )

        return models.Trace(times=result["times"], states=np.transpose(states).tolist())

    spec = specifications.rtamt.parse_dense("always (alt >= 0)", {"alt": 0})
    optimizer = optimizers.UniformRandom[float]()
    options = staliro.TestOptions(
        runs=1,
        tspan=TSPAN,
        iterations=500,
        static_inputs={
            "phi": math.pi / 4 + np.array([-math.pi / 20, math.pi / 30]),
            "theta": -math.pi / 2 * 0.8 + np.array([0, math.pi / 20]),
            "psi": -math.pi / 4 + np.array([-math.pi / 8, math.pi / 8]),
        },
    )

    results = staliro.test(inner, spec, optimizer, options)
    result = results[0]

    return min(e.cost for e in result.evaluations)


if __name__ == "__main__":
    optimizer = optimizers.UniformRandom(min_cost=0.0)
    options = staliro.TestOptions(
        runs=1,
        iterations=100,
        static_inputs={"alt": (2200, 2500)},
    )

    results = staliro.test(outer, optimizer, options)
    result = results[0]
    min_cost = min(e.cost for e in result.evaluations)

    print(f"Minimum cost found: {min_cost}")