A Parallel Decomposition Scheme for Solving Long-Horizon Optimal Control Problems

TL;DR

This paper introduces a parallel temporal decomposition method for long-horizon optimal control problems, ensuring convergence under specific conditions related to boundary perturbation decay and overlap size, with demonstrated effectiveness on quadratic and non-convex problems.

Contribution

It proposes a novel parallel decomposition scheme with convergence guarantees based on boundary perturbation decay, applicable to both convex and non-convex problems.

Findings

Scheme converges if boundary perturbations decay asymptotically.

LQ problems satisfy the decay condition, ensuring convergence.

Overlap size influences convergence rate and scheme performance.

Abstract

We present a temporal decomposition scheme for solving long-horizon optimal control problems. In the proposed scheme, the time domain is decomposed into a set of subdomains with partially overlapping regions. Subproblems associated with the subdomains are solved in parallel to obtain local primal-dual trajectories that are assembled to obtain the global trajectories. We provide a sufficient condition that guarantees convergence of the proposed scheme. This condition states that the effect of perturbations on the boundary conditions (i.e., initial state and terminal dual/adjoint variable) should decay asymptotically as one moves away from the boundaries. This condition also reveals that the scheme converges if the size of the overlap is sufficiently large and that the convergence rate improves with the size of the overlap. We prove that linear quadratic problems satisfy the asymptotic decay condition, and we discuss numerical strategies to determine if the condition holds in more general cases. We draw upon a non-convex optimal control problem to illustrate the performance of the proposed scheme.

Equations

$$\underset{\begin{array}{c}{x}_{M:N}\\ u_{M:N-1}\end{array}}{\min }$$

$$x_M=x_M^{\ast }({\lambda }_M)$$

$$x_i=f(x_{i-1},u_{i-1})({\lambda }_i)\forall i\in {\mathbb{I}}_{M+1:N}$$

$$g(x_i,u_i)\le 0({\mu }_i),\forall i\in {\mathbb{I}}_{M:N-1}$$

$$M_k^{\omega }=\max (M_k-\omega ,M),N_k^{\omega }=\min (N_k+\omega ,N)$$

$$z_{M_k:N_k}^{(t+1)}\leftarrow {\mathcal{P}}_{M_k^{\omega }:N_k^{\omega }}(x_{M_k^{\omega }}^{(t)},{\lambda }_{N_k^{\omega }}^{(t)}),k\in {\mathbb{I}}_{1:K},t\in {\mathbb{Z}}_{>0}.$$

$${\mathcal{L}}_{M:N}(x_{M:N},u_{M:N-1},{\lambda }_{M:N},{\mu }_{M:N-1};x_M^{\ast },{\lambda }_N^{\ast })$$

$$:=\sum _{i=M}^{N-1}\ell (x_i,u_i)+{\lambda }_M^{\top }(x_M-x_M^{\ast })+({\lambda }_N^{\ast }{)}^{\top }{x}_N$$

$$+\sum _{i=M+1}^N{\lambda }_i^{\top }(x_i-f(x_{i-1},u_{i-1}))+\sum _{i=M}^{N-1}{\mu }_i^{\top }g(x_i,u_i)$$

$${\lambda }_i^{\ast }+{\left({\nabla }_x{f}_i^{\ast }\right)}^{\top }{\lambda }_{i+1}^{\ast }+{\nabla }_x{\ell }_i^{\ast }+{\left({\nabla }_x{g}_i^{\ast }\right)}^{\top }{\mu }_i^{\ast }$$

$${\left({\nabla }_u{f}_i^{\ast }\right)}^{\top }{\lambda }_{i+1}^{\ast }+{\nabla }_u{\ell }_i^{\ast }+{\left({\nabla }_u{g}_i^{\ast }\right)}^{\top }{\mu }_i^{\ast }$$

$$x_{i+1}^{\ast }$$

$$g_i^{\ast }\le 0,{\mu }_i^{\ast }\ge 0,\mathrm{diag}({\mu }_i^{\ast })g_i^{\ast }$$

$$\begin{array}{rl}& \Vert z_i^{\ast }-{\tilde{z}}_i{\Vert }_{\infty }\le {\epsilon }_{i-M}\Vert x_M^{\ast }-{\tilde{x}}_M{\Vert }_{\infty }+{\epsilon }_{N-i}\Vert {\lambda }_N^{\ast }-{\tilde{\lambda }}_N{\Vert }_{\infty }\end{array}$$

$$\Vert z_{M:N}^{(t)}-z_{M:N}^{\ast }{\Vert }_{\infty }\le {\left(2{\epsilon }_{\omega }\right)}^t\Vert z_{M:N}^{(0)}-z_{M:N}^{\ast }{\Vert }_{\infty }$$

$$\Vert z_i^{(t+1)}-z_i^{\ast }{\Vert }_{\infty }\le$$

$${\epsilon }_{i-M_k^{\omega }}\Vert x_{M_k^{\omega }}^{(t)}-x_{M_k^{\omega }}^{\ast }{\Vert }_{\infty }+{\epsilon }_{N_k^{\omega }-i}\Vert {\lambda }_{N_k^{\omega }}^{(t)}-{\lambda }_{N_k^{\omega }}^{\ast }{\Vert }_{\infty }$$

$$\Vert z_i^{(t+1)}-z_i^{\ast }{\Vert }_{\infty }$$

$${\Vert z_i^{(t+1)}-z_i^{\ast }\Vert }_{\infty }\le 2{\epsilon }_{\omega }{\Vert (x_{M_k^{\omega }}^{(t)},{\lambda }_{N_k^{\omega }}^{(t)})-(x_{M_k^{\omega }}^{\ast },{\lambda }_{N_k^{\omega }}^{\ast })\Vert }_{\infty }.$$

$$\Vert z^{(t+1)}-z^{\ast }{\Vert }_{\infty }\le 2{\epsilon }_{\omega }\Vert z_{M:N}^{(t)}-z_{M:N}^{\ast }{\Vert }_{\infty }$$

$$=({\nabla }_x{f}_{N_k}^{(t)}{)}^{\top }({\lambda }_{N_k+1}^{(t)}-{\hat{\lambda }}_{N_k+1}^{(t,k)}).$$

$$=({\nabla }_u{f}_{N_k}^{(t)}{)}^{\top }({\lambda }_{N_k+1}^{(t)}-{\hat{\lambda }}_{N_k+1}^{(t,k)}).$$

$$=x_{N_k+1}^{(t)}-{\hat{x}}_{N_k+1}^{(t,k)}$$

$$r^{(t)}:=\underset{k\in {\mathbb{I}}_{1:K-1}}{\max }\Vert {\hat{x}}_{N_k+1}^{(t,k)}-x_{N_k+1}^{(t)}{\Vert }_{\infty }$$

$$s^{(t)}:=\underset{k\in {\mathbb{I}}_{1:K-1}}{\max }\Vert {\hat{\lambda }}_{N_k+1}^{(t,k)}-{\lambda }_{N_k+1}^{(t)}{\Vert }_{\infty },$$

$$r^{(t)}<{ϵ}_{\text{pr,tol}},s^{(t)}<{ϵ}_{\text{du,tol}}.$$

$$\left[\begin{array}{cc}H& G^{\top }\\ G& \end{array}\right]\left[\begin{array}{c}w\\ \lambda \end{array}\right]=\left[\begin{array}{c}{\zeta }^{\ast }\\ {\xi }^{\ast }\end{array}\right],$$

$$H$$

$$G$$

$$Z:=\left[\begin{array}{ccccccc}& & & & & & \\ 1& & & & & & \\ {\overline{x}}_1& & & & & & \\ {\overline{u}}_1& 1& & & & & \\ ⋮& ⋮& & & & & \\ {\overline{x}}_{n_x}& {\overline{x}}_{n_x-1}& \ddots & & & & \\ {\overline{u}}_{n_x}& {\overline{u}}_{n_x-1}& \cdots & 1& & & \\ & {\overline{x}}_{n_x}& \cdots & {\overline{x}}_1& & & \\ & {\overline{u}}_{n_x}& \cdots & {\overline{u}}_1& 1& & \\ & & \ddots & ⋮& ⋮& & \\ & & & {\overline{x}}_{n_x}& {\overline{x}}_{n_x-1}& \ddots & \\ & & & {\overline{u}}_{n_x}& {\overline{u}}_{n_x-1}& \cdots & 1\\ & & & & {\overline{x}}_{n_x}& \cdots & {\overline{x}}_1\end{array}\right]$$

$$Y:=\left[\begin{array}{ccccccc}I& & & & & & \\ U_1& & & & & & \\ X_2& I& & & & & \\ ⋮& ⋮& & & & & \\ U_{n_x}& U_{n_x-1}& \ddots & & & & \\ & X_{n_x}& \cdots & I& & & \\ & U_{n_x}& \cdots & U_1& & & \\ & & & X_1& I& & \\ & & \ddots & ⋮& ⋮& & \\ & & & U_{n_x}& U_{n_x-1}& \ddots & \\ & & & & X_{n_x}& \cdots & I\end{array}\right]$$

$$\underset{p\in {\mathbb{R}}^{N-1}}{\min }(Zp+Y{\xi }^{\ast }{)}^{\top }H(Zp+Y{\xi }^{\ast })-(Z^{\top }{\zeta }^{\ast }{)}^{\top }p.$$

$$p^{\ast }={\left(Z^{\top }HZ\right)}^{-1}\left(Z^{\top }{\zeta }^{\ast }+Z^{\top }HY{\xi }^{\ast }\right).$$

$$w^{\ast }=$$

$${\lambda }^{\ast }=$$

$$w^{\ast }-\tilde{w}=$$

$${\lambda }^{\ast }-\tilde{\lambda }=$$

$$L_H$$

Full text

I Introduction

Long-horizon optimal control problems (OCPs) arise in model predictive control (MPC) applications such as chemical process systems [1], autonomous vehicle steering [2], and battery systems [3]. They also appear in other application domains such as chemical production planning [4] and electricity production planning [5]. Different decomposition techniques have been reported in the literature to improve computational tractability of these problems including dual decomposition [6], alternating direction method of multipliers [7], dual dynamic programming [8], Gauss-Seidel schemes [9, 10], and parallel Newton schemes [11]. Such decomposition techniques allow scalable solutions of long-horizon OCPs by the use of parallel computers.

In this work, we study the convergence properties of a new decomposition paradigm that uses overlapping time domains. This approach is motivated by overlapping Schwartz schemes used for the solution of partial differential equations (PDEs) [12, 13]. In the proposed scheme, the time domain is decomposed into a set of partially overlapping subdomains. Subproblems associated with subdomains are solved in parallel to obtain local primal-dual trajectories by using current primal-dual information from the neighboring subdomains. The subdomain trajectories are then assembled (this can be interpreted as a projection operator) to update the primal-dual trajectories over the entire domain, and the procedure is repeated. We provide a sufficient condition that guarantees convergence of the proposed decomposition algorithm when applied to OCPs. This condition applies to linear and nonlinear problems. Specifically, the condition indicates that convergence of the decomposition scheme is guaranteed provided that the sensitivity of the primal-dual (state-adjoint) trajectories to perturbations in initial states and terminal cost gradients decay asymptotically as one moves away from the boundaries. We call this condition asymptotic decay of sensitivity (ADS). The condition also reveals that the algorithm converges provided that the size of the overlap is sufficiently large and that the convergence rate improves as the size of the overlapping regions increases.

ADS-like properties have been recently explored in the literature. Xu et. al. recently showed that an ADS condition (in the primal space) holds under a time-varying and inequality-constrained linear quadratic (LQ) control setting [15, 16]. To prove this, the authors assumed uniformly complete controllability and exploited the algebraic structure of the Riccati equation. The authors used the primal ADS property to show that trajectories of an overlapping temporal decomposition scheme approximate those of the long-horizon problem and that the approximation error decays as the size of the overlap increases. The authors also showed that receding horizon control provides approximate trajectories and that the error converges as the size of the overlapping regions increase. These works do not provide an algorithmic scheme that delivers optimal trajectories (for a given size of the overlapping region), as we do in this work. Shin et. al. recently established a primal ADS condition for general graph-structured, unconstrained quadratic programs and showed that this condition guarantees the convergence of an iterative overlapping decomposition algorithm [14]. The condition established in this work is specialized to constrained OCPs and operates in the primal-dual space. A primal ADS property has also been established in [17, Lemma 5]. This condition is exploited by the authors to establish the stability of economic MPC.

The paper is organized as follows. In Section II, we present the basic setting and describe the proposed decomposition scheme. In Section III, we propose a primal-dual ADS condition that guarantees convergence of the algorithm. In Section IV, we show that ADS holds for a simplified LQ setting. In Section V, we demonstrate the proposed scheme using a nonlinear economic MPC problem.

II Basic Definitions and Setting

We consider an OCP with a time domain set $\mathbb{I}_{M:N}$, an initial state $x^{*}_{M}\in\mathbb{R}^{n_{x}}$, and a terminal cost gradient $\lambda^{*}_{N}\in\mathbb{R}^{n_{x}}$ of the form:

[TABLE]

We denote this problem as $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$. Here, $x_{i}\in\mathbb{R}^{n_{x}}$ and $u_{i}\in\mathbb{R}^{n_{u}}$ are the state and input variables at time $i$; $\lambda_{i}\in\mathbb{R}^{n_{x}}$ and $\mu_{i}\in\mathbb{R}^{n_{g}}$ are the dual variables associated with (1b)-(1c) and (1d), respectively. Symbol $f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}$ is the dynamic mapping, $\ell:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}$ is the stage cost mapping, and $g:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{g}}$ is the constraint mapping. We denote the set of real numbers and the set of integers as $\mathbb{R}$ and $\mathbb{Z}$, respectively, and we define $\mathbb{I}_{M:N}:=\mathbb{Z}\cap[M,N]$. We use the syntax $(x_{1},x_{2},\cdots,x_{n}):=\begin{bmatrix}x_{1}^{\top}&x_{2}^{\top}&\cdots&x_{n}^{\top}\end{bmatrix}^{\top}$ and $x_{M:N}:=(x_{M},x_{M+1},\cdots,x_{N})$. We denote primal-dual pairs as $z_{i}:=(x_{i},\lambda_{i})$.

Remark 1

From the KKT conditions of (1), it follows that $\lambda_{N}=\lambda^{*}_{N}$. Therefore, incorporating the terminal penalty term $(\lambda^{*}_{N})^{\top}x_{N}$ in the objective function as (1a) essentially constitute the terminal constraint of the dual variable $\lambda$ at $i=N$.

We now describe the proposed decomposition scheme for solving $\mathcal{P}_{M:N}(x_{M}^{*},\lambda_{N}^{*})$. We partition the time domain $\mathbb{I}_{M:N}$ into a collection of $K$ non-overlapping (i.e., disjoint) sets of the form $\mathbb{I}_{M_{1}:N_{1}},\cdots,\mathbb{I}_{M_{K}:N_{K}}$. We also define a collection of overlapping sets $\mathbb{I}_{M^{\omega}_{1}:N^{\omega}_{1}},\cdots,\mathbb{I}_{M^{\omega}_{K}:N^{\omega}_{K}}$ satisfying

[TABLE]

for any $k\in\mathbb{I}_{1:K}$. We call $\omega\in\mathbb{Z}_{>0}$ the size of the overlap. The non-overlapping and overlapping time domains with $\omega=1$ are illustrated in Figure 1.

The scheme starts with an initial guess of the primal $x^{(0)}_{M:N}$ and dual $\lambda^{(0)}_{M:N}$ trajectories with $x^{(0)}_{M}=x_{M}^{*}$ and $\lambda^{(0)}_{N}=\lambda_{N}^{*}$. For each subdomain index $k\in\mathbb{I}_{1:K}$ and iteration counter $t\in\mathbb{Z}_{>0}$, subproblems $\mathcal{P}_{M^{\omega}_{k}:N^{\omega}_{k}}(x^{(t)}_{M_{k}^{\omega}},\lambda^{(t)}_{N_{k}^{\omega}})$ are solved in parallel. Each subproblem requires current primal-dual information from the neighboring subproblems (this information enters in the initial condition and terminal penalty). The solution of subproblem $\mathcal{P}_{M^{\omega}_{k}:N^{\omega}_{k}}(x^{(t)}_{M_{k}^{\omega}},\lambda^{(t)}_{N_{k}^{\omega}})$ yields the local primal-dual trajectories $\hat{z}^{(t+1,k)}_{M^{\omega}_{k}:N^{\omega}_{k}}$. These trajectories are restricted to the non-overlapping subdomains $\mathbb{I}_{M_{k}:N_{k}}$ by using the restriction $z^{(t+1)}_{M_{k}:N_{k}}:=\hat{z}^{(t+1,k)}_{M_{k}:N_{k}}$. The solutions at $\mathbb{I}_{M^{\omega}_{k}:N^{\omega}_{k}}\setminus\mathbb{I}_{M_{k}:N_{k}}$ are discarded. This restriction procedure over $k\in\mathbb{I}_{1:K}$ can be seen as a projection that assembles the entire primal-dual trajectory $z^{(t+1)}_{M:N}$ that is in turn used as the next guess in the algorithm.

The trajectory assembling procedure is sketched in Figure 1. The proposed decomposition scheme is summarized in Algorithm 1 and can be stated compactly as follows:

[TABLE]

III Convergence Results

We now provide a sufficient condition for (1) that guarantees convergence of the decomposition algorithm (2). We begin by making the following existence and uniqueness assumption on solutions of problem (1).

Assumption 1** (Existence and uniqueness of solution)**

There exist sets $X$ and $\Lambda\subseteq\mathbb{R}^{n_{x}}$ such that, for any $M,N\in\mathbb{Z}_{>0}$ with $M<N$ and $(x^{*}_{M},\lambda^{*}_{N})\in X\times\Lambda$, there exists a unique primal-dual solution $z^{*}_{M:N}$ of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ with $z_{i}\in X\times\Lambda$ and for any $i\in\mathbb{I}_{M:N}$.

Remark 2

The solution of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ always satisfies $x_{M}=x^{*}_{M}$ (by (1b)) and $\lambda_{N}=\lambda^{*}_{N}$ (see Remark 1), and thus we can write the solution of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ as $z^{*}_{M:N}=(x^{*}_{M:N},\lambda^{*}_{M:N})$. Because the scheme (2) is mostly concerned with the state and the adjoints (this is the information exchanged between domains), we will not explicitly indicate the controls $u^{*}_{M:N-1}$ and multipliers $\mu^{*}_{M:N-1}$ in the nomenclature.

We now state the following principle of optimality result:

Lemma 1

Assume that problem (1) satisfies Assumption 1 and consider $P,M,N,L\in\mathbb{Z}_{>0}$ with $P\leq M<N\leq L$ and $(x^{*}_{P},\lambda^{*}_{L})\in X\times\Lambda$. Let $z^{*}_{P:L}$ be the solution of $\mathcal{P}_{P:L}(x^{*}_{P},\lambda^{*}_{L})$, then $z^{*}_{M:N}$ is the solution of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$.

Proof:

We define the Lagrangian $\mathcal{L}_{M:N}(\cdot;\cdot)$ of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ as:

[TABLE]

From Assumption 1, the solution of $\mathcal{P}_{1:L}(x^{*}_{1},\lambda^{*}_{L})$ exists and is unique and this implies that the KKT conditions hold at $x^{*}_{M:N}$, $u^{*}_{M:N-1}$, $\lambda^{*}_{M:N}$, $\mu^{*}_{M:N-1}$; that is,

[TABLE]

holds for any $i\in\mathbb{I}_{M:N-1}$, where $f^{*}_{i}:=f(x^{*}_{i},u^{*}_{i})$, $\nabla f^{*}_{i}:=\nabla f(x^{*}_{i},u^{*}_{i})$, $\nabla g^{*}_{i}:=\nabla g(x^{*}_{i},u^{*}_{i})$, $\nabla\ell^{*}_{i}:=\nabla\ell(x^{*}_{i},u^{*}_{i})$.

From (4), we observe that the KKT conditions of $\mathcal{P}_{P:L}(x^{*}_{1},\lambda^{*}_{L})$ with $(x^{*}_{P:L},u^{*}_{P:L-1},\lambda^{*}_{P:L},\mu^{*}_{P:L-1})$ imply that the KKT conditions of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ hold with $(x^{*}_{M:N},u^{*}_{M:N-1},\lambda^{*}_{M:N},\mu^{*}_{M:N-1})$. This consistency is the result of using the duals as a terminal penalty in (1).

By Assumption 1, we have that $(x^{*}_{M},\lambda^{*}_{N})\in X\times\Lambda$. This implies that the solution of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ exists and is unique, and thus the KKT point is the unique solution. Thus, $x^{*}_{M:N-1}$, $u^{*}_{M:N-1}$, $\lambda^{*}_{M:N}$, and $\mu^{*}_{M:N-1}$ form the solutions of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$. ∎

Lemma 1 justifies the structure of (1); specifically, primal and dual information of the neighboring subproblems should be incorporated as initial states and terminal penalties, respectively. Furthermore, it implies that the primal-dual trajectory of the subproblems can be assembled to obtain the optimal primal-dual trajectory of the entire problem if the boundary conditions are set to the optimal values. From this we also see that a receding-horizon control scheme delivers the solution of the long-horizon problem if the terminal cost gradients are obtained from the dual solution of the long-horizon problem. Recently the effect of terminal cost gradient on the solution trajectory of optimal control problems and its connection with dual variables has been investigated in [17], [18]. In particular, it is shown that economic MPC can achieve asymptotic stability without terminal constraints by incorporating the dual of the underlying steady-state problem as a terminal cost gradient. Furthermore, it is shown that LQ economic MPC can be stabilized by adaptively tuning the terminal cost gradient. Such observations align with our observations of Lemma 1. Specifically, the quality of the solutions of OCPs (which in general improves with the horizon length) can also be improved by choosing a suitable terminal cost gradient.

We now state our sufficiency condition for convergence, that we call asymptotic decay of sensitivity (ADS).

Property 1** (Asymptotic decay of sensitivity)**

Assume that (1) satisfies Assumption 1. For given $M$, $N\in\mathbb{Z}_{>0}$ with $M<N$ and $(x^{*}_{M},\lambda^{*}_{N})$, $(\tilde{x}_{M},\tilde{\lambda}_{N})\in X\times\Lambda$, let $z^{*}_{M:N}$ and $\tilde{z}_{M:N}$ be the solutions of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ and $\mathcal{P}_{M:N}(\tilde{x}_{M},\tilde{\lambda}_{N})$, respectively. There exist $\{\varepsilon_{i}\}_{i\in\mathbb{Z}_{\geq 0}}$ with $\varepsilon_{i}\rightarrow 0$ as $i\rightarrow\infty$ such that:

[TABLE]

holds for any $i\in\mathbb{I}_{M:N}$.

Property 1 implies that the solution $z_{i}^{*}=(x^{*}_{i},\lambda^{*}_{i})$ of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ at time $i$ becomes less sensitive to perturbations in the initial state $x^{*}_{M}$ and the terminal cost gradient $\lambda^{*}_{N}$ as the time index $i$ moves away from the boundary.

Remark 3

It is important that the sequence $\{\varepsilon_{i}\}_{i\in\mathbb{Z}_{>0}}$ is a uniform parameter that does not depend on $M$ and $N$ and that does not depend on the choice of the boundary conditions (i.e., $x^{*}_{M}$ and $\lambda^{*}_{N}$). This enables a uniform bound that holds for different subproblems.

Remark 4** (Relation with turnpike properties)**

Property 1 is related but not equivalent to so-called turnpike properties of OCPs (see [19, 20]). Property 1 establishes a relationship between two solution trajectories, while a classical turnpike property compares a single solution trajectory with the steady-state solution. In addition, Property 1 requires asymptotic convergence of the difference between the solutions, while the steady-state turnpike requires a bound on the number of time indexes with $\|x_{i}-x_{\text{s}}\|>\epsilon_{\text{turnpike}}$. An in-depth investigation of the relation between Property 1 and time-varying turnpikes is subject to future work.

We now prove that the ADS property provides a sufficient condition guaranteeing convergence of the proposed decomposition scheme.

Theorem 1** (Convergence)**

Assume that (1) satisfies Assumption 1 and that Property 1 holds. Furthermore, consider scheme (2) with overlap $\omega$ and let $z^{*}_{M:N}$ be the solution of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$. We have that:

[TABLE]

Proof:

From Lemma 1 we have that the solution of $\mathcal{P}_{M:N}(x^{*}_{M},\lambda^{*}_{N})$ on $\mathbb{I}_{M^{\omega}_{k}:N^{\omega}_{k}}$ can be obtained from $\mathcal{P}_{M^{\omega}_{k},N^{\omega}_{k}}(x^{*}_{M^{\omega}_{k}},\lambda^{*}_{N^{\omega}_{k}})$. For each $k\in\mathbb{I}_{1:K}$, applying Property 1 to $\mathcal{P}_{M^{\omega}_{k},N^{\omega}_{k}}(x^{*}_{M^{\omega}_{k}},\lambda^{*}_{N^{\omega}_{k}})$ and $\mathcal{P}_{M^{\omega}_{k},N^{\omega}_{k}}(x^{(t)}_{M^{\omega}_{k}},\lambda^{(t)}_{N^{\omega}_{k}})$ yields

[TABLE]

for any $i\in\mathbb{I}_{M_{k}:N_{k}}$ and $t\in\mathbb{Z}_{>0}$. For $i\in\mathbb{I}_{M_{1}:N_{1}}$, we have that $x^{(t)}_{1}=x^{*}_{1}$ and $i\leq N^{\omega}_{1}-\omega$ yield

[TABLE]

From (7), we have that

[TABLE]

Equation (8) establishes (5). ∎

Theorem 1 indicates that the recursion (2) converges to the solution of the full problem if the size of the overlap $\omega$ is sufficiently large. Furthermore, the convergence rate $2\varepsilon_{\omega}$ converges asymptotically to zero with the size of the overlap $\omega$ (i.e., with a maximal overlap, the iteration converges in one iteration). This reveals a powerful feature of the proposed scheme: one can control the convergence rate by choosing the size of the overlap $\omega$. However, as we increase $\omega$, we also increase the complexity of the subproblem and thus a trade-off exists. Accordingly, the selection of a suitable $\omega$ should consider the convergence rate and the subproblem complexity.

We derive termination criteria for the parallel scheme based on the violation of the KKT conditions. We have that (4d) holds at each iteration. We also have that (4a)-(4c) hold except for $i\in\{N_{1},\cdots,N_{K-1}\}$, which are the boundary indices. For those, the residuals can be evaluated as follows.

[TABLE]

Residuals (9) can be derived by using the fact that (4a)-(4c) hold when $\lambda^{(t+1)}_{i+1}$ and $x^{(t+1)}_{i+1}$ are replaced with $\hat{\lambda}^{(t+1,k)}_{i+1}$ and $\hat{x}^{(t+1,k)}_{i+1}$. Finally, we define the primal-dual residuals:

[TABLE]

and establish the termination criteria:

[TABLE]

IV Linear Quadratic Problems

In this section we show that Property 1 holds for problems with linear and controllable dynamics, convex quadratic objectives, no inequalities, and a single-variable input. Specifically, we make the following assumption. The analysis reveals connections with well-known optimization sensitivity results.

Assumption 2

Consider problem (1) with:

$\ell(x,u):=x^{\top}Qx-f^{\top}x+ru^{2}$* with $Q>0$ and $r>0$* 2. 2.

$f(x,u):=Ax+bu+c$* with $(A,b)$ controllable.* 3. 3.

There are no inequality constraints.

Assumption 2 guarantees that Assumption 1 holds with $X$, $\Lambda=\mathbb{R}^{n_{x}}$. Now we state the main theorem of this section.

Theorem 2

Property 1 holds for OCPs (1) satisfying Assumption 2.

Proof:

Without loss of generality we assume that $M=1$. For any $\mathbb{I}_{1:N}$, $x_{1}^{*}$, $\lambda^{*}_{N}\in\mathbb{R}^{n_{x}}$, the solution of $\mathcal{P}_{1:N}(x^{*}_{1},\lambda^{*}_{N})$ exists and is unique and can be obtained from

[TABLE]

where we define

[TABLE]

Since $(A,b)$ is controllable, the controllability matrix $\mathcal{C}:=\begin{bmatrix}b&Ab&\cdots&A^{n_{x}-1}b\end{bmatrix}$ has full row rank. This implies that $\mathcal{C}u+A^{n_{x}}b=0$ has a unique solution. Using this fact, we can construct $\overline{x}_{1},\cdots,\overline{x}_{n_{x}}\in\mathbb{R}^{n_{x}}$ and $\overline{u}_{1},\cdots,\overline{u}_{n_{x}}\in\mathbb{R}$ with: $\overline{x}_{i+1}=A\overline{x}_{i}+b\overline{u}_{i}$ for any $i\in\mathbb{I}_{0:n_{x}}$, $\overline{x}_{0}=0$, $\overline{u}_{0}=1$, and $\overline{x}_{n_{x}+1}=0$. We construct $Z\in\mathbb{R}^{(Nn_{x}+N-1)\times(N-1)}$ as:

[TABLE]

One can show that $GZ=0$ and that $Z$ has full column rank by using the lower-triangular structure of $Z$. Observe that $G\in\mathbb{R}^{Nn_{x}\times(Nn_{x}+N-1)}$ has full row rank. By the fundamental theorem of linear algebra, the null space of $G$ has dimension $N-1$; consequently, columns of $Z$ span the null space of $G$. Similarly, there exists a unique solution of $\mathcal{C}u+A^{n_{x}}e_{i}=0$ for any $e_{i}$ where $e_{i}$ is the $i$th standard unit vector of $\mathbb{R}^{n_{x}}$. Using this observation, we construct $X_{1},\cdots,X_{n_{x}}\in\mathbb{R}^{n_{x}\times n_{x}}$ and $U_{1},\cdots,U_{n_{x}}\in\mathbb{R}^{1\times n_{x}}$ with $X_{1}=I$, $X_{n_{x}+1}=0$, and $X_{i+1}=AX_{i}+bU_{i}$ for any $i\in\mathbb{I}_{1:n_{x}}$. Now consider

[TABLE]

and observe that $GY=I$ holds. We now apply a null-space projection to (12) using $Z$ and $Y$. This way we obtain the equivalent unconstrained QP

[TABLE]

The solution $p^{*}$ of (15) is given by

[TABLE]

We have that $w^{*}=Zp^{*}+Y\xi^{*}$ and $\lambda^{*}=-Y^{\top}Hw^{*}+Y^{\top}\zeta^{*}$. Using (16), we can write

[TABLE]

where $\overline{H}:=(Z^{\top}HZ)^{-1}$ and $W:=Y^{\top}HZ$. Since $Z$ and $Y$ are independent of the choice of $(x^{*}_{1},\lambda^{*}_{N})$, we obtain the solution of the form (17) with different boundary conditions $(\tilde{x}_{1},\tilde{\lambda}_{N})$. Thus, we have that

[TABLE]

where $\Delta\zeta:=\zeta^{*}-\tilde{\zeta}$ and $\Delta\xi:=\xi^{*}-\tilde{\xi}$. Considering

[TABLE]

and we can see that $\|Z\|_{\infty}\leq L_{Z}$, $\|Z^{\top}\|_{\infty}\leq L_{Z^{\top}}$, $\|Y\|_{\infty}\leq L_{Y}$, $\|Y^{\top}\|_{\infty}\leq L_{Y^{\top}}$, $\|H\|_{\infty}\leq L_{H}$, $\|W\|\leq L_{W}$, and $\|W^{\top}\|\leq L_{W^{\top}}$. Quantities defined in (19) are all uniform (i.e., does not depend on the length of problem and boundary conditions).

The bandwidth $\mathcal{B}(\cdot)$ of a matrix is defined as the smallest integer such that $|i-j|\leq\mathcal{B}(\cdot)$ for any $(\cdot)_{i,j}\neq 0$. Note that the bandwidth of $Z^{\top}HZ$ is not greater than $n_{x}$. The following is a modification of [14, Corollary 1].

Proposition 1

Consider positive definite $\Gamma\in\mathbb{R}^{n\times n}$. Suppose that $\lambda(\Gamma)\in[\lambda_{\min},\lambda_{\max}]$ for some $\lambda_{\min},\lambda_{\max}\in\mathbb{R}_{>0}$. Then the following holds.

[TABLE]

Proof:

The proof is given in [14]. ∎

The following lemma establishes that there exist uniform upper and lower bounds for the eigenvalues of $Z^{\top}HZ$.

Lemma 2

Consider (1), such that Assumption 2 and let $Z$ be from (13). Then there exist uniform parameters $\lambda_{1},\lambda_{2}\in\mathbb{R}_{>0}$ (independent of the choice $N$) such the eigenvalues $\lambda$ of $Z^{\top}HZ$ satisfy $\lambda\in[\lambda_{1},\lambda_{2}]$.

Proof:

The upper bound comes from

[TABLE]

and thus we define $\lambda_{2}:=L^{2}L_{H}$. We can find $\lambda_{1}$ from the lower bound of $p^{\top}(Z^{\top}HZ)p$ for $p\in\mathbb{R}^{N-1}$ with $\|p\|_{2}=1$. We have that

[TABLE]

where $p_{i}=0$ for $i\leq 0$ for convenience and

[TABLE]

Now we show that $\overline{x}_{1},\cdots,\overline{x}_{n_{x}}$ are linearly independent. To establish a contradiction, suppose that $\overline{x}_{1},\cdots,\overline{x}_{n_{x}}$ are linearly dependent; then there exist not all-zero $\alpha_{1},\cdots,\alpha_{n_{x}}$ such that $\alpha_{1}\overline{x}_{1}+\cdots+\alpha_{n_{x}}\overline{x}_{n_{x}}=0$ holds. This implies that

[TABLE]

By rearranging, we obtain

[TABLE]

Since $u_{0}=1$ and $\alpha_{1},\cdots,\alpha_{n_{x}}$ are not all-zero, a non-zero coefficient exists in (21), and thus $b,\cdots,A^{n_{x}-1}b$ are linearly dependent. This contradicts the assumption that $(A,b)$ is controllable. Consequently, $\overline{x}_{1},\cdots,\overline{x}_{n_{x}}$ are linearly independent.

From the linear independence of $\overline{x}_{1},\cdots,\overline{x}_{n_{x}}$, one can show that $\hat{Z}$ has full column rank and thus $\hat{H}$ is positive definite. From (20) one can show that

[TABLE]

Observe that $\hat{H}$ is independent of the length of the problem. Thus, we can set uniform parameter $\lambda_{1}:=\lambda_{\min}(\hat{H})$. ∎

By Proposition 1 and Lemma 2, we have that

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are the uniform upper and lower bound of the eigenvalues of $Z^{\top}HZ$ established in Lemma 2.

Finally, inspecting (18) we observe that $Z$ and $W$ are sparse. In particular, there exists uniform parameter $N_{s}$ such that the following holds for $|j-i|\geq N_{s}$.

[TABLE]

where $\alpha(i):=\mathbb{I}_{(i-1)(n_{x}+1)+1:in_{x}+i-1}$ and $\beta(i):=\mathbb{I}_{(i-1)n_{x}+1:in_{x}}$ (these index sets correspond to the index of $x_{i}$ and $\lambda_{i}$ among the entries of $w$ and $\lambda$, respectively). Using the sparsity structure identified in (24), one can show the following quantities are less than or equal to $L^{2}\|(\overline{H})_{i\pm N_{s},j\pm N_{s}}\|_{\infty}$

[TABLE]

where we use the syntax $i\pm N_{s}:=\mathbb{I}_{i-N_{s}:i+N_{s}}$ and define $(\overline{H})_{i,j}:=0$ if $\{i,j\}\not\subseteq\mathbb{I}_{1:N-1}$ for convenience. From (18) and (25), we have

[TABLE]

By (23) and $\mathbf{1}_{i=N}\leq\rho^{N-i-2N_{s}-1}$, we have that

[TABLE]

for $i\in\mathbb{I}_{M:N}$ and $\rho:=(\lambda_{2}-\lambda_{1})/(\lambda_{2}+\lambda_{1})$. We define

[TABLE]

Note that $\{\varepsilon_{i}\}_{i\in\mathbb{Z}_{\geq 0}}$ is a uniform parameter and (27)-(28) establish Property (1). This concludes the proof. ∎

Theorem 2 implies that the reduced Hessian is positive definite–a key requirement in optimization sensitivity results as the solution must be locally unique and bounded perturbations yield bounded differences in solutions. We thus expect that Theorem 2 can be generalized to LQ OCPs with multiple inputs, time-variant objectives and dynamics, and inequality constraints by exploiting algebraic properties. Here, we focus on a simple setting due to space limitations and to keep the presentation clear. For an even more general setting (nonlinear, inequality-constrained) setting, one can seek to validate the sensitivity property numerically by using simulations. We show how to do this in the next section.

V Numerical Examples

We use a nonlinear OCP to illustrate that Property 1 guarantees convergence of the decomposition algorithm and to highlight that the approach achieves faster solutions than off-the-shelf solvers. We consider a nonlinear economic MPC problem for a chemical reactor [20, 21]. The dynamics of the reactor are given by:

[TABLE]

The control step length is $1$ sec and the differential equations are discretized using an implicit Euler scheme with a step length of $0.25$ sec. We implemented a parallel version of the scheme in Julia. Problems were formulated in the modeling language JuMP [22], and were solved with the nonlinear programming solver IPOPT [23]. The scheme was executed on an Intel Xeon CPU E5-2698 v3 processor running at 2.30GHz.

We use numerical simulations to verify that Property 1 holds. Here, we assess the sensitivity of the primal-dual trajectories against perturbations to the initial state and the terminal cost gradient. The reference problem is formulated with $N=600$ (i.e., $10$ mins) and boundary conditions with $(x^{*}_{1},\lambda^{*}_{N})$. The boundary conditions are perturbed as $(x^{*}_{1},\lambda^{*}_{N})+\delta$, where $\delta$ is sampled from a normal random variable. It is known that the OCP with (29) has a periodic optimal solution when $\rho_{\text{reg}}=0$ and a steady-state solution with $\rho_{\text{reg}}=0.5$ [20]. The solutions of the reference (unperturbed) problem and 30 samples of perturbed problems are shown in Figure 2 ($\rho_{\text{reg}}=0.5$) and Figure 3 ($\rho_{\text{reg}}=0$). When $\rho_{\text{reg}}=0.5$ we see that, for both primal and dual solutions, the distances from the reference trajectories are small in the center of the time domain and grow as we approach the boundary. On the other hand, when $\rho_{\text{reg}}=0$, such convergence is not observed. This indicates that Property 1 holds for $\rho_{\text{reg}}=0.5$ but does not hold for $\rho_{\text{reg}}=0$. This implies that the ADS property strongly depends on the stage cost $\ell$. For high-dimensional systems where it is difficult to graphically assess the ADS property, one may consider assessing the error trajectory based on the norm of deviation from the reference trajectory.

We next verify that the decomposition scheme converges and that the convergence rate improves with the size of the overlap. A problem with $N=4,800$ (i.e., $80$ mins) was decomposed into $8$ subdomains and solved with $\omega=2$, $4$, and $8$ and $\rho_{\text{reg}}=0.5$. Figure 4 (left) shows the evolution of the residuals and confirms that the rate improves with the overlap. The algorithm did not converge when $\rho_{\text{reg}}=0$. These observations reinforce the key role of Property 1.

We also explore the computational efficiency achieved via parallel decomposition. A problem with $N=4,800$ was solved using IPOPT and we compared the solution time against that of the decomposition scheme with an overlap of $\omega=8$ and $\rho_{\text{reg}}=0.5$. The solution found by the decomposition scheme was equal to the solution from IPOPT. Figure 4 (right) clearly illustrates that the solution time decreases as the number of computing cores increases (this also increases the number of subdomains). We can see that, if a sufficient number of cores are used, the proposed scheme becomes faster than IPOPT.

VI Conclusions

We have presented a temporal decomposition scheme for long-horizon OCPs that solves subproblems on overlapping time domains. We have proposed an asymptotic sensitivity decay property that guarantees convergence of the algorithm. This property indicates that the sensitivity of the primal and dual trajectories to perturbations on the boundary conditions (initial state and terminal cost gradient) decays asymptotically as one move away from the boundaries. This property also indicates that the scheme converges if the overlap is sufficiently large and the convergence rate improves with the size of the overlap. We have demonstrated that the solution time of long-horizon OCPs can be improved with the proposed decomposition method.