Embedded nonlinear model predictive control for obstacle avoidance using PANOC

TL;DR

This paper presents a real-time obstacle avoidance method using PANOC, a fast optimization algorithm, with a new modeling framework for various obstacle shapes, demonstrated on an autonomous vehicle.

Contribution

Introducing a novel obstacle modeling framework compatible with PANOC for real-time control in complex environments.

Findings

PANOC achieves fast convergence suitable for embedded systems.

The framework handles nonconvex, nonlinear obstacles like polytopes and ellipsoids.

Successful implementation on a lab-scale autonomous vehicle.

Abstract

We employ the proximal averaged Newton-type method for optimal control (PANOC) to solve obstacle avoidance problems in real time. We introduce a novel modeling framework for obstacle avoidance which allows us to easily account for generic, possibly nonconvex, obstacles involving polytopes, ellipsoids, semialgebraic sets and generic sets described by a set of nonlinear inequalities. PANOC is particularly well-suited for embedded applications as it involves simple steps, its implementation comes with a low memory footprint and its fast convergence meets the tight runtime requirements of fast dynamical systems one encounters in modern mechatronics and robotics. The proposed obstacle avoidance scheme is tested on a lab-scale autonomous vehicle.

Figures

Equations

$$x_{k+1}=f_k(x_k,u_k).$$

$$O_{kj}=\{z\in {\Re }^{n_d}:h_{kj}^i(z)>0,i\in {\mathbb{N}}_{[1,m_{kj}]}\}.$$

$$z_k\notin O_{kj},\text{ for }j\in {\mathbb{N}}_{[1,q_k]}.$$

$$g_k(x_k,u_k)\le 0,$$

$$u_k\in U_k,$$

$$\mathrm{minimize}$$

$$subject\; to$$

$$x_{k+1}=f_k(x_k,u_k),k\in {\mathbb{N}}_{[0,N-1]},$$

$$u_k\in U_k,k\in {\mathbb{N}}_{[0,N-1]},$$

$$z_k\notin O_{kj},j\in {\mathbb{N}}_{[1,q_k]},k\in {\mathbb{N}}_{[0,N]}$$

$$g_k(x_k,u_k)\le 0,k\in {\mathbb{N}}_{[0,N]},$$

$$g_N(x_N)\le 0.$$

$$O=\{z\in {\Re }^{n_d}:h^i(z)>0,i\in {\mathbb{N}}_{[1,m]}\},$$

$$h^{i_0}(z)\le 0,\text{for some}{i}_0\in {\mathbb{N}}_{[1,m]},$$

$${\psi }_O(z)\text{dfn}\frac{1}{2}\prod _{i=1}^m\text{plus}{h^i(z)}^2=0.$$

$$\nabla {\psi }_O(z)=\{\begin{array}{ll}\sum _{i=1}^m{h}^i(z)\prod _{j\ne i}(h^j(z){)}^2\nabla h^i(z),& \text{if}x\in O,\\ 0,& \text{otherwise.}\end{array}$$

$$\min \{t_1,\dots ,t_m\}=0,$$

$$h^i(z)\le t_i,\text{ for }i\in {\mathbb{N}}_{[1,m]}.$$

$${\tilde{h}}_k(z_k)=\sum _{j=1}^{q_k}{\eta }_{kj}\prod _{i=1}^{m_{kj}}\text{plus}{h_{kj}^i(z)}^2,$$

$${\tilde{\ell }}_k(x,u)$$

$${\tilde{\ell }}_N(x)$$

$$\mathrm{minimize}$$

$$subject\; to$$

$$x_{k+1}=f_k(x_k,u_k),k\in {\mathbb{N}}_{[0,N-1]},$$

$$u_k\in U_k,k\in {\mathbb{N}}_{[0,N-1]},$$

$$\underset{u\in U}{\mathrm{minimize}}\ell (u),$$

$$F_0(u)$$

$$F_{k+1}(u)$$

$$\ell (u)\text{dfn}{\tilde{\ell }}_N(F_N(u))+\sum _{k=0}^{N-1}{\tilde{\ell }}_k(F_k(u),u_k).$$

$$u^{\nu +1}=T_{\gamma }(u^{\nu })\text{dfn}{\Pi }_U(u^{\nu }-\gamma \nabla \ell (u^{\nu })),$$

$$u^{\star }=T_{\gamma }(u^{\star }).$$

$$R_{\gamma }(u)=\frac{1}{\gamma }(u-T_{\gamma }(u)).$$

$$u^{\nu +1}=u^{\nu }-H_{\nu }{R}_{\gamma }(u^{\nu }),$$

$$Q_{\gamma }^{\ell }(v;u)=\ell (u)+\nabla \ell (u{)}^{\top }(v-u)+\text{nicefrac}12\gamma \Vert u-v{\Vert }^2.$$

$${\phi }_{\gamma }(u)=\underset{v\in U}{\inf }{Q}_{\gamma }^{\ell }(v;u).$$

$${\phi }_{\gamma }(u)=\ell (u)-\frac{\gamma }{2}\Vert \nabla \ell (u){\Vert }^2+\frac{1}{2\gamma }{\mathrm{dist}}_U^2(u-\gamma \nabla \ell (u)).$$

Full text

I Introduction

I-A Background and Contributions

Autonomous navigation in obstructed environments is a key element in emerging applications such as driverless cars, fleets of automated vehicles in warehouses and aerial robots performing search-and-rescue expeditions. Well-known approaches for finding collision-free motion trajectories are graph-search methods [1], virtual potential field methods [2] or methods using the concept of velocity obstacles [3].

Recent motion planning research focuses on optimization-based strategies. Here, an optimal motion trajectory is sought while collision-avoidance requirements are imposed as constraints thereon. There is a broad range of such formulations. The most straightforward approach constrains the Euclidean distance between the vehicle and obstacle centers [4]. This, however, only allows to separate spheres and ellipsoids as the computation of distances from arbitrary sets is generally not an easy operation.

Other approaches demand the existence of a hyperplane between the vehicle and obstacle at each time instant [5, 6]. This allows the separation of general convex sets. Some methods formulate the collision avoidance requirement by mixed-integer constraints [7]. These types of problems are, however, cumbersome to solve in real-time. For some specific problem types it is possible to transform the system dynamics such that collision-avoidance translates into simple box constraints on the states [8, 9].

Model predictive control (MPC) is a powerful control strategy where control actions are computed by optimizing a cost function which is chosen in order to achieve a control task. Constraints on states, inputs and outputs can be seamlessly incorporated into such a framework. When the system dynamics is linear, the constraints are affine and the cost functions are quadratic, the associated optimization problem is a quadratic program. There exists a mature machinery of convex optimization algorithms [10, 11] which are fast, robust and possess global convergence guarantees which can be used to solve these problems.

Nevertheless, the dynamics of most systems of interest are better modeled by nonlinear equations and constraints are often nonconvex. This situation is very common in obstacle avoidance involving nonconvex optimization problems. These are commonly solved using sequential quadratic programming (SQP) [10] and interior point (IP) methods [12] which are not well-suited for embedded applications with tight runtime requirements. Nonlinear MPC is often performed using the real-time iteration scheme proposed in [13] which trades speed for accuracy and is accompanied by global convergence guarantees under certain assumptions.

For an algorithm to be suitable for an embedded implementation, it needs to involve only simple steps. This deems methods of the forward-backward-splitting (FBS) type [14], such as the proximal gradient method [15, Sec. 2.3], appealing candidates. FBS-type algorithms can be used to solve nonlinear optimal control problems with simple input constraints via first eliminating the state sequence and expressing the cost as a function of the sequence of inputs alone — the so-called single shooting formulation. However, despite its simplicity, FBS, like all first-order methods, can exhibit slow convergence. Its convergence rate is at best Q-linear with a Q-factor close to one for ill-conditioned problems such as most nonlinear MPC problems.

In this work we propose a new modeling framework for generic constraints which can accommodate general nonconvex sets. The proposed methodology assumes that the obstacles are described by a set of nonlinear inequalities and does not require the computation of projections or distances to them. Then, the obstacle avoidance constraints are written as a nonlinear equality constraint involving a smooth function which, in turn, is relaxed using a penalty function.

The resulting problems are solved using PANOC, a proximal averaged Newton-type method for optimal control, which was recently proposed in [16]. Gradients of the cost function can be efficiently computed using automatic differentiation toolboxes such as CasADi. The algorithm is simple to implement, yet robust, since it combines projected gradient iterations with quasi-Newtonian directions to achieve fast convergence.

The proposed framework is tested on a number of simulation scenarios where we show that it is possible to avoid obstacles of complex shape described by nonlinear inequalities. PANOC is compared with SQP and IP methods and is found to be significantly faster. Furthermore, we present experimental results on a lab-scale robotic platform which runs a C implementation of PANOC.

I-B Notation

Let $\N$ be the set of nonnegative integers, $\N_{[k_{1},k_{2}]}$ be the set of integers in the interval $[k_{1},k_{2}]$ and $\barre=\Re\cup\{+\infty\}$ be the set of extended real numbers. For a matrix $A\in\Re^{m\times n}$, we denote its transpose by $A^{\scriptstyle\top}$. For $x\in\Re$, we define the operator $\plus{x}=\max\{x,0\}$. For a nonempty closed convex set $U\subseteq\Re^{n}$, the projection onto $U$ is the operator $\operatorname{\Pi}_{U}(v)=\operatorname*{\arg\!\min}_{u\in U}\|u-v\|$. The distance from the set $U$ is defined as $\operatorname{dist}_{U}(v)=\inf_{u\in U}\|u-v\|$. The class of continuously differentiable functions $f{}:{}\Re^{n}\to\Re$ is denoted by $C^{1}$. The subset of $C^{1}$ of functions with Lipschitz-continuous gradient is denoted as $C^{1{},{}1}$. We use the notation $C^{1{},{}1}_{{L}}$ for $C^{1{},{}1}$ functions with $L$-Lipschitz gradients.

II NMPC for obstacle avoidance

II-A Problem statement

Kinematic equations lead to continuous-time nonlinear dynamical systems of the form $\dot{x}=f_{c}(x,u,t)$ where $x\in\Re^{n_{x}}$ is the system state, typically a vector comprising of position, velocity and orientation data, and $u\in\Re^{n_{u}}$ is the control signal. We assume that the position coordinates $z\in\Re^{n_{d}}$ are part of the state vector. The continuous-time dynamics can be discretized (for instance, using an explicit Runge-Kutta method) leading to a discrete-time dynamical system of the form

[TABLE]

As it is typically the case in practice, we assume that $f_{k}:\Re^{n_{x}}\times\Re^{n_{u}}\to\Re^{n_{x}}$ are smooth mappings.

The objective of the navigation controller is to steer the controlled vehicle from an initial state $x_{0}$ to a target state $x_{\mathrm{ref}}$, typically a position in space together with a desired orientation. At the same time, the vehicle has to avoid certain, possibly moving, obstacles which are described by open sets $O_{kj}\subseteq\Re^{n_{d}}$, $j\in\N_{[1,q_{k}]}$, each described by

[TABLE]

Sets $O_{kj}$ need not be convex. Obstacle avoidance constraints can be concisely written as

[TABLE]

Moreover, the vehicle is only allowed to move in a domain which is described by the inequalities

[TABLE]

where $g_{k}:\Re^{n_{x}}\times\Re^{n_{u}}\to\Re^{n_{c}}$ is a $C^{2}$ mapping and $\leq$ is meant in the element-wise sense. Control actions $u_{k}$ are constrained in a closed compact set $U_{k}$, that is

[TABLE]

on which it is easy to project and hereafter shall be assumed to be convex. Sets $U_{k}$ often represent box constraints of the form $U_{k}=\{u\in\Re^{n_{u}}{}:{}u_{\mathrm{min}}\leq u\leq u_{\mathrm{max}}\}$.

II-B Nonlinear model predictive control

Nonlinear model predictive control problems arising in obstacle avoidance can be written in the following form

[TABLE]

The stage costs $\ell_{k}:\Re^{n_{x}}\times\Re^{n_{u}}\to\Re$ for $k\in\N_{[0,N-1]}$ in (6a) are $C^{1{},{}1}$ functions penalizing deviations of the state from the reference (destination and orientation) and may be taken to be quadratic functions of the form $\ell(x_{k},u_{k})=(x_{k}-x_{\mathrm{ref}})^{{\scriptscriptstyle\top}}Q_{k}(x_{k}-x_{\mathrm{ref}})+(u_{k}-u_{\mathrm{ref}})^{{\scriptscriptstyle\top}}R_{k}(u_{k}-u_{\mathrm{ref}})$. The terminal cost $\ell_{N}:\Re^{n_{x}}\to\Re$ in (6a) is a $C^{1{},{}1}$ function such as $\ell_{N}(x_{N})=(x_{N}-x_{\mathrm{ref}})^{{\scriptscriptstyle\top}}Q_{N}(x_{N}-x_{\mathrm{ref}})$.

II-C Reformulation of obstacle avoidance constraints

Consider an obstacle described as the intersection of a finite number of strict nonlinear inequalities

[TABLE]

where $h^{i}:\Re^{n_{d}}\to\Re$ are $C^{1{},{}1}$ functions. The constraint $z\notin O$ — cf. (6e) — is satisfied if and only if

[TABLE]

or, equivalently, $\plus{h^{i_{0}}(z)}^{2}=0$. This constraint can then be encoded as

[TABLE]

We have expressed the obstacle avoidance constraints as a nonlinear equality constraint. We should remark that, unlike approaches based on the distance-to-set function [17], function $\psi_{O}$ in (9) is a $C^{1}$ function of $z$. Indeed, $\psi_{O}$ is differentiable in $\Re^{n_{d}}\setminus O$ (and equal to [math]), it is differentiable in $O$ with gradient

[TABLE]

Note that $\nabla\psi_{O}$ on the boundary of $O$ vanishes, so it is everywhere continuous. If, additionally, $O$ is bounded, $\psi_{O}$ is $C^{1,1}$.

The formulation of Eq. (9) can be used for obstacles described by quadratic constraints of the form $O=\{z\in\Re^{n_{d}}{}:{}1-(z-c)^{\scriptscriptstyle\top}E(z-c)>0\},$ such as balls and ellipsoids. Then, the associated equality constraint becomes $\plus{1-(z-c)^{\scriptscriptstyle\top}E(z-c)}^{2}=0.$ Polyhedral obstacles of the form $O=\{z\in\Re^{n_{d}}{}:{}b_{i}-a_{i}^{{\scriptscriptstyle\top}}z>0,\in\N_{[1,m]}\},$ with $b_{i}\in\Re$ and $a_{i}\in\Re^{n_{d}}$, can also be accommodated by the constraint $\prod_{i=1}^{m}\plus{b_{i}-a_{i}^{\scriptscriptstyle\top}z}^{2}=0.$

Equation (9) can also be used to describe general nonconvex constraints such as ones described by semi-algebraic sets where $h^{i}(z)$ are polynomials as well as any other obstacle which is available in the aforementioned representation.

Obstacle avoidance constraints (9) are equivalent to the existence of $t\in\Re^{m}$ so that

[TABLE]

This observation reveals a link to vertical complementarity constraints which have been studied extensively in the literature [18, 19].

II-D Relaxation of constraints

Due to the fact that modeling errors and disturbances may lead to the violation of imposed constraints and infeasibility of the MPC optimization problem, it is common practice in MPC to replace state constraints by appropriate penalty functions known as soft constraints. Quadratic penalty functions are often used for this purposes revealing a clear link between this approach and the quadratic penalty method in numerical optimization [15, Sec. 4.2.1].

Equality constraints, such as the ones arising in the reformulation of the obstacle avoidance constraints in Section II-C, can be relaxed by means of soft constraints. Indeed, constraints of the form $\Phi_{k}(z_{k})=0$, $k\in\N_{[1,N]}$, where $\Phi_{k}:\Re^{n_{d}}\to\Re_{+}$ are $C^{2}$ functions, can be relaxed by introducing the penalty functions $\tilde{\Phi}_{k}(z_{k})=\eta_{k}\Phi_{k}(z_{k})$ for some weight factors $\eta_{k}>0$.

That said, constraints like (9) for a set of time-varying obstacles $O_{kj}=\{z\in\Re^{n_{d}}{}:{}h_{kj}^{i}(z)>0,i\in\N_{[1,m_{kj}]}\}$, with $j\in\N_{[1,q_{k}]}$, can be relaxed by the appending the following term in the original cost function

[TABLE]

for some positive weight factors $\eta_{kj}>0$.

Note that the proposed approach for dealing with obstacle avoidance constraints requires only a representation of the obstacles in the generic form (7) and does not call for the computation of distances to the obstacles as in distance-based methods [4, 17], nor does it require the obstacles to be convex sets.

Similarly, inequality constraints of the form $g_{k}(x_{k},u_{k})\leq 0$ and $g_{N}(x_{N})\leq 0$ can be relaxed by introducing the penalty functions $\tilde{g}_{k}(x,u)=\beta_{k}\plus{g_{k}(x,u)}^{2}$ and $\tilde{g}_{N}(x)=\beta_{N}\plus{g_{N}(x)}^{2}$ for positive weights $\beta_{k}>0$, $k\in\N_{[0,N]}$.

We may now relax the state constraints in (6) by defining the modified stage cost and terminal cost functions

[TABLE]

leading to the following relaxed optimization problem without state constraints

[TABLE]

II-E NMPC problem formulation

In this section we cast the nonlinear MPC problem (6) as

[TABLE]

where the optimization is carried out over vectors $u=(u_{0},\ldots,u_{N-1})\in\Re^{n}$, with $n=Nu_{u}$ and $U\dfn U_{0}\times U_{1}\times\cdots\times U_{N-1}$ and $\ell:\Re^{n}\to\Re$ is a real-valued $C^{1{},{}1}_{{L_{\ell}}}$ function.

We introduce the following sequence of functions $F_{k}:\Re^{n}\to\Re^{n_{x}}$ for $k\in\N_{[0,N-1]}$

[TABLE]

Define the smooth function $\ell:\Re^{n}\to\Re$

[TABLE]

The gradient of function $\ell$ in (16) can be computed by means of the reverse mode of automatic differentiation (also known as adjoint method or backpropagation) as shown in Alg. 1 [20].

Problem (14) is in a form that allows the application of the projected gradient iteration

[TABLE]

with $\gamma>0$. In particular, if $\ell\in C^{1{},{}1}_{{L_{\ell}}}$ and $\gamma<\nicefrac{{2}}{{L_{\ell}}}$, then all accumulation points of (17), $u^{\star}$, are fixed points of $T_{\gamma}$ called $\gamma$-critical points, that is [10, Prop. 2.3.2]

[TABLE]

III Fast nonlinear MPC

The problem of finding a fixed point for $T_{\gamma}$ can be reduced to the equivalent problem of finding a zero of the fixed-point residual operator which is defined as the operator

[TABLE]

This motivates the adoption of a Newton-type iterative scheme

[TABLE]

where $H_{\nu}$ are invertible linear operators, appropriately chosen so as to encode first-order information about $R_{\gamma}$. This is done by enforcing the inverse secant condition $s^{\nu}=H_{\nu+1}y^{\nu}$, for $s^{\nu}=u^{\nu+1}-u^{\nu}$ and $y^{\nu}=R_{\gamma}(u^{\nu+1})-R_{\gamma}(u^{\nu})$ and can be obtained by quasi-Newtonian methods such as the limited-memory BFGS (L-BFGS) method [10] which is free from matrix operations, requires only a limited number of inner products and is suitable for embedded implementations. The main weakness of this approach is that convergence is only guaranteed in a neighborhood of a $\gamma$-critical point $u^{\star}$. We shall describe a globalization procedure which hinges on the notion of the forward-backward envelope function.

III-A The Forward-Backward envelope function

The forward-backward envelope (FBE) is an exact, continuous, real-valued merit function for (14) [21, 22, 16, 23, 24]. Function $\ell$ can be approximated at a point $u\in U$ by the quadratic upper bound

[TABLE]

The FBE is then defined as

[TABLE]

This construction is illustrated in Fig. 1. Provided that it is easy to compute the distance to $U$, the FBE can be easily computed by

[TABLE]

Therefore, the computation of the FBE is of the same complexity as that of a forward-backward step.

The FBE possesses several favorable properties, perhaps the most important being that it is real-valued, continuous and for $\gamma\in(0,\nicefrac{{1}}{{L_{\ell}}})$ shares the same (local/strong) minima with the original problem (14). This means that (14) is reduced to an unconstrained minimization problem. Moreover, if $\ell\in C^{2}$, then $\varphi_{\gamma}\in C^{1}$ with $\nabla\varphi_{\gamma}(u)=(I-\gamma\nabla^{2}\ell(u))R_{\gamma}(u)$.

III-B PANOC Algorithm

We employ the proximal averaged Newton-type method for optimal control, for short PANOC, which was recently proposed in [16]. PANOC performs fast Newton-type updates while an FBE-based linesearch endows it with global convergence properties while it uses the same oracle as the projected gradient method.

The iterative scheme, which is presented in Alg. 2, involves the computation of a projected gradient point $\bar{u}^{\nu}$ in step 2 and an L-BFGS direction $d^{\nu}$ in step 4. L-BFGS obviates the need to store or explicitly update matrices $H_{\nu}$ in (20) by storing a number $\mu$ of past values of $s^{\nu}$ and $y^{\nu}$. The computation of $d^{\nu}$ requires only inner products which amount to a maximum of $4\mu n$ scalar multiplications. In step 5, the iterates are updated using a convex combination of the projected gradient update direction $-\gamma R_{\gamma}(u^{\nu})$ and a fast quasi-Newtonian direction $d^{\nu}$. The algorithm is terminated when $\|R_{\gamma}(u^{\nu})\|_{\infty}$ drops below a specified tolerance.

In line 5, the backtracking line search procedure ensures that a sufficient decrease condition is satisfied using the FBE as a merit function. Under mild assumptions, eventually only fast updates are activated and updates reduce to $u^{\nu{}+{}1}=u^{\nu}{}+{}d^{\nu}$. The sequence of fixed-point residuals $\{r^{\nu}\}_{\nu\in\N}$ converges to [math] square summably, while PANOC produces sequences of iterates, $\{u^{\nu}\}_{\nu\in\N}$ and $\{\bar{u}^{\nu}\}_{\nu\in\N}$, whose cluster points are $\gamma$-critical points.

In absence of a Lipschitz constant $L_{\ell}$, the algorithm can be initialized with an estimate, e.g., $L_{\ell}^{0}={\|\nabla_{u}\ell(u^{0}+\delta u)-\nabla_{u}\ell(u^{0})\|}/{\|\delta u\|},$ where $\delta u\in\Re^{n}$ is a small perturbation vector. Then, step 2 in Alg. 2 needs to be replaced by the backtracking procedure of Alg. 3 which updates $L_{\ell}$, $\sigma$ and $\gamma$.

Alg. 3 updates $L_{\ell}$, $\sigma$ and $\gamma$ only a finitely many times, so, it does not affect the convergence properties of the algorithm.

In MPC, we may warm start the algorithm using the optimal solution of the previous MPC instance. Unlike SQP and IP methods, PANOC does not require the solution of linear systems or quadratic programming problems at every iteration and it involves only very simple operations such as vector additions, scalar and inner products. Additionally, PANOC converges globally, that is, from any initial point $u^{0}\in\Re^{n}$.

IV Simulations

The proposed methodology is validated on the control of a mobile robot carrying a trailer shown in Fig. 2. Assuming zero slip of the trailer wheels, the nonlinear kinematics is

[TABLE]

where the state vector $x=(p_{x},p_{y},\theta)$ comprises the coordinates $p_{x}$ and $p_{y}$ of the trailer and the heading angle $\theta$. The input $u=(u_{x},u_{y})$ is a velocity reference which is tracked by a low-level controller. The distance between the center of mass of the trailer and the fulcrum connecting to the towing vehicle is $L=0.5\,\mathrm{m}$. In Fig. 3 we present four obstacle avoidance scenarios involving, among other, obstacles described by polynomial and trigonometric functions. The system dynamics given in (25) is discretized using the fourth-order Runge-Kutta method and we use memory $\mu=10$ for the L-BFGS directions in Alg. 2.

The single shooting formulation of Section II is solved with PANOC, the interior point solver IPOPT, the forward-backward splitting (FBS) implementation of ForBES and the SQP of MATLAB’s fmincon. The problem was also brought in a multiple shooting formulation where obstacle avoidance constraints were imposed as in (9) and it was solved using IPOPT and SQP. A comparison of runtime is shown in Fig. 4. All simulations were executed on an Intel i5-6200U CPU with $12\,\mathrm{GB}$ RAM machine running Ubuntu 14.04. PANOC exhibits very low runtime and outperforms all other solvers by approximately two orders of magnitude.

V Experimental results and discussion

The proposed NMPC algorithm is validated on an in-house mobile robot with a trailer (Fig. 2). The robot has four independent DC motors with Mecanum wheels, which render it holonomic. A low-level microcontroller implements a velocity controller for each motor. Therefore, from a high-level perspective, the robot can be treated as a velocity-steered holonomic device. A trailer is attached to the robot and the angle between robot and trailer is measured with a rotary potentiometer. A ceiling camera detects the robot’s absolute position and orientation using a marker attached on top of the robot. This information is sent to the robot over Wi-Fi and is merged with local encoder measurements to retrieve a fast and accurate estimate of its pose. An Odroid XU4 platform runs the NMPC algorithm on board the robotic platform. PANOC is implemented in C following the C89 standard without external dependencies and can, therefore, be readily executed on embedded systems. The C code for performing Alg. 1 was generated using the AD tool CasADi [25].

The NMPC algorithm is used to steer the vehicle to a desired destination and trailer orientation $x_{\mathrm{ref}}=(3.77,1.40,0.0)$. The system dynamics is discretized using the Euler method. A quadratic cost function is employed with $Q_{k}=Q_{N}=0.1\cdot I_{3}$ and $R_{k}=0.01\cdot I_{2}$. Box constraints are imposed on the inputs with $u_{\mathrm{min}}=-0.8\,\nicefrac{\mathrm{m}}{\mathrm{s}}$ and $u_{\mathrm{max}}=0.8\,\nicefrac{\mathrm{m}}{\mathrm{s}}$. The obstacles are modelled as discussed in Section II-C and are slightly enlarged for safety so as to completely avoid collisions. A horizon length of $N{}={}50$ is used and the NMPC control rate is $10\,\mathrm{Hz}$. Fig. 5 shows the resulting motion of the robot and the predicted position trajectories at four time instants. PANOC is terminated when $\|r^{\nu}\|_{\infty}\leq 10^{-6}$ or if the number of iterations reaches $500$. The number of iterations $\bar{\nu}$, the solving time $t_{\mathrm{s}}$ on the Odroid XU4 and the norm of the fixed point residual at termination at each time instant $k$ are shown in Fig. 6.

Additional material and videos from the experiments are found at https://kul-forbes.github.io/PANOC.

VI Conclusions and research directions

We proposed a novel framework which enables us to encode nonconvex obstacle avoidance constraints as a smooth nonlinear equality constraint. This offers a flexible modeling framework and allows the formulation of nonlinear MPC problems with obstacles of complex nonconvex geometry. The resulting MPC problem is solved with PANOC which has favorable theoretical convergence properties and outperforms state-of-the-art NLP solvers. This framework was tested using a C89 implementation of PANOC on a lab-scale system.

Future work will focus on the development of a proximal Lagrangian framework for the online adaptation of the weight parameters in the obstacle avoidance penalty functions — cf. (12) — so that (predicted) constraint violations, modeled by $\tilde{h}(z_{k})$, are below a desired tolerance. Moreover, the use of semismooth Newton directions in PANOC will lead to quadratic convergence and superior performance [21].