Na, Anitescu, and Kolar

A Fast Temporal Decomposition Procedure for NLDP

A Fast Temporal Decomposition Procedure for Long-horizon Nonlinear Dynamic
Programming

Sen Na ICSI and Department of Statistics, University of California, Berkeley,
<EMAIL_ADDRESS>Mihai Anitescu Mathematics and Computer Science Division,
Argonne National Laboratory<EMAIL_ADDRESS>Mladen Kolar Booth School of
Business, University of Chicago<EMAIL_ADDRESS>
We propose a fast temporal decomposition procedure for solving long-horizon
nonlinear dynamic programs. The core of the procedure is sequential quadratic
programming (SQP) that utilizes a differentiable exact augmented Lagrangian as
the merit function. Within each SQP iteration, we approximately solve the
Newton system using an overlapping temporal decomposition strategy. We show
that the approximate search direction is still a descent direction of the
augmented Lagrangian, provided the overlap size and penalty parameters are
suitably chosen, which allows us to establish the global convergence.
Moreover, we show that a unit stepsize is accepted locally for the approximate
search direction, and further establish a uniform, local linear convergence
over stages. This local convergence rate matches the rate of the recent
Schwarz scheme [38]. However, the Schwarz scheme has to solve nonlinear
subproblems to optimality in each iteration, while we only perform a single
Newton step instead. Numerical experiments validate our theories and
demonstrate the superiority of our method.

nonlinear dynamic programming; temporal decomposition; sequential quadratic
programming; augmented Lagrangian

##  1 Introduction

We consider nonlinear equality-constrained dynamic programs (NLDPs):

$\displaystyle\min_{{\boldsymbol{x}},{\boldsymbol{u}}}\;\;\;$
$\displaystyle\sum_{k=0}^{N-1}g_{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k})+g_{N}({\boldsymbol{x}}_{N}),$
(1a) s.t.
$\displaystyle{\boldsymbol{x}}_{k+1}=f_{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k}),\quad
k=0,1,\ldots,N-1,$ (1b)
$\displaystyle{\boldsymbol{x}}_{0}=\bar{\boldsymbol{x}}_{0},$ (1c)

where ${\boldsymbol{x}}_{k}\in\mathbb{R}^{n_{x}}$ is the state variable,
${\boldsymbol{u}}_{k}\in\mathbb{R}^{n_{u}}$ is the control variable,
$g_{k}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}$
($g_{N}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}$) is the cost function,
$f_{k}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}$
is the dynamical constraint function, $\bar{\boldsymbol{x}}_{0}$ is the given
initial state, and $N$ is the temporal horizon length. In the control
literature, Problem (1) is also called nonlinear optimal control problem. This
paper focuses on solving (1) with a large $N$.

Large-scale nonlinear problems pose significant computational challenges due
to the nonlinearity of the problems and the large number of variables that
need to be optimized. Numerous solvers have been developed to address the
scalability issue from different aspects. For example, IPOPT [60] is a primal-
dual interior point method that employs the logarithmic barrier function with
the filter line-search step. It enjoys both global and local superlinear
convergence [58, 59]. As another example, Knitro [10] integrates the
advantages of two complementary methods—the interior point method and the
active-set method—in a unified package, and achieves the robust performance.
We refer to [43, 13, 61] for the other scalable high-performance solvers for
nonlinear programs.

The long-horizon NLDPs are of special interest as they appear in a variety of
applications including portfolio management [56], autonomous vehicle [23, 68],
and power planning [15]. The aforementioned centralized solvers, while
exploiting the problem-dependent sparse structures to accelerate the
computation, are not particularly designed for NLDPs. As these solvers do not
take advantage of generic properties of dynamical systems, they are deficient
when directly implemented on Problem (1), especially when we have limited
computing resources. The generic properties of dynamical systems, such as
sensitivity and controllability, play a key role in developing various
efficient algorithms [65, 34, 52, 54, 53]. Further, the centralized solvers
are not flexible enough to be implemented on different types of computing
hardware. Their performance heavily relies on a single high-speed processor,
which sometimes is inaccessible. These limitations motivate us to design a
parallel, DP-oriented procedure for (1).

In contrast to solving NLDPs in real time [16, 18, 17, 67], the dominant class
of distributed offline methods is the decomposition-based methods, where the
full problem is decomposed into multiple subproblems, which are then solved in
a parallel environment [63]. For example, temporal decomposition [1, 30],
Lagrangian dual decomposition [31, 2], and alternating direction method of
multipliers (ADMM) [7, 42] are widely adopted in practice. The decomposition-
based methods are preferable when the problem scale (i.e., $N$ in our case) is
so large that solving (1) on a single processor is computationally
prohibitive, but multiple processors can be used to solve subproblems in
parallel.

This paper contributes to the literature on the overlapping temporal
decomposition (OTD) methods [1, 65, 52, 38]. A TD method decomposes the full
temporal horizon $[0,N]$ into several short horizons
$[0,N]\subseteq\cup_{i}[n_{i},n_{i+1}]$; constructs subproblems
$\\{{\mathcal{P}}^{i}\\}_{i}$ associated with short horizons
$\\{[n_{i},n_{i+1}]\\}_{i}$; solves all subproblems in parallel; and retrieves
the full-horizon solution by composing the subproblem solutions sequentially.
Different from the exclusive decomposition in TD, OTD extends each short
horizon $[n_{i},n_{i+1}]$ by $b$ stages on two ends to encourage information
exchange between two adjacent subproblems. The composition is performed by
using only the exclusive part, $[n_{i},n_{i+1}]$, of each subproblem’s
solution. OTD was first empirically studied in [1] on power planning problems,
and rigorously analyzed for linear-quadratic convex DPs in [65]. In this work,
the authors showed under the controllability and boundedness conditions that
$\max_{k}\|(\hat{{\boldsymbol{x}}}_{k}-{\boldsymbol{x}}^{\star}_{k};\hat{{\boldsymbol{u}}}_{k}-{\boldsymbol{u}}^{\star}_{k})\|\leq
C\rho^{b}$ for some constants $C>0$ and $\rho\in(0,1)$. Here,
$(\hat{{\boldsymbol{x}}}_{k},\hat{{\boldsymbol{u}}}_{k})$ is the OTD output at
the stage $k$ and
$({\boldsymbol{x}}^{\star}_{k},{\boldsymbol{u}}^{\star}_{k})$ is the optimal
solution. For general NLDPs, [52, 38] recently proposed an overlapping Schwarz
scheme as an application of OTD on nonlinear problems. We will review the
Schwarz scheme in section 2 but briefly introduce it here to motivate our
study.

Following the same spirit as OTD, the Schwarz scheme decomposes $[0,N]$ by
$[0,N]\subseteq~{}\cup_{i}~{}[m_{1}^{i},m_{2}^{i}]$, where $m_{1}^{i}=n_{i}-b$
and $m_{2}^{i}=n_{i+1}+b$ are two boundaries of the extended intervals. The
$i$-th subproblem ${\mathcal{P}}^{i}={\mathcal{P}}^{i}({\boldsymbol{d}}_{i})$
is parameterized by some boundary variables
${\boldsymbol{d}}_{i}=({\boldsymbol{d}}_{m_{1}^{i}};{\boldsymbol{d}}_{m_{2}^{i}})$.
In the $\tau$-th iteration, one first specifies the boundary variables
${\boldsymbol{d}}_{i}^{\tau}$ by the solutions of adjacent subproblems (noting
that $m_{1}^{i},m_{2}^{i}\notin[n_{i},n_{i+1}]$); then solves the subproblems
$\\{{\mathcal{P}}^{i}({\boldsymbol{d}}_{i}^{\tau})\\}_{i}$ to optimality in
parallel and updates the boundary variables to
${\boldsymbol{d}}_{i}^{\tau+1}$. The solutions to subproblems are finally
composed to derive a full-horizon solution. The Schwarz scheme is demonstrated
in Figure 1. Clearly, the OTD step corresponds to one iteration of the Schwarz
scheme. [38] empirically showed that the Schwarz scheme significantly improves
the efficiency of ADMM; and may be as efficient as the centralized solver
IPOPT, but provide significant flexibility on different computing
environments. By adjusting the overlap size $b$ of the decomposition, the
scheme adapts to centralized or to decentralized environments. Furthermore,
unlike ADMM whose convergence for nonlinear problems is established for some
setups that may not directly apply to (1) (see [27, 62]), the Schwarz scheme
exhibits a uniform local linear convergence. Under the same conditions of OTD
on linear-quadratic convex DPs in [65], [38] showed
$\max_{k}\|({\boldsymbol{x}}_{k}^{\tau}-{\boldsymbol{x}}^{\star}_{k};{\boldsymbol{u}}_{k}^{\tau}-{\boldsymbol{u}}^{\star}_{k})\|\leq(C\rho^{b})^{\tau}\max_{k}\|({\boldsymbol{x}}_{k}^{0}-{\boldsymbol{x}}^{\star}_{k};{\boldsymbol{u}}_{k}^{0}-{\boldsymbol{u}}^{\star}_{k})\|$,
where $({\boldsymbol{x}}_{k}^{\tau},{\boldsymbol{u}}_{k}^{\tau})$ is the
$\tau$-th iterate of the Schwarz scheme. Despite this promising result, the
Schwarz scheme has two limitations.

Figure 1: Demonstration of the Schwarz scheme. The horizontal line is the full
horizon $[0,N]$. The red vertical lines are knots of the exclusive intervals;
the blue vertical lines are boundaries of the extended intervals. The
subproblem $i$ parameterized by
${\boldsymbol{d}}_{i}=({\boldsymbol{d}}_{m_{1}^{i}},{\boldsymbol{d}}_{m_{2}^{i}})$
is solved, but only the exclusive part of the solution is used in the
composition.

First, only local convergence guarantee is established while global
convergence is still unknown. In fact, as we will explain later, the Schwarz
scheme does not converge globally in general. There is no guarantee that the
composed solution of the subproblems is the solution to the full problem, even
if we correctly specify the boundary variables for each subproblem. Thus, we
arise the question: Q1: How to design an OTD-based procedure that converges
globally? Second, the Schwarz scheme is computationally expensive. In each
iteration, the scheme solves all nonlinear subproblems
$\\{{\mathcal{P}}^{i}({\boldsymbol{d}}_{i})\\}_{i}$ to optimality. However, we
want to know: Q2: Is it necessary to solve subproblems to optimality in each
iteration to enjoy (uniform) local linear convergence?

We answer Q1 and Q2 by designing a fast OTD (FOTD) procedure. Our procedure is
inspired by [63], where the author applied the sequential quadratic
programming (SQP) to solve (1), and a parallel block-banded linear solver to
solve the exact Newton system. By using the unit stepsize, [63] conducted a
local analysis. In this paper, we propose a FOTD procedure to integrate global
and local analyses. Similar to [63], FOTD is also built upon an SQP scheme.
However, it utilizes an exact augmented Lagrangian merit function to
adaptively select the stepsize via line search; and adopts an OTD strategy to
approximately solve the Newton system. While there are many distributed
methods for solving the Newton system exactly or approximately, e.g., [39, 40,
30], we specifically focus on OTD in order to build a relation to the Schwarz
scheme. We prove that FOTD converges globally. A key technical step is to show
that the approximate direction is a descent direction of the augmented
Lagrangian, so that the iterates are improved towards the KKT point. This
technical step justifies the choice of the augmented Lagrangian merit
function. Furthermore, we show that the unit stepsize for the approximate
direction is accepted locally, so that FOTD enjoys a uniform local linear
convergence. Such a linear convergence matches the one of the Schwarz scheme
[38]; however, FOTD requires much fewer computations as it does not solve
nonlinear subproblems to optimality. Our experiments validate the theorems and
demonstrate the superiority of FOTD.

Structure of the paper: In section 2, we present the preliminaries of OTD and
review the Schwarz scheme. In section 3, we introduce our FOTD procedure. In
section 4, we study the approximation error of the Newton system. The global
and local convergence results are established in section 5 and section 6,
respectively. Numerical experiments are presented in section 7 followed by
conclusions in section 8. All the proofs are collected in appendices to make
the main paper compact.

Notation: For an integer $n$, $[n]\coloneqq\\{0,1,\ldots,n\\}$. For two
integers $n,m$, we abuse the interval notation and let $[n,m]$, $(n,m)$,
$[n,m)$, $(n,m]$ be the corresponding index sets. All vectors in the paper are
column vectors. For two vectors ${\boldsymbol{a}}$, ${\boldsymbol{b}}$,
$({\boldsymbol{a}};{\boldsymbol{b}})$ denotes the column vector that stacks
${\boldsymbol{a}}$ and ${\boldsymbol{b}}$ sequentially. For a vector-valued
function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, $\nabla
f\in\mathbb{R}^{n\times m}$ is its Jacobian matrix. We let $\|\cdot\|$ denote
the $\ell_{2}$ norm for vectors and the spectral norm for matrices. For a
symmetric matrix $A$, $\lambda_{\min}(A)$ denotes its smallest eigenvalue. We
let $I$ be the identity matrix and ${\boldsymbol{0}}$ be the zero matrix,
whose dimensions can be inferred from the context. We also reserve the
following notation.
${\boldsymbol{x}}={\boldsymbol{x}}_{0:N}=({\boldsymbol{x}}_{0};\ldots;{\boldsymbol{x}}_{N})$
is the state vector;
${\boldsymbol{u}}={\boldsymbol{u}}_{0:N-1}=({\boldsymbol{u}}_{0};\ldots{\boldsymbol{u}}_{N-1})$
is the control vector;
${\boldsymbol{z}}_{k}=({\boldsymbol{x}}_{k};{\boldsymbol{u}}_{k})$
(${\boldsymbol{z}}_{N}={\boldsymbol{x}}_{N}$) is the state-control pair at
stage $k$, and
${\boldsymbol{z}}={\boldsymbol{z}}_{0:N}=({\boldsymbol{z}}_{0};\ldots;{\boldsymbol{z}}_{N})$.
We let $n_{z}=(N+1)n_{x}+Nn_{u}$ and ${\boldsymbol{z}}\in\mathbb{R}^{n_{z}}$.
We may also express ${\boldsymbol{z}}=({\boldsymbol{x}},{\boldsymbol{u}})$
when explicitly specifying the components ${\boldsymbol{x}},{\boldsymbol{u}}$
of ${\boldsymbol{z}}$.

##  2 Preliminaries

We start by setting up an overlapping temporal decomposition (OTD). Given the
full horizon $[0,N]$, we decompose it into $M$ exclusive intervals with knots
$0=n_{0}<n_{1}<\cdots<n_{M}=N$. For example, we can choose evenly spaced knots
$n_{i}=i\cdot N/M$ (suppose $M$ is a divisor of $N$). Each one of these
intervals is then extended by $b\geq 1$ stages on two ends; that is, the
boundaries of the extended intervals are

$m_{1}^{i}=\max\\{n_{i}-b,\;0\\},\quad\quad
m_{2}^{i}=\min\\{n_{i+1}+b,\;N\\},\quad\quad i=0,1,\ldots,M-1.$ (2)

To simplify the notation, we use $m_{1},m_{2}$ to denote the boundaries of a
general extended interval $i$. We note that two successive extended intervals
overlap on $2b$ stages.

For the extended interval $i$, we define the corresponding subproblem as

$\displaystyle{\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i}):\quad\min_{{\tilde{\boldsymbol{x}}}_{i},{\tilde{\boldsymbol{u}}}_{i}}\;\;\;$
$\displaystyle\sum_{k=m_{1}}^{m_{2}-1}g_{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k})+{\tilde{g}}_{m_{2}}({\boldsymbol{x}}_{m_{2}};{\boldsymbol{d}}_{i,2:4}),$
(3a) s.t.
$\displaystyle{\boldsymbol{x}}_{k+1}=f_{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k}),\quad
k\in[m_{1},m_{2}),$ (3b)
$\displaystyle{\boldsymbol{x}}_{m_{1}}={\boldsymbol{d}}_{i,1},$ (3c)

where ${\tilde{\boldsymbol{x}}}_{i}={\boldsymbol{x}}_{m_{1}:m_{2}}$ and
${\tilde{\boldsymbol{u}}}_{i}={\boldsymbol{u}}_{m_{1}:m_{2}-1}$ are the state
and control variables;
${\boldsymbol{d}}_{i}={\boldsymbol{d}}_{i,1:4}=(\bar{\boldsymbol{x}}_{m_{1}};\bar{\boldsymbol{x}}_{m_{2}};\\\
\bar{\boldsymbol{u}}_{m_{2}};\bar{\boldsymbol{\lambda}}_{m_{2}+1})$ are the
boundary variables;
${\tilde{g}}_{m_{2}}({\boldsymbol{x}}_{m_{2}};{\boldsymbol{d}}_{i,2:4})=g_{N}({\boldsymbol{x}}_{N})$
if $i=M-1$ (i.e., the last subproblem), otherwise for $\mu>0$,

${\tilde{g}}_{m_{2}}({\boldsymbol{x}}_{m_{2}};{\boldsymbol{d}}_{i,2:4})=g_{m_{2}}({\boldsymbol{x}}_{m_{2}},\bar{\boldsymbol{u}}_{m_{2}})-\bar{\boldsymbol{\lambda}}_{m_{2}+1}^{T}f_{m_{2}}({\boldsymbol{x}}_{m_{2}},\bar{\boldsymbol{u}}_{m_{2}})+\frac{\mu}{2}\|{\boldsymbol{x}}_{m_{2}}-\bar{\boldsymbol{x}}_{m_{2}}\|^{2}.$
(4)

For the boundary variables
${\boldsymbol{d}}_{i}=({\boldsymbol{d}}_{i,1};{\boldsymbol{d}}_{i,2:4})$,
${\boldsymbol{d}}_{i,1}=\bar{\boldsymbol{x}}_{m_{1}}$ is the initial state
variable;
${\boldsymbol{d}}_{i,2:4}=~{}(\bar{\boldsymbol{x}}_{m_{2}};\bar{\boldsymbol{u}}_{m_{2}};\\\
\bar{\boldsymbol{\lambda}}_{m_{2}+1})$ are the terminal state, control, and
dual variables. The boundary variables
${\boldsymbol{d}}_{i}=({\boldsymbol{d}}_{i,1};{\boldsymbol{d}}_{i,2:4})$ are
given to the subproblem. In (4), $\mu$ is a uniform quadratic penalty
parameter independent from $i$. In what follows, we denote by
${\tilde{\boldsymbol{z}}}_{i}=({\tilde{\boldsymbol{x}}}_{i},{\tilde{\boldsymbol{u}}}_{i})$
the ordered state-control vector of the subproblem $i$.

The subproblem
${\mathcal{P}}^{i}={\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i})$ in (3) is
essentially the truncation of the full problem (1) onto the interval
$[m_{1},m_{2}]$, except that the initial state is fixed at
${\boldsymbol{d}}_{i,1}=\bar{\boldsymbol{x}}_{m_{1}}$ and the terminal cost on
${\boldsymbol{x}}_{m_{2}}$ is adjusted by
${\tilde{g}}_{m_{2}}({\boldsymbol{x}}_{m_{2}};{\boldsymbol{d}}_{i,2:4})$. The
subproblem is parameterized by ${\boldsymbol{d}}_{i}$; we always specify
${\boldsymbol{d}}_{i}$ based on the previous iterate before solving
${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i})$. The formula (4) was first
proposed by [35] for analyzing the real-time model predictive control schemes.
The middle term depending on $\bar{\boldsymbol{\lambda}}_{m_{2}+1}$ and
$f_{m_{2}}$ reduces the KKT residual brought by the horizon truncation, and
the quadratic penalty term convexifies the subproblem. The benefits of (4)
will be clearer later. We clarify two corner cases: for $i=0$ and $m_{1}=0$,
$\bar{\boldsymbol{x}}_{m_{1}}$ is $\bar{\boldsymbol{x}}_{0}$ from Problem (1);
for $i=M-1$ and $m_{2}=N$, the adjustment on the terminal cost is restored. We
mention that there exist other subproblem formulations in some restrictive
setups: if $g_{m_{2}}$, $f_{m_{2}}$ are separable, then
$\bar{\boldsymbol{u}}_{m_{2}}$ is not needed in (4); if $g_{m_{2}}$ is
quadratic and convex and $f_{m_{2}}$ is affine, then we can let $\mu=0$. All
these formulations ensure that the subproblem
${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i})$ is well-defined (e.g., is
lower bounded) given a well-defined full problem.

Before introducing the Schwarz scheme, we need the following notation. We let
${\boldsymbol{\lambda}}={\boldsymbol{\lambda}}_{0:N}$ be the dual variables of
(1), where ${\boldsymbol{\lambda}}_{0}\in\mathbb{R}^{n_{x}}$ is associated
with the constraint in (1c) and
${\boldsymbol{\lambda}}_{k+1}\in\mathbb{R}^{n_{x}}$ for $k\in[N-1]$ is
associated with the $k$-th constraint in (1b). Similarly, we let
${\tilde{{\boldsymbol{\lambda}}}}_{i}={\boldsymbol{\lambda}}_{m_{1}:m_{2}}$ be
the dual variables of (3). For the state variables
${\tilde{\boldsymbol{x}}}_{i}$ (similar for
${\tilde{\boldsymbol{u}}}_{i},{\tilde{\boldsymbol{z}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i}$)
and $k\in[m_{1},m_{2}]$, ${\tilde{\boldsymbol{x}}}_{i,k}$ denotes the variable
at the stage $k$ in the subproblem $i$. For example, $k=n_{i}$ belongs to both
subproblem $i-1$ and subproblem $i$; thus
${\tilde{\boldsymbol{x}}}_{i-1,n_{i}}$ and
${\tilde{\boldsymbol{x}}}_{i,n_{i}}$ refer to different variables at the same
stage. The primal-dual solution of
${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i})$ is denoted by
$({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}))$.
We also define composition and decomposition operators in the next definition.

######  Definition 2.1 (composition and decomposition)

For the subproblem variables
$\\{({\tilde{\boldsymbol{z}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i})\\}_{i=0}^{M-1}$,
we define a composition operator $\mathcal{C}$ as
$\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i})\\}_{i})=({\boldsymbol{z}},{\boldsymbol{\lambda}})$,
where
$({\boldsymbol{z}}_{k},{\boldsymbol{\lambda}}_{k})=({\tilde{\boldsymbol{z}}}_{i,k},{\tilde{{\boldsymbol{\lambda}}}}_{i,k})$
if $k\in[n_{i},n_{i+1})$ for $i\in[M-1]$, and
${\boldsymbol{z}}_{N}={\tilde{\boldsymbol{x}}}_{M-1,N}$ and
${\boldsymbol{\lambda}}_{N}={\tilde{{\boldsymbol{\lambda}}}}_{M-1,N}$.
Conversely, for the full-horizon variable
$({\boldsymbol{z}},{\boldsymbol{\lambda}})$, we define a decomposition
operator $\mathcal{D}$ as
$\mathcal{D}({\boldsymbol{z}},{\boldsymbol{\lambda}})=\\{\mathcal{D}_{i}({\boldsymbol{z}},{\boldsymbol{\lambda}})\\}_{i=0}^{M-1}$,
where
$\mathcal{D}_{i}({\boldsymbol{z}},{\boldsymbol{\lambda}})=({\tilde{\boldsymbol{z}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i})$
with
${\tilde{\boldsymbol{z}}}_{i}=({\tilde{\boldsymbol{x}}}_{i},{\tilde{\boldsymbol{u}}}_{i})=({\boldsymbol{x}}_{m_{1}:m_{2}},{\boldsymbol{u}}_{m_{1}:m_{2}-1})$
and
${\tilde{{\boldsymbol{\lambda}}}}_{i}={\boldsymbol{\lambda}}_{m_{1}:m_{2}}$.

From Definition 2.1, we see that the variables on the overlapping stages are
discarded during the composition. That is,
${\tilde{\boldsymbol{x}}}_{i}={\tilde{\boldsymbol{x}}}_{i,m_{1}:m_{2}}$
(similar for
${\tilde{\boldsymbol{u}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i}$) only
contributes ${\tilde{\boldsymbol{x}}}_{i,n_{i}:n_{i+1}-1}$ to the full vector.

Schwarz scheme. Given the $\tau$-th iterate
$({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$, we specify the
boundary variables by
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{x}}_{m_{1}}^{\tau};\\\
{\boldsymbol{x}}_{m_{2}}^{\tau};{\boldsymbol{u}}_{m_{2}}^{\tau};{\boldsymbol{\lambda}}_{m_{2}+1}^{\tau})$;
solve the subproblems
$\\{{\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i}^{\tau})\\}_{i}$ in parallel
to optimality; and obtain the solutions
$\\{({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i}$.
Then,
$({\boldsymbol{z}}^{\tau+1},{\boldsymbol{\lambda}}^{\tau+1})=\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i})$.
The Schwarz scheme is displayed in Algorithm 1. The convergence of the scheme
is achieved by iteratively updating ${\boldsymbol{d}}_{i}$. One expects that,
as $\tau$ increases, ${\boldsymbol{d}}_{i}^{\tau}$ becomes more precise so
that the solution of ${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i}^{\tau})$
approaches to the truncated full solution. Thus, the composed solution of
subproblems can recover the full solution.

The following result characterizes the relation between the optimality
conditions of the subproblems (3) and the full problem (1).

Algorithm 1 Overlapping Schwarz Decomposition Procedure

1:Input: initial iterate $({\boldsymbol{z}}^{0},{\boldsymbol{\lambda}}^{0})$
with ${\boldsymbol{x}}_{0}^{0}=\bar{\boldsymbol{x}}_{0}$, a scalar $\mu>0$;

2:for $\tau=0,1,2,\ldots$ do

3: for $i=0,1,\ldots,M-1$ (in parallel) do

4: Let
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{x}}_{m_{1}}^{\tau};{\boldsymbol{x}}_{m_{2}}^{\tau};{\boldsymbol{u}}_{m_{2}}^{\tau};{\boldsymbol{\lambda}}_{m_{2}+1}^{\tau})$
(note that ${\boldsymbol{d}}_{M-1}^{\tau}={\boldsymbol{x}}_{m_{1}}^{\tau}$);

5: Solve ${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i}^{\tau})$ to optimality
and obtain the solution
$({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))$;

6: end for

7: Let
$({\boldsymbol{z}}^{\tau+1},{\boldsymbol{\lambda}}^{\tau+1})=\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i})$;

8:end for

######  Theorem 2.2 (relation of KKT conditions)

For any scalar $\mu$, we have two cases:

(i) Suppose $({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$ is a
KKT point of (1), then
$\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$,
$\forall i\in[M-1]$, is a KKT point of
${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i}^{\star})$ where
${\boldsymbol{d}}_{i}^{\star}=({\boldsymbol{x}}^{\star}_{m_{1}};{\boldsymbol{x}}^{\star}_{m_{2}};{\boldsymbol{u}}^{\star}_{m_{2}};{\boldsymbol{\lambda}}^{\star}_{m_{2}+1})$;

(ii) Suppose
$({\tilde{\boldsymbol{z}}}_{i}^{\star},{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star})=({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}))$,
$\forall i\in[M-1]$, is a KKT point of
${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i})$ with any boundary variables
${\boldsymbol{d}}_{i}$ satisfying
${\boldsymbol{d}}_{0,1}=\bar{\boldsymbol{x}}_{0}$, then
$\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i}^{\star},{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star})\\}_{i})$
is a KKT point of (1) if and only if the solutions of any two successive
subproblems are compatible at the common boundaries. That is, for any
$i\in[1,M-1]$ and knot $k=n_{i}$, we have

${\tilde{\boldsymbol{x}}}_{i,k}^{\star}={\tilde{\boldsymbol{x}}}_{i-1,k}^{\star},\quad
A_{k-1}^{T}({\tilde{{\boldsymbol{\lambda}}}}_{i,k}^{\star}-{\tilde{{\boldsymbol{\lambda}}}}_{i-1,k}^{\star})={\boldsymbol{0}},\quad
B_{k-1}^{T}({\tilde{{\boldsymbol{\lambda}}}}_{i,k}^{\star}-{\tilde{{\boldsymbol{\lambda}}}}_{i-1,k}^{\star})={\boldsymbol{0}},$
(5)

where
$A_{k-1}=\nabla_{{\boldsymbol{x}}_{k-1}}^{T}f_{k-1}\in\mathbb{R}^{n_{x}\times
n_{x}}$,
$B_{k-1}=\nabla_{{\boldsymbol{u}}_{k-1}}^{T}f_{k-1}\in\mathbb{R}^{n_{x}\times
n_{u}}$ are the Jacobian matrices of $f_{k-1}$ evaluated at
${\tilde{\boldsymbol{z}}}_{i-1,k-1}^{\star}$.

###### Proof 2.3

See Appendix 9.1.

From Theorem 2.2(i), we can see that the truncated KKT point
$\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$ is
also a KKT point of the subproblem, provided the boundary variables are
correctly specified (i.e.,
${\boldsymbol{d}}_{i}={\boldsymbol{d}}_{i}^{\star}$). However, Theorem 2.2(ii)
provides a negative view: even if one obtains a KKT point for each subproblem
with any boundary variables, composing the solutions together does not
necessarily result in a KKT point of the full problem. Only if the successive
KKT points are compatible on knots $\\{n_{i}\\}_{i=1}^{M-1}$ can we guarantee
to have a full-horizon KKT point. The subtlety lies in the fact that the
solution at stage $n_{i}$ is from subproblem $i$, while the solution at stage
$n_{i}-1$ is from subproblem $i-1$. Thus, (5) is needed to link solutions from
two successive subproblems. By Theorem 2.2(ii), we know the Schwarz scheme,
which simply composes subproblem solutions, may not converge globally in
general.

However, if $({\boldsymbol{z}}^{0},{\boldsymbol{\lambda}}^{0})$ is
sufficiently close to
$({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$, [38] showed that
for $i\in[M-1]$, ${\mathcal{P}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ has a unique
solution
$\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$ in
a neighborhood of ${\boldsymbol{d}}_{i}^{\star}$. Thus, one expects
$({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\rightarrow\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$
as $\tau\rightarrow\infty$, and
$\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i})\rightarrow\mathcal{C}(\\{\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})\\}_{i})=({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$.
Then, the local convergence of the Schwarz scheme is ensured. Specifically,
[38] proved for some $C>0$ and $\rho\in(0,1)$ that

$\max_{k}\|({\boldsymbol{z}}_{k}^{\tau}-{\boldsymbol{z}}^{\star}_{k};{\boldsymbol{\lambda}}_{k}^{\tau}-{\boldsymbol{\lambda}}^{\star}_{k})\|\leq(C\rho^{b})^{\tau}\cdot\max_{k}\|({\boldsymbol{z}}_{k}^{0}-{\boldsymbol{z}}^{\star}_{k};{\boldsymbol{\lambda}}_{k}^{0}-{\boldsymbol{\lambda}}^{\star}_{k})\|.$
(6)

In addition to lacking global convergence, the Schwarz scheme is
computationally intensive, since Line 5 of Algorithm 1 requires finding
optimal solutions to nonlinear subproblems. In the next section, we relax this
computational requirement by using an SQP scheme, and finally show in section
6 that (6) holds even when Line 5 is substituted by one Newton step (cf.
Theorem 6.7).

##  3 Embedding OTD into SQP

We note that the lack of global convergence of the Schwarz scheme is due to
the lack of a coordinator, which can monitor the convergence progress towards
the full solution. This motivates us to solve (1) under an SQP framework. For
each SQP iteration, we apply OTD to approximately solve the Newton system,
which corresponds to a linear-quadratic DP problem. Although other distributed
methods are also applicable, we are particularly interested in OTD since it
reveals a nice relation to the Schwarz scheme. By embedding OTD into SQP, we
are able to establish global convergence, which resolves one of the
limitations of the Schwarz scheme.

We write Problem (1) in a compact form by

$\min_{{\boldsymbol{z}}}\;\;g({\boldsymbol{z}}),\;\;\quad\text{s.t.}\;\;f({\boldsymbol{z}})={\boldsymbol{0}},$

where

$\displaystyle g({\boldsymbol{z}})$
$\displaystyle=\sum_{k=0}^{N}g_{k}({\boldsymbol{z}}_{k})=\sum_{k=0}^{N-1}g_{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k})+g_{N}({\boldsymbol{x}}_{N}),$
(7) $\displaystyle f({\boldsymbol{z}})$
$\displaystyle=({\boldsymbol{x}}_{0}-\bar{\boldsymbol{x}}_{0};{\boldsymbol{x}}_{1}-f_{0}({\boldsymbol{z}}_{0});\ldots;{\boldsymbol{x}}_{N}-f_{N-1}({\boldsymbol{z}}_{N-1}))$
$\displaystyle=({\boldsymbol{x}}_{0}-\bar{\boldsymbol{x}}_{0};{\boldsymbol{x}}_{1}-f_{0}({\boldsymbol{x}}_{0},{\boldsymbol{u}}_{0});\ldots;{\boldsymbol{x}}_{N}-f_{N-1}({\boldsymbol{x}}_{N-1},{\boldsymbol{u}}_{N-1})).$

The Lagrangian function is
$\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})=g({\boldsymbol{z}})+{\boldsymbol{\lambda}}^{T}f({\boldsymbol{z}})$
and the KKT conditions are
$\nabla_{{\boldsymbol{z}}}\mathcal{L}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})={\boldsymbol{0}}$,
$\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})={\boldsymbol{0}}$.
The SQP scheme applies Newton’s method on the KKT system. In particular, given
the $\tau$-th iterate
$({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$, the Newton
direction
$(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})$ is
obtained by solving

$\begin{pmatrix}{\hat{H}}^{\tau}&(G^{\tau})^{T}\\\
G^{\tau}&{\boldsymbol{0}}\end{pmatrix}\begin{pmatrix}\Delta{\boldsymbol{z}}^{\tau}\\\
\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}=-\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}^{\tau}\end{pmatrix},$ (8)

where
$\nabla_{{\boldsymbol{z}}}\mathcal{L}^{\tau}=\nabla_{{\boldsymbol{z}}}\mathcal{L}({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$,
$\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}^{\tau}=\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$,
$G^{\tau}=\nabla_{{\boldsymbol{z}}}^{T}f({\boldsymbol{z}}^{\tau})$, and
${\hat{H}}^{\tau}$ is a modification of the Hessian
$H^{\tau}=\nabla^{2}_{{\boldsymbol{z}}}\mathcal{L}({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$
that preserves the block diagonal structure of $H^{\tau}$. The goal of the
modification is to let ${\hat{H}}^{\tau}$ be positive definite in the null
space $\\{{\boldsymbol{z}}:G^{\tau}{\boldsymbol{z}}={\boldsymbol{0}}\\}$, if
$H^{\tau}$ is not. For example, a simple structure-preserving modification is
the Levenberg-style modification [20]: ${\hat{H}}^{\tau}=H^{\tau}+\gamma I$
for a suitably large $\gamma>0$. Other Hessian modification methods are
referred to in [41]. The explicit formula of the system (8) is displayed in
Problem (14).

Given the Newton direction
$(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})$ from
(8), the SQP iterate is updated as

$\begin{pmatrix}{\boldsymbol{z}}^{\tau+1}\\\
{\boldsymbol{\lambda}}^{\tau+1}\end{pmatrix}=\begin{pmatrix}{\boldsymbol{z}}^{\tau}\\\
{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}+\alpha_{\tau}\begin{pmatrix}\Delta{\boldsymbol{z}}^{\tau}\\\
\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix},$ (9)

with the stepsize $\alpha_{\tau}$ being selected by passing a line search
condition based on a merit function. We employ the following differentiable
exact augmented Lagrangian as the merit function

$\mathcal{L}_{\eta}({\boldsymbol{z}},{\boldsymbol{\lambda}})=\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})+\frac{\eta_{1}}{2}\|\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})\|^{2}+\frac{\eta_{2}}{2}\|\nabla_{{\boldsymbol{z}}}\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})\|^{2},$
(10)

where $\eta=(\eta_{1},\eta_{2})$ are the penalty parameters. The first penalty
biases the feasibility error, while the second penalty biases the optimality
error. The function (10) is called exact augmented Lagrangian, since one can
show that the solution of the unconstrained problem
$\min_{{\boldsymbol{z}},{\boldsymbol{\lambda}}}\mathcal{L}_{\eta}({\boldsymbol{z}},{\boldsymbol{\lambda}})$
is also the solution of (1), provided $\eta_{1}$ is large enough and
$\eta_{2}$ is small enough [3, Proposition 4.15]. We refer to [45, 46, 69, 36,
37] for more studies of (10) on constrained optimization problems. We should
mention that there are many penalty functions that can be used as a merit
function; however, (10) is particularly suitable and important for our
analysis. First, recalling the subproblem terminal cost (4), we need an
approximation of the terminal dual variable, which means that the merit
function has to endure a dual perturbation. Such a requirement rules out the
penalty functions that depend only on the primal variables ${\boldsymbol{z}}$.
Second, for the SQP schemes, the differentiable merit functions such as (10)
can overcome the Maratos effect and locally accept a unit stepsize, which is
critical to have a fast local convergence rate. In contrast, the non-smooth
merit functions suffer from the Maratos effect and require non-trivial
modifications (e.g., the second-order correction) of the SQP schemes to
achieve a fast local rate [41, Chapter 15.5]. We will clearly see the benefits
of the merit function (10) later. In particular, see the discussion below
Theorem 5.3 and see Theorem 6.1 as well.

With (10), the stepsize $\alpha_{\tau}$ is selected to make the following
Armijo condition hold

$\mathcal{L}_{\eta}^{\tau+1}\leq\mathcal{L}_{\eta}^{\tau}+\beta\alpha_{\tau}\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}_{\eta}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}_{\eta}^{\tau}\end{pmatrix}^{T}\begin{pmatrix}\Delta{\boldsymbol{z}}^{\tau}\\\
\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix},$ (11)

where $\beta\in(0,1/2)$ is a prespecified parameter,
$\mathcal{L}_{\eta}^{\tau}=\mathcal{L}_{\eta}({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$
(similar for $\mathcal{L}_{\eta}^{\tau+1}$,
$\nabla\mathcal{L}_{\eta}^{\tau}$), and

$\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}_{\eta}({\boldsymbol{z}},{\boldsymbol{\lambda}})\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}_{\eta}({\boldsymbol{z}},{\boldsymbol{\lambda}})\end{pmatrix}=\begin{pmatrix}I+\eta_{2}H({\boldsymbol{z}},{\boldsymbol{\lambda}})&\eta_{1}G^{T}({\boldsymbol{z}})\\\
\eta_{2}G({\boldsymbol{z}})&I\end{pmatrix}\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})\end{pmatrix}.$
(12)

Among the steps in (8), (9), and (11), solving the Newton system in (8) is
often the most computationally expensive step. Thus, we apply the
decomposition method, OTD, to solve (8) approximately. We obtain an
approximate direction
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$,
and then use it in the steps of (9) and (11). We introduce some extra
notation. We let
$H({\boldsymbol{z}},{\boldsymbol{\lambda}})=\nabla_{{\boldsymbol{z}}}^{2}\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})=\operatorname{diag}(H_{0},\ldots,H_{N})$
be the Hessian of $\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})$ with
respect to ${\boldsymbol{z}}$, where

$H_{k}({\boldsymbol{z}}_{k},{\boldsymbol{\lambda}}_{k+1})=\begin{pmatrix}Q_{k}&S_{k}^{T}\\\
S_{k}&R_{k}\end{pmatrix}=\begin{pmatrix}\nabla_{{\boldsymbol{x}}_{k}}^{2}\mathcal{L}&\nabla_{{\boldsymbol{x}}_{k}{\boldsymbol{u}}_{k}}\mathcal{L}\\\
\nabla_{{\boldsymbol{u}}_{k}{\boldsymbol{x}}_{k}}\mathcal{L}&\nabla_{{\boldsymbol{u}}_{k}}^{2}\mathcal{L}\end{pmatrix},\quad\forall
k\in[N-1],\quad
H_{N}({\boldsymbol{z}}_{N})=\nabla_{{\boldsymbol{x}}_{N}}^{2}\mathcal{L}.$
(13)

Note that $H$ does not depend on ${\boldsymbol{\lambda}}_{0}$. We also let
$A_{k}({\boldsymbol{z}}_{k})=\nabla_{{\boldsymbol{x}}_{k}}^{T}f_{k}({\boldsymbol{z}}_{k})$
and
$B_{k}({\boldsymbol{z}}_{k})=\nabla_{{\boldsymbol{u}}_{k}}^{T}f_{k}({\boldsymbol{z}}_{k})$.
To simplify the notation, we suppress the evaluation points in
$\\{H_{k},A_{k},B_{k}\\}$ and suppress the iteration index $\tau$ to refer to
a general $\tau$-th iteration of (8). With the same block-diagonal structure
as $H$, we have
${\hat{H}}=\operatorname{diag}({\hat{H}}_{0},\ldots,{\hat{H}}_{N})$, and use
${\hat{Q}}_{k},{\hat{S}}_{k},{\hat{R}}_{k}$ to denote the components of
${\hat{H}}_{k}$.

Solving (8) is equivalent to solving the following linear-quadratic DP
problem:

$\displaystyle\min_{{\boldsymbol{p}},{\boldsymbol{q}}}\;\;\;$
$\displaystyle\sum_{k=0}^{N-1}\left\\{\frac{1}{2}\begin{pmatrix}{\boldsymbol{p}}_{k}\\\
{\boldsymbol{q}}_{k}\end{pmatrix}^{T}\begin{pmatrix}{\hat{Q}}_{k}&{\hat{S}}_{k}^{T}\\\
{\hat{S}}_{k}&{\hat{R}}_{k}\end{pmatrix}\begin{pmatrix}{\boldsymbol{p}}_{k}\\\
{\boldsymbol{q}}_{k}\end{pmatrix}+\begin{pmatrix}\nabla_{{\boldsymbol{x}}_{k}}\mathcal{L}\\\
\nabla_{{\boldsymbol{u}}_{k}}\mathcal{L}\end{pmatrix}^{T}\begin{pmatrix}{\boldsymbol{p}}_{k}\\\
{\boldsymbol{q}}_{k}\end{pmatrix}\right\\}+\frac{1}{2}{\boldsymbol{p}}_{N}^{T}{\hat{Q}}_{N}{\boldsymbol{p}}_{N}+\nabla_{{\boldsymbol{x}}_{N}}^{T}\mathcal{L}{\boldsymbol{p}}_{N},$
(14a) s.t.
$\displaystyle{\boldsymbol{p}}_{k+1}=A_{k}{\boldsymbol{p}}_{k}+B_{k}{\boldsymbol{q}}_{k}-\nabla_{{\boldsymbol{\lambda}}_{k+1}}\mathcal{L},\quad
k\in[N-1],$ (14b)
$\displaystyle{\boldsymbol{p}}_{0}=-\nabla_{{\boldsymbol{\lambda}}_{0}}\mathcal{L}.$
(14c)

We let ${\boldsymbol{\zeta}}={\boldsymbol{\zeta}}_{0:N}$ be the dual variables
and
${\boldsymbol{w}}={\boldsymbol{w}}_{0:N}=({\boldsymbol{p}},{\boldsymbol{q}})$
be the primal variables of Problem (14). By the recursive constraints in
(14b)-(14c), we know that the Jacobian $G$ (in any iteration) has a
“staircase” structure: for all $k\in[N]$, the $(k,2k-1)$-block matrix is an
identity matrix. Thus, $G$ has full row rank. Suppose ${\hat{H}}$ is positive
definite in $\\{{\boldsymbol{z}}:G{\boldsymbol{z}}={\boldsymbol{0}}\\}$, then
the KKT matrix in (8) is nonsingular; (14) has a unique global solution
$({\boldsymbol{w}}^{\star},{\boldsymbol{\zeta}}^{\star})$; and
$(\Delta{\boldsymbol{z}},\Delta{\boldsymbol{\lambda}})=({\boldsymbol{w}}^{\star},{\boldsymbol{\zeta}}^{\star})$
[41, Theorem 16.2]. Note that the linear terms in (14a) are given by the
gradients of the Lagrangian. With this formulation, our dual solution
${\boldsymbol{\zeta}}^{\star}$ of (14) is the dual direction
$\Delta{\boldsymbol{\lambda}}$. When the linear terms in (14a) are given by
the gradients of the objective $g$, the dual solution
${\boldsymbol{\zeta}}^{\star}$ would be the dual iterate
${\boldsymbol{\lambda}}+\Delta{\boldsymbol{\lambda}}$.

We now apply OTD on (14). By an analogy to (3), the subproblem $i$,
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$, is defined as

$\displaystyle\min_{{\tilde{\boldsymbol{p}}}_{i},{\tilde{\boldsymbol{q}}}_{i}}\;\;\;$
$\displaystyle\sum_{k=m_{1}}^{m_{2}-1}\left\\{\frac{1}{2}\begin{pmatrix}{\boldsymbol{p}}_{k}\\\
{\boldsymbol{q}}_{k}\end{pmatrix}^{T}\begin{pmatrix}{\hat{Q}}_{k}&{\hat{S}}_{k}^{T}\\\
{\hat{S}}_{k}&{\hat{R}}_{k}\end{pmatrix}\begin{pmatrix}{\boldsymbol{p}}_{k}\\\
{\boldsymbol{q}}_{k}\end{pmatrix}+\begin{pmatrix}\nabla_{{\boldsymbol{x}}_{k}}\mathcal{L}\\\
\nabla_{{\boldsymbol{u}}_{k}}\mathcal{L}\end{pmatrix}^{T}\begin{pmatrix}{\boldsymbol{p}}_{k}\\\
{\boldsymbol{q}}_{k}\end{pmatrix}\right\\}$ (15a)
$\displaystyle\quad+\frac{1}{2}\begin{pmatrix}{\boldsymbol{p}}_{m_{2}}\\\
\bar{\boldsymbol{q}}_{m_{2}}\end{pmatrix}^{T}\begin{pmatrix}{\hat{Q}}_{m_{2}}&{\hat{S}}_{m_{2}}^{T}\\\
{\hat{S}}_{m_{2}}&{\hat{R}}_{m_{2}}\end{pmatrix}\begin{pmatrix}{\boldsymbol{p}}_{m_{2}}\\\
\bar{\boldsymbol{q}}_{m_{2}}\end{pmatrix}+{\boldsymbol{p}}_{m_{2}}^{T}(\nabla_{{\boldsymbol{x}}_{m_{2}}}\mathcal{L}-A_{m_{2}}^{T}\bar{\boldsymbol{\zeta}}_{m_{2}+1})+\frac{\mu}{2}\|{\boldsymbol{p}}_{m_{2}}-\bar{\boldsymbol{p}}_{m_{2}}\|^{2},$
s.t.
$\displaystyle{\boldsymbol{p}}_{k+1}=A_{k}{\boldsymbol{p}}_{k}+B_{k}{\boldsymbol{q}}_{k}-\nabla_{{\boldsymbol{\lambda}}_{k+1}}\mathcal{L},\quad
k\in[m_{1},m_{2}),$ (15b)
$\displaystyle{\boldsymbol{p}}_{m_{1}}=\bar{\boldsymbol{p}}_{m_{1}}.$ (15c)

With slight abuse of notation,
${\boldsymbol{d}}_{i}={\boldsymbol{d}}_{i,1:4}=(\bar{\boldsymbol{p}}_{m_{1}};\bar{\boldsymbol{p}}_{m_{2}};\bar{\boldsymbol{q}}_{m_{2}};\bar{\boldsymbol{\zeta}}_{m_{2}+1})$
are the boundary variables;
${\tilde{\boldsymbol{w}}}_{i}=({\tilde{\boldsymbol{p}}}_{i},{\tilde{\boldsymbol{q}}}_{i})$,
${\tilde{\boldsymbol{\zeta}}}_{i}={\boldsymbol{\zeta}}_{m_{1}:m_{2}}$ are the
primal, dual variables of
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$; and
$({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}))$
is the solution.

FOTD scheme. We set the stage to present the FOTD procedure. FOTD consists of
three steps: given the $\tau$-th iterate
$({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$,

Step 1: Compute Hessian ${\hat{H}}^{\tau}$, Jacobian $G^{\tau}$, and KKT
residual vector $\nabla\mathcal{L}^{\tau}$.

Step 2: Solve subproblems
$\\{{{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i}^{\tau})\\}_{i}$
with
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}})$
and obtain solutions
$\\{({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i}$.
Then,
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})=\mathcal{C}(\\{({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i})$.

Step 3: Update the iterate by (9) with $\alpha_{\tau}$ being selected by a
line search based on (11).

We summarize FOTD in Algorithm 2 and present several remarks.

######  Remark 3.1 (necessity of the quadratic penalty)

Even if the full problem (14) has a unique global solution, this does not
necessarily imply that the subproblem (15) has a unique solution. Consider the
following example. Let $N=2$, $n_{x}=n_{u}=1$,
${\hat{Q}}_{0}={\hat{R}}_{0}=1$, ${\hat{Q}}_{1}=-{\hat{R}}_{1}=-2$,
${\hat{Q}}_{2}=3$, ${\hat{S}}_{0}={\hat{S}}_{1}=0$,
$A_{0}=A_{1}=B_{0}=B_{1}=1$, and all linear terms in (14a) are zero. Then, the
full problem: $\min_{{\boldsymbol{p}},{\boldsymbol{q}}}$
$p_{0}^{2}+q_{0}^{2}-2p_{1}^{2}+2q_{1}^{2}+3p_{2}^{2}$, s.t. $p_{0}=0$,
$p_{1}=p_{0}+q_{0}$, $p_{2}=p_{1}+q_{1}$, has a unique global solution
$({\boldsymbol{p}}^{\star},{\boldsymbol{q}}^{\star})=({\boldsymbol{0}},{\boldsymbol{0}})$.
However, when we truncate onto $[0,1]$, the subproblem has a quadratic
objective with the square matrix $\operatorname{diag}(1,1,-2+\mu)$ and
constraints $p_{0}=0$, $p_{1}=p_{0}+q_{0}$. Thus, by plugging the constraints
into the objective, we can easily obtain that, when $\mu<1$, the subproblem is
unbounded below; when $\mu=1$, the subproblem has infinitely many solutions;
when $\mu>1$, the subproblem has a unique global solution. Thus, to ensure
that the subproblem has a unique global solution, a quadratic penalty with
large enough $\mu$ on the terminal stage is necessary (cf. (15a)).

Algorithm 2 A Fast Overlapping Temporal Decomposition Procedure

1:Input: initial iterate $({\boldsymbol{z}}^{0},{\boldsymbol{\lambda}}^{0})$
with ${\boldsymbol{x}}_{0}^{0}=\bar{\boldsymbol{x}}_{0}$, scalars
$\mu,\eta_{1},\eta_{2}>0$, $\beta\in(0,1/2)$;

2:for $\tau=0,1,2,\ldots$ do

3: Compute ${\hat{H}}^{\tau}$, $G^{\tau}$, and $\nabla\mathcal{L}^{\tau}$;

4: for $i=0,1,\ldots,M-1$ (in parallel) do

5: Let
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}})$
(note that ${\boldsymbol{d}}_{M-1}^{\tau}={\boldsymbol{0}}$);

6: Solve ${{\mathcal{L}}{\mathcal{P}}}^{i}_{\mu}({\boldsymbol{d}}_{i}^{\tau})$
to optimality and obtain the solution
$({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))$;

7: end for

8: Let
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})=\mathcal{C}(\\{({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i})$;

9: Select $\alpha_{\tau}$ by a line search based on (11) and update the
iterate by (9);

10:end for

###### Remark 3.2

We can express the optimality conditions of
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ using a Newton
system like (8). We observe that the KKT matrix depends on
$\\{A_{k},B_{k},{\hat{H}}_{k}\\}_{k=m_{1}}^{m_{2}-1}\cup\\{{\hat{Q}}_{m_{2}}+\mu
I\\}$, but not on ${\boldsymbol{d}}_{i}$. Thus, the uniqueness of the solution
of ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ is
independent from ${\boldsymbol{d}}_{i}$. By Theorem 2.2(i), if
${\boldsymbol{d}}_{i}^{\star}=(\Delta{\boldsymbol{x}}_{m_{1}};\Delta{\boldsymbol{z}}_{m_{2}};\Delta{\boldsymbol{\lambda}}_{m_{2}+1})$,
then
$({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\star}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\star}))=(\Delta{\boldsymbol{x}}_{m_{1}:m_{2}},\Delta{\boldsymbol{u}}_{m_{1}:m_{2}-1},\Delta{\boldsymbol{\lambda}}_{m_{1}:m_{2}})$.
However, since there is no good guess for a search direction, we always set
${\boldsymbol{d}}_{i}=({\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}})$
in Algorithm 2 (Line 5).

###### Remark 3.3

Comparing Line 6 of Algorithm 2 with Line 5 of Algorithm 1, FOTD also solves
subproblems $\\{{{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}\\}_{i}$ to optimality.
However, ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$ is a linear-quadratic DP,
which can be solved efficiently while solving ${\mathcal{P}}_{\mu}^{i}$ is
more expensive. The problem ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$ can be
regarded as a linearization of the problem ${\mathcal{P}}_{\mu}^{i}$; thus
solving ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$ corresponds to computing one
Newton step of ${\mathcal{P}}_{\mu}^{i}$.

###### Remark 3.4

Two SQP components, Hessian modification and line search, are less addressed
in this paper. Performing them efficiently in a parallel environment is of
course desirable from a practical aspect, while we put the implementation
details of this regard aside in the paper but briefly discuss in this remark.
Our paper is mainly concerned with solving the Newton system (8), which always
requires more computations.

For line search, we note that the augmented Lagrangian $\mathcal{L}_{\eta}$
and its gradient $\nabla\mathcal{L}_{\eta}$ are separable in stages. Thus, we
can easily evaluate them in parallel, where each processor computes a short
horizon e.g., $[n_{i},n_{i+1})$, and a coordinator is needed to sum up the
results of all processors to check (11).

For Hessian modification, a desirable modification ${\hat{H}}$ should be close
to the Hessian $H$ whenever $H$ satisfies the second-order sufficient
condition (SOSC, i.e., (16) with ${\hat{H}}$ being replaced by $H$). If $H$
does not satisfy the condition, we regularize $H$, e.g.
${\hat{H}}=H+(\gamma_{RH}+\|H\|)I$, to let ${\hat{H}}$ satisfy the condition,
which ensures that the full Newton system (8) has a unique solution. See
Assumptions 4 and 6. Thus, we have to check the positive definiteness of the
reduced Hessian $Z^{T}HZ$, where the columns of $Z$ span the null space of the
Jacobian $G=\nabla_{{\boldsymbol{z}}}^{T}f$. It is equivalent to testing for
positive definiteness of $H+c\cdot G^{T}G$ with a scalar $c$ sufficiently
large [41]. Note that $H+c\cdot G^{T}G$ is a block-tridiagonal matrix; thus,
we can apply the parallel Cholesky factorization to check its definiteness in
a parallel environment [63, 64, 12]. As analyzed in [12, Section 2.2], the
total flops of the parallel Cholesky of $M$ processors are less than
$7N(n_{x}+n_{u})^{3}$, so the average flop of a single processor is in order
of $O(N(n_{x}+n_{u})^{3}/M)$, which does not grow with $N$ if the ratio $N/M$
is a constant.

######  Remark 3.5 (extensions to inequality constraints)

The FOTD procedure can be generalized to inequality-constrained problems. In
particular, with inequality constraints
$h_{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k})\leq{\boldsymbol{0}}$,
$k\in[N-1]$, the full SQP problem (14) will additionally have the linearized
inequality constraints
$C_{k}{\boldsymbol{p}}_{k}+D_{k}{\boldsymbol{q}}_{k}\leq{\boldsymbol{0}}$,
$k\in[N-1]$, where $C_{k}=\nabla_{{\boldsymbol{x}}_{k}}^{T}h_{k}$ and
$D_{k}=\nabla_{{\boldsymbol{u}}_{k}}^{T}h_{k}$ are the Jacobian matrices of
$h_{k}$. The FOTD scheme still applies OTD on the full problem (14) so that
the subproblems (15) are inequality-constrained quadratic programs (IQPs).
Several methods with warm-start strategies can be applied on IQPs, such as the
interior-point methods and active-set methods. Within the FOTD scheme, the
exact augmented Lagrangian function (10) should also be adapted to the one
that can handle inequality constraints. See [47] for a particular choice. We
should mention that the alternative designs with the same flavor of FOTD are
also available for dealing with inequality constraints. For example, we can
exploit an active-set SQP scheme, where in each iteration we only consider the
inequality constraints that are in the identified active set and regard them
as equalities. Then, the full problem (14) and the subproblem (15) are still
equality-constrained QPs (EQPs), which are easier to solve than IQPs. Since
the design and analysis of inequality constraints are quite involved, we leave
the above extensions of FOTD to the future, and mainly focus on connecting
FOTD with the Schwarz scheme on equality-constrained problems in this paper.

######  Remark 3.6 (relationships to other schemes)

Besides the Schwarz scheme, FOTD is related to several other methods for
optimal control problems. We introduce the connections in this remark. We
emphasize that the two critical components of FOTD, overlapping temporal
decomposition and augmented Lagrangian merit function, have not been
investigated in the following methods. Also, the convergence of FOTD highly
relies on the sensitivity analysis of NLDPs in [34, 38], which differs from
the following methods.

(a) Direct multiple shooting methods. The direct multiple shooting methods
decompose the full horizon into multiple exclusive short intervals, construct
the subproblems associated with short intervals, and compose the solution
trajectories of subproblems under certain matching conditions. See [4, 5] for
more details. The multiple shooting methods are often employed for real-time
nonlinear

model predictive control [50, 17, 29], where the subproblems are solved
sequentially and an approximate but fast control feedback is desired for each
subproblem. Many efficient solvers that exploit the problem sparsity
structures can be used within the methods, such as qpDUNES [24] and HPMPC
[25]. Compared to the aforementioned literature, FOTD solves (1) in an offline
fashion (i.e., the short intervals do not vary with the sampling time) and
applies an overlapping decomposition on the SQP problem instead of on (1).
More importantly, as analyzed in section 5, FOTD does not use any matching
conditions to compose the subproblem solutions, but requires a large overlap
size $b$. Such a design is inspired by the sensitivity analysis of NLDPs. The
exclusive stages are far from the boundaries where the system perturbations
occur, and the sensitivity results state that the effects of boundary
perturbations decay exponentially fast away from the boundaries. Thus, the
solutions at the exclusive middle stages are still accurate enough even
without matching conditions.

(b) Alternating direction method of multipliers (ADMM). ADMM is another
popular method for optimal control problems [42, 7], where one introduces a
set of consensus constraints to split the full problem into multiple
subproblems, and solves the subproblems in each iteration in parallel. We note
that our terminal cost (4) is conceptually similar to the quadratic proximal
term in [42, (3)]. However, FOTD and ADMM have a few key differences. First,
ADMM is designed based on the augmented Lagrangian method, while FOTD is
designed based on SQP. Thus, their primal-dual updating schemes are quite
different. Second, even if we use an augmented Lagrangian merit function in
FOTD, the function is different from the one in ADMM [7] in that it has a
quadratic penalty on the optimality error (i.e., the last term in (10)). Also,
we use the augmented Lagrangian for the stepsize selection, not for the
direction computation, while ADMM combines the two steps together. Third,
similar to the multiple shooting methods, ADMM decomposes the full problem by
introducing extra consensus constraints without overlaps. Thus, ADMM is
sensitive to the parameter $\mu$, while the overlaps in FOTD largely suppress
the boundary perturbations brought by different choices of $\mu$ (and
imprecise boundary variables ${\boldsymbol{d}}_{i}$). See [38, Figure 6] for
empirical evidence.

(c) Iterative linear-quadratic regulator (ILQR). Another method for optimal
control problems is ILQR [32]. In each step, the search direction of control
variables is solved from a QP that is defined by quadratic approximation of
objective (1a) with linear approximations of constraints (1b). A parallel
implementation of ILQR was designed in [30]. In addition to the difference
that FOTD employs an overlapping decomposition for parallelism, the QP in (14)
also differs from the one in ILQR in that its objective is a quadratic
approximation of the Lagrangian. Further, FOTD updates the state, control, and
dual variables with the same stepsize selected by the line search on the
augmented Lagrangian (10), while ILQR updates the control variables with the
line search on the objective (1a), and the state variables are computed by
applying (1b) in a forward pass.

######  Remark 3.7 (QP solvers for the subproblems)

The subproblems (15) of FOTD are QPs. Many efficient QP or linear system
solvers with wart-start strategies can be adopted for FOTD. For example, the
conjugate gradient (CG) and minimal residual (MINRES) are popular methods for
solving the Newton systems, and the popular solvers qpDUNES [24], HPMPC [25],
and FORCES [19] that exploit the sparsity structures of QPs can also be
applied on (15). We note that the above QP solvers generally require a
positive definite Hessian matrix which (15) does not have. Fortunately, [57]
proposed a convexification procedure to resolve this issue.

##  4 Error of the approximate search direction

We study the difference between
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
and $(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})$,
where the former is the OTD approximation (Step 2) and the latter is the exact
Newton direction of (8). We state the assumptions that are required for the
analysis.

[lower bound on the reduced Hessian] For any iteration $\tau\geq 0$, we let
$Z^{\tau}$ be a matrix whose columns are orthonormal vectors and span the null
space $\\{{\boldsymbol{w}}:G^{\tau}{\boldsymbol{w}}={\boldsymbol{0}}\\}$. We
assume that there exists a uniform constant $\gamma_{RH}>0$ that is
independent of $\tau$, such that

$(Z^{\tau})^{T}{\hat{H}}^{\tau}Z^{\tau}\succeq\gamma_{RH}\cdot I.$ (16)

The matrix on the left hand side is called the reduced Hessian.

[controllability] For any stage $k\in[N]$ and an integer $t\in[N-k]$, the
controllability matrix is defined as

$\Xi_{k,t}({\boldsymbol{z}}_{k:k+t-1})=\begin{pmatrix}B_{k+t-1}&A_{k+t-1}B_{k+t-2}&\cdots&(\prod_{l=1}^{t-1}A_{k+l})B_{k}\end{pmatrix}\in\mathbb{R}^{n_{x}\times
tn_{u}}.$

For any iteration $\tau\geq 0$, we assume that there exist a uniform constant
$\gamma_{C}>0$ and an integer $t>0$ that are independent of $\tau$, such that
for any $k\in[N-t]$, there exists $t_{k}^{\tau}\in[1,t]$ so that
$\Xi_{k,t_{k}^{\tau}}^{\tau}(\Xi_{k,t_{k}^{\tau}}^{\tau})^{T}\succeq\gamma_{C}I$,
where
$\Xi_{k,t_{k}^{\tau}}^{\tau}=\Xi_{k,t_{k}^{\tau}}({\boldsymbol{z}}_{k:k+t_{k}^{\tau}-1}^{\tau})$.

[upper boundedness] For any iteration $\tau\geq 0$, there exists a uniform
constant $\Upsilon_{upper}$ that is independent of $\tau$, such that for any
$k$,
$\max\\{\|{\hat{H}}^{\tau}_{k}\|,\;\|A_{k}^{\tau}\|,\;\|B_{k}^{\tau}\|\\}\leq\Upsilon_{upper}$
.

Since $G^{\tau}=\nabla_{{\boldsymbol{z}}}^{T}f({\boldsymbol{z}}^{\tau})$ has
full row rank, Assumption 4 ensures that Problem (14) has a unique global
solution and the KKT matrix in (8) is invertible [41, Lemma 16.1]. This
assumption is standard in the SQP literature [3] and weaker than the linear-
quadratic convex DP setup studied in [65]. Assumption 4 is specifically used
for DP problems [28, 65, 66, 34, 35, 38]. It ensures that the linearized
dynamical system
${\boldsymbol{p}}_{k+1}=A_{k}{\boldsymbol{p}}_{k}+B_{k}{\boldsymbol{q}}_{k}$
is controllable in at most $t$ steps (in each iteration). That is, given any
initial state ${\boldsymbol{p}}_{k}$ and target state
${\boldsymbol{p}}_{k+t_{k}}$, we can always evolve from ${\boldsymbol{p}}_{k}$
to ${\boldsymbol{p}}_{k+t_{k}}$ by specifying a suitable control sequence
$\\{{\boldsymbol{q}}_{j}\\}_{j=k}^{k+t_{k}-1}$. Assumption 4 is standard in
both SQP and DP literature [3, 65, 34]. In section 5, we will impose a
compactness condition on the SQP iterates, which naturally implies the upper
boundedness of the Hessian and Jacobian matrices. We show in Lemma 4.1 that
Assumptions 4 and 4 imply a uniform lower bound on $G^{\tau}(G^{\tau})^{T}$.

###### Lemma 4.1

Suppose $\max\\{\|A_{k}^{\tau}\|,\;\|B_{k}^{\tau}\|\\}\leq\Upsilon_{upper}$,
then Assumption 4 implies that $G^{\tau}(G^{\tau})^{T}\succeq\gamma_{G}I$
where

$\gamma_{G}=\gamma_{G}(\gamma_{C},t,\Upsilon_{upper})\coloneqq\left(\frac{\gamma_{C}}{\gamma_{C}+\frac{\Upsilon_{upper}^{t+1}}{\Upsilon_{upper}-1}}\right)^{2}\cdot\frac{\min\\{1,\gamma_{C}\\}}{(1+\Upsilon_{upper})^{2t}}.$
(17)

###### Proof 4.2

See Appendix 10.1.

The next result shows that
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ has a unique
solution if $\mu$ is large enough.

###### Lemma 4.3

Suppose Assumptions 4-4 hold for Problem (14). Let

$\bar{\mu}=\bar{\mu}(\gamma_{C},t,\Upsilon_{upper})\coloneqq\frac{32\Upsilon_{upper}^{4t+1}}{\gamma_{C}}.$
(18)

If $\mu\geq\bar{\mu}$, then the reduced Hessian of Problem
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ in (15),
defined similarly to (16), is lower bounded by $\gamma_{RH}I$ for any
iteration $\tau\geq 0$ and any boundary variables ${\boldsymbol{d}}_{i}$. This
implies that ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$
has a unique global solution.

###### Proof 4.4

See Appendix 10.2.

Lemma 4.3 shows that
$\\{{{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})\\}_{i=0}^{M-1}$
are solvable for any iteration $\tau\geq 0$ if $\mu\geq\bar{\mu}$, where
$\bar{\mu}$ is independent of the iteration index $\tau$ and the subproblem
index $i$. An immediate consequence is that Assumptions 4-4 hold for the
subproblems ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$ as well.

###### Corollary 4.5

Suppose Assumptions 4-4 hold for Problem (14) and $\mu\geq\bar{\mu}$ with
$\bar{\mu}$ given by (18). Then, the three conditions: lower bound on the
reduced Hessian, controllability, and upper boundedness, hold for the
subproblems
$\\{{{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})\\}_{i\in[M-1]}$
with any ${\boldsymbol{d}}_{i}$. Furthermore, the condition constants are
independent of $i$ and $\tau$. Specifically, the subproblem
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ for $i\in[M-1]$
satisfies

(i) Assumption 4: the reduced Hessian of
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ is lower
bounded by $\gamma_{RH}I$.

(ii) Assumption 4: the controllability condition of
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ is satisfied
with the same constants $(\gamma_{C},t)$.

(iii) Assumption 4: the boundedness condition of
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$ is satisfied
with constant $\Upsilon_{upper}+\mu$.

###### Proof 4.6

See Appendix 10.3.

We are now able to control the error
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau}-\Delta{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau}-\Delta{\boldsymbol{\lambda}}^{\tau})$.
To ease the notation, we suppress the iteration index $\tau$. The study of the
error
$({\tilde{\Delta}}{\boldsymbol{z}}-\Delta{\boldsymbol{z}},{\tilde{\Delta}}{\boldsymbol{\lambda}}-\Delta{\boldsymbol{\lambda}})$
relies on the primal-dual sensitivity analysis of NLDPs in [34, 38], which
requires Assumptions 4-4. We apply the sensitivity results on the subproblems;
thus Corollary 4.5 is critical. It shows that Assumptions 4-4 hold for the
subproblems as long as they hold for the full problem. In principle, the
sensitivity results suggest that, if we perturb the objective and constraints
on one stage, then the perturbation effects on the optimal solution decay
exponentially fast away from that perturbed stage. In the OTD setup, the
perturbations occur at the two boundaries of the extended intervals. Thus, the
composition that uses only the exclusive part $[n_{i},n_{i+1})$ of the
solution preserves all accurate variables.

Let us first bound
$({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}))-({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\prime}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\prime}))$
for any ${\boldsymbol{d}}_{i},{\boldsymbol{d}}_{i}^{\prime}$. Recall that
$({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}))$
denotes the unique global solution of
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i})$. In the
following presentation, we use $C$ and $C^{\prime}$ to denote universal
constants that are independent of the iteration index $\tau$ and the
subproblem index $i$.

###### Theorem 4.7

Let Assumptions 4-4 hold for Problem (14) and $\mu\geq\bar{\mu}$ with
$\bar{\mu}$ given by (18). There exist constants $C^{\prime}>0$ and
$\rho\in(0,1)$ independent of $i$ and $\tau$, such that for all
$k\in[m_{1},m_{2}]$,

$\max\\{\|{\tilde{\boldsymbol{w}}}_{i,k}^{\star}({\boldsymbol{d}}_{i})-{\tilde{\boldsymbol{w}}}_{i,k}^{\star}({\boldsymbol{d}}_{i}^{\prime})\|,\;\|{\tilde{\boldsymbol{\zeta}}}_{i,k}^{\star}({\boldsymbol{d}}_{i})-{\tilde{\boldsymbol{\zeta}}}_{i,k}^{\star}({\boldsymbol{d}}_{i}^{\prime})\|\\}\leq
C^{\prime}(\rho^{k-m_{1}}\|{\boldsymbol{d}}_{i,1}-{\boldsymbol{d}}_{i,1}^{\prime}\|+\rho^{m_{2}-k}\|{\boldsymbol{d}}_{i,2:4}-{\boldsymbol{d}}_{i,2:4}^{\prime}\|)$

for any two boundary variables ${\boldsymbol{d}}_{i}$ and
${\boldsymbol{d}}_{i}^{\prime}$.

###### Proof 4.8

See Appendix 10.4.

The next theorem characterizes the approximation error of the Newton
direction.

######  Theorem 4.9 (error of approximate direction)

Let Assumptions 4-4 hold for Problem (14) and $\mu\geq\bar{\mu}$ with
$\bar{\mu}$ given by (18). There exist constants $C>0$ and $\rho\in(0,1)$
independent of $\tau$, such that

$\|({\tilde{\Delta}}{\boldsymbol{z}}^{\tau}-\Delta{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau}-\Delta{\boldsymbol{\lambda}})\|\leq
C\rho^{b}\|(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})\|.$
(19)

###### Proof 4.10

See Appendix 10.5.

Theorem 4.9 suggests that the error of the approximate direction decays
exponentially fast in terms of the overlap size $b$. We can naturally expect
that
$({\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
contains enough information to decrease the merit function in each iteration,
provided $b$ is large. We study how the inexactness of the direction affects
SQP in the next section.

##  5 Global convergence of FOTD

To enable general distributed methods for solving (8), we suppose that a
decomposition method outputs a direction
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
satisfying

$\|({\tilde{\Delta}}{\boldsymbol{z}}^{\tau}-\Delta{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau}-\Delta{\boldsymbol{\lambda}}^{\tau})\|\leq\delta\cdot\|(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})\|,\quad\text{
for }\delta\in(0,1).$ (20)

Theorem 4.9 shows that FOTD specializes (20) with $\delta=C\rho^{b}$. In this
section, we study the convergence of SQP with the merit function (10) and the
direction
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$.
We require a compactness assumption to strengthen Assumption 4 to hold in a
compact set. We recall that the SQP iterates are generated by (9) and (11).

[compactness] There exists a compact set $\mathcal{Z}\times\Lambda$, where
$\mathcal{Z}=\mathcal{Z}_{0}\times\cdots\times\mathcal{Z}_{N}$ and
$\Lambda=\Lambda_{0}\times\cdots\times\Lambda_{N}$, such that
$({\boldsymbol{z}}^{\tau}+\alpha{\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau}+\alpha{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})\in\mathcal{Z}\times\Lambda$,
$\forall\tau\geq 0,\alpha\in[0,1]$; and we assume that $\\{g_{k},f_{k}\\}$ are
thrice continuously differentiable and

$\max\\{\sup_{\mathcal{Z}_{k}\times\Lambda_{k+1}}\|H_{k}({\boldsymbol{z}}_{k},{\boldsymbol{\lambda}}_{k+1})\|,\;\sup_{\mathcal{Z}_{k}}\|A_{k}({\boldsymbol{z}}_{k})\|,\;\sup_{\mathcal{Z}_{k}}\|B_{k}({\boldsymbol{z}}_{k})\|\\}\leq\Upsilon_{upper},$
(21)

for a constant $\Upsilon_{upper}>0$ independent of $k$.

We use the same constant $\Upsilon_{upper}$ as in Assumption 4 for notational
simplicity. The compactness assumption is standard in the SQP literature [3,
Proposition 4.15]. One equivalent assumption is to assume that the sublevel
set
$\\{({\boldsymbol{z}},{\boldsymbol{\lambda}}):\mathcal{L}_{\eta}({\boldsymbol{z}},{\boldsymbol{\lambda}})\leq\mathcal{L}_{\eta}^{0}\\}$
is contained in a compact set. Furthermore, assuming the existence of the
third derivatives on $\\{g_{k},f_{k}\\}$ is common for the augmented
Lagrangian merit function (10), since $\nabla^{2}\mathcal{L}_{\eta}$ requires
$\nabla^{3}g_{k}$ and $\nabla^{3}f_{k}$ [3, 69, 36, 37]. Note that the third
derivatives are only required in the analysis, and not computed or used in the
algorithm.

The following lemma is an immediate consequence of Assumption 5.

###### Lemma 5.1

Under Assumption 5, there exists a constant $\Upsilon_{HG}>0$ such that

$\max\\{\sup_{\mathcal{Z}\times\Lambda}\|H({\boldsymbol{z}},{\boldsymbol{\lambda}})\|,\sup_{\mathcal{Z}}\|G({\boldsymbol{z}})\|\\}\leq\Upsilon_{HG}.$
(22)

Further, if Assumptions 4-4 hold as well, there exists a constant
$\Upsilon_{KKT}>0$ that is independent of $\tau$, such that

$\left\|\begin{pmatrix}{\hat{H}}^{\tau}&(G^{\tau})^{T}\\\
G^{\tau}&{\boldsymbol{0}}\end{pmatrix}^{-1}\right\|\leq\Upsilon_{KKT},\quad\forall\tau\geq
0.$

###### Proof 5.2

See Appendix 11.1

The bound (22) suggests that $\|G^{\tau}\|\leq\Upsilon_{HG}$, $\forall\tau\geq
0$ (since ${\boldsymbol{z}}^{\tau}\in\mathcal{Z}$). Since (21) implies
$\max\\{\|A_{k}^{\tau}\|,\;\\\ \|B_{k}^{\tau}\|\\}\leq\Upsilon_{upper}$ for
any $k\in[N-1]$, under Assumption 5, we only need
$\|{\hat{H}}^{\tau}\|\leq~{}\Upsilon_{upper}$ from Assumption 4. To simplify
the presentation, we define
$\Upsilon=\max\\{\Upsilon_{upper},\Upsilon_{KKT},\Upsilon_{HG}\\}$.

We now show that
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
is a descent direction of $\mathcal{L}_{\eta}^{\tau}$ provided $\eta_{1}$ is
sufficiently large and $\eta_{2},\delta$ are sufficiently small.

###### Theorem 5.3

Suppose Assumptions 4, 4, 4, 5 hold for the SQP iterates with a search
direction
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
satisfying (20). If

$\eta_{1}\geq\frac{17}{\eta_{2}\gamma_{G}},\quad\eta_{2}\leq\frac{\gamma_{RH}}{12\Upsilon^{2}},\quad\delta\leq\frac{\eta_{2}\gamma_{G}}{9\eta_{1}\Upsilon^{2}},$
(23)

where $\gamma_{G}$ is defined in (17), then

$\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}_{\eta}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}_{\eta}^{\tau}\end{pmatrix}^{T}\begin{pmatrix}{\tilde{\Delta}}{\boldsymbol{z}}^{\tau}\\\
{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}\leq-\frac{\eta_{2}}{2}\left\|\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}^{\tau}\end{pmatrix}\right\|^{2}.$
(24)

###### Proof 5.4

See Appendix 11.2.

By the results of Theorem 5.3 and the mean value theorem, we immediately know
that, under the presented conditions, the stepsize $\alpha_{\tau}$ to satisfy
the Armijo condition (11) can be found by the backtracking line search, and
the SQP iterates with the direction
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
make successful progress towards a stationary point.

Theorem 5.3 also shows the importance of using the exact augmented Lagrangian
merit function (10). In particular, as proved in (67), we rely on a critical
property: for suitably chosen $\eta=(\eta_{1},\eta_{2})$, there exists a small
constant $\kappa(\eta_{2})>0$ (depending on $\eta_{2}$) such that

$\displaystyle\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}_{\eta}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}_{\eta}^{\tau}\end{pmatrix}^{T}\begin{pmatrix}\Delta{\boldsymbol{z}}^{\tau}\\\
\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}\leq-\kappa(\eta_{2})\bigg{(}\underbrace{\left\|\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}^{\tau}\end{pmatrix}\right\|^{2}}_{\text{ensure
convergence}}+\underbrace{\left\|\begin{pmatrix}\Delta{\boldsymbol{z}}^{\tau}\\\
\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}\right\|^{2}}_{\text{allow
approximation}}\bigg{)}.$ (25)

The first term is needed for global convergence when we accumulate the descent
over iterations. The second term allows us to replace
$(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})$ by
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$
because, as proved in (69),

$\left\|\begin{pmatrix}\nabla_{{\boldsymbol{z}}}\mathcal{L}_{\eta}^{\tau}\\\
\nabla_{{\boldsymbol{\lambda}}}\mathcal{L}_{\eta}^{\tau}\end{pmatrix}\right\|\left\|\begin{pmatrix}{\tilde{\Delta}}{\boldsymbol{z}}^{\tau}-\Delta{\boldsymbol{z}}^{\tau}\\\
{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau}-\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}\right\|\lesssim\delta\left\|\begin{pmatrix}\Delta{\boldsymbol{z}}^{\tau}\\\
\Delta{\boldsymbol{\lambda}}^{\tau}\end{pmatrix}\right\|^{2},$

where $\lesssim$ means that the inequality holds up to a constant multiplier.
Therefore, for small enough $\delta$, the margin
$\|(\Delta{\boldsymbol{z}}^{\tau},\Delta{\boldsymbol{\lambda}}^{\tau})\|$ in
(25) allows for an approximation error. Our analysis rules out the exact
penalty merit functions that depend only on the primal variables
${\boldsymbol{z}}$. Such a function $\mathcal{M}({\boldsymbol{z}})$ satisfies
$\nabla_{{\boldsymbol{z}}}^{T}\mathcal{M}^{\tau}\Delta{\boldsymbol{z}}^{\tau}\leq-\kappa\|\Delta{\boldsymbol{z}}^{\tau}\|^{2}$
and
$\|\nabla_{{\boldsymbol{z}}}\mathcal{M}^{\tau}\|\lesssim\|\Delta{\boldsymbol{z}}^{\tau}\|$
(e.g., [22, 26, 49, 44, 3]). If we approximate $\Delta{\boldsymbol{z}}^{\tau}$
by ${\tilde{\Delta}}{\boldsymbol{z}}^{\tau}$, Theorems 4.7, 4.9 suggest that
the approximation error
$\|{\tilde{\Delta}}{\boldsymbol{z}}^{\tau}-\Delta{\boldsymbol{z}}^{\tau}\|$
also depends on $\|\Delta{\boldsymbol{\lambda}}^{\tau}\|$. Thus, the margin
$\|\Delta{\boldsymbol{z}}^{\tau}\|^{2}$ may not be enough to endure the
approximation error. This reveals the benefits of using the primal-dual exact
merit functions based on (10).

We summarize the global convergence results in the next theorem.

######  Theorem 5.5 (global convergence)

Suppose Assumptions 4, 4, 4, 5 hold for the SQP iterates with
$\\{({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})\\}_{\tau}$
satisfying (20) and parameters $(\eta,\delta)$ satisfying (23), then
$\|\nabla\mathcal{L}^{\tau}\|\rightarrow 0$ as $\tau\rightarrow\infty$.
Furthermore, for the FOTD iterates, if $\mu\geq\bar{\mu}$ with $\bar{\mu}$
given by (18) as well as $b$ satisfies

$b\geq\frac{\log\left(9C\eta_{1}\Upsilon^{2}/(\eta_{2}\gamma_{G})\right)}{\log(1/\rho)}$
(26)

with $C>0$ and $\rho\in(0,1)$ from Theorem 4.9, then the FOTD iterates
$\\{({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})\\}_{\tau}$
generated by Algorithm 2 satisfy $\|\nabla\mathcal{L}^{\tau}\|\rightarrow 0$
as $\tau\rightarrow\infty$.

###### Proof 5.6

See Appendix 11.3.

We have shown the global convergence of FOTD. Theorem 5.5 complements the
local convergence result of the Schwarz scheme [38] and answers Q1 raised in
section 1. The next section studies the local convergence of FOTD. We show
that FOTD and the Schwarz have the same local behavior.

To end this section, we discuss an adaptivity extension of our scheme.

######  Remark 5.7 (adaptivity on penalty parameters)

Given the conditions (23) on $(\eta,\delta=C\rho^{b})$, it is possible to
design a scheme that adaptively selects the suitable penalty parameters $\eta$
and overlap size $b$. Since the upper/lower bound constants in (23) do not
depend on $\tau$, we can design a while loop to achieve this goal. In
particular, we let $\nu>1$ be fixed. While (24) does not hold, we let

$\eta_{2}\leftarrow\eta_{2}/\nu,\quad\eta_{1}\leftarrow\eta_{1}\nu^{2},\quad\delta\leftarrow\delta/\nu^{4}.$

Then, we know $\eta_{1}\eta_{2}$ increases by a factor of $\nu$, and
$\eta_{2}$ and $\delta\eta_{1}/\eta_{2}$ decrease by a factor of $1/\nu$.
Thus, for sufficiently large $\tau$, all parameters are stabilized. We also
note that $\delta\leftarrow\delta/\nu^{4}$ is equivalent to letting
$b\leftarrow b+4\log\nu/\log(1/\rho)$, where only $\rho$ has to be tuned
manually.

##  6 Local convergence of FOTD

From now on, we suppose the FOTD iterates satisfy
$({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$
$\rightarrow({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$ as
$\tau\rightarrow\infty$. For two positive sequences $\\{a_{\tau}\\}$ and
$\\{b_{\tau}\\}$, $a_{\tau}=O(b_{\tau})$ if $a_{\tau}/b_{\tau}$ is uniformly
bounded over $\tau$; $a_{\tau}=o(b_{\tau})$ if $a_{\tau}/b_{\tau}\rightarrow
0$ as $\tau\rightarrow\infty$. Our local analysis is divided into three steps:

  1. (a)

we show that $\alpha_{\tau}=1$ is selected for the Armijo condition (11) when
$\tau$ is large.

  2. (b)

we show a relationship between FOTD and the Schwarz scheme: Line 6 in
Algorithm 2 is equivalent to performing a single Newton step for subproblems
in Line 5 in Algorithm 1.

  3. (c)

we prove that FOTD converges linearly, with a linear rate that decays
exponentially in $b$.

We present additional assumptions for local analysis.

[Hessian approximation vanishes] We assume that
$\|H^{\tau}-{\hat{H}}^{\tau}\|=o(1)$, where
$H^{\tau}=\nabla_{{\boldsymbol{z}}}^{2}\mathcal{L}^{\tau}$ is the Lagrangian
Hessian and ${\hat{H}}^{\tau}$ is its approximation.

A vanishing Hessian modification is typical for the local analysis of SQP to
have a superlinear (or quadratic) convergence [6]. As discussed in Remark 3.4,
we can check (16) for Hessian $H^{\tau}$ to decide if a modification is
needed, since (16) holds locally if SOSC is satisfied at
$({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$. Equivalently, we
check the positive definiteness of $H^{\tau}+c(G^{\tau})^{T}G^{\tau}$ for a
constant $c$, which is a block-tridiagonal matrix. A parallel Cholesky
decomposition is applicable in this regard.

[local Lipschitz continuity]

We assume there exists a constant $\Upsilon_{L}$ independent of $k$ such that,
for any two points $({\boldsymbol{z}},{\boldsymbol{\lambda}})$ and
$({\boldsymbol{z}}^{\prime},{\boldsymbol{\lambda}}^{\prime})$ sufficiently
close to $({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$,

$\displaystyle\max\\{\|A_{k}({\boldsymbol{z}}_{k})-A_{k}({\boldsymbol{z}}_{k}^{\prime})\|,\;\|B_{k}({\boldsymbol{z}}_{k})-B_{k}({\boldsymbol{z}}_{k}^{\prime})\|\\}\leq$
$\displaystyle\Upsilon_{L}\|{\boldsymbol{z}}_{k}-{\boldsymbol{z}}_{k}^{\prime}\|,$
$\displaystyle\|H_{k}({\boldsymbol{z}}_{k},{\boldsymbol{\lambda}}_{k+1})-H_{k}({\boldsymbol{z}}_{k}^{\prime},{\boldsymbol{\lambda}}_{k+1}^{\prime})\|\leq$
$\displaystyle\Upsilon_{L}\left\|({\boldsymbol{z}}_{k}-{\boldsymbol{z}}_{k}^{\prime};{\boldsymbol{\lambda}}_{k+1}-{\boldsymbol{\lambda}}_{k+1}^{\prime})\right\|.$

Assumption 6 strengthens the boundedness condition in Assumption 5 to the
(local) Lipschitz continuity. We start the local analysis with Step (a).

Step (a): A unit stepsize is accepted. We have the following theorem.

###### Theorem 6.1

Suppose Assumptions 4, 4, 4, 5, 6, 6 hold for the SQP iterates with search
directions
$\\{({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})\\}_{\tau}$
satisfying (20) and parameters $(\eta,\delta)$ satisfying

$\eta_{1}\geq\frac{17}{\eta_{2}\gamma_{G}},\quad\eta_{2}\leq\frac{\gamma_{RH}}{12\Upsilon^{2}},\quad\delta\leq\frac{1/2-\beta}{3/2-\beta}\cdot\frac{\eta_{2}\gamma_{G}}{9\eta_{1}\Upsilon^{2}},$
(27)

then $\alpha_{\tau}=1$ for all sufficiently large $\tau$.

###### Proof 6.2

See Appendix 12.1.

The condition (27) on $\delta$ is stronger than (23) up to a multiplier
depending on $\beta$. For FOTD, we suppose $\mu\geq\bar{\mu}$ with $\bar{\mu}$
given by (18), and let $\delta=C\rho^{b}$ in (27) to get a condition on $b$ as
in (29).

Step (b): A relation between FOTD and the Schwarz scheme. We establish a
relationship between FOTD and the Schwarz scheme. We will show that, by
specifying
${\boldsymbol{d}}^{\tau}=({\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}})$
in Problem (15), FOTD is equivalent to performing one Newton step for
subproblems (3) with a warm-start initialization in the Schwarz scheme. The
result is summarized in the next theorem.

###### Theorem 6.3

Given the current iterate
$({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$, we consider two
procedures:

(a) Schwarz with a warm-start initialization: for $i\in[M-1]$, we specify
boundary variables
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{x}}_{m_{1}}^{\tau};{\boldsymbol{x}}_{m_{2}}^{\tau};{\boldsymbol{u}}_{m_{2}}^{\tau};{\boldsymbol{\lambda}}_{m_{2}+1}^{\tau})$;
perform one full Newton step for
${\mathcal{P}}_{\mu}^{i}({\boldsymbol{d}}_{i}^{\tau})$ (3) at
$\mathcal{D}_{i}({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$ and
get
$({\tilde{\boldsymbol{z}}}_{i}^{\tau+1},{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\tau+1})$;
then let
$({\boldsymbol{z}}^{\tau+1},{\boldsymbol{\lambda}}^{\tau+1})=\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i}^{\tau+1},{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\tau+1})\\}_{i})$.

(b) FOTD scheme without Hessian approximation: for $i\in[M-1]$, we specify
boundary variables
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}})$;
solve ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}({\boldsymbol{d}}_{i}^{\tau})$ in
(15) with ${\hat{H}}_{k}^{\tau}=H_{k}^{\tau}$, $\forall k\in[m_{1},m_{2}]$;
obtain
$({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))$;
then let
$({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})=\mathcal{C}(\\{({\tilde{\boldsymbol{w}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}),{\tilde{\boldsymbol{\zeta}}}_{i}^{\star}({\boldsymbol{d}}_{i}^{\tau}))\\}_{i})$
and update as
$({\boldsymbol{z}}^{\tau+1},{\boldsymbol{\lambda}}^{\tau+1})=({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})+({\tilde{\Delta}}{\boldsymbol{z}}^{\tau},{\tilde{\Delta}}{\boldsymbol{\lambda}}^{\tau})$.

Then, starting from $({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})$,
both procedures generate the same next iterate
$({\boldsymbol{z}}^{\tau+1},{\boldsymbol{\lambda}}^{\tau+1})$.

###### Proof 6.4

See Appendix 12.2.

Theorem 6.3 reveals a strong relation between the Schwarz scheme and FOTD.
Locally, FOTD can be seen as an improvement of the warm-start Schwarz scheme,
where a single Newton step is performed for the subproblems, instead of
solving the subproblems to optimality as the original Schwarz did. The warm
initialization is recommended by [38] for practical purpose, which avoids the
case where the solutions of the same subproblem in different iterations are
very distinct.

Step (c): local linear convergence of FOTD. We now establish the local
convergence rate for FOTD. For $\tau\geq 0$ we let

$\Psi^{\tau}\coloneqq\max_{k\in[N]}\Psi_{k}^{\tau}\coloneqq\max_{k\in[N]}\|({\boldsymbol{z}}_{k}^{\tau},{\boldsymbol{\lambda}}_{k}^{\tau})-({\boldsymbol{z}}^{\star}_{k},{\boldsymbol{\lambda}}^{\star}_{k})\|.$
(28)

We require the following lemma that shows the one-step error recursion. We use
$C_{1},C_{2}$ to denote generic constants that are ensured to exist, but may
differ from the constant $C$ in Theorem 4.9. The constant $\rho$ is from
Theorem 4.9.

###### Lemma 6.5

Consider the FOTD iterates
$\\{({\boldsymbol{z}}^{\tau},{\boldsymbol{\lambda}}^{\tau})\\}_{\tau}$ under
Assumptions 4, 4, 4, 5, 6, 6 and suppose $\mu\geq\bar{\mu}$ with $\bar{\mu}$
given by $\eqref{equ:mu:cond}$ and $\eta$ satisfies (27). There exist
constants $C_{1}>0$ and $\rho\in(0,1)$ independent of
$\tau,\eta_{1},\eta_{2},\beta$, such that, if $b$ satisfies

$b\geq\frac{\log\left\\{(3/2-\beta)C_{1}\eta_{1}/((1/2-\beta)\eta_{2})\right\\}}{\log(1/\rho)},$
(29)

then, for all sufficiently large $\tau$ (in each case below, $m_{1},m_{2}$
depend on $i$ correspondingly, see (2)),

$\displaystyle\Psi_{k}^{\tau+1}\leq$ $\displaystyle
o(\Psi^{\tau})+C_{1}\left\\{\rho^{k-m_{1}}\|{\boldsymbol{x}}_{m_{1}}^{\tau}-{\boldsymbol{x}}^{\star}_{m_{1}}\|+\rho^{m_{2}-k}\left\|\begin{pmatrix}{\boldsymbol{z}}_{m_{2}}^{\tau}-{\boldsymbol{z}}^{\star}_{m_{2}}\\\
{\boldsymbol{\lambda}}_{m_{2}+1}^{\tau}-{\boldsymbol{\lambda}}^{\star}_{m_{2}+1}\end{pmatrix}\right\|\right\\},\;\;\begin{subarray}{c}\forall
k\in[n_{i},n_{i+1})\\\ i\in[M-2]\end{subarray}\;\;,$ (30a)
$\displaystyle\Psi_{k}^{\tau+1}\leq$ $\displaystyle
o(\Psi^{\tau})+C_{1}\rho^{k-m_{1}}\|{\boldsymbol{x}}_{m_{1}}^{\tau}-{\boldsymbol{x}}^{\star}_{m_{1}}\|,\quad\forall
k\in[n_{M-1},n_{M}].$ (30b)

###### Proof 6.6

See Appendix 12.3.

Lemma 6.5 relies on the decay structure of the KKT matrix inverse, established
in [35, Lemma 2]. In particular, the authors showed that, if Assumptions 4-4
hold for subproblems (15) (verified in Corollary 4.5), the block matrices of
KKT inverse corresponding to each stage have an exponentially decay structure
(see [35, Figure 2]). The core of such a result is the primal-dual sensitivity
analysis of NLDPs [34, 38], which we also made use of in the proof of Theorem
4.7.

Seeing from (30), the first $o(\Psi^{\tau})$ term is the algorithmic
convergence rate, which is superlinear and achieved by the SQP framework; the
second term
$C_{1}\\{\rho^{k-m_{1}}\|{\boldsymbol{x}}_{m_{1}}^{\tau}-{\boldsymbol{x}}^{\star}_{m_{1}}\|+\rho^{m_{2}-k}\|({\boldsymbol{z}}_{m_{2}}^{\tau},{\boldsymbol{\lambda}}_{m_{2}+1}^{\tau})\\\
-({\boldsymbol{z}}^{\star}_{m_{2}},{\boldsymbol{\lambda}}^{\star}_{m_{2}+1})\|\\}$
has a linear convergence rate, which is brought by the horizon truncation in
OTD. Specifically, the perturbations come from the misspecification of
boundary variables: we use
${\boldsymbol{d}}_{i}^{\tau}=({\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}};{\boldsymbol{0}})$
while
${\boldsymbol{d}}_{i}^{\star}=(\Delta{\boldsymbol{x}}_{m_{1}};\Delta{\boldsymbol{x}}_{m_{2}};\Delta{\boldsymbol{u}}_{m_{2}};\Delta{\boldsymbol{\lambda}}_{m_{2}+1})$.

Our result in Lemma 6.5 is different from [35, Theorem 2], where the
algorithmic convergence rate is quadratic as they used the unperturbed Hessian
$H^{\tau}$ in each iteration. The Hessian modification is necessary in our
case as it ensures the global convergence. Our result is also different from
[38, Theorem 7], which has no algorithmic rate as they solved the subproblems
to optimality. Our result reveals that, even if we do not solve the
subproblems to optimality, as long as the algorithmic rate is faster than
linear (superlinear in our case), the perturbation rate will always dominate
for large $\tau$, which is linear. Therefore, a local linear convergence is
still achieved.

We summarize the local convergence guarantee in the next theorem.

######  Theorem 6.7 (uniform linear convergence)

Under the setup of Lemma 6.5, we have for all sufficiently large $\tau$ that
$\Psi^{\tau+1}\leq 3C_{1}\rho^{b}\cdot\Psi^{\tau}$ for constants $C_{1}>0$ and
$\rho\in(0,1)$ from Lemma 6.5.

###### Proof 6.8

We use the facts that
${\boldsymbol{x}}_{0}^{\tau}={\boldsymbol{x}}^{\star}_{0}=\bar{\boldsymbol{x}}_{0}$
(see the proof of Theorem 4.9) and $\min\\{k-m_{1},m_{2}-k\\}\geq b$ for all
$k\in[n_{i},n_{i+1}]$, and apply (30). Then, we obtain $\Psi_{k}^{\tau+1}\leq
3C_{1}\rho^{b}\Psi^{\tau}$ for all $k$. Thus, by (28), we complete the proof.

We finally summarize our convergence analysis of FOTD in the next theorem.

######  Theorem 6.9 (convergence of FOTD in Algorithm 2)

Consider Algorithm 2 under Assumptions 4-4, 5, and suppose $\mu$ satisfies
(18). There exist constants $C_{2}>0,\rho\in(0,1)$ independent of $\tau$ and
algorithmic parameters $\eta_{1},\eta_{2},\beta$, such that if

$\eta_{1}\geq\frac{17}{\eta_{2}\gamma_{G}},\quad\eta_{2}\leq\frac{\gamma_{RH}}{12\Upsilon^{2}},\quad
b\geq\frac{\log\left\\{(3/2-\beta)C_{2}\eta_{1}/((1/2-\beta)\eta_{2})\right\\}}{\log(1/\rho)},$

then $\|\nabla\mathcal{L}^{\tau}\|\rightarrow 0$ as $\tau\rightarrow\infty$.
Moreover, if Assumptions 6, 6 hold locally as well, then for all sufficiently
large $\tau$, $\alpha_{\tau}=1$ and

$\Psi^{\tau+1}\leq(C_{2}\rho^{b})\Psi^{\tau}.$ (31)

###### Proof 6.10

By Theorems 5.5, 6.1, 6.7, Lemma 6.5 and rescaling $C_{2}$ properly, we
complete the proof.

The result in (31) matches the local result (6) proved in [38]; thus, we
answer the Q2 raised in section 1—it is not necessary to solve the subproblems
to optimality for achieving the local linear convergence. As mentioned
earlier, as long as the algorithmic rate is faster than linear (e.g., one
Newton step for each subproblem), the local linear rate, induced by the decay
of the sensitivity of perturbations, is always achieved. We should also
mention that both the Schwarz scheme and FOTD have linear rates of the form
$C\rho^{b}$. Since the analyses of both algorithms only claim the existence of
some constants $C>0$ and $\rho\in(0,1)$ (see Theorem 6.9 and [38, Theorem 8]),
we can always use larger constants between the two algorithms to make their
linear rates identical. In fact, comparing their constant $C$ can be difficult
since it depends on the sharpness of the derivation and various quantities of
the problem (e.g., $\Upsilon,\gamma_{G}$ etc.). However, the constant $\rho$
is the same, and is from [34, Theorem 5.7]. More importantly, both their
linear rates decay exponentially in the overlap size $b$.

##  7 Numerical experiments

We first conduct a numerical experiment on a toy NLDP in [35]:

$\displaystyle\min_{{\boldsymbol{x}},{\boldsymbol{u}}}\;\;$
$\displaystyle\sum_{k=0}^{N-1}\left\\{2\cos(x_{k}-d_{k})^{2}+C_{1}(x_{k}-d_{k})^{2}-C_{2}(u_{k}-d_{k})^{2}\right\\}+C_{1}x_{N}^{2},$
(32a) s.t. $\displaystyle x_{k+1}=x_{k}+u_{k}+d_{k},\quad\forall k\in[N-1],$
(32b) $\displaystyle x_{0}=0.$ (32c)

Here, $n_{x}=n_{u}=1$, and references $\\{d_{k}\\}$ are specified later. As
checked in [35], if $C_{1}-2>4|C_{2}|$, then
$Z^{T}({\boldsymbol{z}})H({\boldsymbol{z}},{\boldsymbol{\lambda}})Z({\boldsymbol{z}})\succeq(C_{1}-2-4|C_{2}|)/4\cdot
I\;\;\forall{\boldsymbol{z}},{\boldsymbol{\lambda}}$. Thus, we simply let
${\hat{H}}^{\tau}=H^{\tau}$ in implementation.

We implement seven methods: one centralized method IPOPT [60] (which is our
baseline), and six parallel methods including the proposed FOTD, the Schwarz
method [38], the direct multiple shooting method [4], the iterative
differential dynamic programming (IDDP) method [55], ILQR [32], and ADMM [42].
The scheme of IDDP is almost the same as ILQR, except that the control
variables of QPs are computed by rolling out the policies along the original
nonlinear dynamics rather than the linearized dynamics [48]. Since the
dynamics (32b) are linear, IDDP and ILQR are identical in this case. We will
study a temperature control problem of thin plates later, where we then have
nonlinear dynamics. We aim to demonstrate three points in this experiment: (i)
FOTD converges globally while the Schwarz may not converge within a reasonable
computational budget for some initializations. (ii) FOTD exhibits at least
linear convergence locally, and the larger overlap size $b$ leads to the
faster convergence. (iii) FOTD is a superior parallel method. It is as
competitive as the centralized solver IPOPT, and robust to the penalty
parameter $\mu$ (cf. (15a)). To illustrate the first point, we generate random
initial iterates that are shared by all methods, and see if each method
converges for all initializations. To illustrate the second point, we plot
$\|\nabla\mathcal{L}^{\tau}\|$ v.s. $\tau$ for FOTD with different $b$, and
see how $\|\nabla\mathcal{L}^{\tau}\|$ behaves on the tail. To illustrate the
third point, we compare the KKT residual $\|\nabla\mathcal{L}^{\tau}\|$ and
running time for all methods, and vary $\mu$ for FOTD to test its robustness.

Simulation setting: We consider three cases in Table 1. For all parallel
methods, we use the same horizon decomposition; that is, we decompose $[0,N]$
evenly with length $M$ for each (exclusive) short interval. Different from the
methods of multiple shooting, ILQR, IDDP, and ADMM, the Schwarz and FOTD have
overlaps between two successive short intervals. We vary the overlap size $b$
from $\\{1,5,25\\}$. For FOTD, we let $\eta_{1}=10$, $\eta_{2}=\beta=0.1$, and
vary $\mu$ in a wide range $\\{1,25,125\\}$. The setup of $\mu$ is shared by
the Schwarz and ADMM, both of which also require the penalty parameter. When
doing the backtracking line search, we decrease the stepsize by a factor of
$0.9$ each time until the line search condition is satisfied. For each case in
Table 1 and each setup of $b$ and $\mu$, we generate 5 initial iterates for
FOTD that are shared by other methods: one is
$({\boldsymbol{z}}^{0},{\boldsymbol{\lambda}}^{0})=({\boldsymbol{0}},{\boldsymbol{0}})$
and the other four are from
$x_{k}^{0},u_{k}^{0},\lambda_{k}^{0}\stackrel{{\scriptstyle
iid}}{{\sim}}\text{Uniform}(-10^{5},10^{5})$ (with $x_{0}^{0}=0$). For those
methods that have to solve nonlinear subproblems (i.e., the Schwarz, multiple
shooting, ADMM), we apply Julia/JuMP package [21] with IPOPT solver [60]. For
all methods except ADMM, we stop the iteration if

$\|\nabla\mathcal{L}^{\tau}\|\leq 10^{-6}\quad\quad\text{ OR
}\quad\quad\|({\boldsymbol{z}}^{\tau+1}-{\boldsymbol{z}}^{\tau};{\boldsymbol{\lambda}}^{\tau+1}-{\boldsymbol{\lambda}}^{\tau})\|\leq
10^{-6}.$ (33)

Since FOTD, ILQR, and IDDP have comparable computations in each iteration
(i.e., they all solve QPs), we regard them as converged if they trigger the
condition (33) within 40 iterations budget (we see from Figure 2 that FOTD
actually needs much less iterations). For the multiple shooting and Schwarz
that solve NLDPs, we reduce the iteration budget a bit to 30 for sake of a
comparable total computation cost. For ADMM, we observe in our experiment that
its KKT sequence has a long flat tail (see [38, Figure 6]). Thus, we prefer to
stop the ADMM iteration early by relaxing the condition (33) to
$\|\nabla\mathcal{L}^{\tau}\|\leq 10^{-6}$ OR
$\|({\boldsymbol{z}}^{\tau+1}-{\boldsymbol{z}}^{\tau};{\boldsymbol{\lambda}}^{\tau+1}-{\boldsymbol{\lambda}}^{\tau})\|\leq
10^{-3}$, and increase its iteration budget to 100. For the methods that do
not compute ${\boldsymbol{\lambda}}$ (i.e., multiple shooting, ADMM, ILQR,
IDDP), we let
${\boldsymbol{\lambda}}^{\tau}=-(G^{\tau}(G^{\tau})^{T})^{-1}G^{\tau}\nabla
g({\boldsymbol{z}}^{\tau})$ when evaluating $\nabla\mathcal{L}^{\tau}$. Note
that ${\boldsymbol{\lambda}}^{\tau}\rightarrow{\boldsymbol{\lambda}}^{\star}$
as ${\boldsymbol{z}}^{\tau}\rightarrow{\boldsymbol{z}}^{\star}$. In addition
to the above settings, we try three linear system solvers for FOTD: one is
sparse LU (which is a default choice and adopted by other methods), and the
other two are generalized minimal residual (GMRES) and induced dimension
reduction (IDR) methods implemented in Julia/IterativeSolvers package.

Table 1: Simulation setups Cases | $N$ | $M$ | $(C_{1},C_{2})$ | $d_{k}$
---|---|---|---|---
Case 1 | $5000$ | $50$ | $(8,1)$ | $1$
Case 2 | $5000$ | $100$ | $(15,3)$ | $100\sin(k)^{2}$
Case 3 | $10000$ | $100$ | $(12,2)$ | $5\sin(k)$

Result summary: First, we investigate whether different methods converge
within the computational budget for all five initializations. We observe that
the multiple shooting, ILQR, and FOTD converge by triggering (33) for all
three cases in Table 1 and for all initializations and setups. The Schwarz
converges for most of cases, but does not converge within the budget for Case
2 with $b=1$. In particular, for this case, there are $2,4,5$ out of $5$
initializations that the Schwarz does not trigger (33) when setting
$\mu=1,25,125$, respectively. For ADMM, it converges for all three cases with
$\mu=1$, but does not converge for all three cases with $\mu=25$ and $125$. We
do not claim that ADMM diverges with large $\mu$ (its KKT is actually below
$10^{-3}$ with $>2000$ iterations), but we clearly see that ADMM is not as
robust as the Schwarz and FOTD to the parameter $\mu$.

Second, we draw the KKT convergence plots for FOTD in Figure 2. We apply the
sparse LU solver for solving QPs and, for Cases 1, 2, and 3, we take $\mu=1$,
$25$, and $125$ as examples respectively. We emphasize that the other setups
of $\mu$ of each case have similar convergence behavior (as revealed by Tables
2-4). From Figure 2, we observe that FOTD converges for all initializations
with different $\mu$ for different cases and exhibits between linear and
superlinear convergence locally. Its performance is robust to $\mu$, and a
larger $b$ generally leads to a faster convergence (cf. Figures 2g-2i). Our
observation is consistent with Theorem 6.9.

(a) Case 1, $b=1$

(b) Case 1, $b=5$

(c) Case 1, $b=25$

(d) Case 2, $b=1$

(e) Case 2, $b=5$

(f) Case 2, $b=25$

(g) Case 3, $b=1$

(h) Case 3, $b=5$

(i) Case 3, $b=25$

Figure 2: Convergence figures. The three plots on each row correspond to the same case in Table 1, while the three plots on each column correspond to the same setup of $b$. Each plot has 5 lines corresponding to 5 initial iterates. The blue line that is separated from the other four lines corresponds to the initial point $({\boldsymbol{z}}^{0},{\boldsymbol{\lambda}}^{0})=({\boldsymbol{0}},{\boldsymbol{0}})$. We observe that FOTD converges for all initializations. It exhibits between linear and superlinear convergence locally, and a larger $b$ generally leads to a faster convergence. Table 2: The KKT residual and running time for all methods implemented on Case 1. For the Schwarz and FOTD, the smallest KKT residual and least running time among different setups of $b$ are highlighted for each setup of $\mu$. Type | Method | KKT residual $(10^{-7})$ | Time (sec.)
---|---|---|---
Centralized | IPOPT | 172.686 | 0.851
Decomposed | MultiShoot | 23.375 | 7.863
ILQR | 29.268 | 7.047
ADMM | 2603.063 | 36.125
|  | $b=1$ | $b=5$ | $b=25$ | $b=1$ | $b=5$ | $b=25$
FOTD (sparse LU) | $\mu=1$ | 13.324 | 0.878 | 0.0979 | 0.474 | 0.455 | 0.881
$\mu=25$ | 4.828 | 1.618 | 0.0980 | 0.460 | 0.449 | 0.885
$\mu=125$ | 7.210 | 0.328 | 0.0977 | 0.459 | 0.453 | 0.881
FOTD (GMRES) | $\mu=1$ | 4.686 | 0.672 | 0.549 | 0.566 | 0.537 | 1.019
$\mu=25$ | 12.841 | 0.897 | 0.132 | 0.597 | 0.559 | 1.026
$\mu=125$ | 4.997 | 0.251 | 0.451 | 0.621 | 0.604 | 1.041
FOTD (IDR) | $\mu=1$ | 4.051 | 0.275 | 0.106 | 0.619 | 0.468 | 0.875
$\mu=25$ | 2.688 | 0.832 | 0.0982 | 0.500 | 0.467 | 0.896
$\mu=125$ | 5.206 | 0.380 | 0.0978 | 0.536 | 0.491 | 0.878
Schwarz | $\mu=1$ | 0.448 | 0.298 | 0.0418 | 1.969 | 2.197 | 2.009
$\mu=25$ | 1.451 | 0.454 | 0.0418 | 2.116 | 2.199 | 2.027
$\mu=125$ | 3.467 | 0.521 | 0.0422 | 2.173 | 2.185 | 2.160
Table 3: The KKT residual and running time for all methods implemented on Case 2. For the Schwarz and FOTD, the smallest KKT residual and least running time among different setups of $b$ are highlighted for each setup of $\mu$. The dash “-” means the scheme does not trigger (33) within the budget. Type | Method | KKT residual $(10^{-7})$ | Time (sec.)
---|---|---|---
Centralized | IPOPT | 31.227 | 1.075
Decomposed | MultiShoot | 15.998 | 10.126
ILQR | 26.685 | 9.334
ADMM | 11841.583 | 35.133
|  | $b=1$ | $b=5$ | $b=25$ | $b=1$ | $b=5$ | $b=25$
FOTD (sparse LU) | $\mu=1$ | 607.625 | 117.596 | 1.108 | 0.708 | 0.599 | 0.863
$\mu=25$ | 364.783 | 106.214 | 0.848 | 0.699 | 0.601 | 0.906
$\mu=125$ | 240.339 | 0.268 | 4.753 | 0.735 | 0.557 | 1.026
FOTD (GMRES) | $\mu=1$ | 798.660 | 115.055 | 56.258 | 1.224 | 1.084 | 1.503
$\mu=25$ | 193.723 | 5.581 | 12.745 | 1.265 | 1.077 | 1.652
$\mu=125$ | 738.168 | 62.946 | 33.445 | 1.310 | 1.029 | 1.520
FOTD (IDR) | $\mu=1$ | 662.446 | 37.652 | 8.290 | 0.828 | 0.695 | 1.098
$\mu=25$ | 310.175 | 48.570 | 7.662 | 0.848 | 0.724 | 1.203
$\mu=125$ | 160.98 | 9.470 | 15.706 | 0.900 | 0.754 | 1.245
Schwarz | $\mu=1$ | 11.675 | 14.911 | 1.916 | 2.427 | 2.316 | 2.055
$\mu=25$ | 14.402 | 14.770 | 1.203 | 2.852 | 2.335 | 2.057
$\mu=125$ | - | 22.714 | 9.227 | - | 2.229 | 2.155
Table 4: The KKT residual and running time for all methods implemented on Case 3. For the Schwarz and FOTD, the smallest KKT residual and least running time among different setups of $b$ are highlighted for each setup of $\mu$.  Type | Method | KKT residual $(10^{-7})$ | Time (sec.)
---|---|---|---
Centralized | IPOPT | 285.911 | 1.405
Decomposed | MultiShoot | 5.652 | 18.571
ILQR | 19.983 | 15.732
ADMM | 28467.487 | 38.571
|  | $b=1$ | $b=5$ | $b=25$ | $b=1$ | $b=5$ | $b=25$
FOTD (sparse LU) | $\mu=1$ | 69.097 | 159.907 | 1.471 | 1.067 | 1.048 | 1.480
$\mu=25$ | 23.867 | 6.988 | 0.714 | 1.161 | 0.942 | 1.480
$\mu=125$ | 4.944 | 23.393 | 1.906 | 1.066 | 1.038 | 1.513
FOTD (GMRES) | $\mu=1$ | 4.293 | 28.472 | 0.0544 | 1.382 | 1.373 | 1.903
$\mu=25$ | 47.298 | 12.676 | 0.624 | 1.490 | 1.284 | 1.877
$\mu=125$ | 5.838 | 11.084 | 17.412 | 1.559 | 1.355 | 1.897
FOTD (IDR) | $\mu=1$ | 144.449 | 5.971 | 3.612 | 1.110 | 0.989 | 1.571
$\mu=25$ | 40.627 | 0.735 | 0.532 | 1.147 | 1.020 | 1.557
$\mu=125$ | 31.889 | 29.349 | 0.0760 | 1.140 | 1.062 | 1.574
Schwarz | $\mu=1$ | 2.363 | 1.691 | 0.0154 | 4.914 | 3.994 | 2.686
$\mu=25$ | 1.383 | 0.134 | 0.0166 | 4.441 | 3.485 | 2.711
$\mu=125$ | 3.127 | 0.182 | 0.0155 | 4.937 | 3.536 | 2.715
Table 5: The KKT residual and running time for all methods implemented on the temperature control problem. For the Schwarz and FOTD, the smallest KKT residual and least running time among different setups of $b$ are highlighted for each setup of $\mu$. For ADMM, the best results among different $\mu$ are highlighted. Type | Method |  KKT residual $(10^{-7})$ | Time (sec.)
---|---|---|---
Centralized | IPOPT | 0.889 | 16.934
Decomposed | MultiShoot | 3937.071 | 104.021
ILQR | 2941.207 | 89.655
IDDP | 2619.498 | 93.832
ADMM | $\mu=1$ | $\mu=25$ | $\mu=125$ | $\mu=1$ | $\mu=25$ | $\mu=125$
1680.350 | 11684.271 | 61455.896 | 76.732 | 224.930 | 479.194
|  | $b=1$ | $b=5$ | $b=25$ | $b=1$ | $b=5$ | $b=25$
FOTD (sparse LU) | $\mu=1$ | 13.621 | 13.488 | 8.050 | 16.914 | 15.134 | 25.531
$\mu=25$ | 28.266 | 8.024 | 5.316 | 19.973 | 18.023 | 28.752
$\mu=125$ | 11.448 | 3.112 | 2.712 | 19.416 | 18.356 | 28.841
FOTD (GMRES) | $\mu=1$ | 17.837 | 15.228 | 2.368 | 16.653 | 16.729 | 21.293
$\mu=25$ | 61.917 | 6.623 | 0.573 | 21.023 | 17.350 | 22.830
$\mu=125$ | 15.052 | 2.937 | 1.704 | 23.336 | 17.871 | 23.571
FOTD (IDR) | $\mu=1$ | 14.934 | 13.515 | 3.998 | 17.106 | 17.425 | 21.393
$\mu=25$ | 33.351 | 6.954 | 8.766 | 18.459 | 16.515 | 21.888
$\mu=125$ | 28.830 | 8.553 | 7.119 | 20.244 | 17.318 | 22.750
Schwarz | $\mu=1$ | 10.995 | 8.750 | 1.303 | 23.873 | 21.155 | 21.724
$\mu=25$ | 12.976 | 4.260 | 2.991 | 24.503 | 22.410 | 23.306
| $\mu=125$ | 5.124 | 3.501 | 1.996 | 23.506 | 20.431 | 22.359

Third, we report the KKT residual and running time for all methods. We average
the results over the convergent runs among five runs (corresponding to five
initializations). Since ADMM converges within the budget for $\mu=1$ only, we
report its results under this setup. The results for Cases 1, 2, and 3 are
summarized in Tables 2, 3, and 4, respectively. From the tables, we have the
following observations. (i) The proposed FOTD and Schwarz outperform other
parallel methods such as the multiple shooting, ILQR, and ADMM, among which
ADMM has the worst performance. We believe the reasons are two folds. First,
the Schwarz [38] and FOTD employ an overlapping decomposition. The overlaps
facilitate the information exchange between subproblems, and effectively
suppress the system perturbations brought by the horizon truncation and
different choices of $\mu$. Second, as also observed in [38, Figure 6], the
ADMM iterates often generate a small stepsize, which is less effective than
the stepsize that is selected by the line search. As the augmented Lagrangian
method, ADMM also suffers when it is initialized with a poor penalty parameter
and/or poor Lagrange multipliers, especially for nonconvex problems [14]. Note
that our objective coefficient for the control variables in (32a) is $-C_{2}$
with $C_{2}>0$; thus, (32) is nonconvex even if we have a linear system in the
constraints. (ii) With three different (but efficient) QP solvers, FOTD
performs equally well, although GMRES has slightly longer running time (within
0.3 sec.) than the other two solvers. Thus, different linear system solvers
can be employed in FOTD to accelerate its computation. (iii) Between Schwarz
and FOTD, the Schwarz tends to attain a smaller KKT residual than FOTD, which
is more significant when the overlap size $b$ is as small as $1$. For the
majority of cases, both methods attain the smallest KKT residual when $b$ is
as large as 25, while attain the largest KKT residual when $b$ is as small as
1. As for the running time, FOTD consistently converges faster than the
Schwarz for different choices of $\mu$ and $b$. The running time of FOTD is
comparable to that of the centralized solver IPOPT. For the three choices of
$b$, the Schwarz tends to converge faster for a large $b$ than for moderate or
small $b$. This is because that the Schwarz performs less iterations when $b$
is large even if solving each subproblem is also more expensive. On the
contrary, FOTD converges faster for a moderate $b$, which reveals the trade-
off between the total number of iterations and the computation cost of a
single iteration. Overall, our experiments demonstrate that both FOTD and
Schwarz are superior parallel methods and robust to the parameter $\mu$. FOTD
is more efficient than the Schwarz, and is as efficient as the popular solver
IPOPT. However, as the parallel method, FOTD (and Schwarz) offers more
flexibility to the computing environments and can be applied when a single
high-speed processor is not accessible.

Thin plate temperature control: We apply FOTD on a thin plate temperature
control problem studied in [38]. We refer to [33, 8, 9] for the physical
background of the problem. In particular, we let $k\in[0,1]$ be the continuous
time index and ${\boldsymbol{\omega}}\in\Omega=[0,1]\times[0,1]$ be the
position in the domain $\Omega$. Then, we consider the following problem

$\displaystyle\min_{x,u}\;\;$
$\displaystyle\int_{0}^{1}\int_{{\boldsymbol{\omega}}\in\Omega}\left\\{({\boldsymbol{x}}({\boldsymbol{\omega}},k)-d({\boldsymbol{\omega}},k))^{2}+u({\boldsymbol{\omega}},k)^{2}\right\\}d{\boldsymbol{\omega}}\;dk,$
(34a) s.t. $\displaystyle\frac{\partial x({\boldsymbol{\omega}},k)}{\partial
k}=\nabla^{2}x({\boldsymbol{\omega}},k)+u({\boldsymbol{\omega}},k)+\frac{2h_{c}}{\kappa_{c}t_{c}}(T_{c}-x({\boldsymbol{\omega}},k))+\frac{2\epsilon_{c}\sigma_{c}}{\kappa_{c}t_{c}}(T_{c}^{4}-x({\boldsymbol{\omega}},k)^{4}),$
(34b) $\displaystyle\hskip 284.52756pt\forall
k\in[0,1],\forall{\boldsymbol{\omega}}\in\Omega,$ $\displaystyle
x({\boldsymbol{\omega}},k)=0,\quad\forall({\boldsymbol{\omega}},k)\in\Omega\times\\{0\\}\text{
or }\partial\Omega\times[0,1],$ (34c)

where $\nabla^{2}$ is the Laplace operator, that is,
$\nabla^{2}x=\frac{\partial^{2}x}{\partial{\boldsymbol{\omega}}_{1}^{2}}+\frac{\partial^{2}x}{\partial{\boldsymbol{\omega}}_{2}^{2}}$,
and $d({\boldsymbol{\omega}},k)$ in the objective (34a) is the prespecified
desired temperature. Here, the PDE constraints in (34b) are governed by a
controlled heat equation (i.e., the first two terms) with extra convection and
radiation terms (i.e., the third and fourth terms). All the symbols with a
subscript “c” are the prespecified constants. In particular, $h_{c}$ is the
convection coefficient; $\kappa_{c}$ is the thermal conductivity;
$\epsilon_{c}$ is the emissivity coefficient; $\sigma_{c}$ is the Stefan-
Boltzmann constant; $T_{c}$ is the ambient temperature; and $t_{c}$ is the
plate thickness. We unify the coefficients of heat equation for simplicity.
The initial and boundary conditions are in (34c).

In our implementation, we discretize $\Omega$ by a $4\times 4$ mesh grid, i.e.
$2\times 2$ mesh grid in the interior. The temporal horizon is decomposed by
5000 evenly spaced knots with $50$ knots for each subproblem. The setups of
all methods are as before, and we follow [33] to set up the problem
parameters. In particular, we let $h_{c}=1$, $\kappa_{c}=400$,
$\epsilon_{c}=0.5$, $\sigma_{c}=5.67\times 10^{-8}$, $T_{c}=300$,
$t_{c}=0.01$, and let $d({\boldsymbol{\omega}},k)=\sin(k)$. The KKT residual
and running time of the methods are summarized in Table 5. From the table, we
again observe that the Schwarz and FOTD outperform other four parallel
methods. The Schwarz attains smaller KKT residual than FOTD, while FOTD
converges faster than the Schwarz. Overall, our experiment shows the
superiority of the overlapping decomposition-based methods.

##  8 Conclusion

This paper proposes a fast overlapping temporal decomposition (FOTD) procedure
for solving long-horizon NLDPs in (1). FOTD relies on the sequential quadratic
programming (SQP) and incorporates SQP with OTD technique. We establish global
convergence and uniform, local linear convergence for FOTD. The local result
matches [38], while FOTD requires fewer computations in each iteration (cf.
Theorem 6.3).

Considering the improvement of the performance of the Schwarz scheme over ADMM
[38], and the relation between FOTD and the Schwarz, we believe the extension
of FOTD is worth studying. One of the drawbacks of FOTD is the separation of
modifying the Hessian matrix and solving the linear-quadratic subproblems
(steps 1 and 2 in section 3), which leads to two separate factorizations for
the subproblem matrices. Such a drawback does not enlarge the flops order of a
single processor, but indeed results in a suboptimal flop multiplier. A more
desirable algorithm should perform the Hessian modification and solve the
subproblem in a single machine jointly, with one factorization, such as
parallel quasi-Newton scheme [11]. Further, we can embed OTD into more
advanced SQP frameworks, such as the trust region-SQP. We can also replace the
line search step by the filter step, and study the behavior of the approximate
direction obtained by OTD on the filter step.

In addition, FOTD can be applied on graph-structured problems, seeing that a
similar exponential decay of sensitivity for graph-structured problems was
established in [51]. Finally, a reasonable conjecture for the improved
performance of the Schwarz scheme and FOTD over ADMM (and other parallel
methods) is the lack of information exchange among subproblems in ADMM. No
overlaps are adopted in ADMM. Thus, whether we can embed OTD into ADMM to
improve the performance of ADMM, and whether the OTD-based ADMM exhibits a
similar convergence rate as FOTD are interesting future research directions.

## Acknowledgments.

This material was based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research (ASCR)
under Contract DE-AC02-06CH11347 and by NSF through award CNS-1545046.

##  9 Proofs of results in section 2

###  9.1 Proof of Theorem 2.2

We start by writing the KKT conditions of Problem (1) and Problem (3). For
$k\in[N-1]$, we let
$A_{k}({\boldsymbol{z}}_{k})=\nabla_{{\boldsymbol{x}}_{k}}^{T}f_{k}({\boldsymbol{z}}_{k})$,
$B_{k}({\boldsymbol{z}}_{k})=\nabla_{{\boldsymbol{u}}_{k}}^{T}f_{k}({\boldsymbol{z}}_{k})$
be the Jacobian matrices. Then the Lagrange function of (1) is

$\mathcal{L}({\boldsymbol{z}},{\boldsymbol{\lambda}})=\sum_{k=0}^{N-1}\\{g_{k}({\boldsymbol{z}}_{k})+{\boldsymbol{\lambda}}_{k}^{T}{\boldsymbol{x}}_{k}-{\boldsymbol{\lambda}}_{k+1}^{T}f_{k}({\boldsymbol{z}}_{k})\\}+\\{g_{N}({\boldsymbol{x}}_{N})+{\boldsymbol{\lambda}}_{N}^{T}{\boldsymbol{x}}_{N}\\}-{\boldsymbol{\lambda}}_{0}^{T}\bar{\boldsymbol{x}}_{0}.$
(35)

Thus, the KKT conditions of (1) are

$\displaystyle\nabla_{{\boldsymbol{x}}_{k}}g_{k}({\boldsymbol{z}}_{k})+{\boldsymbol{\lambda}}_{k}-A_{k}^{T}({\boldsymbol{z}}_{k}){\boldsymbol{\lambda}}_{k+1}=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[N-1],$ (36a)
$\displaystyle\nabla_{{\boldsymbol{u}}_{k}}g_{k}({\boldsymbol{z}}_{k})-B_{k}^{T}({\boldsymbol{z}}_{k}){\boldsymbol{\lambda}}_{k+1}=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[N-1],$ (36b)
$\displaystyle\nabla_{{\boldsymbol{x}}_{N}}g_{N}({\boldsymbol{z}}_{N})+{\boldsymbol{\lambda}}_{N}=\;$
$\displaystyle{\boldsymbol{0}},$ (36c)
$\displaystyle{\boldsymbol{x}}_{k+1}-f_{k}({\boldsymbol{z}}_{k})=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[N-1],$ (36d)
$\displaystyle{\boldsymbol{x}}_{0}-\bar{\boldsymbol{x}}_{0}=\;$
$\displaystyle{\boldsymbol{0}}.$ (36e)

Similarly, the KKT conditions of
${\mathcal{P}}^{i}_{\mu}({\boldsymbol{d}}_{i})$ in (3) are

$\displaystyle\nabla_{{\boldsymbol{x}}_{k}}g_{k}({\tilde{\boldsymbol{z}}}_{i,k})+{\tilde{{\boldsymbol{\lambda}}}}_{i,k}-A_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i,k+1}=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[m_{1},m_{2}),$ (37a)
$\displaystyle\nabla_{{\boldsymbol{u}}_{k}}g_{k}({\tilde{\boldsymbol{z}}}_{i,k})-B_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i,k+1}=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[m_{1},m_{2}),$ (37b)
$\displaystyle\nabla_{{\boldsymbol{x}}_{m_{2}}}g_{m_{2}}({\tilde{\boldsymbol{x}}}_{i,m_{2}},\bar{\boldsymbol{u}}_{m_{2}})+{\tilde{{\boldsymbol{\lambda}}}}_{i,m_{2}}-A_{m_{2}}^{T}({\tilde{\boldsymbol{x}}}_{i,m_{2}},\bar{\boldsymbol{u}}_{m_{2}})\bar{\boldsymbol{\lambda}}_{m_{2}+1}$
$\displaystyle+\mu({\tilde{\boldsymbol{x}}}_{i,m_{2}}-\bar{\boldsymbol{x}}_{m_{2}})=\;$
$\displaystyle{\boldsymbol{0}},$ (37c)
$\displaystyle{\tilde{\boldsymbol{x}}}_{i,k+1}-f_{k}({\tilde{\boldsymbol{z}}}_{i,k})=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[m_{1},m_{2}),$ (37d)
$\displaystyle{\tilde{\boldsymbol{x}}}_{i,m_{1}}-\bar{\boldsymbol{x}}_{m_{1}}=\;$
$\displaystyle{\boldsymbol{0}}.$ (37e)

(i). By Definition 2.1,
$\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})=({\boldsymbol{x}}^{\star}_{m_{1}:m_{2}},{\boldsymbol{u}}^{\star}_{m_{1}:m_{2}-1},{\boldsymbol{\lambda}}^{\star}_{m_{1}:m_{2}})$.
Letting ${\boldsymbol{d}}_{i}={\boldsymbol{d}}_{i}^{\star}$ in (37c), (37e),
we see (37) is a subsystem of (36) so that
$\mathcal{D}_{i}({\boldsymbol{z}}^{\star},{\boldsymbol{\lambda}}^{\star})$
satisfies (37). Thus, the statement holds.

(ii). Suppose
$({\tilde{\boldsymbol{z}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i})=({\tilde{\boldsymbol{z}}}_{i}^{\star}({\boldsymbol{d}}_{i}),{\tilde{{\boldsymbol{\lambda}}}}_{i}^{\star}({\boldsymbol{d}}_{i}))$
satisfies conditions (37). Since $[n_{i},n_{i+1})\subseteq[m_{1},m_{2})$, we
know from (37a), (37b), (37d) that
$({\tilde{\boldsymbol{z}}}_{i,n_{i}:n_{i+1}-1},{\tilde{{\boldsymbol{\lambda}}}}_{i,n_{i}:n_{i+1}-1})$
satisfies

$\displaystyle\nabla_{{\boldsymbol{x}}_{k}}g_{k}({\tilde{\boldsymbol{z}}}_{i,k})+{\tilde{{\boldsymbol{\lambda}}}}_{i,k}-A_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i,k+1}=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[n_{i},n_{i+1}),$ (38a)
$\displaystyle\nabla_{{\boldsymbol{u}}_{k}}g_{k}({\tilde{\boldsymbol{z}}}_{i,k})-B_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i,k+1}=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[n_{i},n_{i+1}),$ (38b)
$\displaystyle{\tilde{\boldsymbol{x}}}_{i,k+1}-f_{k}({\tilde{\boldsymbol{z}}}_{i,k})=\;$
$\displaystyle{\boldsymbol{0}},\quad\forall k\in[n_{i},n_{i+1}),$ (38c)

which is a subset of (36a), (36b), (36d). We consider the composed point
$({\boldsymbol{z}},{\boldsymbol{\lambda}})=\mathcal{C}(\\{({\tilde{\boldsymbol{z}}}_{i},{\tilde{{\boldsymbol{\lambda}}}}_{i})\\}_{i})$.
By Definition 2.1, we know

$\displaystyle{\boldsymbol{z}}_{k}={\tilde{\boldsymbol{z}}}_{i,k},\quad\text{if
}k\in[n_{i},n_{i+1})\text{ for some
}i\in[M-1],\quad\quad{\boldsymbol{z}}_{N}={\tilde{\boldsymbol{x}}}_{M-1,N},$
(39)
$\displaystyle{\boldsymbol{\lambda}}_{k}={\tilde{{\boldsymbol{\lambda}}}}_{i,k},\quad\text{if
}k\in[n_{i},n_{i+1})\text{ for some
}i\in[M-1],\quad\quad{\boldsymbol{\lambda}}_{N}={\tilde{{\boldsymbol{\lambda}}}}_{M-1,N}.$

Thus, (36c) is implied by (37c) for $i=M-1$, and (36e) is implied by (37e) for
$i=0$. For (36b), we use the decomposition
$[N-1]=\cup_{i=0}^{M-1}[n_{i},n_{i+1})$. For any $k\in[N-1]$, we have two
cases.

(a). $k\neq n_{i+1}-1$, $\forall i\in[M-2]$. Then, $k\in[n_{i},n_{i+1}-1)$ for
some $i\in[M-2]$, or $k\in[n_{M-1},n_{M})$. Thus, $k+1\in[n_{i}+1,n_{i+1})$
for $i\in[M-2]$ or $k+1\in[n_{M-1}+1,N]$. By (39), we know for both cases that
${\boldsymbol{z}}_{k}$ and ${\boldsymbol{\lambda}}_{k+1}$ are from the same
subproblem, i.e. ${\boldsymbol{z}}_{k}={\tilde{\boldsymbol{z}}}_{i,k}$ and
${\boldsymbol{\lambda}}_{k+1}={\tilde{{\boldsymbol{\lambda}}}}_{i,k+1}$. Thus,
(36b) is implied by (38b).

(b). $k=n_{i+1}-1$ for $i\in[M-2]$. Then, $k+1=n_{i+1}$ and thus
${\boldsymbol{\lambda}}_{k+1}={\tilde{{\boldsymbol{\lambda}}}}_{i+1,k+1}$.
Comparing (36b) with (38b), we know that

$\displaystyle\nabla_{{\boldsymbol{u}}_{k}}g_{k}({\boldsymbol{z}}_{k})-B_{k}^{T}({\boldsymbol{z}}_{k}){\boldsymbol{\lambda}}_{k+1}={\boldsymbol{0}}$
$\displaystyle\Longleftrightarrow\nabla_{{\boldsymbol{u}}_{k}}g_{k}({\tilde{\boldsymbol{z}}}_{i,k})-B_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i+1,k+1}={\boldsymbol{0}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{pequ:1b}}}{{\Longleftrightarrow}}B_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i+1,k+1}=B_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i,k+1}.$
(40)

Following the same derivation, we consider (36a) and can easily get

$\nabla_{{\boldsymbol{x}}_{k}}g_{k}({\boldsymbol{z}}_{k})+{\boldsymbol{\lambda}}_{k}-A_{k}^{T}({\boldsymbol{z}}_{k}){\boldsymbol{\lambda}}_{k+1}={\boldsymbol{0}}\stackrel{{\scriptstyle\eqref{pequ:1a}}}{{\Longleftrightarrow}}A_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i+1,k+1}=A_{k}^{T}({\tilde{\boldsymbol{z}}}_{i,k}){\tilde{{\boldsymbol{\lambda}}}}_{i,k+1},$
(41)

if $k=n_{i+1}-1$ for $i\in[M-2]$. Finally, we consider (36d), where we have
two cases.

(a). $k\neq n_{i+1}-1$, $\forall i\in[M-2]$. As before, ${\boldsymbol{z}}_{k}$
and ${\boldsymbol{x}}_{k+1}$ are from the same subproblem, so that (36d) is
implied by (38c).

(b). $k=n_{i+1}-1$ for $i\in[M-2]$. Then,
${\boldsymbol{x}}_{k+1}={\tilde{\boldsymbol{x}}}_{i+1,k+1}$. Comparing (36d)
with (38c), we know

${\boldsymbol{x}}_{k+1}=f_{k}({\boldsymbol{z}}_{k})\Longleftrightarrow{\tilde{\boldsymbol{x}}}_{i+1,k+1}=f_{k}({\tilde{\boldsymbol{z}}}_{i,k})\stackrel{{\scriptstyle\eqref{pequ:1c}}}{{\Longleftrightarrow}}{\tilde{\boldsymbol{x}}}_{i+1,k+1}={\tilde{\boldsymbol{x}}}_{i,k+1}.$
(42)

Combining (9.1), (41), and (42), we complete the proof.

##  10 Proofs of results in section 4

###  10.1 Proof of Lemma 4.1

We suppress the iteration index $\tau$ since the result holds uniformly over
$\tau$. By Assumption 4, we know that for any integer $k$, there exists an
integer $t_{k}\in[1,t]$ such that
$\Xi_{k,t_{k}}\Xi_{k,t_{k}}^{T}\succeq\gamma_{C}I$. Let us define the knots
recursively by $k_{j+1}=k_{j}+t_{k_{j}}$, $j=0,1,\ldots,J-1$ with $k_{0}=0$.
We suppose $J$ is large enough such that $k_{J}>N-t$. Thus, we derive a
horizon decomposition
$[N-1]=\cup_{j=0}^{J-1}[k_{j},k_{j+1}-1]\cup[k_{J},N-1]$. Since
$G=\nabla^{T}_{{\boldsymbol{z}}}f$, we can apply the definition of $f$ in (7)
and write $GG^{T}$ explicitly. We have

$GG^{T}=\left(\begin{smallmatrix}I&-A_{0}^{T}\\\
-A_{0}&\begin{subarray}{c}I+A_{0}A_{0}^{T}+B_{0}B_{0}^{T}\end{subarray}&\ddots\\\
&\ddots&\ddots&-A_{N-1}^{T}\\\
&&-A_{N-1}&\begin{subarray}{c}I+A_{N-1}A_{N-1}^{T}+B_{N-1}B_{N-1}^{T}\end{subarray}\end{smallmatrix}\right).$
(43)

We let $k_{j}^{\prime}=k_{j+1}-1=k_{j}+t_{k_{j}}-1$, and define matrices
$\\{{\mathcal{T}}_{j}\\}_{j\in[J]}$ corresponding to each interval
$[k_{j},k_{j}^{\prime}]$ as (${\mathcal{T}}_{J}$ is defined similarly, except
that the bottom corner block is $I+A_{N-1}A_{N-1}^{T}+B_{N-1}B_{N-1}^{T}$)

${\mathcal{T}}_{j}=\left(\begin{smallmatrix}I&-A_{k_{j}}^{T}\\\
-A_{k_{j}}&\begin{subarray}{c}I+A_{k_{j}}A_{k_{j}}^{T}+B_{k_{j}}B_{k_{j}}^{T}\end{subarray}&\ddots\\\
&\ddots&\ddots&-A_{k_{j}^{\prime}}^{T}\\\
&&-A_{k_{j}^{\prime}}&\begin{subarray}{c}A_{k_{j}^{\prime}}A_{k_{j}^{\prime}}^{T}+B_{k_{j}^{\prime}}B_{k_{j}^{\prime}}^{T}\end{subarray}\end{smallmatrix}\right),\quad\forall
j\in[J-1].$

By the expression in (43), it suffices to show that each ${\mathcal{T}}_{j}$
is lower bounded away from zero. Then, we have
$\lambda_{\min}(GG^{T})\geq\min_{j\in[J]}\lambda_{\min}({\mathcal{T}}_{j})$.
Let us consider ${\mathcal{T}}_{J}$ first. We let
$A_{m}^{n}=A_{m}A_{m-1}\cdots A_{n}$ and define a matrix
${\mathcal{T}}_{J}^{1/2}$ as

$\displaystyle{\mathcal{T}}_{J}^{1/2}\coloneqq$
$\displaystyle\left(\begin{smallmatrix}I\\\ -A_{k_{J}}&I&&&&-B_{k_{J}}\\\
&\ddots&\ddots&&&&\ddots\\\ &&&I&&&&-B_{N-2}\\\
&&&-A_{N-1}&I&&&&-B_{N-1}\end{smallmatrix}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{smallmatrix}I\\\ -A_{k_{J}}&I\\\ &&\ddots\\\
&&&I\\\
&&&&I\end{smallmatrix}\right)\times\cdots\times\left(\begin{smallmatrix}I\\\
&I\\\ &&\ddots\\\ &&&I\\\
&&&-A_{N-1}&I\end{smallmatrix}\right)\times\left(\begin{smallmatrix}I\\\
&I&&&&-B_{k_{J}}\\\ &&\ddots&&&\vdots&\ddots\\\
&&&I&&-A_{N-2}^{k_{J}+1}B_{k_{J}}&&-B_{N-2}\\\
&&&&I&-A_{N-1}^{k_{J}+1}B_{k_{J}}&&-A_{N-1}B_{N-2}&-B_{N-1}\end{smallmatrix}\right)$
$\displaystyle\eqqcolon$ $\displaystyle P_{k_{J}}\cdots P_{N-1}Q_{J}.$

Then, we use the fact that $Q_{J}Q_{J}^{T}\succeq I$ and have

$\displaystyle{\mathcal{T}}_{J}={\mathcal{T}}_{J}^{1/2}({\mathcal{T}}_{J}^{1/2})^{T}=P_{k_{J}}\cdots
P_{N-1}Q_{J}Q_{J}^{T}P_{N-1}^{T}\cdots P_{k_{J}}^{T}\succeq P_{k_{J}}\cdots
P_{N-1}P_{N-1}^{T}\cdots P_{k_{J}}^{T}.$

Since $\|A_{k}\|\leq\Upsilon_{upper}$ for $k\in[k_{J},N-1]$, we know
$\|P_{k}^{-1}\|\leq 1+\Upsilon_{upper}$. Thus, $P_{k}P_{k}^{T}\succeq
1/(1+\Upsilon_{upper})^{2}I$. By the above display and using $N-k_{J}\leq t$,
we have $\lambda_{\min}({\mathcal{T}}_{J})\succeq
1/(1+\Upsilon_{upper})^{2t}$. We then consider ${\mathcal{T}}_{j}$ for
$j\in[J-1]$ similarly. By the same matrix multiplication, we have

${\mathcal{T}}_{j}^{1/2}\coloneqq\left(\begin{smallmatrix}I\\\
-A_{k_{j}}&I&&&-B_{k_{j}}\\\ &\ddots&\ddots&&&\ddots\\\
&&&I&&&-B_{k_{j}^{\prime}-1}\\\
&&&-A_{k_{j}^{\prime}}&&&&-B_{k_{j}^{\prime}}\end{smallmatrix}\right)=P_{k_{j}}\cdots
P_{k_{j}^{\prime}}Q_{j},$

where (slightly different from $Q_{J}$)

$Q_{j}=\left(\begin{smallmatrix}I\\\ &I&&&-B_{k_{j}}\\\
&&\ddots&&\vdots&\ddots\\\
&&&I&-A_{k_{j}^{\prime}-1}^{k_{j}+1}B_{k_{j}}&&-B_{k_{j}^{\prime}-1}\\\
&&&&-A_{k_{j}^{\prime}}^{k_{j}+1}B_{k_{j}}&&-A_{k_{j}^{\prime}}B_{k_{j}^{\prime}-1}&-B_{k_{j}^{\prime}}\end{smallmatrix}\right)\eqqcolon\begin{pmatrix}I&Q_{j,1}\\\
{\boldsymbol{0}}&Q_{j,2}\end{pmatrix}.$

Since $Q_{j,2}$ is just the controllability matrix $\Xi_{k_{j},t_{k_{j}}}$
(except that the components are in a reverse order), Assumption 4 implies that
$Q_{j,2}Q_{j,2}^{T}\succeq\gamma_{C}I$. Furthermore, since
$\max\\{\|A_{k}\|,\;\|B_{k}\|\\}\leq\Upsilon_{upper}$, $\forall
k\in[k_{j},k_{j}^{\prime}-1]$, we get

$\max\\{\|Q_{j,2}\|,\;\|Q_{j,1}\|\\}\leq\Upsilon_{upper}+\Upsilon_{upper}^{2}+\cdots\Upsilon_{upper}^{k_{j}^{\prime}-k_{j}+1}\leq\sum_{i=1}^{t}\Upsilon_{upper}^{i}\leq\frac{\Upsilon_{upper}^{t+1}}{\Upsilon_{upper}-1},$

where the second inequality is due to $t_{k_{j}}\leq t$. Finally, noting that

$Q_{j}Q_{j}^{T}=\begin{pmatrix}I+Q_{j,1}Q_{j,1}^{T}&Q_{j,1}Q_{j,2}^{T}\\\
Q_{j,2}Q_{j,1}^{T}&Q_{j,2}Q_{j,2}^{T}\end{pmatrix},$

we apply the algebra result in [34, Lemma 4.8(ii)] and obtain

$\displaystyle
Q_{j}Q_{j}^{T}\succeq\left(\frac{\gamma_{C}}{\gamma_{C}+\frac{\Upsilon_{upper}^{t+1}}{\Upsilon_{upper}-1}}\right)^{2}\cdot\min\\{1,\;\gamma_{C}\\}\cdot
I\eqqcolon\gamma_{Q}I.$

Thus, using $P_{k}P_{k}^{T}\succeq 1/(1+\Upsilon_{upper})^{2}I$, we obtain

$\displaystyle{\mathcal{T}}_{j}={\mathcal{T}}_{j}^{1/2}({\mathcal{T}}_{j}^{1/2})^{T}=P_{k_{j}}\cdots
P_{k_{j}^{\prime}}Q_{j}Q_{j}^{T}P_{k_{j}^{\prime}}^{T}\cdots
P_{k_{j}}^{T}\succeq\frac{\gamma_{Q}}{(1+\Upsilon_{upper})^{2t}}I.$

Since $\gamma_{Q}<1$, we let $\gamma_{G}=\gamma_{Q}/(1+\Upsilon_{upper})^{2t}$
and have $\lambda_{\min}({\mathcal{T}}_{j})\succeq\gamma_{G}$, $\forall
j\in[J]$. This finishes the proof.

###  10.2 Proof of Lemma 4.3

We adapt the proof of [35, Lemma 1 and Theorem 1]. We suppress the iteration
index $\tau$. The reduced Hessian of ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$
is independent from ${\boldsymbol{d}}_{i}$, which only affects linear terms.
From (15b)-(15c), we know that the Jacobian matrix of the constraints in
${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$ has full row rank. Thus, it suffices
to show that the reduced Hessian of ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$
is lower bounded by $\gamma_{RH}I$. Let

${\hat{H}}^{(i)}=\operatorname{diag}({\hat{H}}_{m_{1}},{\hat{H}}_{m_{1}+1},\ldots,{\hat{H}}_{m_{2}-1},{\hat{Q}}_{m_{2}}+\mu
I)$

be the Hessian of ${{\mathcal{L}}{\mathcal{P}}}_{\mu}^{i}$. We only need to
show that
${\tilde{\boldsymbol{w}}}_{i}^{T}{\hat{H}}^{(i)}{\tilde{\boldsymbol{w}}}_{i}\geq\gamma_{RH}\|{\tilde{\boldsymbol{w}}}_{i}\|^{2}$
for any
${\tilde{\boldsymbol{w}}}_{i}=({\tilde{\boldsymbol{p}}}_{i},{\tilde{\boldsymbol{q}}}_{i})\neq{\boldsymbol{0}}$
satisfying

$\displaystyle{\tilde{\boldsymbol{p}}}_{i,k+1}=$ $\displaystyle
A_{k}{\tilde{\boldsymbol{p}}}_{i,k}+B_{k}{\tilde{\boldsymbol{q}}}_{i,k},\quad
k\in[m_{1},m_{2}),$ (44a) $\displaystyle{\tilde{\boldsymbol{p}}}_{i,m_{1}}=$
$\displaystyle{\boldsymbol{0}}.$ (44b)

We take the last subproblem (i.e., $i=M-1$) as an example to illustrate the
proof idea. For $i=M-1$, we have $\mu=0$. For any
${\tilde{\boldsymbol{w}}}_{M-1}=({\tilde{\boldsymbol{p}}}_{M-1},{\tilde{\boldsymbol{q}}}_{M-1})\neq{\boldsymbol{0}}$
satisfying (44), we extend ${\tilde{\boldsymbol{w}}}_{M-1}$ forward by filling
with ${\boldsymbol{0}}$ to get a full-horizon vector ${\boldsymbol{w}}$. That
is ${\boldsymbol{w}}=({\boldsymbol{0}};{\tilde{\boldsymbol{w}}}_{M-1})$. We
can verify that
${\boldsymbol{w}}=({\boldsymbol{p}},{\boldsymbol{q}})\in\\{{\boldsymbol{w}}:G{\boldsymbol{w}}={\boldsymbol{0}},{\boldsymbol{w}}\neq{\boldsymbol{0}}\\}$.
Therefore,

${\tilde{\boldsymbol{w}}}_{M-1}^{T}{\hat{H}}^{(M-1)}{\tilde{\boldsymbol{w}}}_{M-1}={\boldsymbol{w}}^{T}{\hat{H}}{\boldsymbol{w}}\geq\gamma_{RH}\|{\boldsymbol{w}}\|^{2}=\gamma_{RH}\|{\tilde{\boldsymbol{w}}}_{M-1}\|^{2},$

where the first and last equalities are due to the construction of
${\boldsymbol{w}}$, and the middle inequality is due to Assumption 4. For
$i\in[M-2]$, we consider the following two cases.

Case 1: $m_{2}\geq N-t$. Given
${\tilde{\boldsymbol{w}}}_{i}=({\tilde{\boldsymbol{p}}}_{i},{\tilde{\boldsymbol{q}}}_{i})\neq{\boldsymbol{0}}$
satisfying (44), we can still extend it forward by filling with
${\boldsymbol{0}}$. However, extending backward with ${\boldsymbol{0}}$ will
make the full vector ${\boldsymbol{w}}$ outside the space
$\\{{\boldsymbol{w}}:G{\boldsymbol{w}}={\boldsymbol{0}}\\}$. Instead, we
construct the following extension. We let
${\boldsymbol{q}}_{m_{2}:N-1}={\boldsymbol{0}}$,
${\boldsymbol{p}}_{k+1}=A_{k}{\boldsymbol{p}}_{k}$ for $k=[m_{2},N)$. Thus,

${\boldsymbol{w}}=({\boldsymbol{0}};{\tilde{\boldsymbol{w}}}_{i};{\boldsymbol{0}};{\boldsymbol{p}}_{m_{2}+1};{\boldsymbol{0}};{\boldsymbol{p}}_{m_{2}+2}\;\ldots;{\boldsymbol{0}};{\boldsymbol{p}}_{N})\in\\{{\boldsymbol{w}}:G{\boldsymbol{w}}={\boldsymbol{0}},{\boldsymbol{w}}\neq{\boldsymbol{0}}\\}.$
(45)

Moreover, by Assumption 4,

$\|{\boldsymbol{p}}_{m_{2}+1:N}\|^{2}=\sum_{k=m_{2}+1}^{N}\|{\boldsymbol{p}}_{k}\|^{2}\leq\sum_{j=1}^{N-m_{2}}\Upsilon_{upper}^{2j}\|{\boldsymbol{p}}_{m_{2}}\|^{2}\\\
=\frac{\Upsilon_{upper}^{2}(\Upsilon_{upper}^{2(N-m_{2})}-1)}{\Upsilon_{upper}^{2}-1}\|{\boldsymbol{p}}_{m_{2}}\|^{2}\leq\frac{\Upsilon_{upper}^{2t+2}-\Upsilon_{upper}^{2}}{\Upsilon_{upper}^{2}-1}\|{\boldsymbol{p}}_{m_{2}}\|^{2}.$
(46)

The last inequality uses $N-m_{2}\leq t$. Further,

$\displaystyle{\boldsymbol{w}}^{T}{\hat{H}}{\boldsymbol{w}}$