Solving Control Problems

Given an STL formula and a dynamical system, the STL synthesis problem is to discover a system trajectory (inputs, states, outputs) such that the output signal satisfies the specification.

This page documents our implementations of several state-of-the-art synthesis algorithms. They are categorized by the software used to define and solve the underlying optimization problem.

Quick summary:

  • DrakeMICPSolver: the standard mixed-integer convex programming (MICP) approach. Only linear systems and linear predicates are supported. Finds a globally optimal solution. This tends to be the fastest method for short-horizon specifications.

  • DrakeSos1Solver: an improved MICP approach which uses fewer binary variables. Tends to outperform the standard MICP on long-horizon specifications. Only linear systems and linear predicates are supported. Finds a globally optimal solution.

  • DrakeSmoothSolver: optimizes over a smooth approximation of the STL robustness measure. Finds a locally optimal solution. Works with nonlinear systems and predicates.

  • GurobiMICPSolver: identical to DrakeMICPSolver, but uses Gurobi’s python bindings instead of Drake’s.

  • ScipyGradientSolver: the simplest (and slowest) method. Optimizes over the (non-smooth) STL robustness measure directly using scipy.minimize. Finds a locally optimal solution. Works with nonlinear systems and predicates.

Drake

Drake is a modeling, simulation, and control toolbox for robotic systems. It includes a convienient interface for fomulating and solving optimization problems and provides bindings to numerous specialized solvers, including Gurobi, Mosek, SNOPT, Ipopt, and many others.

A complete list of supported solvers and the types of optimization problems they address can be found here.

Installation instructions for Drake can be found here.

DrakeMICPSolver

class stlpy.solvers.DrakeMICPSolver(spec, sys, x0, T, M=1000, robustness_cost=True, solver='gurobi', presolve=True, verbose=True)

Bases: DrakeSTLSolver

Given an STLFormula \(\varphi\) and a LinearSystem, solve the optimization problem

\[ \begin{align}\begin{aligned}\min & -\rho^{\varphi}(y_0,y_1,\dots,y_T) + \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\\\text{s.t. } & x_0 \text{ fixed}\\& x_{t+1} = A x_t + B u_t\\& y_{t} = C x_t + D u_t\\& \rho^{\varphi}(y_0,y_1,\dots,y_T) \geq 0\end{aligned}\end{align} \]

using mixed-integer convex programming. This gives a globally optimal solution, but may be computationally expensive for long and complex specifications.

Note

This class implements the algorithm described in

Belta C, et al. Formal methods for control synthesis: an optimization perspective. Annual Review of Control, Robotics, and Autonomous Systems, 2019. https://dx.doi.org/10.1146/annurev-control-053018-023717.

Warning

Drake must be compiled from source to support the Gurobi MICP solver. See https://drake.mit.edu/from_source.html for more details.

Drake’s naive branch-and-bound solver does not require Gurobi or Mosek, and can be used with the bnb solver option, but this tends to be very slow.

Parameters

  • spec – An STLFormula describing the specification.

  • sys – A LinearSystem describing the system dynamics.

  • x0 – A (n,1) numpy matrix describing the initial state.

  • T – A positive integer fixing the total number of timesteps \(T\).

  • M – (optional) A large positive scalar used to rewrite min and max as mixed-integer constraints. Default is 1000.

  • robustness_cost – (optional) Boolean flag for adding a linear cost to maximize the robustness measure. Default is True.

  • solver – (optional) String describing the solver to use. Must be one of ‘gurobi’, ‘mosek’, or ‘bnb’.

  • presolve – (optional) A boolean indicating whether to use gurobi’s presolve routines. Only affects the gurobi solver. Default is True.

  • verbose – (optional) A boolean indicating whether to print detailed solver info. Default is True.

AddControlBounds(u_min, u_max)

Add upper and lower bounds on the control inputs \(u_t\) to the optimization problem:

\[u_{min} \leq u_t \leq u_{max} \quad \forall t\]
Parameters

  • u_min – A (m,) numpy array specifying the minimum control input

  • u_max – A (m,) numpy array specifying the maximum control input

AddQuadraticCost(Q, R)

Add a quadratic running cost to the optimization problem:

\[\min \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\]
Parameters

  • Q – A (n,n) numpy array representing the state penalty matrix \(Q\)

  • R – A (m,m) numpy array representing the control penalty matrix \(R\)

AddRobustnessConstraint(rho_min=0.0)

Add a constraint on the robustness measure to the optimization problem:

\[\rho^{\varphi}(y_0,y_1,\dots,y_T) \geq \rho_{min}\]
Parameters

rho_min – (optional) Minimum robustness measure \(\rho_{min}\). Defaults to 0, which enforces STL satisfaction.

AddStateBounds(x_min, x_max)

Add upper and lower bounds on the state variables \(x_t\) to the optimization problem:

\[x_{min} \leq x_t \leq x_{max} \quad \forall t\]
Parameters

  • x_min – A (n,) numpy array specifying \(x_{min}\)

  • x_max – A (n,) numpy array specifying \(x_{max}\)

Solve()

Solve the STL syntheis optimization problem and return an optimal trajectory.

Return x

A (n,T) numpy array containing the optimal state \(x_t\) for each timestep.

Return u

A (m,T) numpy array containing the optimal control \(x_t\) for each timestep.

Return rho

A scalar indicating the optimal robustness value.

Return solve_time

The time it took the solver to find a solution, in seconds.

Note

x and u are returned as None if the optimization problem is infeasible or the solver is unable to find a solution.

DrakeSos1Solver

class stlpy.solvers.DrakeSos1Solver(spec, sys, x0, T, M=1000, robustness_cost=True, solver='gurobi', presolve=True, verbose=True)

Bases: DrakeMICPSolver

Given an STLFormula \(\varphi\) and a LinearSystem, solve the optimization problem

\[ \begin{align}\begin{aligned}\min & -\rho^{\varphi}(y_0,y_1,\dots,y_T) + \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\\\text{s.t. } & x_0 \text{ fixed}\\& x_{t+1} = A x_t + B u_t\\& y_{t} = C x_t + D u_t\\& \rho^{\varphi}(y_0,y_1,\dots,y_T) \geq 0\end{aligned}\end{align} \]

using mixed-integer convex programming. This method uses fewer binary variables by encoding disjunction with a Special Ordered Set of Type 1 (SOS1) constraint.

Note

This class implements the encoding described in

Kurtz V, et al. Mixed-Integer Programming for Signal Temporal Logic with Fewer Binary Variables. IEEE Control Systems Letters, 2022. https://arxiv.org/abs/2204.06367.

Warning

Drake must be compiled from source to support the Gurobi MICP solver. See https://drake.mit.edu/from_source.html for more details.

Drake’s naive branch-and-bound solver does not require Gurobi or Mosek, and can be used with the bnb solver option, but this tends to be very slow.

Parameters

  • spec – An STLFormula describing the specification.

  • sys – A LinearSystem describing the system dynamics.

  • x0 – A (n,1) numpy matrix describing the initial state.

  • T – A positive integer fixing the total number of timesteps \(T\).

  • M – (optional) A large positive scalar used to rewrite min and max as mixed-integer constraints. Default is 1000.

  • robustness_cost – (optional) Boolean flag for adding a linear cost to maximize the robustness measure. Default is True.

  • solver – (optional) String describing the solver to use. Must be one of ‘gurobi’, ‘mosek’, or ‘bnb’.

  • presolve – (optional) A boolean indicating whether to use gurobi’s presolve routines. Only affects the gurobi solver. Default is True.

  • verbose – (optional) A boolean indicating whether to print detailed solver info. Default is True.

AddControlBounds(u_min, u_max)

Add upper and lower bounds on the control inputs \(u_t\) to the optimization problem:

\[u_{min} \leq u_t \leq u_{max} \quad \forall t\]
Parameters

  • u_min – A (m,) numpy array specifying the minimum control input

  • u_max – A (m,) numpy array specifying the maximum control input

AddQuadraticCost(Q, R)

Add a quadratic running cost to the optimization problem:

\[\min \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\]
Parameters

  • Q – A (n,n) numpy array representing the state penalty matrix \(Q\)

  • R – A (m,m) numpy array representing the control penalty matrix \(R\)

AddRobustnessConstraint(rho_min=0.0)

Add a constraint on the robustness measure to the optimization problem:

\[\rho^{\varphi}(y_0,y_1,\dots,y_T) \geq \rho_{min}\]
Parameters

rho_min – (optional) Minimum robustness measure \(\rho_{min}\). Defaults to 0, which enforces STL satisfaction.

AddStateBounds(x_min, x_max)

Add upper and lower bounds on the state variables \(x_t\) to the optimization problem:

\[x_{min} \leq x_t \leq x_{max} \quad \forall t\]
Parameters

  • x_min – A (n,) numpy array specifying \(x_{min}\)

  • x_max – A (n,) numpy array specifying \(x_{max}\)

Solve()

Solve the STL syntheis optimization problem and return an optimal trajectory.

Return x

A (n,T) numpy array containing the optimal state \(x_t\) for each timestep.

Return u

A (m,T) numpy array containing the optimal control \(x_t\) for each timestep.

Return rho

A scalar indicating the optimal robustness value.

Return solve_time

The time it took the solver to find a solution, in seconds.

Note

x and u are returned as None if the optimization problem is infeasible or the solver is unable to find a solution.

DrakeSmoothSolver

class stlpy.solvers.DrakeSmoothSolver(spec, sys, x0, T, k=2.0, verbose=True)

Bases: DrakeSTLSolver

Given an STLFormula \(\varphi\) and a LinearSystem, solve the optimization problem

\[ \begin{align}\begin{aligned}\max ~& \rho^{\varphi}(y_0,y_1,\dots,y_T)\\\text{s.t. } & x_0 \text{ fixed}\\& x_{t+1} = f(x_t, u_t)\\& y_{t} = g(x_t, u_t)\end{aligned}\end{align} \]

where \(\rho^{\varphi}\) is defined using a smooth approximation of min and max. We then use one of the general nonlinear solvers availible in Drake to find a locally optimal solution.

Note

This class implements the algorithm described in

Gilpin, Y, et al. A Smooth Robustness Measure of Signal Temporal Logic for Symbolic Control. IEEE Control Systems Letters, 2021. https://arxiv.org/abs/2006.05239.

Note

This method is most effective when used in conjunction with the SNOPT sparse SQP solver. SNOPT is included with the binary version of Drake, but requires a license when built from source.

Parameters

  • spec – An STLFormula describing the specification.

  • sys – A LinearSystem describing the system dynamics.

  • x0 – A (n,1) numpy matrix describing the initial state.

  • T – A positive integer fixing the total number of timesteps \(T\).

  • k – (optional). A smoothing parameter characterizing the tightness of the smooth approximation. Larger values give a tighter approximation.

  • verbose – (optional) A boolean indicating whether to print detailed solver info. Default is True.

AddControlBounds(u_min, u_max)

Add upper and lower bounds on the control inputs \(u_t\) to the optimization problem:

\[u_{min} \leq u_t \leq u_{max} \quad \forall t\]
Parameters

  • u_min – A (m,) numpy array specifying the minimum control input

  • u_max – A (m,) numpy array specifying the maximum control input

AddQuadraticCost(Q, R)

Add a quadratic running cost to the optimization problem:

\[\min \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\]
Parameters

  • Q – A (n,n) numpy array representing the state penalty matrix \(Q\)

  • R – A (m,m) numpy array representing the control penalty matrix \(R\)

AddRobustnessConstraint(rho_min=0.0)

Add a constraint on the robustness measure to the optimization problem:

\[\rho^{\varphi}(y_0,y_1,\dots,y_T) \geq \rho_{min}\]
Parameters

rho_min – (optional) Minimum robustness measure \(\rho_{min}\). Defaults to 0, which enforces STL satisfaction.

AddStateBounds(x_min, x_max)

Add upper and lower bounds on the state variables \(x_t\) to the optimization problem:

\[x_{min} \leq x_t \leq x_{max} \quad \forall t\]
Parameters

  • x_min – A (n,) numpy array specifying \(x_{min}\)

  • x_max – A (n,) numpy array specifying \(x_{max}\)

Solve()

Solve the STL syntheis optimization problem and return an optimal trajectory.

Return x

A (n,T) numpy array containing the optimal state \(x_t\) for each timestep.

Return u

A (m,T) numpy array containing the optimal control \(x_t\) for each timestep.

Return rho

A scalar indicating the optimal robustness value.

Return solve_time

The time it took the solver to find a solution, in seconds.

Note

x and u are returned as None if the optimization problem is infeasible or the solver is unable to find a solution.

Gurobi

Solvers in this section use the python bindings of Gurobi, a commercial optimizer that handles a variety of problem classes, including convex programs, mixed-integer programs, and non-convex quadratic programs with quadratic constraints. Free academic licenses are available.

GurobiMICPSolver

class stlpy.solvers.GurobiMICPSolver(spec, sys, x0, T, M=1000, robustness_cost=True, presolve=True, verbose=True)

Bases: STLSolver

Given an STLFormula \(\varphi\) and a LinearSystem, solve the optimization problem

\[ \begin{align}\begin{aligned}\min & -\rho^{\varphi}(y_0,y_1,\dots,y_T) + \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\\\text{s.t. } & x_0 \text{ fixed}\\& x_{t+1} = A x_t + B u_t\\& y_{t} = C x_t + D u_t\\& \rho^{\varphi}(y_0,y_1,\dots,y_T) \geq 0\end{aligned}\end{align} \]

with Gurobi using mixed-integer convex programming. This gives a globally optimal solution, but may be computationally expensive for long and complex specifications.

Note

This class implements the algorithm described in

Belta C, et al. Formal methods for control synthesis: an optimization perspective. Annual Review of Control, Robotics, and Autonomous Systems, 2019. https://dx.doi.org/10.1146/annurev-control-053018-023717.

Parameters

  • spec – An STLFormula describing the specification.

  • sys – A LinearSystem describing the system dynamics.

  • x0 – A (n,1) numpy matrix describing the initial state.

  • T – A positive integer fixing the total number of timesteps \(T\).

  • M – (optional) A large positive scalar used to rewrite min and max as mixed-integer constraints. Default is 1000.

  • robustness_cost – (optional) Boolean flag for adding a linear cost to maximize the robustness measure. Default is True.

  • presolve – (optional) A boolean indicating whether to use Gurobi’s presolve routines. Default is True.

  • verbose – (optional) A boolean indicating whether to print detailed solver info. Default is True.

AddControlBounds(u_min, u_max)

Add upper and lower bounds on the control inputs \(u_t\) to the optimization problem:

\[u_{min} \leq u_t \leq u_{max} \quad \forall t\]
Parameters

  • u_min – A (m,) numpy array specifying the minimum control input

  • u_max – A (m,) numpy array specifying the maximum control input

AddQuadraticCost(Q, R)

Add a quadratic running cost to the optimization problem:

\[\min \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\]
Parameters

  • Q – A (n,n) numpy array representing the state penalty matrix \(Q\)

  • R – A (m,m) numpy array representing the control penalty matrix \(R\)

AddRobustnessConstraint(rho_min=0.0)

Add a constraint on the robustness measure to the optimization problem:

\[\rho^{\varphi}(y_0,y_1,\dots,y_T) \geq \rho_{min}\]
Parameters

rho_min – (optional) Minimum robustness measure \(\rho_{min}\). Defaults to 0, which enforces STL satisfaction.

AddStateBounds(x_min, x_max)

Add upper and lower bounds on the state variables \(x_t\) to the optimization problem:

\[x_{min} \leq x_t \leq x_{max} \quad \forall t\]
Parameters

  • x_min – A (n,) numpy array specifying \(x_{min}\)

  • x_max – A (n,) numpy array specifying \(x_{max}\)

Solve()

Solve the STL syntheis optimization problem and return an optimal trajectory.

Return x

A (n,T) numpy array containing the optimal state \(x_t\) for each timestep.

Return u

A (m,T) numpy array containing the optimal control \(x_t\) for each timestep.

Return rho

A scalar indicating the optimal robustness value.

Return solve_time

The time it took the solver to find a solution, in seconds.

Note

x and u are returned as None if the optimization problem is infeasible or the solver is unable to find a solution.

Scipy

Solvers in this section are based on the (relatively simple) scipy minimize optimizer, which can be installed with pip:

pip install scipy

Alternative installation instructions can be found here.

ScipyGradientSolver

class stlpy.solvers.ScipyGradientSolver(spec, sys, x0, T, method='slsqp', verbose=True)

Bases: STLSolver

Given an STLFormula \(\varphi\) and a NonlinearSystem, solve the optimization problem

\[ \begin{align}\begin{aligned}\min & - \rho^{\varphi}(y_0,y_1,\dots,y_T) + \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\\\text{s.t. } & x_0 \text{ fixed}\\& x_{t+1} = f(x_t, u_t)\\& y_{t} = g(x_t, u_t)\end{aligned}\end{align} \]

using a shooting method and the scipy.optimize solver.

Warning

This solver uses finite-differences to approximate the gradient of the (non-smooth) cost. As such, this method is likely to scale extremely poorly.

Parameters

  • spec – An STLFormula describing the specification.

  • sys – A NonlinearSystem describing the system dynamics.

  • x0 – A (n,1) numpy matrix describing the initial state.

  • T – A positive integer fixing the total number of timesteps \(T\).

  • method – (optional) String characterizing the optimization algorithm to use. See the scipy docs for more details. Default is Sequential Least Squares ("slsqp").

  • verbose – (optional) A boolean indicating whether to print detailed solver info. Default is True.

AddQuadraticCost(Q, R)

Add a quadratic running cost to the optimization problem:

\[\min \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\]
Parameters

  • Q – A (n,n) numpy array representing the state penalty matrix \(Q\)

  • R – A (m,m) numpy array representing the control penalty matrix \(R\)

Solve()

Solve the STL syntheis optimization problem and return an optimal trajectory.

Return x

A (n,T) numpy array containing the optimal state \(x_t\) for each timestep.

Return u

A (m,T) numpy array containing the optimal control \(x_t\) for each timestep.

Return rho

A scalar indicating the optimal robustness value.

Return solve_time

The time it took the solver to find a solution, in seconds.

Note

x and u are returned as None if the optimization problem is infeasible or the solver is unable to find a solution.

Write Your Own Solver

All the solvers described above inherit from the STLSolver class. To implement your own optimization-based solver, all you need to do is create a new class that inherits from STLSolver.

class stlpy.solvers.base.STLSolver(spec, sys, x0, T, verbose)

Bases: ABC

A simple abstract base class defining a common solver interface for different optimization-based STL synthesis methods.

This class considers variations on the trajectory synthesis problem

\[ \begin{align}\begin{aligned}\min & \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\\\text{s.t. } & x_0 \text{ fixed }\\& x_{t+1} = f(x_t, u_t)\\& y_t = g(x_t, u_t)\\& \rho^{\varphi}(y_0,y_1,\dots,y_T) \geq 0\end{aligned}\end{align} \]

where \(Q \succeq 0\) and \(R \succeq 0\) are cost weights, \(f\) and \(g\) define the system dynamics, and \(\rho\) is the robustness measure associated with a given STL specification \(\varphi\).

Possible variations include using the robustness measure \(\rho^\varphi\) as a cost, dropping the quadratic running cost, and removing the hard satisfaction constriant \(\rho^{\varphi}\geq 0\).

Parameters

  • spec – An STLFormula describing the specification.

  • sys – An NonlinearSystem characterizing the system dynamics.

  • x0 – A (n,1) numpy array representing the initial state \(x_0\).

  • T – A positive integer fixing the total number of timesteps \(T\).

  • verbose – A boolean specifying whether to print detailed solver info.

abstract AddControlBounds(u_min, u_max)

Add upper and lower bounds on the control inputs \(u_t\) to the optimization problem:

\[u_{min} \leq u_t \leq u_{max} \quad \forall t\]
Parameters

  • u_min – A (m,) numpy array specifying the minimum control input

  • u_max – A (m,) numpy array specifying the maximum control input

abstract AddDynamicsConstraints()

Add the dynamics constraints

\[ \begin{align}\begin{aligned}& x_0 \text{ fixed }\\& x_{t+1} = f(x_t, u_t)\\& y_t = g(x_t, u_t)\end{aligned}\end{align} \]

to the optimization problem.

abstract AddQuadraticCost(Q, R)

Add a quadratic running cost to the optimization problem:

\[\min \sum_{t=0}^T x_t^TQx_t + u_t^TRu_t\]
Parameters

  • Q – A (n,n) numpy array representing the state penalty matrix \(Q\)

  • R – A (m,m) numpy array representing the control penalty matrix \(R\)

abstract AddRobustnessConstraint(rho_min=0.0)

Add a constraint on the robustness measure to the optimization problem:

\[\rho^{\varphi}(y_0,y_1,\dots,y_T) \geq \rho_{min}\]
Parameters

rho_min – (optional) Minimum robustness measure \(\rho_{min}\). Defaults to 0, which enforces STL satisfaction.

abstract AddRobustnessCost()

Add the robustness measure as a (linear) cost to the optimization problem:

\[\min -\rho^{\varphi}(y_0,y_1,\dots,y_T).\]
abstract AddSTLConstraints()

Add constraints to the optimization problem to define the robustness measure

\[\rho^{\varphi}(y_0,y_1,\dots,y_T).\]
abstract AddStateBounds(x_min, x_max)

Add upper and lower bounds on the state variables \(x_t\) to the optimization problem:

\[x_{min} \leq x_t \leq x_{max} \quad \forall t\]
Parameters

  • x_min – A (n,) numpy array specifying \(x_{min}\)

  • x_max – A (n,) numpy array specifying \(x_{max}\)

abstract Solve()

Solve the STL syntheis optimization problem and return an optimal trajectory.

Return x

A (n,T) numpy array containing the optimal state \(x_t\) for each timestep.

Return u

A (m,T) numpy array containing the optimal control \(x_t\) for each timestep.

Return rho

A scalar indicating the optimal robustness value.

Return solve_time

The time it took the solver to find a solution, in seconds.

Note

x and u are returned as None if the optimization problem is infeasible or the solver is unable to find a solution.