The Control Toolbox - An Open-Source C++ Library for Robotics, Optimal and Model Predictive Control

01/12/2018 ∙ by Markus Giftthaler, et al. ∙ ETH Zurich 0

We introduce the Control Toolbox (CT), an open-source C++ library for efficient modelling, control, estimation, trajectory optimization and model predictive control. The CT is applicable to a broad class of dynamic systems, but features additional modelling tools specially designed for robotics. This paper outlines its general concept, its major building blocks and highlights selected application examples. The CT was designed for intuitive modelling of systems governed by ordinary differential- or difference equations. It supports rapid prototyping of cost functions and constraints and provides common interfaces for different optimal control solvers. To date, we support Single Shooting, the iterative Linear-Quadratic Regulator, Gauss-Newton Multiple Shooting and classical Direct Multiple Shooting. We provide interfaces to different NLP and linear-quadratic solvers, such as IPOPT, SNOPT, HPIPM, or a custom Riccati solver. The CT was designed with performance for online control in mind and allows to solve large-scale optimal control problems highly efficiently. Some of the key features enabling fast run-time performance are full support for Automatic Differentiation, derivative code generation and thorough multi-threading. For robotics problems, the we offer an interface to a fully auto-differentiable rigid-body dynamics modelling engine. In combination with derivative code generation, this allows for an unprecedented performance in solving optimal control problems for complex articulated robotic systems.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

I-a What is the Control Toolbox?

A common task for robotics researchers and practitioners is to model systems, implement equations of motion and design model-based controllers, estimators, planning algorithms, etc. Sooner or later, one is confronted with questions of efficient implementation, computing derivatives, formulating cost functions and constraints or running controllers in a model-predictive control fashion.

The Control Toolbox is specially designed for these tasks. It is written entirely in C++ and has a strong focus on highly efficient code that can be run online (in the loop) on robots or other actuated hardware. A significant contribution of the CT is its implementation of optimal control algorithms, spanning a range from simple LQR reference implementations to constrained Model Predictive Control. The CT supports Automatic Differentiation (Auto-Diff) and allows to generate derivative code for arbitrary scalar and vector-valued functions. The toolbox was designed with usability in mind, allowing users to apply advanced concepts such as Nonlinear Model Predictive Control (NMPC) or numerical optimal control quickly and with minimal effort. In contrast to highly integrated frameworks, the CT follows a modular library approach. Thus, it can be easily interfaced with external modeling and solver frameworks. Additionally, building blocks such as Automatic-Differentiation are usable in several applications including optimal control, classical feedback control (e.g. LQRs), Kalman filtering or sensitivity analysis. In summary, several elements of a planning and control pipeline for an application such as joint space PID control, inverse dynamics control and motion planning can be developed using a single library.

The CT has been designed to provide the tools needed for fast development and evaluation of control methods while being optimized for efficiency and allowing for online operation. While the emphasis lies on control, the provided tools can also be used for simulation, estimation or other optimization applications.

There are four key components of the control toolbox: Modelling of continuous and discrete-time dynamical systems, Automatic Differentiation, optimal control algorithms and Rigid Body Dynamics algorithms. For many of these components, we have carefully selected existing implementations and provide a seamless integration between them without creating a rigid framework. Instead, CT offers easy to use tools for fast prototyping, such as numerical integrators or LQR design. More complex approaches are provided in the form of reference implementations, such as generating a whole-body Nonlinear Model Predictive Control setup for a robot based solely on a semantic description.

I-B Related Work

When looking at software tools in robotics, the library that matches the scope of CT the closest is Drake [1]. However, Drake started out as a Matlab implementation, which made it unsuitable for hard-realtime and online control. While its codebase is gradually moving towards C++, not all features have been ported yet and Auto-Diff support is limited. Another popular software library is MuJoCo [2], which excels at simulation but follows a closed-source policy.

There are many optimal control toolboxes outside of the robotics community which focus on transcribing and solving nonlinear optimal control problems and influenced the development of CT. Notable examples are ACADO [3], its successor ACADOS [4], PSOPT [5], the closed-source toolbox MUSCOD [6] as well as the commercial tools GPOPS [7] and ForcesPro [8]. Broadly speaking, there are three categories of solvers frequently used in optimal control: linear-quadratic optimal control solvers such as HPIPM [9], general quadratic programming solvers such as qpOases [10], and general-purpose nonlinear-programming packages, such as IPOPT [11], SNOPT [12] or NLopt [13]. Developing such solvers is a research field by itself and not the primary scope of the CT. Therefore, we provide an interface to SNOPT, IPOPT, and HPIPM as well as custom implementations of iLQR [14] and Gauss-Newton Multiple Shooting [15]. This gives the user the opportunity to evaluate different solver types and implementations.

For Automatic Differentiation, we rely on CppAD paired with the code-generation framework CppADCodeGen [16]. For a more detailed review on Auto-Diff frameworks as well as why symbolic differentiation as e.g. available in Matlab [17], Mathematica [18] or MapleSim [19] compares unfavorably to Auto-Diff regarding speed, we refer to [20].

For modelling rigid body dynamics, there exists a variety of mature libraries, with notable examples being RBDL [21], DART [22], Pinocchio [23], iDynTree [24] and RobCoGen [25]. Due to its modularity, CT can be interfaced with either of these libraries. We decided to provide a reference interface to RobCoGen due to its wide Auto-Diff support and suitability for online control. Support and interfaces to other rigid body dynamics libraries are going to be added to CT once their Auto-Diff support improves.

I-C Scope

For control, especially numerical optimal control in a robotics context, there are many individual libraries available that provide key ingredients such as modeling frameworks, Auto-Diff, and optimal control solvers. However, due to different data representations, modeling assumptions, basic conventions and the lack of reference implementations in a robotics context, integrating these components is a tedious, time-consuming and error-prone process. Also, for researchers entering the field, the correct choice of modeling framework or solver is difficult to make since the solvers’ scope and functionality differ strongly. Furthermore, their performance highly depends on the specific modeling implementation.

For this reason, the CT aims at providing users with the tools to quickly implement well-established control methods and combine their problem with different solvers, transcription methods, and modeling approaches. Since the CT is open-source, it is easy to adapt individual components to specific use cases or integrate custom modeling or solver frameworks. On top of this integration, the CT follows a holistic approach to robot control: the individual components are not only useful for numerical optimal control but can also be employed for classical feedback control, inverse dynamics control, estimation and planning. The CT provides features such as Kalman filtering or contact constraint projection for Rigid Body dynamics [26]. Thus, the CT can be used for several aspects of a robotics control and planning toolchain: from low-level realtime closed-loop control to kinematic planning or dynamic whole-body trajectory optimization.

I-D Structure of This Paper

This paper is structured as follows. In Section II, we present an overview of the CT’s design and implementation and give an outline of its structure. The different main modules of the CT are highlighted in Sections III to VI. Selected application examples are given in Section VII. For real-time applications, optimizing runtime performance is an important issue, on which we comment in Section VIII. The paper is concluded by important links and licence information (Section IX) and acknowledgements to additional contributors in Section X.

Ii Overview

Ii-a Fundamental Dependencies

The CT is written in C++ and has been tested under Ubuntu 14.04 and 16.04 with library versions as provided in the package sources. Building the CT requires a C++ compiler with C++11 support. Since the CT is designed as a toolbox rather than an integrated application, we tried to provide maximum flexibility to the users. Therefore, it is not tied to a specific middleware such as ROS and dependencies are kept at a minimum. The two essential dependencies for CT are Eigen [27] and kindr [28] (which is based on Eigen). Eigen is a popular library for linear algebra in C++ and provides efficient implementations of standard matrix operations as well as more advanced linear algebra methods. Kindr is a header only kinematics library which builds on top of it and provides data types for different rotation representations such as quaternions, Euler angles or rotation matrices.

Ii-B Structure and Modules of the CT

The Control Toolbox consists of three main modules. The core module (ct_core), the optimal control module (ct_optcon) and the rigid body dynamics module (ct_rbd). There is a clear hierarchy between the modules, which means the modules depend on each other in this order. For example, one can use the core module without ct_optcon or ct_rbd.

  • ct_core provides general type definitions and mathematical tools. For example, it contains most data type definitions, definitions for systems and controllers, as well as basic functionality such as numerical integrators for differential equations.

  • ct_optcon builds on top of the ‘core’ module and adds infrastructure for defining and solving optimal control problems. It contains the functionality for defining cost functions, constraints, solver backends and a generic NMPC wrapper.

  • ct_rbd provides tools for modelling rigid body dynamics systems and interfaces with ct_core and ct_optcon.

For testing as well as for giving examples, we provide a fourth module: the ‘models’ module (ct_models) contains various robot models including a quadruped, a robotic arm, a normal quadrotor and a quadrotor with a slung load. These four different modules are detailed in Sections III-VI.

Iii Core Module

Iii-a Basic System Definitions

The core module defines basic data types and interfaces to describe non-linear system dynamics of the forms


The right-hand side (RHS) of Equation (1) only depends on the time and state (core::StateVector) and is called a core::System. The RHS of Equation (2) additionally depends on the control input (core::ControlVector) and is named core::ControlledSystem. The dynamics equations can be implemented by the user in any desired way but are currently restricted to ordinary ODEs and difference equations. For the remainder of this section, we limit the scope to a continuous-time perspective. Note that for modeling robotic systems in continuous-time, the rigid body dynamics module provides a variety of tools, which are detailed in Section V.

As the name suggests, the core::ControlledSystem provides the interface for closing a feedback control loop. Every controlled system can take a pointer to a control law deriving from core::Controller. Full flexibility for implementing a policy of general form is given to the user. This includes special cases where the control is merely constant, depending on neither nor , time-varying, only depending on , or a general feedback controller, depending on both and t.

We provide a set of pre-defined control laws, which includes a core::ConstantController with fixed , a classical PID controller (core::PIDController), or a full time-varying core::StateFeedbackController with feedforward term of form .

Iii-B Integration and Simulation

The CT provides different numerical integrators (core::Integrator). We offer own implementations and integrators based on ‘boost odeint’ [29].

The CT currently features fixed-step integrators like Euler and fourth-order Runge-Kutta as well as different (error controlled) variable step integrators. Additionally, for symplectic systems (core::SymplecticSystem) a semi-implicit Euler integrator (core::SymplecticIntegrator) is available, which can help with stiff systems. All integrators take a pointer to a system and return trajectories (core::DiscreteTrajectory), i.e. timed series of states and control inputs (core::StateTrajectory and core::ControlTrajectory

) respectively. These trajectories can be either equidistant in time or unevenly sampled. In both cases, an interpolation strategy can be applied to obtain states and inputs at a specific time which is not directly stored.

For rapid prototyping and testing of control loops, we provide a core::ControlSimulator which allows running controllers and system integration in parallel and in real-time. Please note, however, that the CT cannot replace a high-fidelity physics simulator. For such purposes, we refer for example to [30].

Iii-C Computing Derivatives

Derivative Numerical Computation Setup Error
method Accuracy Speed Time Safety
Analytic Deriv.
Symbolic Engine
Auto-diff Codegen
TABLE I: Comparing different options to obtain derivatives

The CT can be used to compute derivatives of arbitrary vector-valued smooth nonlinear functions . For computing first order derivatives (Jacobians) , the most-widespread methods are

  1. Numerical differentiation, e.g. by the method of finite-differences,

  2. Analytical derivation, e.g. performed manually,

  3. Symbolic math engines,

  4. Automatic Differentiation, also known as Algorithmic Differentiation, with an optional source code generation step.

The different approaches are compared in Table I. Automatic Differentiation allows to conveniently obtain derivative information: it relieves the user from computing analytical derivatives manually or symbolically111Auto-Diff uses graph structures to compute derivatives. Hence, it is inherently different from symbolic engines such as Maple or Maxima., which may be intractable for complex systems. However, it is as accurate and fast as analytic derivatives and outperforms numerical differentiation in terms of accuracy and speed while providing a similar level of convenience. Combining Automatic Differentiation with source code generation (Auto-Diff Codegen) results in the best runtime of the differentiation methods supported by CT. For a detailed review and numerical examples, the interested reader is referred to [20].

Iii-D Linearizing Dynamic Systems

The CT defines the structure of a linear system (core::LinearSystem) as


where and are the Jacobians of a non-linear, time-varying system evaluated at desired setpoints for and . In order to compute this linearization for a non-linear system, CT provides two different helper classes. The core::SystemLinearizer takes a core::ControlledSystem and applies numerical differentiation to compute the Jacobians. Alternatively, the core::AutoDiffLinearizer can be used to apply Auto-Differentiation for the Jacobians which is more accurate than numerical differentiation. Finally, Auto-Differentiation is combined with code generation in core::ADCodegenLinearizer which is as accurate as analytical derivatives and typically fast to evaluate. The code-gen linearizer employs a technique called just-in time compilation (JIT), which generates the derivative code at runtime. Since this can take a few seconds, the derivative code can be stored to file and compiled in separate libraries. Examples for this approach are given in ct_models, see Section VI.

Iii-E Computing Approximated and Exact Sensitivities

Many control algorithms, for example the direct approaches to optimal control shown in Section IV, require a discrete-time approximation of the nonlinear system dynamics of form , where we call and ‘sensitivities’. In many cases it may suffice to approximate these matrices based on the continuous-time counterparts , and a simple Forward-Euler, Backward-Euler or Tustin discretization scheme. The CT provides the core::SensitivityApproximation class, which can be used to compute such low-order approximations in a straight-forward way.

However, especially when aiming at a coarse time-discretization while dealing with a highly nonlinear dynamic systems, it can be beneficial to use higher-order integration schemes to compute , . The core::SensitivityIntegrator solves the integrals

for a given starting time and time-step by means of integrating a Sensitivity ODE. Special cases for obtaining exact sensitivities for symplectic integration schemes are included, too. Exact sensitivities can help to robustify and improve the convergence behavior of many optimal control algorithms in the CT, which are summarized below.

Iv Optimal Control Module

A broad variety of model-based optimal control tasks can be formulated as continuous-time optimal control problems. From a robotics perspective, this includes tasks such as agile flight, reaching an object in a cluttered scene, moving a mobile manipulator or quadrupedal locomotion. In direct optimal control, the continuous-time optimal control problem is first transcribed into a numerically tractable discrete problem. Two possible ways to complete this step are:

  1. Transcribing the problem into a nonlinear program (NLP) using multiple-shooting, single shooting or direct collocation and subsequently solving it using standard NLP solvers such as IPOPT or SNOPT.

  2. Using iterative Riccati-based shooting methods derived from the Principle of Optimality such as DDP [31], their Gauss-Newton counterparts, iLQR [14] or Gaus-Newton Multiple Shooting (GNMS) [15]. These methods are popular due to their overall efficiency and linear time complexity.

The package ct_optcon covers both classical off-the-shelf NLP solvers and custom Riccati-based solutions, paired with different flavors of Single and Multiple Shooting. An important design feature is the CT’s modularity, which allows combining different cost functions, dynamics, constraints, and solvers in an almost arbitrary way and therefore allows for rapid prototyping of optimal control setups, including Nonlinear Model Predictive Control.

Iv-a Cost Functions

The cost function package provides means of quickly prototyping objective functions based on a highly modular approach. A CT cost function is assumed to consist of a sum of elementary cost function building blocks, which are called ‘terms’. Each term evaluates to a scalar as a function of the current time, control input, and state and derives from optcon::TermBase.

The overall cost function is designed such that it holds intermediate terms and final terms, which can be assigned individually. The intermediate and final costs are then given as the sums over the evaluations of all intermediate and final terms. Equivalently, the intermediate and final derivatives result as the sums of the individual intermediate and final term gradients. The cost function package supports both analytic derivatives for terms as well as Automatic Differentiation and just-in-time compilation (JIT) up to second order derivatives.

We offer a selection of frequently used standard cost function terms, which penalize the deviations from given control and state reference points, including a purely quadratic term (optcon::TermQuadratic), a cross-term (optcon::TermMixed) and a purely linear term (optcon::TermLinear). Furthermore there are terms for tracking reference trajectories in state and control (optcon::TermQuadTracking) and terms which formulate soft constraints on state and control variables (optcon::TermStateBarrier).

All existing terms can be automatically constructed from text-files, in which the cost function weights and parameters can be structured in a simple manner. For custom terms, reading from a file is simple to implement thanks to a pre-specified set of loading methods. Additionally, all terms can be made-time-varying using time-activation functions, which can be used to introduce way-point costs, for example.

Iv-B Constraints

The constraint package generalizes the modular idea presented for cost functions in Section IV-A to vector-valued functions. The corresponding elementary building blocks derive from optcon::ConstraintBase and again support both analytic derivatives, Automatic-Differentiation and Auto-Diff with JIT. For constraints, the terms are not summarized but stacked in a so-called ‘constraint container’. Every container additionally features an upper and a lower bound. For constraints, we currently only support first-order derivatives (optcon::LinearConstraintContainer). To date, the predefined terms include simple linear path inequality constraints and box constraints on states and controls.

Iv-C Optimal Control Problem Containers

A optcon::OptConProblem is a unified container for nonlinear controlled system dynamics, Equation (2), nonlinear cost functions, nonlinear constraints, a time horizon variable and an initial state. It serves as the main interface between a user and the different implementations of optimal control algorithms and NMPC.

Similarly, the container optcon::LQOCProblem is dedicated to constrained linear-quadratic optimal control problems. However, this container is designed to directly store the linearized dynamics, the Jacobians and Hessians of the cost function and the constraint Jacobians in matrix representation.

Iv-D LQR and Linear Quadratic Solvers

The CT provides C++ code for different variants of the classical Linear Quadratic Regulator. We provide direct and iterative solvers for the continuous-time Algebraic Riccati Equation (optcon::CARE), and iterative solvers for the discrete-time Algebraic Riccati Equation (optcon::DARE). Those can be used to design infinite-horizon LQR controllers and state- and disturbance estimators in both continuous- and discrete time. Furthermore, there is a time-varying, finite-horizon discrete-time LQR version available (optcon::FHDTLQR).

For unconstrained linear-quadratic optimal control problems, the CT offers a custom Riccati solver, optcon::GNRiccatiSolver, which achieves high efficiency using advanced options such as fixed Hessian regularization.

For constrained LQ optimal control problems the CT includes an interface to the interior point solver HPIPM [9], which is a competitive solver for constrained optimal control problems: it features linear time complexity thanks to a Riccati factorization and uses a linear algebra package with CPU architecture specific optimization [32].

Iv-E NLP Problems and Solvers

A unified, Eigen-based interface for formulating nonlinear programming problems (optcon::Nlp) and solving them (optcon::NlpSolver) is part of the CT. To date, we provide interfaces to the free interior-point solver IPOPT [11] and the commercial SQP-solver SNOPT [12].

Iv-F Gauss-Newton Shooting Algorithms with Riccati solvers

The CT implements a family of Gauss-Newton Multiple Shooting algorithms in both unconstrained and constrained fashion [15]. This family of algorithms performs Sequential Quadratic Programming on the original nonlinear optimal control problem, uses appropriate Riccati solvers to solve linear-quadratic sub-problems efficiently, and utilizes a line-search over a merit function for globalization. The algorithms employ a piece-wise constant control parameterization. A famous limit case of the family of algorithms is the iterative Linear Quadratic Regulator (iLQR). The details of these algorithms have been extensively covered elsewhere [14, 15]. However, we note that the CT shows how to integrate these algorithms in a single framework at almost identical computational cost. These algorithms are particularly powerful for unconstrained problems with long time horizons or very fine control discretizations. Additionally, at every iteration, they design a time-varying state-feedback control law, which generalizes the policy in the vicinity of the optimal solution.

Iv-G Classical Direct Multiple Shooting

Complementary to GNMS, the CT also implements the original Direct Multiple Shooting (DMS) method by Bock and Plitt [33], which we solve using a classical NLP solver (see Section IV-E. We provide this method separately since it complements the other algorithms in several aspects. While GNMS currently only supports a constant control parameterization, DMS also supports linear interpolation. DMS in combination with IPOPT can furthermore leverage exact Hessians or other Hessian approximations. DMS furthermore supports adaptive step-size integration. Lastly, DMS can make use of more advanced globalization techniques as employed by the NLP solvers, such as complex filter schemes [34]. However, for problems with long time horizons, the DMS implementation cannot compete with GNMS or iLQR at runtime, due to computational limitations of the currently available off-the-shelf NLP solvers.

Iv-H Nonlinear Model Predictive Control

Thanks to a dedicated design of interfaces between solvers and the optimal control problem definition, the CT optimal control problem solvers can be automatically run in Nonlinear Model Predictive Control fashion using the class optcon::MPC. The latter offers options like automatic warm-starting, pre-integration for delay-compensation, different modes to handle time horizons (e.g. receding horizon, fixed time horizon) and offers explicit support for real-time iteration schemes [35]. For a detailed example of NMPC using a GNMS nonlinear optimal control solver, the reader is referred to the ct_optcon online tutorial.

V Rigid Body Dynamics Module

Generally speaking, the main task of the rigid body dynamics module ct_rbd is to provide wrappers that map specialized RBD code into a general ordinary differential equation of form (

2), cost functions, and constraints.

The rigid body dynamics module currently relies on RobCoGen [25], a code-generation framework for rigid body dynamics and kinematics. To model a new robot in the Control Toolbox, an additional code-generation step is required to create the dynamics and kinematics equations based on a user-provided semantic robot description. To date, RobCoGen is the only C++ rigid body dynamics engine that supports Automatic Differentiation, a major ingredient for our NMPC applications, see e.g. [36].

Generating a new robot model based on the code-generation output of RobCoGen is straight-forward and amounts to creating a single header file with only a few lines of code. Essentially, one needs to specify kinematic branches and end-effector locations. In the background, ct_rbd creates containers and wrappers which allow convenient access to the generated robot dynamics and kinematics functions as well as force-transforms and Jacobians. For fixed-base systems, the dynamics container is the class rbd::FixBaseFDSystem, for floating-base systems it is rbd::FloatingBaseFDSystem. The floating-base state is

where and define base orientation and position expressed in the inertial (‘world’) frame. and represent local angular and linear velocity expressed in a body fixed frame. Joint angles and velocities are represented by and , respectively.

For a straight-forward application of nonlinear optimal control to robotic systems, ct_rbd offers wrapper classes which allow running nonlinear optimal control algorithms for any rigid-body dynamics model. As an easy way to handle contacts on arbitrary reference frames, we currently support a soft spring-damper contact model, which is described in detail in [20]. Alternatively, contact forces can be chosen as additional control inputs. Furthermore, the CT allows to

  • generate operational-space models from the generated dynamics equations,

  • augment rigid-body dynamic systems with arbitrary user-defined actuator dynamics models,

  • use a number of pre-defined standard controllers such as joint position controllers plus inverse dynamics,

  • use predefined cost function terms that are specific to robotic systems, e.g. we define auto-differentiable cost function terms for end-effector task-space positioning.

Lastly, we provide an interface for solving inverse kinematics problems using IKFast [37].

Vi Models Module

The models module, ct_models, contains a collection of fix- and floating base robot models which serve as examples of how to include systems in different ways:

  • the quadrotor is a floating-base system which is modeled independent from ct_rbd and can serve as an example of how to implement a system which derives directly from core::ControlledSystem.

  • the inverted pendulum is the simplest system to be modeled using RobCoGen: a fix-base robot with 1 DoF.

  • ‘HyA’ models the fix-base, 6 DoF robot arm from [38].

  • ‘HyQ’ [39] is a quadrupedal robot with 18 DoF.

  • the quadrotor with slung-load is modeled with RobCoGen. It is an example of how to adapt rbd::FloatingBaseFDSystem for robots with unusual actuation.

For systems modeled using RobCoGen, ct_models contains the generated dynamics code. ct_models also gives examples of how to compile derivative code for forward and inverse dynamics into a separately loadable library.

Vii Application Examples

The CT has been validated in a number of projects, including many hardware experiments, demonstrations, and academic publications. The following presents a compact summary. For details, we refer the interested reader to the referenced papers. CT application examples include

  • NMPC for a hexrotor flying through a window222 [40],

  • a quadrotor with rotor failure performing a go-to task333,

  • Trajectory Optimization and full-body Nonlinear Model Predictive Control on different quadruped robots, including performing agile squat jumps444 [41, 36],

  • online trajectory optimization with collision avoidance [20] on a 6 DoF industrial robot arm,

  • the computation of derivatives of constraints and cost functions including complex kinematic chains was demonstrated in hardware experiments in 555 [42].

  • pick-and-place arm motions for mobile manipulator were demonstrated in [43].

In many of the above examples, our solvers reason about full rigid body dynamics models which are not simplified or altered by heuristics. Even for the most complex systems, the quadrupedal robots with 36 states and 12 control inputs, we can run our solvers in nonlinear MPC-fashion at rates higher than 150 Hz. These frequencies can be achieved even for long time horizons over 500 ms and for complicated locomotion tasks without pre-specified contact sequences, locations or timings.

Viii Performance Optimization

The Control Toolbox is optimized for performance and, if used correctly, constitutes one of the fastest implementations for many state-of-the-art control approaches. This section gives an outline of important steps to achieve the best performance. To achieve best runtime performance, the CT can make of two main techniques:

Viii-1 Multithreading

Thorough multithreading can increase the performance of many optimal control algorithms. While some parts of the optimal control algorithms in the CT are strictly sequential (for instance the backward propagation of the Riccati equations), other parts can be entirely parallelized (e.g. the forward integration on separate multiple shooting intervals in DMS and GNMS and computing linear-quadratic approximations about solution candidates). When employing multi-threading, the required computation time decreases approximately linearly with the number of available cores. In practice, a trade-off needs to be achieved between single-core computation power (for the sequential algorithmic parts) and the overall number of cores (for the simultaneous parts). For experimental results on performance gains through multi-threading, the interested reader is referred to [20].

Viii-2 Vectorization

To achieve the best runtime in every core, one can employ the processor’s vectorization capabilities, which are Single Instruction Multiple Data (SIMD) implementations. SIMD is well-known to be particularly profitable an efficient execution of linear algebra operations, such as matrix-vector multiplications. To date, the authors recommend employing AVX instructions [44], as the register sizes of AVX are continuously growing in modern CPUs.

Ix Further Information

The Control Toolbox is released under the Apache Licence, version 2.0. More detailed documentation and a tutorial are available online, The source-code is available at

X Acknowledgements

Developing and maintaining a large software framework is a team effort. We gratefully acknowledge the contributions of our current and former colleagues Farbod Farshidian, Diego Pardo, Timothy Sandy, Jan Carius, Ruben Grandia and Hamza Merzic.

This research has been funded through a Swiss National Science Foundation Professorship award to Jonas Buchli, the NCCR Digital Fabrication and the NCCR Robotics.


  • [1] R. Tedrake, “Drake: A planning, control, and analysis toolbox for nonlinear dynamical systems.”, 2014. retrieved: 01-01-2018.
  • [2] E. Todorov, T. Erez, and Y. Tassa, “MuJoCo: A physics engine for model-based control,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5026–5033, 2012.
  • [3] B. Houska, H. J. Ferreau, and M. Diehl, “ACADO toolkit-An open-source framework for automatic control and dynamic optimization,” Optimal Control Applications and Methods, vol. 32, no. 3, pp. 298–312, 2011.
  • [4] ACADOS, “ACADOS: Fast and Embedded Optimal Control Problem Solver.” retrieved: 2018-03-25.
  • [5] V. M. Becerra, “Solving complex optimal control problems at no cost with psopt,” in Computer-Aided Control System Design (CACSD), 2010 IEEE International Symposium on, pp. 1391–1396, IEEE, 2010.
  • [6] M. Diehl, D. Leineweber, and A. Schäfer, “MUSCOD-II users’ manual,” tech. rep., Interdisciplinary Center for Scientific Computing (IWR), University of Heidelberg, Germany, 2001.
  • [7] M. A. Patterson and A. V. Rao, “Gpops-ii: A matlab software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming,” ACM Trans. Math. Softw., vol. 41, pp. 1:1–1:37, Oct. 2014.
  • [8] Embotech, “FORCES Pro.” retrieved: 2018-01-01.
  • [9] G. Frison, Algorithms and Methods for Fast Model Predictive Control. PhD thesis, Technical University of Denmark, 2015.
  • [10] H. Ferreau, C. Kirches, A. Potschka, H. Bock, and M. Diehl, “qpOASES: A parametric active-set algorithm for quadratic programming,” Mathematical Programming Computation, vol. 6, no. 4, pp. 327–363, 2014.
  • [11] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical programming, vol. 106, no. 1, pp. 25–57, 2006.
  • [12] P. E. Gill, W. Murray, and M. A. Saunders, “SNOPT: An SQP algorithm for large-scale constrained optimization,” SIAM review, vol. 47, no. 1, pp. 99–131, 2005.
  • [13] S. G. Johnson, “The NLopt nonlinear-optimization package.” retrieved: 2018-01-01.
  • [14] E. Todorov and W. Li, “A generalized iterative LQG method for locally-optimal feedback control of constrained nonlinear stochastic systems,” in American Control Conference, 2005. Proceedings of the 2005, pp. 300–306, IEEE, 2005.
  • [15] M. Giftthaler, M. Neunert, M. Stäuble, J. Buchli, and M. Diehl, “A Family of Iterative Gauss-Newton Shooting Methods for Nonlinear Optimal Control.” arXiv:1711.11006 [cs.SY], 2017.
  • [16] B. M. Bell, “CppAD: a package for C++ algorithmic differentiation,” Computational Infrastructure for Operations Research, vol. 57, 2012.
  • [17] MathWorks Inc., “Matlab & Simulink software.”, 2017. [Online, retrieved: 2017-12-23].
  • [18] Wolfram Research, “Mathematica Software.”, 2018. [Online, retrieved: 2018-01-03].
  • [19] Maplesim, “Maplesoft.”, 2018. [Online, retrieved: 2018-03-25].
  • [20] M. Giftthaler, M. Neunert, M. Stäuble, M. Frigerio, C. Semini, and J. Buchli, “Automatic differentiation of rigid body dynamics for optimal control and estimation,” Advanced Robotics, November 2017.
  • [21] M. L. Felis, “RBDL: an efficient rigid-body dynamics library using recursive algorithms,” Autonomous Robots, vol. 14, no. 2, p. 495–511, 2016.
  • [22] K. Liu, “Dynamic Animation and Robotics Toolkit (DART).” [Online, retrieved 02-02-2018].
  • [23] J. Carpentier, F. Valenza, N. Mansard, et al., “Pinocchio: fast forward and inverse dynamics for poly-articulated systems.”, 2015–2018.
  • [24] F. Nori, S. Traversaro, J. Eljaik, F. Romano, A. Del Prete, and D. Pucci, “icub whole-body control through force regulation on rigid noncoplanar contacts,” Frontiers in Robotics and AI, vol. 2, no. 6, 2015.
  • [25] M. Frigerio, J. Buchli, D. G. Caldwell, and C. Semini, “RobCoGen: a code generator for efficient kinematics and dynamics of articulated robots, based on Domain Specific Languages,” Journal of Software Engineering for Robotics (JOSER), vol. 7, no. 1, pp. 36–54, 2016.
  • [26] D. Pardo, M. Neunert, A. Winkler, R. Grandia, and J. Buchli, “Hybrid direct collocation and control in the constraint-consistent subspace for dynamic legged robot locomotion,” in Proceedings of Robotics: Science and Systems, (Cambridge, Massachusetts), July 2017.
  • [27] Eigen Linear Algebra Library, 2017., retrieved: 2017-08-25.
  • [28] C. D. Bellicoso, M. Blösch, R. Diethelm, P. Fankhauser, P. Furgale, C. Gehring, and H. Sommer, “Kindr - Kinematics and Dynamics for Robotics.” accessed 01-01-2018.
  • [29] K. Ahnert and M. Mulansky, “Odeint – Solving Ordinary Differential Equations in C++,” in AIP Conference Proceedings, pp. 1586–1589, 2011.
  • [30] M. A. Sherman, A. Seth, and S. L. Delp, “Simbody: multibody dynamics for biomedical research,” Procedia IUTAM, vol. 2, pp. 241 – 261, 2011. IUTAM Symposium on Human Body Dynamics.
  • [31] D. Mayne, “A second-order gradient method for determining optimal trajectories of non-linear discrete-time systems,” International Journal of Control, vol. 3, no. 1, pp. 85–95, 1966.
  • [32] G. Frison, D. Kouzoupis, T. Sartor, A. Zanelli, and M. Diehl, “BLASFEO: basic linear algebra subroutines for embedded optimization,” 2017. arXiv:1704.02457 [cs.MS].
  • [33] H. Bock and K. Plitt, “A multiple shooting algorithm for direct solution of optimal control problems,” IFAC World Congress Proceedings Volumes, vol. 17, no. 2, pp. 1603 – 1608, 1984.
  • [34] J. Nocedal and S. J. Wright, Numerical Optimization. Springer, 1999.
  • [35] M. Diehl, R. Findeisen, F. Allgower, H. G. Bock, and J. P. Schloder, “Nominal stability of real-time iteration scheme for nonlinear model predictive control,” IEE Proceedings - Control Theory and Applications, vol. 152, pp. 296–308, May 2005.
  • [36] M. Neunert, M. Stäuble, M. Giftthaler, C. D. Bellicoso, J. Carius, C. Gehring, M. Hutter, and J. Buchli, “Whole-Body Nonlinear Model Predictive Control Through Contacts for Quadrupeds,” 2018. IEEE Robotics and Automation Letters.
  • [37] R. Diankov, “ Ikfast: The robot kinematics compiler.” retrieved: 2018-01-01.
  • [38] B. U. Rehman, Design and Control of a Compact Hydraulic Manipulator for Quadruped Robots. PhD thesis, Istituto Italiano di Tecnologia (IIT) and University of Genova, 2016.
  • [39] C. Semini, N. G. Tsagarakis, E. Guglielmino, M. Focchi, F. Cannella, and D. G. Caldwell, “Design of HyQ–a hydraulically and electrically actuated quadruped robot,” Institution of Mechanical Engineers, Journal of Systems and Control Engineering, vol. 225, pp. 831–849, 2011.
  • [40] M. Neunert, C. de Crousaz, F. Furrer, M. Kamel, F. Farshidian, R. Siegwart, and J. Buchli, “Fast nonlinear model predictive control for unified trajectory optimization and tracking,” in IEEE International Conference on Robotics and Automation, 2016.
  • [41] M. Neunert, F. Farshidian, A. W. Winkler, and J. Buchli, “Trajectory optimization through contacts and automatic gait discovery for quadrupeds,” IEEE Robotics and Automation Letters (RA-L), 2017.
  • [42] M. Giftthaler, F. Farshidian, T. Sandy, L. Stadelmann, and J. Buchli, “Efficient Kinematic Planning for Mobile Manipulators with Non-holonomic Constraints Using Optimal Control,” in IEEE International Conference on Robotics and Automation, pp. 3411–3417, May 2017.
  • [43] T. Sandy, M. Giftthaler, K. Dörfler, M. Kohler, and J. Buchli, “Autonomous repositioning and localization of an In Situ Fabricator,” in IEEE International Conference on Robotics and Automation, pp. 2852–2858, May 2016.
  • [44] N. Firasta, M. Buxton, P. Jinbo, K. Nasri, and S. Kuo, “Intel avx: New frontiers in performance improvements and energy efficiency,” Intel white paper, vol. 19, p. 20, 2008.