\subsection{Proof of Lemma \ref{lem:assume-guarantee}} \label{pf:lem:assume-guarantee}

\def\Pi^{0:t}{\Pi^{0:t}}

We require several  definitions. Define the projection operators
\begin{align*}
\Pi^{0:t}_{\mathcal{C}} &:= \bmat{I_{n_x^{\mathcal{C}} (t+1)} & 0_{n_x^{\mathcal{C}} (t+1) \times n_x^{\mathcal{C}} (T - t)}} \Pi_\couple
\end{align*}
for  $t = 0, .., T$.
It follows $x_{\mathcal{C}}^{0:t} = \Pi^{0:t}_{\mathcal{C}} x$ for each time $t \in \{0, .., T\}$. Additionally, we define the function $f_\phi: {\bf R}^{N_x} \times {\bf R}^{N_x^{\mathcal{C}}} \to {\bf R}^{N_x}$ as
\begin{align}
f_\phi (w, v_{\mathcal{C}} ) := {\widetilde{B}} \phi (w,   v_{\mathcal{C}}) + {\widetilde{L}} w + {\widetilde{H}}  v_{\mathcal{C}}. \label{eq:f_phi}
\end{align}
\rev{The function $f_\phi$ maps  $(w,v_{\mathcal{C}})$ to the surrogate system state $\widetilde{x}$  under the control policy $\phi$.}
Lemma \ref{lem:f_phi} shows that the function $f_\phi$ is \emph{strictly causal} in $v_{\mathcal{C}}$, and that $x = f_\phi (w, x_{\mathcal{C}})$.
We state it without proof, as it  follows from the  structure of the  surrogate state equation \eqref{eq:surr_state_i}.

\begin{lemm} \label{lem:f_phi}
Let $\phi \in \Phi$ and $x = B\phi (w, x_{\mathcal{C}}) + L w$. The following properties hold.
\begin{enumerate}[(i)]\setlength{\itemsep}{.05in}
\item Let $t \in \{0, \dots, T\}$ and  $w \in {\bf R}^{N_x}$. It holds that
$$\Pi^{0:t-1}_{\mathcal{C}} (v - v') = 0 \ \Longrightarrow \ \Pi^{0:t}_{\mathcal{C}} f_\phi (w, v_{\mathcal{C}}) = \Pi^{0:t}_{\mathcal{C}} f_\phi (w, v_{\mathcal{C}}').$$
\item The state trajectory satisfies  $x = f_\phi (w, x_{\mathcal{C}})$ for all $w \in {\bf R}^{N_x}$.
\end{enumerate}
\end{lemm}

\vspace{.03in}

Now, fix $w \in {\mathcal{W}}$ and let $x = B\phi (w, x_{\mathcal{C}}) + L w$. \rev{To complete the proof, it suffices to show that $x_{\mathcal{C}}^{0:t} \in \Pi^{0:t}_{\mathcal{C}} \Pi_\couple^\top {\mathcal{V}}_{\mathcal{C}}$ for all $t \in \{0, \dots, T\}.$  In particular, the satisfaction of this condition for $t = T$ implies that $x_{\mathcal{C}} \in \Pi_\couple \Pi_\couple^\top {\mathcal{V}}_{\mathcal{C}} = {\mathcal{V}}_{\mathcal{C}}$. We prove this by induction in $t$.}

\emph{(Base step).} \, We have that $x_{\mathcal{C}} (0) = \widetilde{x}_{\mathcal{C}} (0) \in \Pi_{\mathcal{C}}^{0:0} \Pi_\couple^\top {\mathcal{V}}_{\mathcal{C}}$, where the first equality follows from the initial condition of the surrogate state equation \eqref{eq:surr_state_i}.

\emph{(Induction step).} \, Assume that $x^{0:t-1}_{\mathcal{C}} \in \Pi^{0:t-1}_{\mathcal{C}} \Pi_\couple^\top {\mathcal{V}}_{\mathcal{C}}$. We complete the proof by showing that $x^{0:t}_{\mathcal{C}} \in \Pi^{0:t}_{\mathcal{C}} \Pi_\couple^\top {\mathcal{V}}_{\mathcal{C}}$. Fix a $v \in {\bf R}^{N_x}$ that satisfies $\Pi_\couple v \in {\mathcal{V}}_{\mathcal{C}}$ and $x^{0:t-1}_{\mathcal{C}} = \Pi_{\mathcal{C}}^{0:t-1} v$ (which is guaranteed to exist by our induction hypothesis).
We have that
\begin{align*}
x_{\mathcal{C}}^{0:t} &=  \Pi^{0:t}_{\mathcal{C}}   x =  \Pi^{0:t}_{\mathcal{C}}   f_\phi (w, x_{\mathcal{C}}) = \Pi^{0:t}_{\mathcal{C}}   f_\phi (w, v_{\mathcal{C}})\\
&= \Pi^{0:t}_{\mathcal{C}} \Pi_\couple^\top \Pi_\couple  f_\phi (w, v_{\mathcal{C}}) \in \Pi^{0:t}_{\mathcal{C}} \Pi_\couple^\top {\mathcal{V}}_{\mathcal{C}}.
\end{align*}
The second equality follows from property (ii) in Lemma \ref{lem:f_phi}; the third equality follows from a combination of property (i) in Lemma \ref{lem:f_phi} and the equality $x^{0:t-1}_{\mathcal{C}} = \Pi_{\mathcal{C}}^{0:t-1} v$; the fourth equality follows from the identity that $\Pi^{0:t}_{\mathcal{C}} \Pi_\couple^\top \Pi_\couple = \Pi^{0:t}_{\mathcal{C}}$ for each time $t$; and the final inclusion  follows from the assumption that the policy $\phi$ satisfies the assume-guarantee contract specified by the contract set ${\mathcal{V}}_{\mathcal{C}}$.

\subsection{Proof of Lemma \ref{lem:matrix_inequality}} \label{pf:lem:matrix_inequality}

It suffices to show that the matrix inequality \eqref{eq:matrix_inequality1} is satisfied if the LMI \eqref{eq:LMI} is satisfied.
Define $\beta := \alpha \lambda$, and divide both sides of the matrix inequality \eqref{eq:matrix_inequality1} by $\lambda$. It follows that the matrix inequality \eqref{eq:matrix_inequality1} is satisfied if and only if there exists $\beta \in [0, \lambda]$, such that
\begin{align}
&\nonumber \beta^{-1} \Pi_\couple {P^w} \Sigma {P^w}^\top \Pi_\couple^\top + (\lambda - \beta)^{-1} \Pi_\couple {P^\xi} \Sigma {P^\xi}^\top \Pi_\couple^\top  \\
& \qquad \preceq \Pi_\couple (\lambda \Sigma - Y \Sigma - \Sigma Y^\top + \lambda^{-1} Y \Sigma Y^\top) \Pi_\couple^\top. \label{eq:matrix_inequality2}
\end{align}
Since $\lambda^{-1} Y \Sigma Y^\top \succeq 0$, the matrix inequality \eqref{eq:matrix_inequality2} is satisfied if
\begin{align*}
    \beta^{-1} \Pi_\couple {P^w} \Sigma {P^w}^\top \Pi_\couple^\top + (\lambda - \beta)^{-1} \Pi_\couple {P^\xi} \Sigma {P^\xi}^\top \Pi_\couple^\top  \preceq \Pi_\couple \widetilde{\Sigma} \Pi_\couple^\top,
\end{align*}
where $\widetilde{\Sigma} := \lambda \Sigma - Y \Sigma - \Sigma Y^\top$. It follows from Schur's Lemma that the above matrix inequality is satisfied if and only if the LMI \eqref{eq:LMI} is satisfied. This completes the proof.

\section{Decentralized Control Design via Contracts} \label{sec:feas_control}

In this section, we construct a convex inner approximation of the decentralized control design problem \eqref{opt:decent_dFeedback}
by introducing a surrogate information structure in which the information-coupling states are modeled as \emph{fictitious disturbances}, which are  ``assumed" to take values in a given ``contract'' set.
To ``guarantee'' the satisfaction of this assumption, we impose a contractual constraint on the control policy requiring that the actual information-coupling states induced by the control policy belong to the contract set.
Given a fixed contract set, the resulting problem is a convex disturbance-feedback control design problem, whose feasible policies are shown  to be feasible for original problem \eqref{opt:decent_dFeedback}.

\subsection{Surrogate Information}
We associate a \emph{fictitious disturbance}  $v_i (t) \in {\bf R}^{n_x^i}$ with each subsystem $i \in V$ and time period $t = 0, \dots, T$. We let $v \in {\bf R}^{N_x}$ denote the corresponding fictitious disturbance trajectory induced by these individual elements, which we model as a random vector \rtwo{whose distribution we take as given in this section.}\footnote{\rtwo{In Sec. \ref{sec:joint_opt}, we treat the probability distribution of the fictitious disturbance trajectory $v$ as a design choice. In Section \ref{sec:param_fict}, we propose a  parameterization of  the fictitious disturbance trajectory that facilitates the co-optimization of its distribution with the control policy via semidefinite programming.}}  The support of the fictitious disturbance trajectory, which we denote by ${\mathcal{V}} \subset {\bf R}^{N_x}$, is assumed to be a convex and compact set. We assume that the support of the pair of random vectors $(w, v)$ is given by  set ${\mathcal{W}} \times {\mathcal{V}}$.

Treating the fictitious disturbances  as a surrogates for the information-coupling states, we define the \emph{surrogate local information} associated with subsystem $i$ as
\begin{align*}
\widetilde{\zeta}_i (t) :=  \{v_j^{0:t} | j \in {\mathcal{C}} (i) \} \cup \{w_j^{-1: t-1}  | j \in {\mathcal{N}} (i) \}.
\end{align*}
Given a decentralized control policy $\phi \in \Phi$, the surrogate local information induces a \emph{surrogate control input} for each subsystem $i$ given by
$\widetilde{u}_i (t) := \phi_i ( \widetilde{\zeta}_i (t), t).$
We denote the corresponding \emph{surrogate input trajectory}
by
\begin{align*}
\widetilde{u} := \phi (w, v_{\mathcal{C}}),
\end{align*}
where $v_{\mathcal{C}} := \Pi_\couple v$ is the projection of the fictitious disturbance $v$ onto the subspace associated with the information-coupling states.

\subsection{Surrogate Dynamics}

We replace the information-coupling states with  fictitious disturbances. This induces a \emph{surrogate system state} that evolves according to the following surrogate state equation:
\begin{align}
\nonumber \widetilde{x}_i (t+1) = & \sum_{j \in V \setminus {\mathcal{C}}(i) } A_{ij} (t) \widetilde{x}_j (t) + \sum_{j \in {\mathcal{C}}  (i)} A_{ij} (t) v_j (t)  \\
&  + \sum_{j=1}^N B_{ij} (t) \widetilde{u}_j (t) + w_i (t), \label{eq:surr_state_i}
\end{align}
where $\widetilde{x}_i (t)$ denotes the surrogate state of subsystem $i$ at time $t$. We require that the initial condition of the surrogate system equal that of the true system states agree, i.e.,  $\widetilde{x}_i (0) = x_i (0)$ for each subsystem $i$.
It will be convenient to express the  surrogate state dynamics in terms of the system trajectories according to
\begin{equation}
\widetilde{x} =  {\widetilde{B}} \widetilde{u} +{\widetilde{L}} w + {\widetilde{H}} v_{\mathcal{C}}. \label{eq:surrogate_state}
\end{equation}
\rtwo{It is straightforward to construct the lower block triangular matrices ${\widetilde{B}}$, ${\widetilde{L}}$, and ${\widetilde{H}}$ given the problem data in \eqref{eq:surr_state_i}}

\subsection{Assume-Guarantee Contracts} \label{sec:assume-guarantee}

Under the proposed  surrogate information structure, the information-coupling states are assumed to behave as disturbances with bounded support.
In the following definition, we introduce a class of constraints on the family of admissible control policies, which guarantee that the information-coupling states behave in accordance with this assumption.

\begin{definitio}[Assume-Guarantee Contract] \label{def:assume-guarantee}
A control policy $\phi \in \Phi$ is said to satisfy the \emph{assume-guarantee contract} specified in terms of a \emph{contract set} ${\mathcal{V}}_{\mathcal{C}} \subseteq {\bf R}^{N_x^{\mathcal{C}}}$  if
\ $\Pi_\couple \widetilde{x} \in {\mathcal{V}}_{\mathcal{C}}$   for all  $(w , v_{\mathcal{C}}  ) \in {\mathcal{W}} \times {\mathcal{V}}_{\mathcal{C}}$, \
where $ \widetilde{x} =  {\widetilde{B}} \phi (w,   v_{\mathcal{C}}) + {\widetilde{L}} w + {\widetilde{H}}  v_{\mathcal{C}} $.
\end{definitio}

Here, the set ${\mathcal{V}}_{\mathcal{C}}$ is referred to as a \emph{contract set}, as it specifies the set to which the information-coupling states are both assumed and required to belong.
The satisfaction of the assume-guarantee contract ensures that the \emph{surrogate} information-coupling states  $\widetilde{x}_{\mathcal{C}} := \Pi_\couple \widetilde{x}$ belong to the given contract set.
In the following lemma, we show that the \emph{actual} information-coupling states induced by the policy $u = \phi (w, x_{\mathcal{C}})$ are also guaranteed to belong to the contract set if the assume-guarantee contract is satisfied.

\begin{lemm} \label{lem:assume-guarantee}
Let $\phi \in \Phi$ be a control policy that satisfies the assume-guarantee contract defined by the contract set ${\mathcal{V}}_{\mathcal{C}} \subseteq {\bf R}^{N_x^{\mathcal{C}}}$. It follows that $\Pi_{\mathcal{C}} x\in {\mathcal{V}}_{\mathcal{C}}$ for all $w \in {\mathcal{W}}$, where $x = B \phi (w, x_{\mathcal{C}}) + L w$.
\end{lemm}

Appendix \ref{pf:lem:assume-guarantee} contains the proof of Lemma \ref{lem:assume-guarantee}.
Prop. \ref{prop:feas} provides an inner approximation of the decentralized control design problem \eqref{opt:decent_dFeedback} by introducing an assume-guaranatee contractual constraint.

\begin{propositio} \label{prop:feas}
Let $\phi \in \Phi$ be a feasible control policy for the following problem:
\begin{align}
\begin{alignedat}{8}
&\text{\rm minimize} \quad && {\bf E} \left[ \widetilde{x}^\top R_x \widetilde{x} + \widetilde{u} ^\top R_u \widetilde{u} \right] \\
& \text{\rm subject to} \quad && \phi \in \Phi \\
&&& \hspace{-.105in}\left. \begin{array}{l}
\widetilde{u} = \phi (w,  v_{\mathcal{C}}) \\
\Pi_\couple \widetilde{x} \in {\mathcal{V}}_{\mathcal{C}} \\
\widetilde{x} =  {\widetilde{B}} \widetilde{u} +{\widetilde{L}} w + {\widetilde{H}} v_{\mathcal{C}} \\
F_x \widetilde{x} + F_u \widetilde{u} + F_w w \leq g
\end{array}
\hspace{-0.04in}\right\}  \forall (w, v_{\mathcal{C}}) \in {\mathcal{W}} \times {\mathcal{V}}_{\mathcal{C}}.
\end{alignedat} \label{opt:contract}
\end{align}
It follows that $\phi$ is also  feasible for problem \eqref{opt:decent_dFeedback}.
\end{propositio}

We omit the proof of Prop. \ref{prop:feas}, as it is a direct consequence of Lemma \ref{lem:assume-guarantee}.
Given a fixed contract set $ {\mathcal{V}}_{\mathcal{C}}$, problem \eqref{opt:contract} is a convex disturbance feedback control design problem. The choice of the contract set  and fictitious disturbance distribution does, however, play an important role in determining the performance of the control policies that are optimal solutions to problem \eqref{opt:contract}.
In Section \ref{sec:joint_opt}, we develop a systematic approach to enable the joint optimization of the contract set with the control policy using semidefinite programming.

\section{Problem Formulation} \label{sec:formulation}

\subsection{System Model}

Consider a discrete-time, linear time-varying system consisting of $N$ coupled subsystems whose dynamics are described by
\begin{align} \label{eq:x_i}
x_i(t+1) = \sum_{j=1}^N  \left( A_{ij}(t) x_j(t)   + B_{ij} (t) u_j (t) \right) +  w_i (t),
\end{align}
for $i = 1, \dots , N$.
We denote the \emph{local state}, \emph{local input}, and \emph{local disturbance} associated with  each subsystem $i$ at time $t$ by $x_i(t) \in {\bf R}^{n_x^i}$, $u_i(t) \in {\bf R}^{n_u^i}$, and $w_i (t) \in {\bf R}^{n_x^i}$, respectively.
The system is assumed to evolve  over a finite time horizon $T$, and  the initial condition is assumed to be a random vector with known probability distribution.
Eq. \eqref{eq:x_i} can be presented more compactly as
\begin{align} \label{eq:full_sys}
x(t+1) = A(t) x(t) + B(t) u (t) + w(t).
\end{align}
Here, we denote by $x(t) := (x_1(t), .., x_N(t)) \in {\bf R}^{n_x}$, $u(t) := (u_1(t), .., u_N(t)) \in {\bf R}^{n_u}$, and $w (t) := (w_1(t), .., w_N(t)) \in {\bf R}^{n_x}$ the \emph{full system state, input,} and \emph{disturbance} at time $t$. The dimensions of the system state and input are given by $n_x := \sum_{i=1}^N n_x^i$ and $n_u := \sum_{i=1}^N n_u^i$, respectively.
The input and disturbance trajectories are related to the state trajectory according to
\begin{align}
x = B u + L w, \label{eq:trajectory}
\end{align}
where $x$, $u$, and $w$ denote the system state, input, and disturbance trajectories, respectively. They are defined according to
\begin{align}
x & := (x(0), \dots, x(T))  \in {\bf R}^{N_x},  \label{eq:x_traj_def}\\
u & := (u(0), \dots, u(T-1))  \in {\bf R}^{N_u}, \label{eq:u_traj_def}\\
w & : = (w(-1), w(0), \dots, w(T-1))  \in {\bf R}^{N_x}, \label{eq:w_traj_def}
\end{align}
where the corresponding dimensions are given by
$N_x := n_x(T+1)$ and $N_u := n_u T$.
To simplify the specification of disturbance-feedback affine control policies in the sequel, we adopt the notational convention $w (-1 ) := x(0)$. \rtwo{The lower block triangular
matrices $B$ and $L$  are straightforward to  construct from the problem data in  \eqref{eq:full_sys}.}

\subsection{Disturbance Model}

We assume that disturbance trajectory $w$ is a zero-mean random vector. We denote its second moment matrix  by $M := {\bf E} [w w^\top].$
We let  ${\mathcal{W}} \subset {\bf R}^{N_x}$ denote the support of $w$, which we assume is a convex, compact set  with a non-empty interior. This assumption   ensures that the matrix $M$ is positive definite and  finite valued.

\subsection{System Constraints}

We consider a general family of polyhedral constraints on the state and input trajectories of the form
\begin{align}
F_x x + F_u u + F_w w \leq g \qquad  \forall w \in {\mathcal{W}}, \label{eq:robust_con}
\end{align}
where $F_x \in {\bf R}^{m \times N_x}$, $F_u \in {\bf R}^{m \times N_u}$, $F_w \in {\bf R}^{m \times N_x}$, $g \in {\bf R}^{m}$ are assumed to be given. Note that, in general, such constraints may couple states and inputs across different subsystems and time periods.

\subsection{Information Structure}

We encode the pattern according to which information is shared between subsystems with a directed graph $G_I = (V,E_I)$, which we refer to as the \emph{information graph} of the system.
The vertex set $V = \{1, \dots , N \}$ indexes the subsystems being controlled, and the edge set $E_I \subseteq V \times V$ includes the directed edge
$(i,j)$ if and only if  subsystem $j$ has access to subsystem $i$'s local state measurement at each time $t$.
We let $V_I^-(i)$ denote the in-neighborhood of each subsystem $i \in V$ in the information graph $G_I$.

Each subsystem is assumed to have access to the entire history of its local information up to and including time $t$.
Formally, the \emph{local information} available to each subsystem $i$ at time $t$ is defined as
\begin{align} \label{eq:pol1}
z_i(t) : = \{ x_j^{0:t} \ | \ (j, i) \in E_I \}.
\end{align}
The local control input to each subsystem $i$ is restricted to be a causal function of its local information. That is, the local input to subsystem $i$ at time $t$ is of the form
\begin{align} \label{eq:pol2}
u_i(t)  = \gamma_i(z_i(t), t),
\end{align}
where $\gamma_i(\cdot, t)$ is a measurable function of the local information $z_i (t)$.
We  define the \emph{local control policy} for subsystem $i$ as $\gamma_i : = (\gamma_i(\cdot,0), \dots, \gamma_i(\cdot,T-1))$. We refer to the  collection of local control policies  $\gamma : = (\gamma_1, \dots, \gamma_N)$ as the \emph{decentralized control policy}, which relates the state trajectory $x$ to the input trajectory $u$ according to $u = \gamma (x)$.
Finally, we let $\Gamma$ denote the set of all decentralized control policies respecting the information constraints encoded in Eq. \eqref{eq:pol2}.

\subsection{Decentralized Control Design}
We consider the following family of constrained decentralized control design problems:
\begin{align}
\begin{alignedat}{8}
&\text{minimize} \quad && {\bf E} \left[ x^\top R_x x + u ^\top R_u u \right] \\
& \text{subject to} \quad && \gamma \in \Gamma  \\
&&& \hspace{-.105in}\left. \begin{array}{l}
u = \gamma (x) \\
x = B u + L w  \\
F_x x + F_u u + F_w w \leq g
\end{array}
\hspace{-0.04in}\right\} \forall w \in {\mathcal{W}}.
\end{alignedat} \label{opt:decent}
\end{align}
Here, the cost matrices $R_x \in {\bf R}^{N_x \times N_x}$ and $R_u \in {\bf R}^{N_u \times N_u}$ are  assumed to be symmetric and positive semidefinite.
The tractability of the decentralized control design problem \eqref{opt:decent} depends  on its information structure.
In particular, if the information structure is \emph{partially nested}, then problem \eqref{opt:decent} can be equivalently reformulated (via the Youla parameterization) as a convex program in the space of disturbance feedback policies \cite{Ho1972}.
If the information structure is \emph{nonclassical} (i.e., not partially nested), then problem \eqref{opt:decent} is known to be computationally intractable, in general \cite{Tsitsiklis1985, Sandell1978, Mahajan2012_survey}.

\section{Introduction}

We investigate the design of affine decentralized control policies for
stochastic discrete-time, linear systems that evolve over a finite horizon, and are subject to polyhedral constraints on the state and input trajectories. The computational tractability of such problems depends in part on their information structures  \cite{sandell1974solution, Tsitsiklis1985}.
A decentralized control problem is said to have a \emph{nonclassical} information structure if the information available to one controller can be affected by the control actions of another that it cannot access or reconstruct.
Under such information structures, the calculation of optimal decentralized control policies is known to be computationally intractable, because of the so called signalling incentive for controllers to communicate with each other via the actions they undertake \cite{Witsenhausen1968, sandell1974solution, Tsitsiklis1985}.
To complicate matters further, there may be hard constraints  coupling the local actions and states of different controllers that must be jointly enforced without explicit communication.
In this paper, we  address these challenges by relaxing the requirement that decentralized controllers be optimal with respect to the broad family of all causal policies, and instead search for suboptimal decentralized controllers that can be efficiently computed via convex optimization.

\emph{Related Literature: \ }
\rev{Many decentralized control methods based on  model predictive control (MPC) tools have been proposed \cite{venkat2005stability, Camponogara2002, Jia2001, Bemporad2010, Keviczky2006, dunbar2007distributed, Dunbar2006, lucia2015contract, farina2012distributed, richards2004decentralized, richards2007robust, riverso2012tube, trodden2006robust, trodden2010distributed, hernandez2017distributed, trodden2014cooperative}.}
These approaches typically rely on a decomposition of the decentralized control problem into a collection of decoupled  local control problems, where each subproblem is solved using traditional MPC methods.
A particular class of methods---commonly referred to as ``tube-based'' decentralized MPC---treat the coupling states and inputs affecting each subsystem as independent and bounded disturbances that are assumed to take values in the given state and input constraint sets  \cite{dunbar2007distributed, Dunbar2006, lucia2015contract, farina2012distributed, richards2004decentralized, richards2007robust, riverso2012tube, trodden2006robust, trodden2010distributed, hernandez2017distributed, trodden2014cooperative}.
Decentralized controllers  computed in this manner  may result in overly conservative behaviors for a
number of reasons.
\emph{First}, the treatment of the coupling states and inputs as independent disturbances ignores the potential dynamical coupling between these variables.
\emph{Second}, the over approximation of the coupling state and input trajectory sets by their corresponding state and input constraint sets will likely be  loose for many problem instances. More importantly, the over approximation of the coupling state and input trajectory sets in this manner ignores the fact that these sets depend on the control policy being used to regulate the system, and, therefore, neglects the possibility of co-optimizing their specification with the control policy.

\emph{Contribution: \ }
We provide a computationally tractable method to calculate control policies that are guaranteed to be feasible for constrained decentralized control problems with nonclassical information structures.
Loosely speaking, the proposed approach  eliminates the informational coupling between subsystems by treating the information-coupling states as  disturbances whose trajectories are ``assumed" to take values in  a  ``contract'' set that we design. To ensure the satisfaction of this assumption, we impose a contractual constraint on the control policy that ``guarantees'' that the information-coupling states that it induces indeed  belong to the contract set.
Naturally, this approach yields an inner approximation of the original decentralized control design problem, where the conservatism of the resulting approximation depends on the specification of the contract set.
To limit the suboptimality of such approximations, we formulate a semi-infinite program to co-optimize the decentralized control policy with the location, scale, and orientation of an ellipsoidal contract set. \rtwo{The ability to adjust both the scale and orientation of the contract set enables the construction of contract sets that accurately capture the spatial and inter-temporal coupling between different subsystem states.} We establish a condition on the set of allowable contracts that facilitates the joint optimization of the control policy and the contract set via  semidefinite programming.

Assume-guarantee contracts  have  been utilized in a  variety of control applications \cite{kim2015compositional, chen2019compositional,nuzzo2013contract,sangiovanni2012taming}.
Closer to the setting considered in this paper, there are related papers in the literature that investigate an approach to decentralized control design via the co-optimization of the control policy and the contract sets that govern the interactions between the different subsystems in the network \cite{trodden2017distributed, darivianakis2018decentralized}. \rtwo{The techniques developed in these papers have two primary limitations.
First, they treat all states that result in a physical coupling between different subsystems as  disturbances, while our approach only treats states that induce an informational coupling between subsystems as disturbances.
Second, these techniques only permit the scaling and translation of a base contract set when co-optimizing it with the control policy. To the best of our knowledge, the method proposed in this paper provides the first computationally tractable approach to co-optimizing the control policy with the location, scale, and \emph{orientation} of the contract set,  expanding substantially the family of contracts that can be efficiently optimized over.
We provide an example that illustrates  the advantages of the proposed method.}

\emph{Notation: \ }   Let ${\bf R}$ and ${\bf R}_+$ denote the sets of real and non-negative real numbers, respectively.
Given a collection of vectors $x_1, .., x_N$ and an index set $J \subseteq \{1, .., N\}$,
we let $x_J$ denote the concatenation of the vectors $x_j$  for $j \in J$ in ascending order of their indices. Given a sequence $\{x(t)\}$ and time indices $s \leq t$, we denote the sequence of elements from time $s$ to $t$ by  $x^{s:t} := (x(s), x(s+1), .., x(t) )$.
 Given a block matrix $A$, we let $[A]_{ij}$ denote its $(i,j)$-th block. We denote the trace of a square matrix $A$ by $\trace{A}$.
We denote the Minkowski sum of two sets ${\mathcal{S}}, {\mathcal{T}} \subseteq {\bf R}^n$  by ${\mathcal{S}} \oplus {\mathcal{T}} := \{x + y \,  | \, x \in {\mathcal{S}}, \  y \in {\mathcal{T}} \}$.

\section{Policy-Contract Optimization} \label{sec:joint_opt}

In this section, we provide a semidefinite programming-based method to co-optimize the design of the decentralized control policy together with the contract set that constrains its design.
As part of the proposed approach, we consider a restricted family of control policies that are affinely parameterized in both the disturbance and fictitious disturbance histories. We also parameterize the fictitious disturbance process as a causal affine function of a given (primitive) disturbance process---an approach that is similar in nature to the class of parameterizations that have been recently studied in the context of robust optimization with adjustable uncertainty sets \cite{zhang2017robust}. As one of our primary results in this section, we identify a structural condition on the family of allowable contract sets that facilitates the inner approximation of the resulting policy-contract optimization problem as a semidefinite program.

\subsection{Affine Decentralized Control Policies}
We consider affine disturbance-feedback  policies of the form
\begin{align}
\nonumber \widetilde{u}_i (t) =  u^o_i (t) &+ \sum_{ j \in {\mathcal{N}}  (i) }   \sum_{s=-1}^{t-1} Q_{ij}^w (t,s+1) w_j (s)  \\
&+ \sum_{ j \in {\mathcal{C}} (i)} \sum_{s=0}^t  Q_{ij}^v (t,s) v_j (s), \label{eq:aff1}
\end{align}
for $t = 0, .., T-1$ and $i = 1, .., N$.
We refer to  $u^o_i (t) \in {\bf R}^{n_u^i}$ as the open-loop control input, and the matrices $Q_{ij}^w (t,s+1) \in {\bf R}^{n_u^i \times n_x^j}$ and $ Q_{ij}^v (t,s) \in {\bf R}^{n_u^i \times n_x^j}$ as the feedback control gains.
The affine control policy \eqref{eq:aff1} can be expressed in terms of trajectories as
\begin{align}
\widetilde{u} = u^o +  Q^w  w + Q^v v, \label{eq:aff2}
\end{align}
where the  matrices $Q^w \in {\bf R}^{N_u \times N_x}$ and $Q^v \in {\bf R}^{N_u \times N_x}$ are both $T \times (T+ 1)$ block matrices, whose $(t, s)$-th blocks are defined as
\begin{align}
[Q^w (t,s) ]_{ij} &= \begin{cases}
Q_{ij}^w (t,s) & \text{if } j \in {\mathcal{N}} (i) , \ t \geq s , \\
0 & \text{otherwise},
\end{cases}  \label{eq:Qw}\\
[Q^v (t,s)] _{ij} &= \begin{cases}
Q_{ij}^v (t,s) & \text{if } j \in {\mathcal{C}} (i)  , \  t \geq s , \\
0 & \text{otherwise}.
\end{cases}  \label{eq:Qfic}
\end{align}
for $i = 1, \dots, N$ and  $j = 1, \dots, N$.
We let ${\mathcal{Q}}_N$ and ${\mathcal{Q}}_C$ denote the matrix subspaces  respecting the block sparsity patterns specified according to Eqs. \eqref{eq:Qw} and \eqref{eq:Qfic}, respectively.

\subsection{Affine Parameterization of the Fictitious Disturbance} \label{sec:param_fict}

We focus our analysis on fictitious disturbances that are expressed according to an affine transformation of a \emph{primitive disturbance}. Such a parameterization yields contract sets that have adjustable location, scale, and orientation.
Specifically, we let the random vector $\xi \in {\bf R}^{N_x}$ denote the the \emph{primitive disturbance trajectory}, whose support $\Xi \subseteq {\bf R}^{N_x}$ is assumed to be convex and compact. \rtwo{The distribution of the primitive disturbance trajectory $\xi$ is a design choice. In Assumption \ref{ass:ellip}, we propose a particular distribution that enables the co-optimization of the control policy and the contract set via semidefinite programming.}

We parameterize the fictitious disturbance trajectory affinely in the primitive disturbance as
\begin{align}
v := \overline{v} + Z \xi. \label{eq:zeta}
\end{align}
The parameters  $\overline{v} \in {\bf R}^{N_x}$ and $Z \in {\bf R}^{N_x \times N_x}$  control the shape of the resulting contract set, which is given by ${\mathcal{V}}_{\mathcal{C}} = \Pi_\couple \left( \overline{v} \oplus Z \Xi \right).$
Throughout the paper, we will restrict our attention to transformations \eqref{eq:zeta} in which the matrix parameter $Z$ is both lower triangular and invertible. We denote the set of all such matrices by  ${\mathcal{Z}} \subset {\bf R}^{N_x \times N_x}$.

The specification of the fictitious disturbance according to Eq. \eqref{eq:zeta} induces the following structure in the surrogate control input:
\begin{align}
\widetilde{u} = u^o + Q^v \overline{v} +  Q^w  w + Q^v Z \xi. \label{eq:aff3}
\end{align}
We eliminate the bilinear terms in Eq. \eqref{eq:aff3} through the following  the change of variables:
\begin{align}
 \overline{u}:= u^o + Q^v \overline{v}   \ \text{ and } \ Q^\xi := Q^v Z. \label{eq:change_of_var}
\end{align}
This change of variables gives rise to a reparameterization of the surrogate input trajectory as
\begin{align}
\widetilde{u} = \overline{u} +  Q^w  w + Q^\xi \xi, \label{eq:aff4}
\end{align}
where the matrix $Q^\xi \in {\bf R}^{N_u \times N_x}$ must satisfy the sparsity constraint
$Q^\xi Z^{-1} \in {\mathcal{Q}}_C$
in order to ensure the satisfaction of the original sparsity constraint that $Q^v \in {\mathcal{Q}}_C$.
The affine parameterization of the control policy and the fictitious disturbance  according to   \eqref{eq:aff2} and \eqref{eq:zeta}, respectively, facilitates the co-optimization of the control policy and contract set as follows:
\begin{align}
\begin{alignedat}{8}
&\text{minimize} \quad && {\bf E} \left[ \widetilde{x}^\top R_x \widetilde{x} + \widetilde{u} ^\top R_u \widetilde{u} \right] \\
& \text{subject to} \quad &&Q^w \in {\mathcal{Q}}_N, \,  Q^\xi\in {\bf R}^{N_u \times N_x} , \,  Z \in {\mathcal{Z}} \\
&&& \overline{u} \in {\bf R}^{N_u}, \, \overline{v} \in {\bf R}^{N_x} , \\
&&& Q^\xi Z^{-1} \in {\mathcal{Q}}_C \\
&&& \hspace{-.105in} \left. \begin{array}{l}
v  =  \overline{v}  + Z \xi\\
\widetilde{u} = \overline{u} +  Q^w  w + Q^\xi \xi \\
\widetilde{x} =  {\widetilde{B}} \widetilde{u} +{\widetilde{L}} w + {\widetilde{H}} v_{\mathcal{C}}  \\
\Pi_\couple \widetilde{x} \in \Pi_\couple  \left( \overline{v} \oplus Z {\mathcal{W}} \right) \\
F_x \widetilde{x} + F_u \widetilde{u} + F_w w \leq g
\end{array}
\hspace{-0.04in}\right\} \forall (w, \xi) \in {\mathcal{W}} \times \Xi.
\end{alignedat} \label{opt:inner0_joint}
\end{align}

In contrast to problem \eqref{opt:contract}, where the contract set is considered fixed, problem \eqref{opt:inner0_joint}
allows for the co-optimization of the control  policy  together  with  the  contract set.
This increased flexibility in control design does,  however, result in the loss of problem convexity. Specifically,  the nonconvexity in  problem \eqref{opt:inner0_joint} stems from the assume-guarantee contractual constraint on the affine control policy, and the bilinear sparsity constraint  $Q^\xi Z^{-1} \in {\mathcal{Q}}_C$. In Section \ref{sec:restrict}, we construct explicit convex inner approximations for each of these constraints, yielding a convex inner approximation of problem \eqref{opt:inner0_joint} in the form of a semidefinite program.

\subsection{Restricting the Contract Set}
\label{sec:restrict}

In what follows, we introduce additional restrictions on the set of allowable of matrix parameters $\mathcal{Z}$ that permit the reformulation of the nonconvex bilinear constraint $Q^\xi Z^{-1} \in {\mathcal{Q}}_C$ as the linear constraint $Q^\xi \in {\mathcal{Q}}_C$.
\sam{We do so by imposing an additional subspace constraint on $Z$ that guarantees the invariance of the subspace ${\mathcal{Q}}_C$ under multiplication with such matrices $Z$.}
First, we impose a structural constraint the matrix parameter $Z$ of the form
\begin{align}
Z = \lambda I - Y. \label{eq:fictransform}
\end{align}
Here, $\lambda \geq 1$ is scalar parameter and $Y \in {\bf R}^{N_x \times N_x}$ is a $(T+1) \times (T+1)$ strictly block lower triangular matrix of the form
\begin{align}
 Y = \bmat{0 \\ Y(1,0) & 0 \\
 \vdots &   \ddots  & \ddots\\
 Y (T, 0 ) &  \cdots & Y(T, T-1) & 0}. \label{eq:Y}
 \end{align}
In addition, each block of the matrix $Y$ is required to be a $N \times N$ block matrix, whose $(i, j)$-th block has dimension $n_x^i \times n_x^j$.   We impose an additional restriction on the structure of the matrix $Y$ in the form of sparsity constraints on each of its blocks that reflect the pattern of informational coupling between subsystems.
Specifically, we encode the pattern of informational coupling between subsystems according to a directed graph $\gr_C := (V, \edges_C)$, where
\begin{align*}
\edges_C := \{(j, i) \in E_I \, | \, j \in {\mathcal{C}} (i) \}.
\end{align*}
We let $V_C^+ (i)$ denote the out-neighborhood of a node $i \in V $ in the \emph{information coupling graph} $\gr_C$.
Using this graph, we impose a sparsity constraint  on each  block of the matrix $Y$ of the form:
\begin{align}
[Y (t, s) ]_{ij} = 0 \quad \text{if } \ V_C^+ (i) \nsubseteq V_C^+ (j)  \label{eq:subspace_Y}
\end{align}
for all $i,j \in \{1, \dots, N\}$, and $t, s \in  \{0,\dots, T\}$. We denote the subspace of all matrices $Y$ that respect these sparsity constraints by
\begin{align*}
    {\mathcal{Y}} (\gr_C) := \{ Y \in {\bf R}^{N_x \times N_x} \ | \ Y \ \text{satisfies  \eqref{eq:Y} \& \eqref{eq:subspace_Y}} \}
\end{align*}

Lemma \ref{lem:QY} establishes the invariance of the subspace ${\mathcal{Q}}_C$ under multiplication by matrices
$Y \in {\mathcal{Y}} (\gr_C)$.
\rtwo{The result can be verified by direct calculation. The proof  is omitted due to space limitations.}

\begin{lemm} \label{lem:QY} If $Q \in {\mathcal{Q}}_C$ and $Y \in {\mathcal{Y}} (\gr_C)$, then $Q Y \in {\mathcal{Q}}_C$.
\end{lemm}

The following useful result is a consequence of Lemma \ref{lem:QY}.

\begin{lemm} \label{lem:change_of_variable}
Let  $Y \in {\mathcal{Y}} (\gr_C)$  and $\lambda \in [1, \infty)$. It follows that $${\mathcal{Q}}_C = \big\{ Q^\xi ( \lambda I - Y)^{-1} \, | \, Q^\xi \in {\mathcal{Q}}_C  \big\}.$$
\end{lemm}
\def F{F}
\prf{Fix $\lambda \in [1, \infty)$ and $Y \in {\mathcal{Y}} (\gr_C)$.
Define the finite-dimensional linear map $F: {\bf R}^{N_u \times N_x} \to {\bf R}^{N_u \times N_x}$ according to
$F (Q) := Q (\lambda I - Y)^{-1}.$
We show that $F ({\mathcal{Q}}_C  ) = {\mathcal{Q}}_C$.
The invertibility of the linear map $F$, in combination with the finite dimensionality of the subspace ${\mathcal{Q}}_C$, implies that $F ({\mathcal{Q}}_C  ) = {\mathcal{Q}}_C$ if and only if $F^{-1} ({\mathcal{Q}}_C ) \subseteq {\mathcal{Q}}_C$. Note that
$F^{-1} (Q) = Q (\lambda I - Y).$
The desired result follows, as Lemma \ref{lem:QY} implies that $QY \in {\mathcal{Q}}_C$ for each $Q \in {\mathcal{Q}}_C$ and $Y \in {\mathcal{Y}} (\gr_C)$.}

It follows from Lemma \ref{lem:change_of_variable} that if we impose the additional  constraint $$Z \in \{ \lambda I - Y \ | \ \lambda \in [1, \infty) \ \text{\&} \ Y \in {\mathcal{Y}} (\gr_C) \},$$ then the bilinear constraint $Q^\xi Z^{-1} \in {\mathcal{Q}}_C$ can be
rewritten as an equivalent linear constraint $Q^\xi \in {\mathcal{Q}}_C$.

\section{Illustrative Example} \label{sec:case_study}

\rtwo{In this section, we apply the control design methodology proposed in this paper to a linear time-invariant system composed of $N = 3$ subsystems, where the local state and input dimensions  of each subsystem are $n_x^i = n_u^i = 1$. The control horizon is set to $T = 15$.
The system matrices are given by
\begin{align*}
A(t) &= \bmat{0.5  & 0 & 0 \\
0.5 & 0.5 & 0 \\
0.5 & 1.2 & -1.2  }  \quad \text{and} \quad
B(t) = \mathop{\mathrm{diag}} (0.1, \,  1, \,  1)
\end{align*}
for all $t = 0, \dots, T-1$.
The system disturbance trajectory $w$ is assumed to be a uniformly distributed random vector over the hypersphere given by ${\mathcal{W}} := \{w \in {\bf R}^{3 (T+1)} \, \left| \, \lVert w \rVert_2 \leq 1  \right. \}$.
The system state and input trajectories are required to satisfy the constraints
\begin{align*}
\lVert x \rVert_{\infty} \leq x_{\max} \quad \text{and} \quad \lVert u \rVert_{\infty} \leq 2.5 \quad  \forall w \in {\mathcal{W}},
\end{align*}
where  $x_{\max} \in {\bf R}_+$ is constraint parameter that we will vary in our numerical studies. The cost matrices are specified as
\begin{align*}
R_x = I_{T+1} \otimes \mathop{\mathrm{diag}} (0.1, \, 0.1, \, 2) \quad \text{and} \quad
R_u = I_{T} \otimes \mathop{\mathrm{diag}} (5,\, 5, \, 1),
\end{align*}
where $\otimes$ denotes the Kronecker product operator.}

\rtwo{In this case study, we evaluate the performance our proposed control design method using two different information graphs depicted in Fig. \ref{fig:info_graph}. The information graph $G_{I_1}$ in Fig. \ref{fig:info_graph}(a) induces a  decentralized control design problem with a \emph{partially nested} information structure, while  the information graph $G_{I_2}$ in Fig. \ref{fig:info_graph}(b) results in a \emph{nonclassical} information structure. We  also compare our approach to a  closely related contract-based decentralized control design method  from  \cite{darivianakis2018decentralized}.}

\input{info_graph}

\vspace{-.1in}

\subsection{Numerical Results and Discussion}

\rtwo{Focusing initially on the class of decentralized control problems induced by the information graph  $G_{I_1}$ in Fig. \ref{fig:info_graph}(a), we compare the performance of our control design approach against that of the method proposed in \cite{darivianakis2018decentralized}. In Fig. \ref{fig:parametric}, we plot the cost incurred by controllers computed based on our method (solid blue line)  and  the cost incurred by controllers computed according to the method in  \cite{darivianakis2018decentralized} (solid red line)
for different values of the state constraint parameter $x_{\max}$ ranging between 1 and 2.5.  We also plot a lower bound (dashed black line) on the optimal value of each problem instance \eqref{opt:decent} using a recently proposed convex relaxation \cite{lin2019convexinfo}.
Notice that our control design technique results in  controllers that  strictly outperform the controllers generated by the method from \cite{darivianakis2018decentralized} for each value of the constraint parameter $x_{\max}$.    Moreover, our method results in controllers that appear to be globally optimal for constraint parameter values $x_{\max} \in  [1.3, \, 2.5]$, as they achieve the lower bound on the optimal cost in that parameter regime.  The cost curves are not plotted for values of the constraint parameter $x_{\max}$ where the underlying control design method fails to return a feasible controller.  Notice that the method from this paper is more successful in generating  feasible controllers for values of the constraint parameter $x_{\max}$ close to one.

\begin{figure}[http]
\centering
\includegraphics[width = 1 \linewidth]{parametric_TAC.pdf}
\caption{\footnotesize \rtwo{For different values of the constraint parameter $x_{\max}$ ranging between 1 and 2.5, we plot the cost incurred by controllers generated by: (i) our method using the information graph $G_{I_1}$ (blue solid line); (ii) our method using the information graph $G_{I_2}$ (green solid line); (iii) the method in \cite{darivianakis2018decentralized} using the information graph $G_{I_1}$ (red solid line). The cost curves are not plotted for values of the constraint parameter $x_{\max}$ where the underlying control design method fails to return a feasible controller. The black dashed line depicts a lower bound on the optimal value of each problem instance using the convex relaxation method from \cite{lin2019convexinfo}.}}

\label{fig:parametric}
\end{figure}

In Fig. \ref{fig:parametric}, we also  plot the cost incurred by controllers generated by our method (solid green line) using the smaller information graph $G_{I_2}$ depicted in Fig. \ref{fig:info_graph}(b).
Despite using less information, our approach still results in controllers that substantially outperform those generated by \cite{darivianakis2018decentralized} for a wide range of constraint parameter values.

There are two possible explanations for these observed differences in performance between the two methods. First, the approximation technique proposed in \cite{darivianakis2018decentralized}  treats all neighboring states to a subsystem as fictitious disturbances, while our approach only treats neighboring states that result in  `information coupling' as fictitious disturbances. Second, the method from \cite{darivianakis2018decentralized} does not permit the rotation of the primitive contract set when co-optimizing its specification with the decentralized control policy.  Because of this limitation, the method in \cite{darivianakis2018decentralized} may fail to generate contract sets that accurately capture the spatial and inter-temporal correlation between different states, resulting in more conservative approximations.
We illustrate this limitation in Fig. \ref{fig:support}, where we plot the state-support sets and contract sets generated by each method for $x_3(t)$ at two consecutive time periods $t=9$ and $t=10$. Notice that our method generates an ellipsoidal contract set that more accurately captures
the correlation between $x_3(9)$ and $x_3(10)$.
}

\begin{figure}[h]
    \centering
    \includegraphics[width = 0.95 \linewidth]{support.pdf}
    \caption{\footnotesize \rtwo{Plots of the state-support sets (solid boundaries) and the contract sets (dashed boundaries) associated with states $(x_3 (9), x_3 (10))$. The sets  generated by the method in \cite{darivianakis2018decentralized} using the  information graph $G_{I_1}$ are depicted in red. The sets  generated by our method using the  information graph $G_{I_2}$ are depicted in blue. The state constraint set is depicted in black.}}
    \label{fig:support}
\end{figure}

\section{Information Decomposition} \label{sec:preliminaries}
The primary difficulty in solving decentralized control design problems stems from  the informational coupling that emerges when a subsystem's local information is affected by prior control actions that it cannot access or reconstruct. With the aim of isolating the effects of these actions on the  information available to each subsystem, we propose an   information  decomposition that partitions the local information available to each subsystem into a partially nested  subset (i.e., an information subset that is unaffected by control actions previously applied to the system) and its complement.
This  decomposition enables an equivalent reformulation of the decentralized control design problem where the control policy is expressed as an explicit function of the system disturbance and the so called \emph{information-coupling} states.
This reformulation will serve as the basis for the contract-based approach to decentralized control design  in Sec. \ref{sec:feas_control}.

\subsection{Decomposition of Local Information}
We decompose the local information available to each subsystem according to a partition of  its in-neighbors in the  information graph $G_I$.  Specifically, for each subsystem $i \in V$, we let ${\mathcal{N}}(i) \subseteq V_I^{-}(i)$ denote the set of in-neighboring subsystems such that the information conveyed by their local state measurements is unaffected by the prior control actions of any subsystem.
This requirement is satisfied if the local information of subsystem $i$  permits the reconstruction of all states and control actions directly affecting the local states of all subsystems belonging to ${\mathcal{N}}(i)$.
We denote the  complement of this set by ${\mathcal{C}}(i) : = V_I^{-}(i) \setminus {\mathcal{N}}(i)$ for each subsystem $i \in V$.

With the goal of providing an explicit characterization of this partition, we first provide a characterization of the physical coupling between different subsystems as reflected by the block sparsity patterns of the system matrices $A$ and $B$. We describe this coupling in terms of a pair of directed graphs, $G_A := (V, E_A)$ and $G_B : = (V, E_B)$, whose edge sets  are defined according to
\begin{align*}
E_A := \{ (j,i) \in V \times V \ | \  \exists t = 0, \dots, T-1 \text{ s.t. } A_{ij} (t)\neq 0\},\\
E_B := \{ (j,i) \in V \times V \ | \  \exists t = 0, \dots, T-1 \text{ s.t. } B_{ij} (t)\neq 0\}.
\end{align*}
We let  $V_A^- (i)$ and $V_B^- (i)$ denote the in-neighborhoods associated with each node $i \in V$ in  $G_A$ and $G_B$, respectively.

Building on this graphical representation of the physical coupling between subsystems,  the following definition  formalizes the class of information decompositions considered in this paper. For each subsystem $i \in V$, define the set
\begin{align*}
{\mathcal{N}}(i) := \{j \in  V_I^- (i) \ | \  \text{\eqref{eq:cond1}, \eqref{eq:cond2} are satisfied}   \},
\end{align*}
where the above conditions  are given by
\begin{align}
V_A^- (j) &\subseteq V_I^- (i) , \label{eq:cond1} \\
\bigcup_{k \in V_B^- (j)} \hspace{-.05in} V_I^- (k) &\subseteq V_I^- (i). \label{eq:cond2}
\end{align}
Condition \eqref{eq:cond1} requires that subsystem $i$ possess access to all states that directly affect subsystem $j$'s state through the system dynamics. Condition \eqref{eq:cond2} requires that subsystem $i$ have access to the local information of each subsystem whose control actions directly affect subsystem $j$'s state. This  ensures that subsystem $i$ is able to reconstruct all control actions that directly affect subsystem $j$'s state.
Collectively, conditions \eqref{eq:cond1} and \eqref{eq:cond2} can be interpreted as a requirement on the \emph{local nesting of information}, in the sense that if $j \in {\mathcal{N}}(i)$, then  subsystem $i$ is assumed to  have access to all states and control actions that directly affect subsystem $j$'s state through the state equation. As a result, subsystem $i$ can explicitly reconstruct the local disturbance $w_j (t)$ acting on any subsystem $j \in {\mathcal{N}}(i)$ based only on its local information $z_i (t+1)$ as follows:
\begin{align*}
\nonumber w_j (t) = x_j (t+1)  - \hspace{-.1in} \sum_{k \in V_A^- (j)} A_{jk}(t) x_k(t)  - \hspace{-.1in}  \sum_{k \in V_B^- (j)} B_{jk} (t) u_k (t) .
\end{align*}

The local states of subsystems not belonging to ${\mathcal{N}}(i)$ may contain information that can be influenced by prior control actions. We refer to these states as the \emph{information-coupling states} associated with subsystem $i$ at stage $t$, and  denote them by $x_{{\mathcal{C}}(i)} (t)$ where  $${\mathcal{C}} (i) := V_I^- (i) \setminus {\mathcal{N}}(i).$$
The collection of information-coupling states across all subsystems are denoted by the  $x_{\mathcal{C}} (t) \in {\bf R}^{n_x^{\mathcal{C}}}$, where ${\mathcal{C}} := \bigcup_{i \in V} \, {\mathcal{C}}(i).$
The \emph{trajectory of information-coupling states} is denoted by
\begin{align*}
x_{\mathcal{C}} := (x_{\mathcal{C}} (0), \dots, x_{\mathcal{C}} (T) ) \in {\bf R}^{N_x^{\mathcal{C}}},
\end{align*}
where $N_x^{\mathcal{C}} := n_x^{\mathcal{C}} (T+1).$
It will be notationally convenient to express the mapping from the state trajectory $x$ to its subvector $x_{\mathcal{C}}$ in terms of the projection operator  $\Pi_\couple: {\bf R}^{N_x} \to {\bf R}^{N_x^{\mathcal{C}}}$, where $x_{\mathcal{C}} = \Pi_{\mathcal{C}} x$.

\begin{remar}[Partially Nested Information] It can be shown that the given information structure is \emph{partially nested} if and only if the set of  information coupling states is empty, i.e., ${\mathcal{C}} = \emptyset$.
\end{remar}

\subsection{Control Input Reparameterization}

The proposed information decomposition $ V_I^- (i) = {\mathcal{N}}(i) \cup {\mathcal{C}}(i)$ suggests  a natural reparameterization of the control policy in terms of the following equivalent information set.

\begin{lemm}[Equivalent Information Sets] \label{lem:equi_info}
Define the information set $\zeta_i (t) := \{x_j^{0:t}  \, | \, j \in {\mathcal{C}} (i) \} \cup \{w_j^{-1: t-1} \, | \, j \in {\mathcal{N}} (i) \}$.
The sets $z_i (t)$ and $\zeta_i (t)$ are functions of each other for each subsystem $i$ and time $t$.
\end{lemm}

We omit the proof of Lemma \ref{lem:equi_info}, as it mirrors the proof of a closely related result in \cite[Lemma 1]{Lin2016}.
Lemma \ref{lem:equi_info} suggests the following equivalent parameterization of the local control input:
\begin{align}
u_i (t) = \phi_i (\zeta_i (t), t), \label{eq:pol_equiv}
\end{align}
where $\phi_i (\cdot ,t)$ is a measurable function of its arguments.
We let $\phi_i := ( \phi_i (\cdot, 0), .., \phi_i (\cdot, T-1) )$ and $\phi := (\phi_1, .., \phi_N)$ denote the reparameterized control policy for each subsystem $i \in V$ and the full system, respectively.
With a slight abuse of notation, we express the input trajectory induced by the reparameterized  policy $\phi$ as $$u = \phi (w, x_{\mathcal{C}}).$$
Finally, we denote by $\Phi$ the set of reparameterized decentralized control policies that respect the information constraints implied by Eq. \eqref{eq:pol_equiv}.
The reparameterization of the control policy  according to Eq. \eqref{eq:pol_equiv}  results in the following  equivalent reformulation of the original decentralized control problem \eqref{opt:decent}:
\begin{align}
\begin{alignedat}{8}
&\text{minimize} \quad && {\bf E} \left[ x^\top R_x x + u ^\top R_u u \right] \\
& \text{subject to} \quad && \phi \in \Phi  \\
&&& \hspace{-.105in}\left. \begin{array}{l}
u = \phi (w, x_{\mathcal{C}}) \\
x = B u + L w  \\
F_x x + F_u u + F_w w \leq g
\end{array}
\hspace{-0.04in}\right\} \forall w \in {\mathcal{W}}.
\end{alignedat} \label{opt:decent_dFeedback}
\end{align}
Problem \eqref{opt:decent_dFeedback} is nonconvex, in general, if the set of information-coupling subsystems is nonempty, i.e.,  ${\mathcal{C}} \neq \emptyset$. In Section \ref{sec:feas_control}, we construct a convex inner approximation of problem \eqref{opt:decent_dFeedback},  where the information-coupling states are assumed to behave as disturbances with bounded support, and the control policy is constrained in a manner that ensures consistency between  the assumed and actual behaviors of the information-coupling states.

\subsection{Semidefinite Programming Approximation} \label{sec:sdp}

\rev{
In this section, we introduce a series of convex restrictions that culminate in the conservative approximation of
problem \eqref{opt:inner0_joint} as a semidefinite program. Our results rely on the following assumption concerning the distributions of the disturbance trajectory $w$ and primitive disturbance trajectory $\xi$.

\begin{assumptio} \label{ass:ellip}
Let $\Sigma \succ 0$. The disturbance trajectory $w$ is assumed to have an ellipsoidal support set given by
\begin{align*}
{\mathcal{W}}: = \big\{ z \in {\bf R}^{N_x} \, \left|  \, z^\top \Sigma^{-1} z \leq 1\right. \big\}.
\end{align*}
\rtwo{We also  assume that the primitive disturbance trajectory $\xi$ is an independent and identically distributed  copy of the disturbance trajectory $w$. In particular, this implies that $\Xi = {\mathcal{W}}$.}
\end{assumptio}

In what follows, we employ  these assumptions to  reformulate the robust linear inequality constraints specified in problem \eqref{opt:inner0_joint} as second order cone constraints. We also construct a conservative approximation of the assume-guarantee contractual constraint in problem \eqref{opt:inner0_joint} according to a finite collection of linear matrix inequalities.
}

To lighten notation, we write the surrogate state trajectory $\widetilde{x}$ as
\begin{align*}
\widetilde{x} = \overline{x} +  {P^w}  w +  {P^\xi} \xi ,
\end{align*}
where $\overline{x}  := {\widetilde{B}} \overline{u} + {\widetilde{H}} \Pi_\couple \overline{v}$, ${P^w}:= {\widetilde{B}} Q^w + {\widetilde{L}}$,  and ${P^\xi}:= {\widetilde{B}} Q^\xi + {\widetilde{H}} \Pi_\couple (\lambda I - Y)$.

The following result provides an equivalent reformulation of the robust linear inequality constraints specified in problem \eqref{opt:inner0_joint} as second-order cone constraints. The proof of Lemma \ref{lem:robust_linear} follows directly from the identity $\sup_{w \in {\mathcal{W}}} c^\top w = \lVert \Sigma^{1/2} c \rVert_2$ for all $c \in {\bf R}^{N_x}$.
\begin{lemm} \label{lem:robust_linear} Let Assumption \ref{ass:ellip} hold.
The  robust  constraints
$$F_x \widetilde{x} + F_u \widetilde{u} + F_w w \leq g \quad \forall  (w,  \xi) \in {\mathcal{W}} \times \Xi$$ are satisfied if and only if
\begin{align}
\nonumber & \left\lVert \Sigma^{1/2} e_i^\top  ( F_x {P^w} +  F_u Q^w + F_w ) \right\rVert_2 + \left\lVert \Sigma^{1/2} e_i^\top ( F_x {P^\xi} + F_u Q^{\xi}) \right\rVert_2 \\
& \qquad \leq e_i^\top ( g - F_x \overline{x} - F_u \overline{u} )  \label{eq:robust_linear}
\end{align}
for $i = 1, \dots, m$, where $e_i \in {\bf R}^m$ is the $i$th standard basis vector.
\end{lemm}

We now address the nonconvexity that stems from the assume-guarantee contractual constraint specified in problem \eqref{opt:inner0_joint}:
\begin{align} \label{eq:assgar}
    \Pi_\couple \widetilde{x} \in \Pi_\couple  \left( \overline{v} \oplus Z {\mathcal{W}} \right) \quad \forall  (w,  \xi) \in {\mathcal{W}} \times \Xi.
\end{align}
First, notice  that, under the stated assumptions, the contractual constraint \eqref{eq:assgar} is equivalent to the following set containment constraint:
\begin{align}
\Pi_\couple \left( \overline{x} \oplus {P^w} {\mathcal{W}} \oplus  {P^\xi} {\mathcal{W}} \right)    \subseteq \Pi_\couple  \left( \overline{v} \oplus  Z {\mathcal{W}}  \right). \label{eq:containment}
\end{align}
The set containment constraint \eqref{eq:containment} requires that the Minkowski sum of two ellipsoids be contained within another ellipsoid.   We leverage on the following known result  to conservatively approximate this set containment constraint by a quadratic matrix inequality.

\begin{lemm}[Theorem 4.2 in \cite{durieu2001multi}] \label{lem:robust_quadratic}
Let Assumption \ref{ass:ellip} hold and let
 $L_1, L_2, L_3 \in {\bf R}^{m \times N_x}$, where it is assumed that $m \leq N_x$. If there exists a scalar $\alpha \in [0, 1]$ such that
\begin{align}
\alpha^{-1} L_1 \Sigma L_1^\top +  (1 - \alpha)^{-1} L_2 \Sigma L_2^\top \preceq  L_3 \Sigma L_3^\top. \label{eq:MI_lemma0}
\end{align}
then $L_1 {\mathcal{W}} \oplus L_2 {\mathcal{W}} \subseteq L_3 {\mathcal{W}}$.
\end{lemm}

\rev{
A direct application of Lemma \ref{lem:robust_quadratic} reveals that the set containment constraint \eqref{eq:containment} is satisfied if there exists a scalar $\alpha \in [0, 1]$ such that
\begin{align}
&\Pi_\couple (\overline{v} - \overline{x} ) = 0 \\
&\nonumber \Pi_\couple \left( \alpha^{-1}  {P^w} \Sigma {P^w}^\top  + (1 - \alpha)^{-1} {P^\xi} \Sigma {P^\xi}^\top \right) \Pi_\couple^\top  \\
& \qquad \preceq \Pi_\couple (\lambda I - Y) \Sigma (\lambda I - Y)^\top \Pi_\couple^\top. \label{eq:matrix_inequality1}
\end{align}
The quadratic matrix inequality \eqref{eq:matrix_inequality1} is  nonconvex in the decision variables $\lambda$ and $Y$. The following result provides a conservative approximation of the quadratic matrix inequality \eqref{eq:matrix_inequality1} as a linear matrix inequality.}  Appendix \ref{pf:lem:matrix_inequality} contains a proof of Lemma \ref{lem:matrix_inequality}.
\begin{lemm} \label{lem:matrix_inequality}
Let Assumption \ref{ass:ellip} hold. The set containment constraint \eqref{eq:containment} is satisfied if there exists a scalar $\beta \in [0, \lambda]$ such that
\begin{align}
&\Pi_\couple \left(   \overline{x} - \overline{v} \right) =0 , \label{eq:LMI0}\\
&\bmat{\Pi_\couple \widetilde{\Sigma} \Pi_\couple^\top & \Pi_\couple {P^w} & \Pi_\couple {P^\xi} \\
{P^w}^\top \Pi_\couple^\top & \beta \Sigma^{-1} & 0 \\
{P^\xi} ^\top \Pi_\couple^\top & 0 & (\lambda - \beta) \Sigma^{-1} } \succeq 0, \label{eq:LMI}
\end{align}
where  $\widetilde{\Sigma} = \lambda \Sigma - Y \Sigma - \Sigma Y^\top.$
\end{lemm}

The combined application of Lemmas \ref{lem:change_of_variable}, \ref{lem:robust_linear}, and \ref{lem:matrix_inequality}
enables the conservative approximation of the nonconvex decentralized control design problem  \eqref{opt:inner0_joint} as a finite-dimensional semidefinite program. We state the resulting approximation in the following proposition.

\begin{propositio} \label{prop:final}
The  following semidefinite program is  a conservative approximation of problem \eqref{opt:inner0_joint}:
\begin{align}  \label{opt:inner4_joint}
\begin{alignedat}{8}
&\text{\rm minimize} \quad && {\rm Tr} \left( {P^\xi}^\top R_x  {P^\xi} M   +  {P^w}^\top R_x  {P^w} M \right) + \overline{x}^\top R_x  \overline{x}\\
&&& +  {\rm Tr} \left({Q^w}^\top R_u Q^w M  +   {Q^\xi}^\top R_u Q^\xi M  \right)+  \overline{u}^\top R_u  \overline{u} \\[3pt]
& \text{\rm subject to} \quad &&Q^w \in {\mathcal{Q}}_N,  \, Q^\xi \in {\mathcal{Q}}_C , \,  Y \in {\mathcal{Y}} (\gr_C) , \, \overline{u} \in {\bf R}^{N_u}\\
&&& {P^w}, \, {P^\xi}  \in {\bf R}^{N_x \times N_x}, \ \overline{v}  , \,  \overline{x} \in {\bf R}^{N_x} ,  \ \lambda , \,  \beta \in {\bf R}_+ \\
&&& \lambda \geq \max \{ 1, \beta \} \\
&&& \overline{x} =  {\widetilde{B}} \overline{u} + {\widetilde{H}} \Pi_\couple \overline{v} \\
&&& {P^w} = {\widetilde{B}} Q^w + {\widetilde{L}} \\
&&& {P^\xi} = {\widetilde{B}} Q^\xi + {\widetilde{H}} \Pi_\couple (\lambda I - Y) \\
&&& \eqref{eq:robust_linear}, \ \eqref{eq:LMI0}, \ \eqref{eq:LMI}.
\end{alignedat}
\end{align}
The decision variables for problem \eqref{opt:inner4_joint} are the matrices $Q^w$, $Q^\xi$, $Y$, ${P^w}$, ${P^\xi}$, the vectors $\overline{u}$, $\overline{v}$, $\overline{x}$, and the scalars $\lambda$ and $\beta$.\footnote{We note that the optimization variables $\overline{x}$, ${P^w}$,  and ${P^\xi}$ have been introduced  to simplify the statement of problem \eqref{opt:inner4_joint}. These additional  variables can be eliminated via direct substitution. }
\end{propositio}

It follows from Prop. \ref{prop:final} that  any feasible solution to problem \eqref{opt:inner4_joint} can be mapped to an affine control policy that is guaranteed to be feasible for the original decentralized control design problem \eqref{opt:decent_dFeedback}, using the change of variables specified in  \eqref{eq:change_of_var}. Specifically, given any feasible solution to problem \eqref{opt:inner4_joint}, the following  affine control policy is guaranteed to be a feasible solution to problem \eqref{opt:decent_dFeedback}:
\begin{align}  \label{eq:final_aff}
u = \overline{u} \, + \, Q^w w \,+ \, Q^v(x- \overline{v}),
\end{align}
where $Q^v = Q^\xi (\lambda I - Y)^{-1}$. Here, $\overline{u} - Q^v \overline{v} \in {\bf R}^{N_u}$ is the open-loop component of the control policy, and $Q^w \in {\mathcal{Q}}_N$ and  $Q^v \in {\mathcal{Q}}_C$  are the feedback control gains.  The satisfaction of the sparsity constraint $Q^v \in {\mathcal{Q}}_C$ is guaranteed by  Lemma \ref{lem:change_of_variable}.


\section{Introduction}
The general problem of learning dynamical models from a time series of observations has a long history spanning many fields
\cite{Nelles_Book_Nonlinear_Sys_Identification_2013,
Sjoeberg_Black_Box_Nonlinear_Dyn_1995,
Chiuso_System_Identification_Review_2019,
Hong_System_Identification_Review_2008}
including in dynamical systems
\cite{Kutz_Brunton_book_ch_ROMs_2019,
Sjoeberg_Black_Box_Nonlinear_Dyn_1995,
Mallet_Coifman_Manifold_Learning_LVM_Dyn_Sys_2015,
Kutz_Lusch_Nonlinear_Embeddings_2018,
Mezic_Koopman_Review_2013,
Ohlberger_Redcuced_Basis_Review_2016,Hesthaven2016,
Crutchfield_Dynamics_Symbolic_1987,
DeVore_Reduced_Basis_Methods_Ch_2017},
control~\cite{Brunton_Kutz_SINDy_2016,
Nelles_Book_Nonlinear_Sys_Identification_2013,
Kutz_DNN_Time_Step_Constraints_2018,DMD_Schmid_2010},
statistics
\cite{Archer_Variational_State_Space_2015,
Jordan_Nonlinear_Dynamics_Theory_2020,
Ghahramani_EM_Nonlinear_Dynamics_1998},
and
machine learning
\cite{Chiuso_System_Identification_Review_2019,
Hong_System_Identification_Review_2008,
Carlberg_Lee_Nonlinear_Dynamics_AE_2020,
Karniadakis_Raissi_Hidden_Phys_PDEs_2018,
Bertozzi_CDMD_2019,
Perdikaris_Phys_Informed_Gen_Models_2018}.  Referred to as system identification in control and engineering, many approaches have been developed starting with linear dynamical systems (LDS).
These includes the Kalman Filter and extensions~\cite{Kalman1960,
DelMoral1997,
Godsill2019,
VanDerMerwe_Kalman_Unscented_2000,
Wan_Kalman_Nonlinear_2000}, Principle Orthogonal Decomposition (POD)~\cite{POD_Intro_Chatterjee_2000,Mendez2018}, and more recently Dynamic Mode Decomposition (DMD)
\cite{DMD_Schmid_2010,
DMD_Kutz_Brunton_Book_2016,
DMD_Theory_and_App_Kutz_2014}
and Koopman Operator approaches
\cite{Mezic_Koopman_Review_2013,Das2019,
Putinar_Mezic_cont_spectrum_Koopman_2020}.
These successful and widely-used approaches rely on assumptions on the model structure, most commonly, that a time-invariant LDS provides a good local approximation or that noise is Gaussian.

There also has been research on
more general nonlinear system identification
\cite{Archer_Variational_State_Space_2015,
Ljung_Nonlinear_Sys_Iden_Review_2019,
Chiuso_System_Identification_Review_2019,
Hong_System_Identification_Review_2008,
Schoen_ML_State_Space_Est_2011,
Kutz_Lusch_Nonlinear_Embeddings_2018,
Jordan_Nonlinear_Dynamics_Theory_2020,
Nelles_Book_Nonlinear_Sys_Identification_2013}.
Nonlinear systems pose many open challenges and fewer unified approaches given the rich behaviors of nonlinear dynamics.  For classes of systems and specific application domains, methods have been developed which make different levels of assumptions about the underlying structure of the dynamics.  Methods for learning nonlinear dynamics include the NARAX and NOE approaches with function approximators based on neural networks and other models classes
\cite{Nelles_Book_Nonlinear_Sys_Identification_2013,
Sjoeberg_Black_Box_Nonlinear_Dyn_1995},
sparse symbolic dictionary methods that are linear-in-parameters such as SINDy
\cite{Brunton_Kutz_SINDy_2016,
Schmidt2009,
Sjoeberg_Black_Box_Nonlinear_Dyn_1995},
and dynamic Bayesian networks (DBNs), such as
Hidden Markov Chains (HMMs) and Hidden-Physics Models
~\cite{Karniadakis_Raissi_Hidden_Phys_PDEs_2018,
Pawar2020,Saul2020,Baum_HMM_1966,Krishnan_GMS_2017,
Ghahramani_EM_Nonlinear_Dynamics_1998}.

A central challenge in learning non-linear dynamics is to obtain representations not only capable of reproducing similar outputs as observed directly in the training dataset but to infer structures that can provide stable more long-term extrapolation capabilities over multiple future steps and input states.
In this work, we develop learning methods aiming to obtain robust non-linear models by providing ways to incorporate more structure and information about the underlying system related to smoothness, periodicity, topology, and other constraints.  We focus particularly on developing Probabilistic Autoencoders (PAE) that incorporate noise-based regularization and priors to learn lower dimensional representations from observations.  This provides the basis of nonlinear state space models for prediction.  We develop methods for incorporating into such representations geometric and topological information about the system.  This facilitates capturing qualitative features of the dynamics to enhance robustness and to aid in interpretability of results.  We demonstrate and perform investigations of our methods to obtain models for reductions of parameterized PDEs and for constrained mechanical systems.

\section*{Learning Nonlinear Dynamics with \\ Variational Autoencoders (VAEs)}
\label{sec_VAR}

We develop data-driven approaches based on a Variational Autoencoder (VAE) framework~\cite{KingmaWellingVAE2014}.  We learn from observation data a set of lower dimensional representations that are used to make predictions for the dynamics.  In practice, data can include experimental measurements, large-scale computational simulations, or solutions of complicated dynamical systems for which we seek reduced models.  Reductions aid in gaining insights for a class of inputs or physical regimes into the underlying mechanisms generating the observed behaviors.  Reduced descriptions are also helpful in many optimization problems in design and in development of controllers
\cite{Nelles_Book_Nonlinear_Sys_Identification_2013}.

Standard autoencoders can result in encodings that yield
unstructured scattered disconnected coding points for system features $\mb{z}$.  VAEs provide probabilistic encoders and decoders where noise provides regularizations that promote more connected encodings, smoother dependence on inputs, and more disentangled feature components~\cite{KingmaWellingVAE2014}.  As we shall discuss, we also introduce other regularizations into our methods to help aid in interpretation of the learned latent representations.

\begin{figure}[ht]
\centerline{\includegraphics[width=0.99\columnwidth]{./fig/fig_encoder_decoder4.png}}
\caption{\textbf{Learning Nonlinear Dynamics.}  Data-driven methods are developed for learning robust models to predict from  $u(x,t)$ the non-linear evolution to $u(x,t + \tau)$ for PDEs and other dynamical systems.  Probabilistic Autoencoders (PAEs) are utilized to learn representations $z$ of $u(x,t)$ in low dimensional latent spaces with prescribed geometric and topological properties.  The model makes predictions using learnable maps that (i) encode an input $u(x,t) \in \mathcal{U}$ as $z(t)$ in latent space \textit{(top)}, (ii) evolve the representation $z(t) \rightarrow z(t + \tau)$ \textit{(top-right)}, (iii) decode the representation $z(t + \tau)$ to predict $\hat{u}(x,t + \tau)$ \textit{(bottom-right)}.
}
\label{fig:vae_learning_dynamics}
\end{figure}

We learn VAE predictors using a Maximum Likelihood Estimation (MLE) approach for the Log Likelihood (LL) $\mathcal{L}_{LL} = \log(p_\theta(\mb{X},\mb{x}))$.  For dynamics of $u(s)$, let $\mb{X} = u(t)$ and $\mb{x} = u(t + \tau)$.  We base $p_\theta$ on the autoencoder framework in Figure~\ref{fig:vae_learning_dynamics} and~\ref{VAE_schematic}.  We use variational inference to approximate the LL by the Evidence Lower Bound (ELBO)~\cite{Blei_VI_Review_2017} to train a model with parameters $\theta$ using encoders and decoders based on minimizing the loss function
\begin{eqnarray}
\nonumber
\theta^* &=& \arg\min_{\theta_e,\theta_d} -\mathcal{L}^B(\theta_e,\theta_d,\theta_\ell;\mb{X}^{(i)},\mb{x}^{(i)}), \\
\label{equ:vae_loss}
\mathcal{L}^B &=&  \mathcal{L}_{RE} + \mathcal{L}_{KL} + \mathcal{L}_{RR}, \\
\nonumber
\mathcal{L}_{RE} &=&
E_{\mathfrak{q}_{\theta_e}(\mb{z}|\mb{X}^{(i)})}\left\lbrack \log \mathfrak{p}_{\theta_d}(\mb{x}^{(i)} | \mb{z}') \right\rbrack  \\
\nonumber
\mathcal{L}_{KL} &=&
-\beta\mathcal{D}_{KL}\left(\mathfrak{q}_{\theta_e}(\mb{z}|\mb{X}^{(i)}) \, \| \,
\tilde{\mathfrak{p}}_{\theta_d}(\mb{z})\right) \\
\nonumber
\mathcal{L}_{RR} &=&
\gamma
E_{\mathfrak{q}_{\theta_e}(\mb{z}'|\mb{x}^{(i)})}\left\lbrack \log \mathfrak{p}_{\theta_d}(\mb{x}^{(i)} | \mb{z}') \right\rbrack.
\end{eqnarray}
The $\mathfrak{q}_{\theta_e}$ denotes the encoding probability distribution and $\mathfrak{p}_{\theta_d}$ the decoding probability distribution.  The loss $\ell = -\mathcal{L}^B$ provides a regularized form of MLE.

\begin{figure}[ht]
\centerline{\includegraphics[width=0.99\columnwidth]{./fig/vae_schematic5.png}}
\caption{\textbf{Variational Autoencoder (VAE).}  VAEs~\cite{KingmaWellingVAE2014} are used to learn representations of the nonlinear dynamics.  Deep Neural Networks (DNNs) are trained (i) to serve as feature extractors to represent functions $u(x,t)$ and their evolution in a low dimensional latent space as $z(t)$  (encoder $\sim \mathfrak{q}_{\theta_e}$), and (ii) to serve as approximators that can construct predictions $u(x,t + \tau)$ using features $z(t+ \tau)$ (decoder $\sim \mathfrak{p}_{\theta_d}$).
}
\label{VAE_schematic}
\end{figure}

The terms $\mathcal{L}_{RE}$ and $\mathcal{L}_{KL}$ arise from the ELBO variational bound $\mathcal{L}_{LL} \geq \mathcal{L}_{RE} + \mathcal{L}_{KL}$ when $\beta = 1$,~\cite{Blei_VI_Review_2017}. This provides a way to estimate the log likelihood that the encoder-decoder reproduce the observed data sample pairs $(\mb{X}^{(i)},\mb{x}^{(i)})$ using the codes $\mb{z'}$ and $\mb{z}$.  Here, we include a latent-space mapping $\mb{z}' = f_{\theta_\ell}(\mb{z})$ parameterized by $\theta_\ell$, which we can use to characterize the evolution of the system or further processing of features. The $\mb{X}^{(i)}$ is the input and $\mb{x}^{(i)}$ is the output prediction.  For the case of dynamical systems, we take $\mb{X}^{(i)} \sim u^i(t)$ a sample of the initial state function $u^i(t)$ and the output $\mb{x}^{(i)} \sim u^i(t + \tau)$ the predicted state function $u^i(t + \tau)$.  We discuss the specific distributions used in more detail below.

The $\mathcal{L}_{KL}$ term involves the Kullback-Leibler Divergence
\cite{Kullback1951,Book_Cover_Thomas_2006} acting similar to a Bayesian prior on latent space to regularize the encoder conditional probability distribution so that for each sample this distribution is similar to $p_{\theta_d}$.  We take $p_{\theta_d} = \eta(0,\sigma_0^2)$ a multi-variate Gaussian with independent components.  This serves (i) to disentangle the features from each other to promote independence, (ii) provide a reference scale and localization for the encodings $\mb{z}$, and (iii) promote parsimonious codes utilizing smaller dimensions than $d$ when possible.

The $\mathcal{L}_{RR}$ term gives a regularization that promotes retaining information in $z$ so the encoder-decoder pair can reconstruct functions.  As we shall discuss, this also promotes organization of the latent space for consistency over multi-step predictions and aids in model interpretability.

We use for the specific encoder probability distributions conditional Gaussians
$\mb{z} \sim \mathfrak{q}_{\theta_e}(\mb{z}|\mb{x}^{(i)}) = \mathfrak{a}(\mb{X}^{(i)},\mb{x}^{(i)}) + \eta(0,\sigma_e^2)$ where
$\eta$ is a Gaussian with variance $\sigma_e^2$, (i.e. $\mathbb{E}^{\mb{X}^i}[z] = \mathfrak{a}$,
$\mbox{Var}^{\mb{X}^i}[z] = \sigma_e^2$). One can think of the learned mean function  $\mathfrak{a}$ in the VAE as corresponding to a typical
encoder
$\mathfrak{a}(\mb{X}^{(i)},\mb{x}^{(i)};\theta_e) = \mathfrak{a}(\mb{X}^{(i)};\theta_e) = \mb{z}^{(i)}$
and the variance function $\sigma_e^2 = \sigma_e^2(\theta_e)$ as providing control of a noise source to further regularize the encoding.  Among other properties, this promotes connectedness of the ensemble of latent space codes.  For the VAE decoder distribution, we take $\mb{x} \sim \mathfrak{p}_{\theta_d}(\mb{x}|\mb{z}^{(i)}) = \mathfrak{b}(\mb{z}^{(i)}) + \eta(0,\sigma_d^2)$.  The learned mean function $\mathfrak{b}(\mb{z}^{(i)};\theta_e)$
corresponds to a typical decoder and the variance function $\sigma_e^2 = \sigma_e^2(\theta_d)$ controls the source of regularizing noise.

The terms to be learned in the VAE framework are $(\mathfrak{a},\sigma_e,f_{\theta_\ell},\mathfrak{b},\sigma_d)$ which are parameterized by $\theta = (\theta_e,\theta_d,\theta_\ell)$.  In practice, it is useful to treat variances $\sigma_{(\cdot)}$ initially as hyper-parameters.
We learn predictors for the dynamics by training over samples of evolution pairs $\{(u_n^i,u_{n+1}^i)\}_{i=1}^m$, where $i$ denotes the sample index and
$u_n^i = u^i(t_n)$ with $t_n = t_0 + n\tau$ for a time-scale $\tau$.

To make predictions, the learned models use the following stages: (i) extract from $u(t)$ the features $z(t)$, (ii) evolve $z(t) \rightarrow z(t + \tau)$, (iii) predict
using $z(t + \tau)$ the $\hat{u}(t + \tau)$, summarized in Figure~\ref{fig:vae_learning_dynamics}.  By composition of the latent evolution map the model makes multi-step predictions of the dynamics.

\section{Learning with Manifold Latent Spaces \\ Roles of Non-Euclidean Geometry and Topology}
For many systems, parsimonious representations can be obtained by working with non-euclidean manifold latent spaces, such as a torus for doubly periodic systems or even non-orientable manifolds, such as a klein bottle as arises in imaging and perception studies~\cite{Carlsson2008}.  For this purpose, we learn encoders $\mathcal{E}$ over a family of mappings to a prescribed manifold $\mathcal{M}$ of the form
$$\mb{z} = \mathcal{E}_\phi(\mb{x}) = \Lambda(\tilde{\mathcal{E}}_{\phi}(\mb{x})) = \Lambda(\mb{w}), \; \; \; \mb{w} = \tilde{\mathcal{E}}_{\phi}(\mb{x}).$$  We take the map $\tilde{\mathcal{E}}_{\phi}(\mb{x}): \mb{x} \rightarrow \mb{w}$, where we represent a smooth closed manifold $\mathcal{M}$ of dimension $m$ in $\mathbb{R}^{2m}$, as supported by the Whitney Embedding Theorem
\cite{Whitney1944}. The $\Lambda$ maps (projects) points $\mb{w} \in \mathbb{R}^{2m}$ to the manifold representation $\mb{z} \in \mathcal{M} \subset \mathbb{R}^{2m}$.  In practice, we accomplish this two ways: (i) we provide an analytic mapping $\Lambda$ to $\mathcal{M}$, (ii) we provide a high resolution point-cloud representation of the target manifold along with local gradients and use for $\Lambda$ a quantized mapping to the nearest point on $\mathcal{M}$.  We provide more details in Appendix A.

This allows us to learn VAEs with latent spaces for $\mb{z}$ with general specified topologies and controllable geometric structures.  The topologies of sphere, torus, klein bottle are intrinsically different than $\mathbb{R}^n$.  This allows for new types of priors such as uniform on compact manifolds or distributions with more symmetry.
As we shall discuss, additional latent space structure also helps in learning more robust representations less sensitive to noise since we can unburden the encoder and decoder from having to learn the embedding geometry and avoid the potential for them making erroneous use of  extra latent space dimensions.  We also have statistical gains since the decoder now only needs to learn a mapping from the manifold $\mathcal{M}$ for reconstructions of $\mb{x}$.  These more parsimonious representations also aid identifiability and interpretability of models.

\section{Related Work}
Many variants of autoencoders have been developed for making predictions of sequential data, including those based on
Recurrent Neural Networks (RNNs) with LSTMs and GRUs~\cite{Schmidhuber_LSTM_1997,
Goodfellow2016,Cho_Bengio_GRU_2014}.  While RNNs provide a rich approximation class for sequential data, they pose for dynamical systems challenges for interpretability and for training to obtain predictions stable over many steps with robustness against noise in the training dataset.  Autoencoders have also been combined with symbolic dictionary learning
for latent dynamics in \cite{Kutz_Discovery_Coordinates_2019} providing some advantages for interpretability and robustness, but require specification in advance of a sufficiently expressive dictionary.  Neural networks incorporating physical information have also been developed that impose stability conditions during training
\cite{Carlberg_RNN_Dynamics_2020,
Carlberg_Lee_Nonlinear_Dynamics_AE_2020,
Erichson_AE_Lyapunov_Stable_Flow_2019}.  The work of \cite{Bengio_Recurrent_LVM_2015}
investigates combining RNNs with VAEs to obtain more robust models for sequential data and considered tasks related to processing speech and handwriting.

In our work we learn dynamical models making use of VAEs to obtain probabilistic encoders and decoders between euclidean and non-euclidean latent spaces to provide additional regularizations to help promote parsimoniousness, disentanglement of features, robustness, and interpretability.  Prior VAE methods used for dynamical systems include
\cite{Hernandez_VAE_Dynamics_2018,
Pearce_Latent_Dynamics_GP_2020,
Girin_Dynamical_VAE_Review_2020,
Chen_VAE_Dynamic_Motion_Primitives_2016,
Pearce_Latent_Dynamics_GP_2020,
Roeder_VAE_Dynamics_Hierarchical_2019}.
These works use primarily euclidean latent spaces and consider applications including human motion capture and ODE systems.
Approaches for incorporating topological information into latent variable representations include the early works by Kohonen on Self-Organizing Maps (SOMs)
\cite{Kohonen_Original_Paper_Self_Organizing_Maps_1982}
and Bishop on Generative Topographical Maps (GTMs)
based on density networks providing a generative approach
\cite{Bishop_Early_Generative_Topographical_Map_1996}.  More recently, VAE methods using non-euclidean latent spaces include
\cite{Jensen_Manifold_Latent_Model_2020,
Kalatzis_VAE_B_Motion_Priors_2020,
Falorsi_Homeomorphic_VAE_2018,
Chen_VAE_Learn_Flat_Manifolds_2020,
Kipf_VAE_Hyperspherical_2018,
Arvanitidis_Noneuclidean_Latent_Space_2018}.
These incorporate the role of geometry by augmenting the prior distribution $\tilde{\mathfrak{p}}_{\theta_d}(z)$ on latent space to bias toward a manifold.  In the recent work
\cite{Rey_Diffusion_VAE_2020}, an explicit projection procedure is introduced, but in the special case of a few manifolds having an analytic projection map.

In our work we develop further methods for more general latent space representations, including non-orientable manifolds, and applications to parameterized PDEs and constrained mechanical systems.  We introduce more general methods for non-euclidean latent spaces in terms of point-cloud representations of the manifold along with local gradient information that can be utilized within general back-propogation frameworks, see Appendix A.  This also allows for the case of manifolds that are non-orientable and having complex shapes.  Our methods provide flexible ways to design and control both the topology and the geometry of the latent space by merging or subtracting shapes or stretching and contracting regions.  We also consider additional types of regularizations for learning dynamical models facilitating multi-step predictions and more interpretable state space models.  In our work, we also consider reduced models for non-linear PDEs, such as Burgers Equations, and learning representations for more general constrained mechanical systems.  We also investigate the role of non-linearities making comparisons with other data-driven models.

\section*{Results}
\label{sec_Results}

\subsection*{Burgers' Equation of Fluid Mechanics: Learning Nonlinear PDE Dynamics}

We consider the nonlinear viscous Burgers' equation
\begin{eqnarray}
u_t = -uu_x + \nu u_{xx},
\label{eqn:burgers}
\end{eqnarray}
where $\nu$ is the viscosity~\cite{Bateman1915,Hopf1950}.  We consider periodic boundary conditions on $\Omega = [0,1]$.  Burgers equation is motivated as a mechanistic model for the fluid mechanics of advective transport and shocks, and serves as a widely used benchmark for analysis and computational methods.

The nonlinear Cole-Hopf Transform $\mathcal{CH}$ can be used to relate Burgers equation to the linear Diffusion equation $\phi_t = \nu \phi_{xx}$~\cite{Hopf1950}.  This provides a representation of the solution $u$
\begin{eqnarray}
\nonumber
\phi(x,t) & = & \mathcal{CH}[u] = \exp \left( -\frac{1}{2\nu}\int_0^x u(x',t)dx'\right) \\
u(x,t) & = & \mathcal{CH}^{-1}[\phi] = -2\nu \frac{\partial}{\partial x} \ln \phi(x,t).
\end{eqnarray}
This can be represented by the Fourier expansion
\begin{eqnarray}
\nonumber
&& \phi(x,t) = \sum_{k=-\infty}^{\infty} \hat{\phi}_k(0) \exp(-4\pi^2 k^2 \nu t) \cdot \exp(i2\pi kx).
\end{eqnarray}
The $\hat{\phi}_k(0) = \mathcal{F}_k[\phi(x,0)]$ and $\phi(x,t) = \mathcal{F}^{-1}[\{\hat{\phi}_k(0) \exp(-4\pi^2 k^2 \nu t)\}]$ with $\mathcal{F}$ the Fourier transform. This provides an analytic representation of the solution of the viscous Burgers equation $u(x,t) = \mathcal{CH}^{-1}[\phi(x,t)]$ where $\hat{\phi}(0) = \mathcal{F}[\mathcal{CH}[u(x,0)]]$.  In general, for nonlinear PDEs with initial conditions within a class of functions $\mathcal{U}$, we aim to learn models that provide predictions $u(t + \tau) = \mathcal{S}_\tau u(t)$ approximating the evolution operator $\mathcal{S}_\tau$ over time-scale $\tau$.  For the Burgers equation, the $\mathcal{CH}$ provides an analytic way to obtain a reduced order model by truncating the Fourier expansion to $|k| \leq n_f/2$.  This provides for the Burgers equation a benchmark model against which to compare our learned models. For general PDEs comparable analytic representations are not usually available, motivating development of data-driven approaches.

We develop VAE methods for learning reduced order models for the responses of nonlinear Burgers Equation when the initial conditions are from a collection of functions $\mathcal{U}$.  We learn VAE models that extract from
$u(x,t)$ latent variables $z(t)$ to predict $u(x,t + \tau)$.  Given the non-uniqueness of representations and to promote interpretability of the model, we introduce the inductive bias that the evolution dynamics in latent space for $z$ is linear of the form $\dot{z} = -\lambda_0 z$, giving exponential decay rate $\lambda_0$.  For discrete times, we take $z_{n+1} = f_{\theta_\ell}(z_n) = \exp(-\lambda_0 \tau) \cdot z_n$, where $\theta_\ell = (\lambda_0)$.  We still consider general nonlinear mappings for the encoders and decoders which are represented by deep neural networks.  We train the model on the pairs $(u(x,t),u(x,t + \tau))$ by drawing $m$ samples of $u^{i}(x,t_i) \in \mathcal{S}_{t_i} \mathcal{U}$ which generates the evolved state under Burgers equation $u^i(x,t_i + \tau)$ over time-scale $\tau$.
We perform VAE studies with parameters $\nu = 2 \times 10^{-2}$, $\tau = 2.5 \times 10^{-1}$ with VAE Deep Neural Networks (DNNs) with layer sizes (in)-400-400-(out), ReLU activations, and $\gamma = 0.5$, $\beta = 1$, and initial standard deviations $\sigma_d = \sigma_e = 4 \times 10^{-3}$.
We show results of our VAE model predictions in Figure~\ref{fig:Burgers_compare_VAE_DMD_etc} and Table~\ref{table:Burgers_compare_VAE_DMD_etc}.

\begin{figure}[ht]
\centerline{\includegraphics[width=0.99\columnwidth]{./fig/comparison_4.png}}
\caption{\textbf{Burgers' Equation: Prediction of Dynamics.}  We consider responses for $\mathcal{U}_1 = \{ u\, | \, u(x,t;\alpha) = \alpha \sin(2\pi x) + (1-\alpha) \cos^3(2 \pi x)\}$.  Predictions are made for the evolution $u$ over the time-scale $\tau$ satisfying equation~\ref{eqn:burgers} with initial conditions in $\mathcal{U}_1$.  We find our nonlinear VAE methods are able to learn with $2$ latent dimensions the dynamics with errors $< 1\%$. Methods such as DMD~\cite{DMD_Schmid_2010,DMD_Theory_and_App_Kutz_2014} with $3$ modes which are only able to use a single linear space to approximate the initial conditions and prediction encounter challenges in approximating the nonlinear evolution.  We find our linear VAE method with $2$ modes provides some improvements, by allowing for using different linear spaces for representing the input and output functions, but at the cost of additional computations.  Results are summarized in Table~\ref{table:Burgers_compare_VAE_DMD_etc}.}
\label{fig:Burgers_compare_VAE_DMD_etc}
\end{figure}

\begin{figure}[ht]
\centerline{\includegraphics[width=0.99\columnwidth]{./fig/manifold_gamma_vs_no_gamma_6.png}}
\caption{\textbf{Burgers' Equation: Latent Space Representations and Extrapolation Predictions.}  We show the latent space representation $z$ of the dynamics for the input functions $u(\cdot,t;\alpha) \in \mathcal{U}_1$.  VAE organizes for $u$ the learned representations $z(\alpha,t)$  in parameter $\alpha$ \textit{(blue-green)} into circular arcs that are concentric in the time parameter $t$, \textit{(yellow-orange)} \textit{(left)}.  The reconstruction regularization with $\gamma$ aligns subsequent time-steps of the dynamics in latent space facilitating multi-step predictions.  The learned VAE model exhibits a level of extrapolation to predict dynamics even for some inputs $u \not \in \mathcal{U}_1$ beyond the training dataset \textit{(right)}.
}
\label{fig:burgers_latent_space}
\label{fig:Burgers_compare_VAE_DMD_etc_outside_training}
\end{figure}

We show the importance of the non-linear approximation properties of our VAE methods in capturing system behaviors by making comparisons with Dynamic Mode Decomposition (DMD)~\cite{DMD_Schmid_2010,DMD_Theory_and_App_Kutz_2014}, Principle Orthogonal Decomposition (POD)~\cite{POD_Intro_Chatterjee_2000}, and a linear variant of our VAE approach.  Recent CNN-AEs have also studied related advantages of non-linear approximations~\cite{Carlberg_Lee_Nonlinear_Dynamics_AE_2020}.  Some distinctions in our work is the use of VAEs to further regularize AEs and using topological latent spaces to facilitate further capturing of structure.  The DMD and POD are widely used and successful approaches that aim to find an optimal linear space on which to project the dynamics and learn a linear evolution law for system behaviors.  DMD and POD have been successful in obtaining models for many applications, including steady-state fluid mechanics and transport problems
\cite{DMD_Theory_and_App_Kutz_2014,
DMD_Schmid_2010}.  However, given their inherent linear approximations they can encounter well-known challenges related to translational and rotational invariances, as arise in advective phenomena and other settings~\cite{Kutz_Brunton_book_ch_ROMs_2019}.    Our comparison studies can be found in Table~\ref{table:Burgers_compare_VAE_DMD_etc}.

\begin{table}[ht]
\begin{center}
{
\fontsize{8.3}{10}\selectfont
\selectfont
\begin{tabular}{l|c|cccc}
\hline
\rowcolor{atz_table1}
\textbf{Method} & \textbf{Dim} & \textbf{0.25s} & \textbf{0.50s} & \textbf{0.75s} & \textbf{1.00s} \\
\cline{1-6}
VAE Nonlinear  & 2 & \textbf{4.44e-3} & \textbf{5.54e-3} & \textbf{6.30e-3} & \textbf{7.26e-3} \\
VAE Linear  & 2 & 9.79e-2 & 1.21e-1 & 1.17e-1 & 1.23e-1\\
DMD & 3 & 2.21e-1 & 1.79e-1 & 1.56e-1 & 1.49e-1 \\
POD & 3 & 3.24e-1 & 4.28e-1 & 4.87e-1 & 5.41e-1\\
\cline{1-6}
Cole-Hopf-2 & 2 &5.18e-1 & 4.17e-1 &3.40e-1 & 1.33e-1\\
Cole-Hopf-4 & 4 & 5.78e-1 & 6.33e-2 & 9.14e-3 & 1.58e-3\\
Cole-Hopf-6 & 6 &1.48e-1 & 2.55e-3 & 9.25e-5 & 7.47e-6\\
\hline
\end{tabular} \\
}
\vspace{0.2cm}
{
\fontsize{7.6}{10}\selectfont
\begin{tabular}{c|ccccc}
\hline
\rowcolor{atz_table1}
$\bsy{\gamma}$ & \textbf{0.00s} & \textbf{0.25s} & \textbf{0.50s} & \textbf{0.75s} & \textbf{1.00s} \\
\cline{1-6}
0.00 & 1.600e-01 & 6.906e-03 & 1.715e-01 & 3.566e-01 & 5.551e-01\\
0.50 & 1.383e-02 & 1.209e-02 & 1.013e-02 & 9.756e-03 & 1.070e-02\\
2.00 & 1.337e-02 & 1.303e-02 & 9.202e-03 & 8.878e-03 & 1.118e-02\\
\hline
\end{tabular}
} \\
\vspace{0.2cm}
{
\fontsize{7.6}{10}\selectfont
\begin{tabular}{c|ccccc}
\hline
\rowcolor{atz_table1}
$\bsy{\beta}$ & \textbf{0.00s} & \textbf{0.25s} & \textbf{0.50s} & \textbf{0.75s} & \textbf{1.00s} \\
\cline{1-6}
0.00 & 1.292e-02 & 1.173e-02 & 1.073e-02 & 1.062e-02 & 1.114e-02\\
0.50 & 1.190e-02 & 1.126e-02 & 1.072e-02 & 1.153e-02 & 1.274e-02\\
1.00 & 1.289e-02 & 1.193e-02 & 7.903e-03 & 7.883e-03 & 9.705e-03\\
4.00 & 1.836e-02 & 1.677e-02 & 8.987e-03 & 8.395e-03 & 8.894e-03\\
\hline
\end{tabular}
}
\end{center}
\caption{\textbf{Burgers' Equation: Prediction Accuracy.}
The reconstruction $L^1$-relative errors in predicting $u(x,t)$ for our VAE methods, Dynamic Model Decomposition (DMD), and Principle Orthogonal Decomposition (POD), and reduction by Cole-Hopf (CH), over multiple-steps and number of latent dimensions (Dim)
\textit{(top)}.  Results when varying the strength of the reconstruction regularization $\gamma$ and prior $\beta$ \textit{(bottom)}.
}
\label{table:Burgers_compare_VAE_DMD_etc}
\end{table}

We also considered how our VAE methods performed when adjusting the parameters $\beta$ for the strength of the prior $\tilde{\mathfrak{p}}$ as in $\beta$-VAEs~\cite{Higgins2017BetaVAE} and
$\gamma$ for the strength of the reconstruction regularization.  The reconstruction regularization has a significant influence
on how the VAE organizes representations in latent space and the accuracy of predictions of the dynamics, especially over multiple steps, see Figure~\ref{fig:Burgers_compare_VAE_DMD_etc_outside_training} and Table~\ref{table:Burgers_compare_VAE_DMD_etc}.     The regularization serves to align representations consistently in latent space facilitating multi-step compositions.  We also found our VAE learned representions capable of some level of extrapolation beyond the training dataset.
When varying $\beta$, we found that larger values improved the multiple step accuracy whereas small values improved the single step accuracy, see Table~\ref{table:Burgers_compare_VAE_DMD_etc}.

\subsection*{Constrained Mechanics: Learning with Non-Euclidean Latent Spaces}

\begin{figure}[ht]
\centerline{\includegraphics[width=0.99\columnwidth]{./fig/arm_schematic_torus_klein_combine2.png}}
\caption{\textbf{VAE Representations of Motions using Manifold Latent Spaces.}  We learn from observations representations for constrained mechanical systems using general non-euclidean manifolds latent spaces $\mathcal{M}$.  The arm mechanism has configurations $\mb{x} = (\mb{x}_1,\mb{x}_2) \in \mathbb{R}^4$.  For rigid segments, the motions are constrained to be on a manifold (torus) $\mathcal{M} \subset \mathbb{R}^4$.  For extendable segments, we can also consider more exotic constraints, such as requiring $\mb{x}_1,\mb{x}_2$ to be on a klein bottle in $\mathbb{R}^4$ \textit{(top)}.  Results of our VAE methods for learned representations for motions under these constraints are shown.  VAE learns the segment length constraint and two nearly decoupled coordinates for the torus dataset that mimic the roles of angles.  VAE learns for the klein bottle dataset two segment motions to generate configurations \textit{(middle and bottom)}.
}
\label{fig:manifold_lvms}
\end{figure}

\begin{table}[ht]
\begin{center}
\fontsize{7.3}{10}
\selectfont
\begin{tabular}{l|cccc}
\hline
\rowcolor{atz_table3}
\multicolumn{1}{l}{\textbf{Torus}} &
\multicolumn{4}{|l}{epoch} \\
\hline
\rowcolor{atz_table1}
method & 1000 & 2000 & 3000 & final\\
\hline
VAE 2-Manifold & \textbf{6.6087e-02} & \textbf{6.6564e-02} & \textbf{6.6465e-02} &
\textbf{6.6015e-02}
\\
VAE $\mathbb{R}^2$ & 1.6540e-01 & 1.2931e-01 & 9.9903e-02 & 8.0648e-02 \\
VAE $\mathbb{R}^4$ & 8.0006e-02 & 7.6302e-02 & 7.5875e-02 & 7.5626e-02\\
VAE $\mathbb{R}^{10}$ & 8.3411e-02 & 8.4569e-02 & 8.4673e-02 & 8.4143e-02\\
\hline
\rowcolor{atz_table1}
with noise $\sigma$ & 0.01 & 0.05 & 0.1 & 0.5 \\
\hline
VAE 2-Manifold & \textbf{6.7099e-02} & \textbf{8.0608e-02} & \textbf{1.1198e-01} & \textbf{4.1988e-01} \\
VAE $\mathbb{R}^2$ & 8.5879e-02 & 9.7220e-02 & 1.2867e-01 & 4.5063e-01 \\
VAE $\mathbb{R}^4$ & 7.6347e-02 & 9.0536e-02 & 1.2649e-01 & 4.9187e-01 \\
VAE $\mathbb{R}^{10}$ & 8.4780e-02 & 1.0094e-01 & 1.3946e-01 & 5.2050e-01 \\
\hline
\rowcolor{atz_table3}
\multicolumn{1}{l}{\textbf{Klein Bottle}} &
\multicolumn{4}{|l}{epoch} \\
\hline
\rowcolor{atz_table1}
method & 1000 & 2000 & 3000 & final\\
\cline{1-5}
VAE 2-Manifold & \textbf{5.7734e-02} & \textbf{5.7559e-02} & \textbf{5.7469e-02} & \textbf{5.7435e-02} \\
VAE $\mathbb{R}^2$ & 1.1802e-01 & 9.0728e-02 & 8.0578e-02 & 7.1026e-02 \\
VAE $\mathbb{R}^4$ & 6.9057e-02 & 6.5593e-02 & 6.4047e-02 & 6.3771e-02 \\
VAE $\mathbb{R}^{10}$ & 6.8899e-02 & 6.9802e-02 & 7.0953e-02 & 6.8871e-02 \\
\cline{1-5}
\hline
\rowcolor{atz_table1}
with noise $\sigma$ & 0.01 & 0.05 & 0.1 & 0.5 \\
\cline{1-5}
VAE 2-Manifold & \textbf{5.9816e-02} & \textbf{6.9934e-02} & \textbf{9.6493e-02} & \textbf{4.0121e-01} \\
VAE $\mathbb{R}^2$ & 1.0120e-01 & 1.0932e-01 & 1.3154e-01 & 4.8837e-01 \\
VAE $\mathbb{R}^4$ & 6.3885e-02 & 7.6096e-02 & 1.0354e-01 & 4.5769e-01 \\
VAE $\mathbb{R}^{10}$ & 7.4587e-02 & 8.8233e-02 & 1.2082e-01 & 4.8182e-01 \\
\end{tabular}
\end{center}
\caption{\textbf{Manifold Latent Variable Model: VAE Reconstruction Errors} The $L^2$-relative errors of reconstruction for our VAE methods.  The final is the lowest value during training.  The manifold latent spaces show improved learning.  When an incompatible topology is used, such as $\mathbb{R}^2$, this can result in
deterioration in learned representations. With noise in the input $\tilde{X} = X + \sigma \eta(0,1)$ and reconstructing the target $X$, the manifold latent spaces also show improvements for learning.
}
\label{table:VAE_Manifold_Reconstruction}
\end{table}

To learn more parsimonous and robust representations of physical systems, we develop methods for latent spaces having geometries and topologies more general than euclidean space.  This is helpful in capturing inherent structure such as periodicities or other symmetries.
We consider physical systems with constrained mechanics, such as the arm mechanism for reaching for objects in figure~\ref{fig:manifold_lvms}.  The observations are taken to be the two locations $\mb{x}_1, \mb{x}_2 \in \mathbb{R}^2$ giving $\mb{x} = (\mb{x}_1,\mb{x}_2) \in \mathbb{R}^4$.  When the segments are rigidly constrained these configurations lie on a manifold (torus).  We can also allow the segments to extend and consider more exotic constraints such as the two points $\mb{x}_1, \mb{x}_2$ must be on a klein bottle in $\mathbb{R}^4$.  Related situations arise in other areas of imaging and mechanics, such as in pose estimation and in studies of visual perception
~\cite{Carlsson_Klein_Bottle_Images_2014,
Carlsson2008,
Sarafianos_Review_Pose_Estimation_2016}.  For the arm mechanics, we can use this prior knowledge to construct a torus latent space represented by the product space of two circles $S^1 \times S^1$.
To obtain a learnable class of manifold encoders, we use the family of maps
$\mathcal{E}_\theta = \Lambda(\tilde{\mathcal{E}}_\theta(x))$,
with $\tilde{\mathcal{E}}_\theta(x)$ into $\mathbb{R}^4$
and
$\Lambda(\mb{w}) = \Lambda(w_1,w_2,w_3,w_4) = (z_1,z_2,z_3,z_4) = \mb{z}$,
where $(z_1,z_2) = (w_1,w_2)/\|(w_1,w_2)\|$, $(z_3,z_4) = (w_3,w_4)/\|(w_3,w_4)\|$,
see VAE Section and Appendix A.  For the case of klein bottle constraints, we use our point-cloud representation of the non-orientable manifold with the parameterized embedding in $\mathbb{R}^4$
\begin{eqnarray*}
\begin{array}{ll}
z_1 = (a + b\cos(u_2))\cos(u_1) &
z_2 = (a + b\cos(u_2))\sin(u_1) \\
z_3 = b\sin(u_2)\cos\left(\frac{u_1}{2}\right) &
z_4 = b\sin(u_2)\sin\left(\frac{u_1}{2}\right),
\end{array}
\end{eqnarray*}
with $u_1,u_2 \in [0,2\pi]$.  The $\Lambda(\mb{w})$ is taken to be the map to the nearest point of the manifold $\mathcal{M}$, which we compute numerically along with the needed gradients for backpropogation as discussed in Appendix A.

Our VAE methods are trained with encoder and decoder DNN's having layers of sizes (in)-100-500-100-(out)
with Leaky-ReLU activations with s = 1e-6 with results reported in Figure~\ref{fig:manifold_lvms} and Table~\ref{table:VAE_Manifold_Reconstruction}.  We find
learning representations is improved by use of the manifold latent spaces, in these trials even showing a slight edge over $\mathbb{R}^4$.  When the wrong topology is used, such as in $\mathbb{R}^2$, we find in both cases a significant deterioration in the reconstruction accuracy, see Table~\ref{table:VAE_Manifold_Reconstruction}.  This arises since the encoder must be continuous and hedge against the noise regularizations.  This results in an incurred penalty for a subset of configurations.  The encoder exhibits non-injectivity and a rapid doubling back over the space to accommodate the decoder by lining up nearby configurations in the topology of the input space manifold to handle noise perturbations in $z$ from the probabilistic nature of the encoding.  We also studied robustness when training with noise for $\tilde{X} = X + \sigma\eta(0,1)$ and measuring accuracy for reconstruction relative to target $X$.  As the noise increases, we see that the manifold latent spaces improve reconstruction accuracy acting as a filter through restricting the representation.  The probabilistic decoder will tend to learn to estimate the mean over samples of a common underlying configuration and with the manifold latent space restrictions is more likely to use a common latent representation.  For $\mathbb{R}^d$ with $d > 2$, the extraneous dimensions in the latent space can result in overfitting of the encoder to the noise.  We see as $d$ becomes larger the reconstruction accuracy decreases, see Table~\ref{table:VAE_Manifold_Reconstruction}.  These results demonstrate how geometric priors can aid learning in constrained mechanical systems.

\section{Conclusions}
We developed VAE's for learning robustly nonlinear dynamics of physical systems by introducing methods for latent representations utilizing general geometric and topological structures.  We demonstrated our methods for learning the non-linear dynamics of PDEs and constrained mechanical systems.  We expect our methods can also be used in other physics-related tasks and problems to leverage prior geometric and topological knowledge for improving learning for nonlinear systems.

\section{Acknowledgments}
{ \fontsize{9.0}{10} \selectfont
Authors research supported by grants DOE Grant ASCR PHILMS DE-SC0019246 and NSF Grant DMS-1616353. Also to R.N.L. support by a donor to UCSB CCS SURF program.  Authors also acknowledge UCSB Center for Scientific Computing NSF MRSEC (DMR1121053) and UCSB MRL NSF CNS-1725797.  P.J.A. would also like to acknowledge a hardware grant from Nvidia.
}

\newpage
\clearpage