Preconditioned warm-started Newton-Krylov methods for MPC with discontinuous control

TL;DR

This paper introduces preconditioned Newton-Krylov methods tailored for solving optimal control problems with discontinuous controls in model predictive control, demonstrating efficient computation and the benefits of warm-start strategies.

Contribution

The work develops and tests Newton-Krylov methods with preconditioning and warm-start techniques for discontinuous control problems in MPC, achieving near-optimal computational complexity.

Findings

GMRES with preconditioning achieves O(N) complexity.

Warm-starting along the horizon improves convergence.

Method successfully applied to a double integrator minimum-time problem.

Abstract

We present Newton-Krylov methods for efficient numerical solution of optimal control problems arising in model predictive control, where the optimal control is discontinuous. As in our earlier work, preconditioned GMRES practically results in an optimal $O(N)$ complexity, where $N$ is a discrete horizon length. Effects of a warm-start, shifting along the predictive horizon, are numerically investigated. The~method is tested on a classical double integrator example of a minimum-time problem with a known bang-bang optimal control.

Figures

Equations

$$\frac{d}{d\tau }x(\tau )=f(\tau ,x(\tau ),u(\tau ),p),0<\tau <1,$$

$$x(0)=x_0,\psi (\tau ,x(\tau ),p){|}_{\tau =1}=0,$$

$$C(\tau ,x(\tau ),u(\tau ),u_d(\tau ),p)=0,0\le \tau \le 1,$$

$$J=\varphi ({\tau }_N,x_N,p)+\sum _{i=0}^{N-1}L({\tau }_i,x_i,u_i,u_{d,i},p)\Delta \tau ,$$

$$\underset{u_i,p}{\min }J,$$

$$x_{i+1}=x_i+f({\tau }_i,x_i,u_i,p)\Delta \tau ,i=0,1,\dots ,N-1,$$

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

$$\psi ({\tau }_N,x_N,p)=0.$$

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

$$X=[x_0^T,\dots ,x_N^T,{\lambda }_0^T,\dots ,{\lambda }_N^T{]}^T$$

$$U=[u_0^T,u_{d,0}^T,{\mu }_0^T\dots ,u_i^T,u_{d,i}^T,{\mu }_i^T,\dots ,$$

$$u_{N-1}^T,u_{d,N-1}^T,{\mu }_{N-1}^T,{\nu }^T,p^T{]}^T.$$

$$\frac{\partial {\mathcal{L}}^T}{\partial X}(X,U)=0\text{mbox}and\frac{\partial {\mathcal{L}}^T}{\partial U}(X,U)=0.$$

$$H(t,x,\lambda ,u,u_d,\mu ,p)=L(t,x,u,u_d,p)$$

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

$$x_{i+1}=x_i+f({\tau }_i,x_i,u_i,p)\Delta \tau ,i=0,\dots ,N-1.$$

$${\lambda }_N=\frac{\partial {\varphi }^T}{\partial x}({\tau }_N,x_N,p)+\frac{\partial {\psi }^T}{\partial x}({\tau }_N,x_N,p)\nu$$

$${\lambda }_i={\lambda }_{i+1}+\frac{\partial H^T}{\partial x}({\tau }_i,x_i,{\lambda }_{i+1},u_i,u_{d,i},{\mu }_i,p)\Delta \tau .$$

$$a(V)=(F[U^{(0)}(t)+hV,x(t),t]-F[U^{(0)}(t),x(t),t]).$$

$$ha(\Delta U/h)=b,\text{mbox}whereb=-F[U^{(0)}(t),x(t),t].$$

$$a(\Delta U)=b\text{mbox}whereb=-F[U^{(0)}(t),x(t),t]$$

$$U(t)=[u_0^T,u_{d,0}^T,{\mu }_0^T,\dots ,\dots ,{\nu }^T,p^T{]}^T,$$

$$a_j(V)=\left(F[U_{j-1}+hV,x_j,t_j]-F[U_{j-1},x_j,t_j]\right)/h,$$

$$A_j{e}_k=a_j(e_k),k=1,\dots ,m,$$

$$a_j(\Delta U_j)=b_j,\text{mbox}where{b}_j=-F[U_{j-1},x_j,t_j].$$

$$\underset{u,p}{\min }p\left[1-w_d\sum _{i=0}^{N-1}\Delta \tau u_d(i)\right],$$

$$\left[\begin{array}[]{c}x(i+1)\\ y(i+1)\end{array}\right]=\left[\begin{array}[]{c}x(i)\\ y(i)\end{array}\right]+p\Delta\tau\left[\begin{array}[]{c}y(i)\\ u(i)\end{array}\right],$$

$$x(0)=x_t,y(0)=y_t,x(N)=x_f,y(N)=y_f,$$

$$u^2(i)+u_d^2(i)=1.$$

$$\underset{u,p}{\min }\frac{{\alpha }_1}{2}(x(N)-x_f{)}^T(x(N)-x_f)$$

$$\left[\begin{array}[]{c}x(i+1)\\ y(i+1)\end{array}\right]=\left[\begin{array}[]{c}x(i)\\ y(i)\end{array}\right]+p\Delta\tau\left[\begin{array}[]{c}y(i)\\ u(i)\end{array}\right],$$

$$x(0)=x_t,y(0)=y_t,$$

Full text

1 Introduction.

We present and study numerically an approach to solution of minimum-time problems by the Model Predictive Control (MPC) method. The approach has been designed in our previous papers [9]-[13]. We apply it here to the double integrator problem, which is a classical example of the minimum-time optimal control problem. One of the main difficulties for numerical solution of this problem is caused by that the optimal control input is a function with large discontinuities, which makes the continuation procedures along the trajectories harder for accurate computation.

The minimum-time optimal control problems have been well studied; see e.g., [1, 2, 3, 6, 28]. In simple cases such as a double integrator with constraints on the control input [15], a closed-form of the feedback law can be found analytically. In general cases, obtaining a closed-form solution is a challenging task and a minimum-time open-loop trajectory is generated numerically by applying either direct methods, such as the direct shooting method, penalty function method, and SQP, or indirect methods such as the multiple shooting method and collocation method [4, 7, 14, 17, 18, 20],[23]-[27],[29]. To provide robustness to unmeasured uncertainties, the open-loop control can be augmented with a feedback stabilizer to the computed open-loop trajectory.

Applying the MPC philosophy to the minimum-time control involves recomputing the open-loop state and control trajectory subject to pointwise-in-time state and control constraints, terminal state constraint, and the current state as the initial conditions. The computed control trajectory is applied at the next time step; see, e.g., [16, 21, 30].

We use a time scaling linear transformation to obtain a fixed-end time problem, with the terminal time appearing as a multiplicative parameter in the dynamic equations. We then convert the model to a discrete time and formulate a prediction problem, where both the control input sequence and the horizon length, now appearing as a parameter, have to be optimized.

We use the continuation method for advancing the predictive horizon in time. Since the optimal control input is a piecewise continuous function and the discontinuities are severe, the perturbations of the unknown vectors during the continuation can be sufficiently large. There are two simple remedies for minimizing these perturbations: (1) application of the Newton-like iterations several times and (2) the use of suitable interpolations for current data when reusing them as a warm-start in the next horizon interval. In our numerical experiments, satisfactory results are obtained by the aid of both techniques. In some situations, the linear interpolation shows impressive improvements, in other situations, multiple application of the Newton-like iterations is the only tool for computations without failure.

In [13], we propose a new preconditioner technique for the Newton-like iterations. It works in the case, when the problem has only few state constraints, e.g., terminal state constraints. The double integrator problem satisfies these restrictions, and we present a construction of such preconditioner for this problem. The preconditioned iterations have the arithmetical cost $O(N)$, where $N$ is the number of grid points in the predictive horizon.

2 Formulation of the prediction problem

We want to develop a numerical method for solving the double integrator problem [2], where a bounded control input $|u|\leq 1$ drives the dynamic system $\ddot{\chi}=u$ from an initial state $\chi(t_{0}),\dot{\chi}(t_{0})$ to a final state $\chi(t_{f}),\dot{\chi}(t_{f})$ in minimum time $t_{f}-t_{0}$. We introduce the state variables $x_{1}=\chi$, $x_{2}=\dot{\chi}$ and convert the problem with the free terminal time $t_{f}$ to a problem with a fixed terminal time by the linear time scaling $\tau=(t-t_{0})(t_{f}-t_{0})$ such that $0\leq\tau\leq 1$. The system dynamics after scaling becomes \frac{d}{d\tau}\left[\begin{array}[]{c}x_{1}\\ x_{2}\end{array}\right]=(t_{f}-t_{0})\left[\begin{array}[]{c}x_{2}\\ u\end{array}\right]. We also replace the inequality $|u(\tau)|\leq 1$ by the equality constraint $u^{2}+u_{d}^{2}-1=0$, $0\leq\tau<1$, where $u_{d}(\tau)$ is a dummy (or slack) argument. Our goal is to minimize the performance index $J=p=t_{f}-t_{0}$ by choosing a suitable control input $u(\tau)$.

A general continuous-time model for the minimum time problem after time scaling reads as follows: given a performance index $J=\phi(\tau,x(\tau),p)|_{\tau=1}+\int_{0}^{1}L(\tau,x,u,u_{d},p)d\tau$, find $\min_{p,u}J$ subject to the equality constraints

[TABLE]

where $x\in R^{n_{x}}$ and $u\in R^{n_{u}}$ are the state and control input, respectively. The vector function $C$ determines equality constraints, where possible dummy arguments $u_{d}$ correspond to inequality constraints. For instance, an inequality constraint $c_{j}(t,x,u)\leq 0$ can be replaced by the equality constraint $c_{j}(t,x,u)+u_{d}^{2}=0$ or by $c_{j}(t,x,u)+exp(u_{d})=0$, and then a suitable penalty term, e.g., $-w_{d}p\int_{0}^{1}u_{d}d\tau$ must be added to the performance index $J$, where $w_{d}>0$ is a small constant required by the interior point method [29].

Let us consider a discretization of the continuous-time model on the uniform grid $\tau_{i}=i\Delta\tau$, where $\Delta\tau=1/N$. The adjustable variables are the control time sequence $u(0)$, $u(1)$, …, $u(N-1)$ and the parameter $p$ that need to be simultaneously optimized. The discrete-time optimal control problem can be recast in the following form: given a performance index

[TABLE]

find

[TABLE]

subject to a fixed initial value $x_{0}$ and the equality constraints

[TABLE]

[TABLE]

[TABLE]

From several choices for the functions $\phi$, $L$ and $\psi$ in the double integrator problem. we have tried two: (1) $\phi(\tau_{N},x_{N},p)=p$, $\psi(\tau_{N},x_{N},p)=x(\tau_{N})-x_{f}$, $L(\tau_{i},x_{i},u_{i},u_{d,i},p)=-w_{d}u_{d,i}p$, and (2) $\phi(\tau_{N},x_{N},p)=p+\frac{\beta}{2}x_{N}^{2}$, $L(\tau_{i},x_{i},u_{i},u_{d,i},p)=-w_{d}u_{d,i}p+\alpha u_{i}^{2}p$ for a small $\alpha>0$ and large $\beta>0$.

To solve the discrete-time optimal control problem by the classical Lagrange multiplier method, we introduce the costate $\lambda$ for the dynamics equations, the Lagrange multiplier $\mu$ for the equality constraints, and Lagrange multiplier $\nu$ for the terminal constraint $\psi$ and derive necessary optimality conditions are by means of the discrete Lagrangian

[TABLE]

where

[TABLE]

and

[TABLE]

The optimality conditions are given by the stationarity conditions

[TABLE]

Further derivations use the Hamiltonian function

[TABLE]

The optimality conditions are reformulated below as a system of nonlinear equations

[TABLE]

with respect to the unknown vector $U(t)$. We recall that in the MPC method $x_{t}$ denotes the current measured state vector, which serves as the initial value $x_{0}$ in the prediction problem, and $t$ is the current system time. The vector function $F(U,x,t)$ is obtained from the optimality conditions after elimination of all states $x_{i}$ and costates $\lambda_{i}$ by the following procedure.

Having the current measured state $x_{t}$, set $x_{0}=x_{t}$ and compute $x_{i}$, $i=1,2\ldots,N$, by the forward recursion

[TABLE]

Then starting with the value

[TABLE]

compute the costate $\lambda_{i}$, $i=N\!-\!1,\ldots,0$, by the backward recursion

[TABLE] 2. 2.

Using the values $x_{i}$ and $\lambda_{i}$, $i=0,1\ldots,N$, computed by the forward and backward recursions, compute the vector function $F[U,x,t]$ as the vector $\frac{\partial\mathcal{L}^{T}}{\partial U}(X,U)$. Explicit formulas for the components of $F[U,x,t]$ can be found in [13].

3 The Newton-Krylov method

The system of nonlinear equations (2.1) with respect to $U(t)$ must be solved in real time at each instance $t$ of the discrete system time. Assume that $U^{(0)}(t)$ is an approximation to the exact solution $U(t)$. Then $F[U(t),x(t),t]-F[U^{(0)}(t),x(t),t]=-F[U^{(0)}(t),x(t),t]$. Let us denote $b=-F[U^{(0)}(t),x(t),t]$ and $\Delta U=U(t)-U^{(0)}(t)$. For a sufficiently small scalar $h>0$, e.g. $h=10^{-8}$, which may be different from the system time step $\Delta t$ and from the horizon time step $\Delta\tau$, we introduce the finite difference operator

[TABLE]

Equation (2.1) is then equivalent to the operator equation

[TABLE]

We suppose that $h$ is sufficiently small such that the operator $a(V)$ is almost linear in the sense that $\|ha(\Delta U/h)-a(\Delta U)\|\leq\delta\|\Delta U\|$ with a tiny constant $\delta>0$. Thus, solution of (2.1) is reduced to solving the operator equation

[TABLE]

for a given approximation $U^{(0)}(t)$ and setting $U(t)=U^{(0)}(t)+\Delta U$.

Given formulas for computing the vector function $F[U,x_{t},t]$, we want to solver equation (2.1) at the points of the uniform grid $t_{j}=t_{0}+j\Delta t$, $i=0,1,\ldots$ by using the above proposed approach.

At the initial state $x_{0}=x(t_{0})$, we find an approximate solution $U_{0}$ to the equation $F[U_{0},x_{0},t_{0}]=0$ by a suitable optimization procedure. The dimension of the vector $u(t)$ is denoted by $n_{u}$. Since

[TABLE]

the first block entry of $U_{0}$, formed from the first $n_{u}$ elements of $U_{0}$, is taken as the control $u_{0}$ at the state $x_{0}$. The next state $x_{1}=x(t_{1})$ is either measured by a sensor or computed by the formula $x_{1}=x_{0}+\Delta tf(t_{0},x_{0},u_{0})$.

The computation of $U_{j}$ for $j>0$ is as follows. By $e_{k}$, we denote the $k$-th column of the $m\times m$ identity matrix, where $m$ is the dimension of the vector $U_{j}$. At time $t_{j}$, we arrive with the state $x_{j}$ and the vector $U_{j-1}$. The vector $U_{j-1}$ serves as the approximation $U^{(0)}_{j}$, and such an approximation is often called the warm start. The difference operator

[TABLE]

implicitly determines an $m\times m$ matrix $A_{j}$ with the columns

[TABLE]

which approximates the symmetric Jacobian matrix $F_{U}[U_{j-1},x_{j},t_{j}]$ so that $a_{j}(V)=A_{j}V+O(h)$. At the current time $t_{j}$, our goal is to solve the equation

[TABLE]

Then we set $U_{j}=U_{j-1}+\Delta U_{j}$ and choose the first $n_{u}$ components of $U_{j}$ as the control $u_{j}$. The next state $x_{j+1}=x(t_{j+1})$ either comes from a sensor, estimated, or computed by the formula $x_{j+1}=x_{j}+\Delta tf(t_{j},x_{j},u_{j})$.

Equation (3.2) can be solved approximately by generating the matrix $A_{j}$ and then solving the system of linear equations $A_{j}\Delta U_{j}=b_{j}$ using a direct method, e.g., the Gaussian elimination.

An alternative method of solving (3.2) is by a Krylov subspace iteration, e.g., by GMRES method [22], for which we do not generate the matrix $A_{j}$ explicitly. Namely, we simply use the operator $a_{j}(V)$ instead of computing the matrix-vector product $A_{j}V$, for arbitrary vectors $V$; cf., [8]. The arithmetic cost of the alternative method can be often less than that of the direct method because the matrix $A_{j}$ is not generated and the $O(N^{3})$ cost due to the Gaussian elimination is avoided.

We refer to the Newton-like refinement by solving equation (3.2) with the GMRES method as the Newton-Krylov method.

Remark 1.

The above numerical procedure uses the so-called warm start. Specifically, the operator $a_{j}(V)$ and the right-hand-side $b_{j}$ are computed at time $t_{j}$ using the value $U_{j-1}$ from the previous time $t_{j-1}$ as an approximation $U^{(0)}_{j}$. This approximation does not take advantage of the structure of the underlying problem. In our case, the vector $U(t)$ contains the components $u(\tau_{i})$, $u_{d}(\tau_{i})$, $\mu(\tau_{i})$, $\nu_{1}$, $\nu_{2}$ and $p$ computed in the prediction horizon $[0,1]$. It is possible to apply, e.g., the linear interpolation of the variables $u(\tau_{i})$, $u_{d}(\tau_{i})$, $\mu(\tau_{i})$ from the previous horizon interval $[t_{j-1},t_{j-1}+p_{j-1}]$ to the current horizon interval $[t_{j},t_{j}+p_{j}]$, where $p_{j-1}$ equals $p$ from $U_{j-1}$ and $p_{j}=p_{j-1}-\Delta t$. Such an interpolation is referred to as the shifting along the predictive horizon; see a related discussion in [4].

Remark 2.

The shifting along the predictive horizon, described in Remark 1, is a light improvement, which may help in some situations. A more cardinal improvement of the warm start is solving the operator equation $a(\Delta U)=b$ and subsequent correction $U+\Delta U$ several times, say, 2-3 Newton-like iterations instead of one.

4 Sparse preconditioner

Convergence of the GMRES iteration can be rather slow if the matrix $A_{j}$ is ill-conditioned. However, the convergence can be improved by preconditioning; see e.g. [22]. A matrix $T$ that approximates a matrix $A^{-1}$ and such that computing the product $Tr$ for an arbitrary vector $r$ is relatively easy, is referred to as a preconditioner. The preconditioning for the system of linear equations $Ax=b$ formally replaces the original system $Ax=b$ with the equivalent preconditioned linear system $TAx=Tb$. If the condition number $\kappa(TA)=\|TA\|\|A^{-1}T^{-1}\|$ of the matrix $TA$ is small, convergence of iterative solvers for the preconditioned system become fast. Instead of $T$, its inverse $M=T^{-1}$ is often called a preconditioner.

Our specific optimization problem over a predictive horizon admits construction of efficient preconditioners $M_{j}$ with a sparse structure. These preconditioners have been introduced in [13], and we refer to this paper for a more detailed description.

The construction of $M_{j}$ is based on several observations. The first one is the fact that the Jacobi matrix $F_{U}$ is the Schur complement of the Jacobi matrix for the Lagrangian $\mathcal{L}$. If we denote the result of the forward and backward recursions by $X=g(U)$, then most entries of $F_{U}$ are approximated by the matrix $\mathcal{L}_{UU}(g(U),U)$ due to the second fact, which affirms that the states $x_{i}$, computed by the forward recursion, and the costates $\lambda_{i}$, computed by the subsequent backward recursion, satisfy the following property: $\partial x_{i_{1}}/\partial u_{i_{2}}=O(\Delta\tau)$, $\partial\lambda_{i_{1}}/\partial u_{i_{2}}=O(\Delta\tau)$, $\partial x_{i_{1}}/\partial\mu_{i_{2}}=0$ and $\partial\lambda_{i_{1}}/\partial\mu_{i_{2}}=O(\Delta\tau)$. The second fact is a corollary of theorems about the derivatives of solutions of ordinary differential equations with respect to a parameter; see, e.g., [19, 5]. The matrix $\mathcal{L}_{UU}(g(U),U)$ is very sparse and easily computable. We have to generate explicitly by the formula $A(e_{k})$ only the last $l$ columns and rows of $F_{U}$, where the integer $l$ denotes the sum of dimensions of $\psi$ and $p$.

As a result, the setup of $M_{j}$, computation of its LU factorization, and application of the preconditioner all cost $O(N)$ floating point operations. The memory requirements are also of order $O(N)$.

5 Simulations of a Double Integrator System

To evaluate the effectiveness of the Newton-Krylov method, we compute a solution to the double integrator system with control input constraints: $\min_{u(t)}T$ subject to $\ddot{\chi}=u$, $\chi(0)=x_{0}$, $\dot{\chi}(0)=y_{0}$, $\chi(T)=x_{f}$, $\dot{\chi}(T)=y_{f}$, $|u|\leq 1$. We consider two discrete-time models of this problem and solve both numerically by the MPC approach, where the predictive horizon is the interval $[t,T]$, $t$ being the current system time.

Model 1 of the predictive problem on the scaled horizon has the form

[TABLE]

subject to

[TABLE]

[TABLE]

[TABLE]

Model 2 of the predictive problem on the scaled horizon has the form

[TABLE]

subject to

[TABLE]

[TABLE]

[TABLE]

We choose $\alpha_{1}=10^{3}$ and $\alpha_{2}=0.1$ in Model 2. The initial state is $x_{0}=-1$, $y_{0}=0$, and the terminal state is $x_{f}=0$, $y_{f}=0$ in both models. The interior point penalty parameter $w_{d}=0.005$ and is fixed as suggested, e.g., in paper [29].

Below we give detailed formulas for an implementation of Model 1. Similar formulas can be easily derived for Model 2.

The components of the discretized problem on the predictive scaled horizon $[0,1]$ are as follows:

• •

$\Delta\tau=1/N$, $\tau_{i}=i\Delta\tau$;

• •

the unknown variables are the state \left[\begin{array}[]{c}x_{i}\\ y_{i}\end{array}\right], the costate \left[\begin{array}[]{c}\lambda_{1,i}\\ \lambda_{2,i}\end{array}\right], the control $u_{i}$, the dummy variable $u_{d,i}$, the Lagrange multipliers $\mu_{i}$ and \left[\begin{array}[]{c}\nu_{1}\\ \nu_{2}\end{array}\right], the parameter $p$;

• •

the state is determined by the forward recursion

[TABLE]

where $i=0,1,\ldots,N-1$;

• •

the costate is determined by the backward recursion ($\lambda_{1,N}=\nu_{1}$, $\lambda_{2,N}=\nu_{2}$)

[TABLE]

where $i=N-1,N-2,\ldots,0$;

• •

the system of equations $F(U,x_{t},y_{t},t)=0$, where

[TABLE]

has the following rows from the top to bottom:

[TABLE]

[TABLE]

[TABLE]

The sparse preconditioner $M_{j}$ is a symmetric matrix with the sparsity structure

[TABLE]

where $M_{11}$ is a block diagonal matrix with the nonsingular diagonal $3\times 3$-blocks, $i=0,\ldots,N-1$,

[TABLE]

The matrix block \left[\begin{array}[]{c}M_{12}\\ M_{22}\end{array}\right] consists of 3 columns and is computed by means of the operator $a_{j}(V)$ as

[TABLE]

where $e_{3N+1},e_{3N+2},e_{3N+3}$ are the columns $3N+1$, $3N+2$ and $3N+3$ of the unit matrix of order $3N+3$. The block $M_{21}$ is transpose to $M_{12}$, $M_{21}=M_{12}^{T}$.

The LU factorization of $M_{j}$ can be computed with $O(N)$ floating point operations, e.g., as

[TABLE]

where $S_{22}=M_{22}-M_{21}M_{11}^{-1}M_{12}$. The setup and application of the preconditioner also cost $O(N)$ operations.

6 Numerical results

All the plots are computed by the GMRES method with the absolute tolerance $10^{-6}$.

Figures 1, 2, 3, 4 show the numerical results for Model 1, where the horizon length $N$ equals 20. The exact and computed time to destination is $t_{f}=2$. The system time between $t_{0}=0$ and $t_{f}$ is uniformly partitioned into 500 subintervals with the time step $t_{j+1}-t_{j}=2/500$. Our MATLAB implementation fails when only 1 Newton’s refinement is used at each time instance $t_{j}$, and it works with 2 Newton’s refinements per step. The failure is not cured by the shifting along the predictive horizon described in Remark 1. When the interval $[0,t_{f}]$ is partitioned into 1000 subintervals, our MATLAB program works with 1 Newton’s refinement at each step $t_{j}$.

Figure 3 shows the Euclidean norm of the residual $F[U_{j-1},x_{j},t_{j}]$ denoted by $\|b\|_{2}$ and the Euclidean norm $\|F\|_{2}$ of the residual $F[U_{j},x_{j},t_{j}]$ after all Newton’s refinemens at time $t_{j}$ (2 refinements in Figure 3).

Figure 4 displays the total number of GMRES iterations at each time $t_{j}$ when the GMRES method is applied without preconditioner. The average number of GMRES iterations per step is about 77. The GMRES method with preconditioning uses only 2 iterations per step thus demonstrating more than the 35x speedup. We recall that our preconditioning provides the $O(N)$ complexity of the prediction problem.

The interpolation by shifting along the predictive horizon from Remark 1 does not give an essential improvement in accuracy for Model 1 with the chosen time steps in the horizon and the system time. The multiple application of Newton’s refinements, however, provides an essential improvement in accuracy and computation without failure.

Figures 5–9 show the numerical results for Model 2. Again, the computed time to destination equals $t_{f}=2$, but the system time between $t_{0}=0$ and $t_{f}$ is partitioned into 200 subintervals. The predictive horizon has the length $N=70$. Our MATLAB program succeeds in computing the control input and trajectory with only 1 Newton’s refinement per step by using the shifting along the predictive horizon. The program fails without the shifting interpolation. The number of GMRES iterations without preconditioning is shown in Fig. 9 and its average value per step is about 35. The number of GMRES iterations with preconditioning shown in Fig. 8 is about 2.5 per step. Hence the speedup is about 14x.

We have run our program for Model 2, when $N=20$, the time interval $[0,t_{f}]$ is partitioned into 200 subintervals, and only 1 Newton’s refinement is used. The program fails to compute without the shifting interpolation and works with the interpolation. Figure 10 displays the computed control input.

Models 1 and 2 demonstrate different behavior in our numerical experiments. Model 2 produces the control input with a smoother approximation of the discontinuity than Model 1 in spite of that its predictive horizon contains 70 grid points against 20 grid points in Model 1. Model 1 requires application of several Newton’s iterations per time step and the shifting interpolation is not enough for the work without failure. Model 2 works perfectly with the shifting interpolation and fails with it. The multiple Newton refinements help for Model 2 but the number of required refinements is not small. Thus, Model 1 with the chosen set of parameters requires multiple Newton iterations, and Model 2 with its own set of parameters requires the shifting interpolation along the horizon.

7 Conclusion

Our numerical experiments demonstrate that the Newton-Krylov method successfully works even in the cases, when the control input is discontinuous. Both the multiple Newton refinement and the interpolation by shifting along the predictive horizon can essentially improve accuracy of the computed solution and avoid failure during computations. The sparse preconditioner gives very good acceleration, e.g., a speedup about 15-30 times. A future research may include more advanced numerical examples and deeper study of the interpolation by shifting along the predictive horizon.

Acknowledgements.

We are very grateful to Ilya Lashuk for his guidance and help in using Model 2.