Least squares dynamics in Newton-Krylov Model Predictive Control

TL;DR

This paper introduces a least squares approach within Newton-Krylov methods for nonlinear Model Predictive Control, offering an alternative to Ohtsuka's C/GMRES that handles inconsistent constraints effectively.

Contribution

It proposes a novel least squares formulation for Newton-Krylov methods in MPC, improving robustness when constraints are inconsistent.

Findings

Numerical tests confirm the effectiveness of the proposed method.

The approach handles constraint inconsistencies better than existing methods.

Modified Newton-Krylov methods show promising results in simulations.

Abstract

Newton-Krylov methods for nonlinear Model Predictive Control are pioneered by T. Ohtsuka under the name "C/GMRES". Ohtsuka eliminates a system state over the horizon from Karush-Kuhn-Tucker stationarity conditions of a Lagrangian using equations of system dynamics. We propose instead using least squares to fit the state to the dynamics and some constraints on the state, if they are inconsistent. Correspondingly modified Newton-Krylov methods are described. Numerical tests demonstrate workability of our modification.

Equations

$$\underset{u,p}{\min }J(u,p),$$

$$J(u,p)=\varphi (x(\tau ),p){|}_{\tau =t+T}+{\int }_t^{t+T}L(\tau ,x(\tau ),u(\tau ),p)d\tau$$

$$\frac{dx}{d\tau }=f(\tau ,x(\tau ),u(\tau ),p)+{\eta }_f,\tau \in [t,t+T],$$

$$g(\tau ,x(\tau ),u(\tau ),p)+{\eta }_g=0,\tau \in [t,t+T],$$

$$x(\tau ){|}_{\tau =t}=x(t),$$

$$C(\tau ,x(\tau ),u(\tau ),p){|}_{\tau \in [t,t+T]}=0,$$

$$\psi (x(\tau ),p){|}_{\tau =t+T}=0.$$

$$\underset{x}{\min }S(x),$$

$$S(x)={\Vert f(\tau ,x(\tau ),u(\tau ),p)-dx/d\tau \Vert }_f^2$$

$$+{\Vert g(\tau ,x(\tau ),u(\tau ),p)\Vert }_g^2,$$

$$\underset{u,p,x}{\min }J(u,p)+S(x),$$

$$\underset{x_1,x_2,\dots x_N}{\min }\sum _{i=0}^{N-1}\Vert (x_{i+1}-x_i)/\Delta \tau -f({\tau }_i,x_i,u_i,p){\Vert }_{W_f^{-1}}^2$$

$$+\Vert g({\tau }_{i+1},x_{i+1},u_{i+1},p){\Vert }_{W_g^{-1}}^2,$$

$$B=\frac{1}{\Delta \tau }\left(\begin{array}{cccc}I& & & \\ -I& I& & \\ & \ddots & \ddots & \\ & & -I& I\end{array}\right)$$

$$R=B\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right]-\left[\begin{array}{c}f({\tau }_0,x_0,u_0,p)+x_0/\Delta \tau \\ f({\tau }_1,x_1,u_1,p)\\ ⋮\\ f({\tau }_{N-1},x_{N-1},u_{N-1},p)\end{array}\right],$$

$$G=\left[\begin{array}{c}g({\tau }_1,x_1,u_1,p)\\ g({\tau }_2,x_2,u_2,p)\\ ⋮\\ g({\tau }_N,x_N,u_N,p)\end{array}\right].$$

$$\underset{x}{\min }{R}^T{W}_f^{-1}R+G{W}_g^{-1}G.$$

$$\nabla G=\left[\begin{array}{ccc}{\nabla }_{x}g({\tau }_1,x_1,u_1,p)& & \\ & {\nabla }_{x}g({\tau }_2,x_2,u_2,p)& \\ & \ddots & \\ & & {\nabla }_{x}g({\tau }_N,x_N,u_N,p)\end{array}\right],$$

$$\nabla R=B-\left[\begin{array}{cccc}0& & & \\ {\nabla }_{x}f({\tau }_1,x_1,u_1,p)& 0& & \\ & {\nabla }_{x}f({\tau }_2,x_2,u_2,p)& 0& \\ & \ddots & & \ddots \\ & & {\nabla }_{x}f({\tau }_{N-1},x_{N-1},u_{N-1},p)& 0\end{array}\right].$$

$$(\nabla R{)}^T{W}_f^{-1}R+(\nabla G{)}^T{W}_g^{-1}G=0.$$

$$\underset{u_i,p}{\min }\left[\varphi (x_N,p)+\sum _{i=0}^{N-1}L({\tau }_i,x_i,u_i,p)\Delta \tau \right],$$

$$C({\tau }_i,x_i,u_i,p)=0,i=0,1,\dots ,N-1,$$

$$\psi (x_N,p)=0.$$

$$\mathcal{L}(X,U)=\varphi (x_N,p)+\sum _{i=0}^{N-1}L({\tau }_i,x_i,u_i,p)\Delta \tau$$

$$+{\lambda }_0^T[x_0-x(t)]$$

$$+[{\lambda }_1^T\dots {\lambda }_N^T][(\nabla R{)}^T{W}_f^{-1}R+(\nabla G{)}^T{W}_g^{-1}G]\Delta \tau$$

$$+\sum _{i=0}^{N-1}{\mu }_i^{T}C({\tau }_i,x_i,u_i,p)\Delta \tau +{\nu }^T\psi (x_N,p),$$

$$U(t)=[u_0^T,\dots ,u_{N-1}^T,{\mu }_0^T,\dots ,{\mu }_{N-1}^T,{\nu }^T,p^T{]}^T.$$

$$\frac{\partial \mathcal{L}}{\partial x}(X,U)=0.$$

$$\displaystyle F[U,x,t]=\left[\begin{array}[]{c}\begin{array}[]{c}\frac{\partial\cal L}{\partial u_{0}}(X,U)\\ \vdots\\ \frac{\partial\cal L}{\partial u_{i}}(X,U)\\ \vdots\\ \frac{\partial\cal L}{\partial u_{N-1}}(X,U)\end{array}\\ \\ \begin{array}[]{c}C(\tau_{0},x_{0},u_{0},p)\Delta\tau\\ \vdots\\ C(\tau_{i},x_{i},u_{i},p)\Delta\tau\\ \vdots\\ C(\tau_{N-1},x_{N-1},u_{N-1},p)\Delta\tau\end{array}\\ \\ \psi(x_{N},p)\\[8.61108pt] \frac{\partial\cal L}{\partial p}(X,U)\end{array}\right].$$

$$F[U(t),x(t),t]=0$$

$$F[U_j,x_j,t_j]-F[U_{j-1},x_j,t_j]=b_j,$$

$$b_j=-F[U_{j-1},x_j,t_j].$$

Full text

I Introduction

The paper is concerned with Model Predictive Control (MPC), see, e.g., [1, 2, 3, 4], for cases, where a model of a state of a system dynamics contradicts to some state constrains, making MPC infeasible. The contradictions may appear, e.g., from uncertainties and inaccuracies. We propose using least squares to fit the state to both the dynamics and the contradicting constraints, within a framework of Newton-Krylov methods for nonlinear MPC (NMPC), pioneered by Ohtsuka [5] for numerical solution of the MPC problems.

As an illustrating example (also used in our numerical tests), let us consider a continuous dynamical system where the state automatically satisfies an equality constraint, e.g., the state is on a smooth manifold, e.g., a sphere in [6]. Discretized dynamical models, used for state prediction over a finite MPC horizon, approximate the continuous model and may not exactly satisfy the state equality constraint of the continuous case. When the predictive horizon is long, the predicted trajectories may deviate far from the manifold determined by the equality constraint. In [6], we propose solving optimal control problems over smooth manifolds by using the so-called “structure preserving integration methods” [7] within the Ohtsuka’s method [5, 8, 9].

In the present work, we introduce a new prediction technology, aimed at removing the inconsistency of the state dynamics with some equality constraint on the state, by means of the least squares. In our inner-outer approach, the inner layer is the least squares fit of the state to the dynamics and the inconsistent constraints on the state, while the outer layer is NMPC solved by the Newton-Krylov methods of Ohtsuka.

We formulate a theoretical framework of two-level MPC, develop a numerical method similar to Ohtsuka’s method, and show numerical results for a test minimum-time problem describing motion on a unit sphere with constrained controls.

II Least squares dynamics in continuous MPC

MPC determines a control input $u(t)$ by solving a prediction model on a finite horizon $[t,t+T]$. We consider a modified, using unknown disturbance vectors $\eta_{f}$ and $\eta_{g}$, variant of the prediction model from [6], where the control $u(\tau)$ and a parameter vector $p$ minimize

[TABLE]

the performance index

[TABLE]

subject to uncertain model dynamics

[TABLE]

uncertain constraint on the state

[TABLE]

and the following certain constraints

[TABLE]

[TABLE]

[TABLE]

The initial value $x(\tau)|_{\tau=t}$ for the time-dependent differential equation (2) is the current state vector $x(t)$ of the dynamic system. The control vector $u=u(\tau)$, which solves the prediction problem, is used as an input to control the dynamic system at time $t$. The components of the vector $p(t)$ are parameters of the system.

The generally nonlinear equation (2) exactly describes the model system dynamics, while the generally nonlinear constraint (3) is also exact. But the disturbance vectors $\eta_{f}$ and $\eta_{g}$ are unknown and, if dropped, may result in inconsistency for arbitrary $u$ and $p$, thus leading to an infeasible MPC problem. Assuming the vector-functions $u$ and $p$ fixed, the disturbance vectors $\eta_{f}$ and $\eta_{g}$ can be minimized with respect to the function $x$ over the horizon via least squares, i.e.

[TABLE]

where

[TABLE]

and $\|\cdot\|_{f}$ and $\|\cdot\|_{g}$ are functional norms, e.g., based on the weighted $L_{2}$ norm of a function $h$ as $\|h\|_{W^{-1}}^{2}=\int h^{T}(\tau)W^{-1}(\tau)h(\tau)d\tau$ with the weight matrix $W^{-1}(\tau)$. We note that the solution $x(\tau)$ over the horizon $\tau\in[t,t+T]$ has the given fixed initial, when $\tau=t$, value $x(t)$.

III Relaxed dynamics alternating minimization

Our discussion in §II motivates relaxing (2) and (3) by simply adding the term $S(x)$ minimized in (7) to the performance index $J(u,p)$ to be minimized in (1), i.e.

[TABLE]

subject to only the certain constraints, i.e. (4), (5), and (6).

Explicitly adding $x$ to the set of minimization variables in (8), may add computation costs to perform minimization, compared to the original setup (1). A well known idea of alternating minimization, see, e.g., [10], may reduce computations by iteratively minimizing $J(u,p)+S(x)$ alternatively and separately with respect to $u,p$ and with respect to $x$.

One can interpret such an alternating minimization as inner-outer approach, where the inner layer is the least squares fit (7) of the state to the dynamics and the inconsistent constraints on the state, while the outer layer is NMPC minimization (1), solved iteratively. Newton-Krylov methods of Ohtsuka [5, 9] are examples of interest of iterative minimization (1) of the performance index. In the next section, we describe in the discrete case, how the original setup from [5], where (2) and (3) are treated as exact certain constraints, can be modified to substitute the least squares minimization (7) for (2) and (3), formulating the discrete Karush-Kuhn-Tucker (KKT) necessary conditions of (1) with the relaxed dynamics.

IV Least squares discrete dynamics in KKT

Continuous formulation of the finite horizon prediction problem stated above can be discretized on a uniform time grid over the horizon $[t,t+T]$ partitioned into $N$ equal time steps of size $\Delta\tau$, and the time-continuous vector functions $x(\tau)$ and $u(\tau)$ are sampled at the grid points $\tau_{i}$, $i=0,1,\ldots,N$ and denoted by the indexed values $x_{i}$ and $u_{i}$ respectively. The integral of the performance cost $J$ over the horizon is approximated by means of the rectangular quadrature rule. The time derivative of the state vector is approximated by the forward difference formula.

Before deriving the Euler equations for the NMPC formulation, we discretize $x$ in the least squares minimization (7),

[TABLE]

keeping the first component $x_{0}=x(t)$ fixed, where $\|\cdot\|_{W^{-1}}$ denote weighted, using a matrix $W$, $2$-norms of vectors.

When the disturbances $\eta_{f}$ and $\eta_{g}$ are of random nature, the covariance matrices $W_{f}$ and $W_{g}$ may be available. In our test examples in §VI, we use the covariance matrices of the form $W_{f}=\alpha^{-1}I$ and $W_{g}=\beta^{-1}I$ with $\alpha=1$ and a suitable scalar $\beta>0$, with $I$ being the identity matrix.

For convenience, we introduce the block bidiagonal matrix

[TABLE]

and the vectors

[TABLE]

[TABLE]

In this notation, the discrete version of the least squares minimization (7) takes the following form,

[TABLE]

The gradients with respect to $x$ of the vectors $G$ and $R$ equal

[TABLE]

[TABLE]

Hence the solution $x_{i}$, $i=1,\ldots N$, of the discrete least squares minimization satisfies the equation

[TABLE]

The discretized optimal control problem NMPC is then formulated as follows:

[TABLE]

subject to the system (9) for $x_{i}$ and the equality constraints

[TABLE]

[TABLE]

Necessary optimality conditions for the discretized finite horizon problem can be derived by means of the discrete Lagrangian function

[TABLE]

where we gather the variables into vectors $X=[x_{i}\;\lambda_{i}]^{T}$, $i=0,1,\ldots,N$, and $U=[u_{i}\;\mu_{i}\;\nu\;p]^{T}$, $i=0,1,\ldots,N-1$. Here, $\lambda=[\lambda_{1}^{T}\ldots\lambda^{T}_{N}]^{T}$ is the costate vector, and $\mu$ is the Lagrange multiplier vector associated with the constraint (10). The terminal constraint (11) is relaxed by the aid of the Lagrange multiplier $\nu$.

Calculating the derivatives of the Lagrangian $\mathcal{L}$ we obtain the necessary optimality KKT conditions, $\mathcal{L}_{\lambda_{i}}=0$, $\mathcal{L}_{x_{i}}=0$, $i=0,1,\ldots,N$, $\mathcal{L}_{u_{j}}=0$, $\mathcal{L}_{\mu_{j}}=0$, $i=0,1,\ldots,N-1$, $\mathcal{L}_{\nu_{k}}=0$, $\mathcal{L}_{p_{l}}=0$.

We further convert the KKT conditions into a nonlinear equation $F[U,x,t]=0$, where the vector $U$ combines the control input $u$, the Lagrange multiplier $\mu$, the Lagrange multiplier $\nu$, and the parameter $p$, all in one vector:

[TABLE]

The vector argument $x$ in $F[U,x,t]$ denotes the current measured or estimated state vector, which serves as the initial vector $x_{0}$ in the following procedure, which eliminates the state variables $x_{i}$ and costate variables $\lambda_{i}$.

Having the current state $x_{0}$, measured or estimated, we compute $x_{i}$, $i=1,2\ldots,N$, by solving least squares equations (9) instead of the forward Euler method $x_{i+1}=x_{i}+f(\tau_{i},x_{i},u_{i},p)\Delta\tau$ of [5].

Then compute the costates $\lambda_{i}$, $i=N,N\!-\!1,\ldots,1$, from the system of linear equations

[TABLE]

The value $\lambda_{N}$ is defined by the differentiation of the term $\nu^{T}\psi(x_{N},p)$ with respect $x$. 2. 2.

Calculate $F[U,x,t]$, using just obtained $x_{i}$ and $\lambda_{i}$, as

[TABLE]

The equation with respect to the unknown vector $U(t)$

[TABLE]

gives the required necessary optimality conditions.

Ohtsuka in [5] proposes solving (15) using Newton-Krylov methods applied a forward-difference approximation to the Jacobian $F_{U}$ as described in [11]. In the next section, we repeat the necessary details, following [5, 9, 12], only slightly modified to take into account the relaxed dynamics.

V Newton-Krylov methods to solve KKT

Let us assume that the dynamic system, which is controlled with the MPC approach, is sampled on a uniform time grid $t_{j}=j\Delta t$, $j=0,1,\ldots$ and denote $x_{j}=x(t_{j})$. Equation (15) must be solved at each time step $t_{j}$ online on the controller board, which is the most computationally challenging part of an NMPC implementation for systems with fast dynamics.

The nonlinear equation $F[U_{j},x_{j},t_{j}]=0$ with respect to the unknown variables $U_{j}$ approximating $U(t_{j})$ is equivalent to the following equation

[TABLE]

where

[TABLE]

Using a small scalar $h>0$, which is, in general, different from the time steps $\Delta t$ and $\Delta\tau$, we introduce, as, e.g., in [11], the forward difference operator

[TABLE]

approximating the derivative $F_{U}[U_{j-1},x_{j},t_{j}](V)$ along the direction $V$. We remark that the equation $F[U_{j},x_{j},t_{j}]=0$ is equivalent to the operator equation $a_{j}(\Delta U_{j}/h)=b_{j}/h$, where $\Delta U_{j}=U_{j}-U_{j-1}$.

Let us introduce an $m\times m$ matrix $A_{j}$ with the columns $A_{j}e_{k}$, $k=1,\ldots,m$, defined by the formula $A_{j}e_{k}=a_{j}(e_{k})$, where $m$ is the dimension of the vector $U$ and $e_{k}$ denotes the $k$-th column of the $m\times m$ identity matrix. The matrix $A_{j}$ is an $O(h)$ approximation of the Jacobian matrix $F_{U}[U_{j-1},x_{j},t_{j}]$, which is symmetric by Theorem 1.

Theorem 1

The Jacobian matrix $F_{U}[U,x,t]$ is symmetric.

Proof:

The equation $\mathcal{L}_{X}(X,U)=0$ is solvable with respect to $X$ due to the solvability of the least squares minimization for $x_{i}$ and a system of linear equations for $\lambda_{i}$. The rest of the proof is identical to that in [12] for the case of the exact dynamics and is provided here for completeness.

Let us denote the solution to $\mathcal{L}_{X}(X,U)=0$ by $X=g(U)$. Then $F[U]=\mathcal{L}_{U}(g(U),U)$ and

[TABLE]

Differentiation of the identity $\mathcal{L}_{U}(g(U),U)=0$ with respect to $U$ gives the identity

[TABLE]

Differentiation of the identity $\mathcal{L}_{X}(g(U),U)=0$ with respect to $U$ gives the identity

[TABLE]

Hence $g_{U}=-\mathcal{L}_{XX}^{-1}(g(U),U)\mathcal{L}_{XU}(g(U),U)$ and

[TABLE]

which is called the Schur complement of the symmetric Hessian matrix of $\mathcal{L}$ at the point $(X,U)=(g(U),U)$. The Schur complement of any symmetric matrix is symmetric. ∎

Suppose that an approximate solution $U_{0}$ to the equation $F[U_{0},x_{0},t_{0}]=0$ is available. Finding sufficiently accurate approximation $U_{0}$ is crucial for success of Newton-like methods and search for it is usually a challenging operation. However, we omit descriptions of suitable methods for finding the starting value $U_{0}$ here because it is unrelated to, although needed for, the “warm-start” procedure described below.

The first block entry of $U_{0}$ is taken as the input control $u_{0}$ at the state $x_{0}$. The next state $x_{1}=x(t_{1})$ is measured by sensors or estimated.

At the time $t_{j}$, $j>1$, we have the state $x_{j}$ and the vector $U_{j-1}$ from the previous time $t_{j-1}$. Our goal is to solve the following equation with respect to $V$:

[TABLE]

Then we set $\Delta U_{j}=hV$, $U_{j}=U_{j-1}+\Delta U_{j}$ and choose the first block component of $U_{j}$ as the control $u_{j}$. The next system state $x_{j+1}=x(t_{j+1})$ is measured by sensors or estimated.

A direct way to solve the operator equation (19) is forming the matrix $A_{j}$ explicitly and then solving the system of linear equations $A_{j}\Delta U_{j}=b_{j}$; e.g., by the Gaussian elimination.

A faster alternative is solving (19) by Krylov iterative methods (such as GMRES [5, 11], or MINRES [13], possibly with preconditioning [12]), where the operator $a_{j}(V)$ is used without explicit construction of the matrix $A_{j}$; cf., [5, 11]. Krylov methods, applied to a finite difference approximation (17) of a Jacobian, are call “Newton-Krylov methods” in [11].

VI Proof of concept numerical example

We numerically simulate a minimum-time motion from an initial state $x_{0}$ to a terminal state $x_{f}$ over the unit two-dimensional sphere in $R^{3}$. The system dynamics is governed by the system of ordinary differential equations

[TABLE]

where the control input $u$ is subject to the inequality constraint $|u-c|\leq r$, which we relax with the equality constraint

[TABLE]

The variable $u_{d}$ is fictitious and controlled by the scalar $w_{d}$ introduced below.

The cost function is $J=p-\int_{t}^{t_{f}}w_{d}u_{d}$, where $p=t_{f}-t$ is the time to destination, and $w_{d}$ is a small positive constant.

We choose the receding horizon coinciding with the interval $[t,t_{f}]$. The horizon is parameterized by the dimensionless time $\tau\in[0,1]$ by means of the linear mapping $\tau\to t+\tau p$. The normalized interval $[0,1]$ is partitioned uniformly into the grid $\tau_{i}=i\Delta\tau$, $i=0,1,\ldots,N$, with the step size $\Delta\tau=1/N$. The discretized variables include the state $x_{i}$ and costate $\lambda_{i}$, the control input $u_{i}$ and slack variable $u_{d,i}$, the Lagrange multipliers $\mu_{i}$ and $\nu$, the parameter $p$.

The uncertain predictive model of the dynamical system on the receding horizon is the forward Euler method

[TABLE]

where

[TABLE]

The truncation error of the Euler methods is the disturbance $\eta_{f}$ in (2). We remark that $\eta_{f}$ is not random here and highly correlated with the state function $x(\tau)$.

It is directly verified that the continuous system dynamics $\dot{x}=A(u)x$ satisfies the equality constraint on the state $x_{i}^{T}x_{i}-1=0$, $i=1,\ldots,N$. Hence the constraint (5) has $g(x_{i})=x_{i}^{T}x_{i}-1$ and $\eta_{g}=0$. The goal of the least squares minimization is to satisfy the constraint (5) “softly.” We note that for this test problem it is possible to satisfy the state constraint $x_{i}^{T}x_{i}-1=0$ exactly by projecting $x_{i+1}$ onto the unit sphere after every step of (20); see, e.g., [6].

Yet another way in this example to satisfy the equality constraint $x_{i}^{T}x_{i}-1=0$ is to use the so-called exponential integrator $x_{j+1}=\exp\left(A(u_{j})x_{j}\right)$, which preserves the norm $\|x_{j}\|_{2}$. We use this exponential integrator for numerical simulation of the system dynamics replacing measurements.

The discretized cost function is

[TABLE]

We choose our least squares approximation of the state $x_{i}$, with the fixed initial value $x_{0}$ and a scalar parameter $\beta\geq 0$,

[TABLE]

The parameter $\beta$ determines the force of satisfying the equality constraint $x_{i}^{T}x_{i}-1=0$: the larger the constant $\beta$ the larger the enforcement.

The least squares minimization problem is equivalent to the system of nonlinear equations

[TABLE]

[TABLE]

where

[TABLE]

The corresponding discrete Lagrangian function then has the following form

[TABLE]

The costate $\lambda$ satisfies the formula

[TABLE]

where $C$ is the block diagonal matrix given by

[TABLE]

The function $F(U,x_{0},t)$, where

[TABLE]

has the following rows from the top to bottom:

[TABLE]

The example is chosen here for historical reasons—one of tests from our prior work [6]. It is not the most beneficial one illustrating effectiveness of the proposed least squares fit of the state with uncertain dynamics and constraints over the horizon, because in this example the state constraint to the sphere can in practice be actually certain, and satisfied with high accuracy by other means; e.g., [6]. The role of this example is a proof of concept.

Remark. A very important circumstance arises in the problems with the state constraints derived from the system dynamics. The number of terminal constraints must be reduced to the dimension of the smooth manifold determined by the equality constraint on the state. In our case, the dimension of the sphere equals 2, and, therefore, the Lagrange multiplier $\nu$ must contain only 2 components instead of 3. In our MATLAB implementation, we keep the components of $\nu$ corresponding to the $x$ and $y$ coordinates of the terminal state, but the last component of the right-hand side in the equation for the costate is set to zero. If the above described reduction of the terminal constraint is not fulfilled, then the subsequent computations lead to singular Jacobians in the Newton-type iterations.

VII Numerical results

We perform several preliminary numerical experiments in MATLAB with the test problem from § VI. Problem (21) is solved by the MATLAB function lsqnonlin for nonlinear least squares problems. The operator equation (19) is solved by the gmres function of MATLAB. The relative error tolerance for the GMRES iterations is $tol=10^{-5}$. The number of grid points on the horizon is $N=10$, the sampling time of the simulation is $\Delta t=1/200$, and $h=10^{-8}$.

Other constants are as follows: $c=0.5$, $r=0.1$, $w_{d}=0.005$, $\beta=10$.

The initial value for $U_{0}$ is computed by the MATLAB function fsolve, which finds a solution to nonlinear equations by Newton-type methods. We note that finding good initial approximation for fsolve may be non-trivial.

The trajectory satisfying the system dynamics $\dot{x}=A(u)x$ has been computed by the simple exponential integrator $x_{j+1}=\exp\left(A(u_{j})x_{j}\right)$ substituting the measurements.

Figure 1 shows the computed trajectory $(x,y,z)$ on the sphere for the test example. Figure 2 left panel shows the $(x,y)$-projection of the computed trajectory. Figure 2 right panel shows the input control variable with the constraints. Figure 3 displays the number of GMRES iterations at the grid points. Finally, Figure 4 displays the 2-norm of the residual function $F[U]$ that is supposed to vanish.

Since our implementation is not optimized, we do not provide the timing. However, we note that computations by the function lsqnonlin are relatively time consuming. We also observe numerically that a successful execution of Newton-type iterations requires solution of the least squares minimization problem with sufficiently high accuracy.

Comparison with other relevant numerical methods based on, e.g., multiple shooting [8] is a subject for future research.

Conclusions

A novel concept of least squares relaxation of state dynamics and some constraints in NMPC calculations of control over the horizon is proposed via alternating minimization and implemented in Newton-Krylov iterative methods. Numerical results for a proof of concept example demonstrate feasibility of computer implementations of the proposed technology. Future research is needed to design numerically efficient algorithms and to test our techniques for large uncertainties.