# Event-Based Control for Synchronization of Stochastic Linear Systems with
Application to Distributed Estimation

Jiaqi Yan  Yilin Mo  Hideaki Ishii  Department of Computer Science, Tokyo
Institute of Technology, Japan. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS> Department of Automation and BNRist, Tsinghua
University, Beijing, China. E-mail<EMAIL_ADDRESS>
###### Abstract

This paper studies the synchronization of stochastic linear systems which are
subject to a general class of noises, in the sense that the noises are bounded
in covariance but might be correlated with the states of agents and among each
other. We propose an event-based control protocol for achieving the
synchronization among agents in the mean square sense and theoretically
analyze the performance of it by using a stochastic Lyapunov function, where
the stability of $c$-martingales is particularly developed to handle the
challenges brought by the general model of noises and the event-triggering
mechanism. The proposed event-based synchronization algorithm is then applied
to solve the problem of distributed estimation in sensor network.
Specifically, by losslessly decomposing the optimal Kalman filter, it is shown
that the problem of distributed estimation can be resolved by using the
algorithms designed for achieving the synchronization of stochastic linear
systems. As such, an event-based distributed estimation algorithm is
developed, where each sensor performs local filtering solely using its own
measurement, together with the proposed event-based synchronization algorithm
to fuse the local estimates of neighboring nodes. With the reduced
communication frequency, the designed estimator is proved to be stable under
the minimal requirements of network connectivity and collective system
observability.

##  1 Introduction

Because of the widespread applications, the past decades have witnessed ever-
growing interests in multi-agent systems (MASs). As one of its fundamental
focuses, the synchronization problem for linear MASs has been studied in the
literature. Given a group of agents modeled by general linear dynamics, the
distributed controllers are proposed therein, towards the end of achieving the
asymptotic consensus on the states of agents. For instance, in the continuous-
time domain, Ma and Zhang (2010) find that minimal requirement for
consensusability is that the agent dynamics has to be stabilizable and the
communication graph must contain a spanning tree for unstable agent dynamics.
On the other hand, You and Xie (2011) consider the MASs with discrete-time
linear dynamics. Under a certain relation between the system matrix and the
Laplacian matrix, the authors show that the synchronization among agents is
able to be achieved with the static local feedback control. Since then,
various control protocols are developed within this field of research, such as
Gu et al. (2011); Xu et al. (2019); Yan et al. (2021).

On the other hand, the autonomous agents are often equipped with embedded
microprocessors and onboard communication that are powered by energy-finite
batteries and thus have limited energy resources. As such, event-triggering
has recently been popular for their capabilities of improving the resource
utilization efficiency (Dimarogonas et al. (2011); Yi et al. (2018); Kadowaki
and Ishii (2014); Nowzari et al. (2019); Yan et al. (2020)). Since the early
works Garcia et al. (2014) and Zhu et al. (2014), remarkable efforts have also
been devoted to developing the event-based control protocols for solving the
synchronization problem in deterministic linear MASs. In these years,
stochastic linear systems with event-triggered schemes have received
particular research attention. For example, Ma et al. (2016) focus on the
stochastic MASs, where the dynamics of each agent is subject to mutually
uncorrelated zero-mean Gaussian white noises. Using linear matrix
inequalities, they provide an event-based controller which facilitates the
synchronization among agents in the mean square sense while decreasing their
communication frequency. Considering state-dependent noises, Ding et al.
(2015) analyze control performance leveraging the theory of input-to-state
stability in probability, and derive sufficient conditions under which the
consensus in probability is reached by using an event-triggered control
protocol. However, despite these works, the results in this area have been
scattered in literature.

In this work, we also focus on the event-based control for stochastic linear
system synchronization. Different from the existing works, this paper
considers a more general class of noises. It includes the independent Gaussian
white noise model in Ma et al. (2016) and the state-dependent noise model in
Ding et al. (2015), since the noises are only assumed to be bounded in
covariance while they might be correlated with the states of agents along time
and among agents. Because of its generality, this problem reveals various
applications in both theoretical and engineering fields (Ding et al. (2015);
Qi et al. (2020)).

Here, we follow the approach of Yang et al. (2021), which finds that, by
performing the local decomposition of the Kalman filter, the problem of
distributed state estimation can be reformulated to that of the
synchronization of stochastic linear MASs. We would therefore show how the
event-based control protocol proposed in this work can be applied to develop
stable distributed estimators as well. In particular, the aim is to estimate
the state of an $n$-dimensional linear time invariant system by using the
measurements from $m$ sensors, where the Kalman filter offers the optimal
solution in a centralized manner (Anderson and Moore (2012)). This paper
presents a novel framework for the event-based distributed implementation of
the centralized Kalman filter, where each sensor node performs local filtering
with its own measurement based on a local decomposition of Kalman filter, and
information fusion by communicating with neighbors through the proposed event-
based synchronization algorithm. The designed estimator is proved to be stable
(i.e., the state estimation error is mean-square bounded at each node side),
which extends the results in Yang et al. (2021) for the full transmission
case.

##  2 Problem Formulation

Let us consider the synchronization of a group of $m$ agents. The dynamics of
each agent $i$ is described by the following stochastic linear system:

$\eta_{i}(k+1)=S\eta_{i}(k)+Bu_{i}(k)+L_{i}z_{i}(k),$ (1)

where $\eta_{i}(k)\in\mathbb{R}^{n}$ and $u_{i}(k)\in\mathbb{R}$ are
respectively the state and control input, and $z_{i}(k)$ is the system noise
with zero mean and bounded covariance, i.e., $\mathbb{E}(z_{i}(k))=0$ and
$\operatorname*{cov}(z_{i}(k))$ is bounded at any time. By assuming so, we
consider the most general class of noises which might be correlated along time
and among agents.

In this paper, we would like to design the control input for each node such
that they are able to achieve synchronization in the presence of the noises
$\\{z_{i}(k)\\}$. To this end, we make the following definitions:

$\displaystyle\bar{\eta}(k)\triangleq\frac{1}{m}\sum_{i=1}^{m}\eta_{i}(k),\;\bar{L}_{z}(k)\triangleq\frac{1}{m}\sum_{i=1}^{m}L_{i}z_{i}(k).$
(2)

Namely, $\bar{\eta}(k)$ and $\bar{z}(k)$ respectively denote the averages of
local states and noises at current time. Then the network of agents is said to
reach synchronization in the mean square sense, if it holds at any $k\geq 0$
that:

  1. 1.

Consistency: The average of local states keeps consistent throughout the
execution, i.e.,

$\bar{\eta}(k+1)=S\bar{\eta}(k)+\bar{L}_{z}(k).$ (3)

  2. 2.

Consensus: The agents can achieve consensus with bounded error covariance,
i.e., there exists a bounded $\Xi$ such that

$\operatorname*{cov}[(\eta_{i}(k)-\bar{\eta}(k))]\leq\Xi.$ (4)

We elucidate these statements as below. First, the consistency condition
claims that the dynamics of average state $\bar{\eta}(k)$ is governed by a
linear system, the input to which is the external noises only. Therefore, the
interaction among agents will only affect the evolution of local states but
not $\bar{\eta}(k)$. Second, due to the presence of uncertain noises, it is
unable for the agents to achieve exact agreement, i.e., to perfectly track the
average state $\bar{\eta}(k)$. As a compromise, in (4), we aim to track
$\bar{\eta}(k)$ with bounded error covariance. We should note that in the
absence of noises, these conditions become identical to the standard ones as
in You and Xie (2011); Xu et al. (2019); Gu et al. (2011).

For a stable system matrix $S$, it is straightforward that simply holding the
zero input, $u_{i}(k)=0$, can realize the synchronization among agents. To
avoid this trivial case, this paper focuses on the unstable system matrix:

###### Assumption 1

At least one eigenvalue of $S$ lies on or outside the unit circle. Moreover,
$(S,B)$ is controllable.

The communication network is modeled by a connected undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{A})$, where
$\mathcal{V}=\\{1,...,m\\}$ is the set of agents,
$\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ is the set of edges.
Moreover, $\mathcal{A}=\left[a_{ij}\right]$ is the weighted adjacency matrix,
where $a_{ij}\geq 0$ and $a_{ij}=a_{ji},\forall i,j\in\mathcal{V}$. Notice
that $(i,j)\in\mathcal{E}$ if and only if $a_{ij}>0$. Let us arrange the
eigenvalues of Laplacian matrix $\mathcal{L}_{\mathcal{G}}$ as

$0=\mu_{1}<\mu_{2}\leq\cdots\leq\mu_{m}.$ (5)

###  2.1 Event-triggered communication and control strategy

This paper aims to design an event-driven strategy, under which
synchronization among agents is reached with low communication frequency. In
particular, for any agent $i$, it broadcasts to neighbors only when a certain
event at its local side is triggered. Based on the received knowledge, agent
$i$ updates its state according to

$u_{i}(k)=\Gamma\sum_{j=1}^{m}a_{ij}(\hat{\eta}_{j}(k)-\hat{\eta}_{i}(k)),$
(6)

where $\Gamma\in\mathbb{R}^{1\times n}$ is the matrix to be determined, and

$\hat{\eta}_{i}(k)\triangleq
S^{(k-k_{s}^{i})}\eta_{i}(k_{s}^{i}),\;\;k\in[k_{s}^{i},k_{s+1}^{i}),$ (7)

such that $\eta_{i}(k_{s}^{i})$ denotes the information most recently
broadcast by agent $i$.

To determine the triggering instants $k_{s}^{i}$ for agent $i$, let us define
the triggering function $f_{i}(k)$ as

$\begin{split}f_{i}(k)&=||\epsilon_{i}(k)||^{2}-(c_{0}+c_{1}\varrho^{k}),\\\
\end{split}$ (8)

where $c_{0},c_{1}$ are positive constants, $\varrho\in(0,1)$, and

$\epsilon_{i}(k)\triangleq\hat{\eta}_{i}(k)-\eta_{i}(k).$ (9)

At the initial phase, communications at all agents are carried out. After
that, each agent updates its local state based on (6) until the triggering
function (8) exceeds $0$. To be specific, once $f_{i}(k)\geq 0$, agent $i$
will be triggered. It then broadcasts its local state $\eta_{i}(k)$ to
neighbors, yielding that $\epsilon_{i}(k)$ is reset to zero. Therefore, the
sequence of triggering instants is determined recursively as
$k_{s+1}^{i}\triangleq\min\left\\{k>k_{s}^{i}\mid f_{i}(k)\geq
0\right\\},\;k_{0}^{i}=0.$

##  3 Performance Analysis

We shall next study the performance of the proposed algorithm. Particularly,
in order to reach the synchronization among agents, it is desired to establish
both consistency and consensus conditions.

Notice that the presence of noises prevents us from directly applying
approaches of Lyapunov stability for deterministic systems to the analysis. We
therefore resort to a stochastic analogue of it. To this end, we will prove an
important preliminary result for our purpose, i.e., $c$-martingale convergence
lemma, which would be used later in the stability analysis:

######  Lemma 1 ($c$-martingale convergence lemma)

Let $\\{\mathcal{F}(t)\\}$ be a filtration in a probability space
$(\Omega,\mathcal{F},\mathcal{P})$. Moreover, let $\\{V(k)\\}$ be a sequence
of non-negative random variables such that

$\mathbb{E}[\Delta V(k)]\leq-\rho V(k)+c(k),$ (10)

where $\mathbb{E}[\Delta
V(k)]\triangleq\mathbb{E}[V(k+1)|\mathcal{F}(k)]-V(k),$ $\rho>0$, and
$\mathbb{E}[c(k)]\leq c<\infty$. Then it holds for $k\geq 0$ that
$\mathbb{E}\left[V(k)\right]$ is bounded.

It follows from (10) that

$0\leq\mathbb{E}[V(k+1)|\mathcal{F}(k)]\leq(1-\rho)V(k)+c(k).$ (11)

The result is thus obvious by taking expectation on both sides of (11), which
yields that

$\displaystyle 0$
$\displaystyle\leq\mathbb{E}\left[V(k+1)\right]\leq\left(1-\rho\right)\mathbb{E}\left[V(k)\right]+c$
(12)
$\displaystyle\leq\left(1-\rho\right)^{k+1}\mathbb{E}\left[V(0)\right]+c\sum_{t=0}^{k+1}\left(1-\rho\right)^{t}.$

The proof is thus finished. $\square$

We are now ready to present the main result on achieving synchronization of
stochastic linear systems:

###### Theorem 1

Suppose that the Mahler measure111The Mahler measure of a matrix is defined as
the absolute product of unstable eigenvalues of it. of $S$ meets the condition

$\prod_{j}|\lambda_{j}^{u}(S)|<\frac{1+\mu_{2}/\mu_{m}}{1-\mu_{2}/\mu_{m}},$
(13)

where $\lambda_{j}^{u}(S)$ denotes the $j$th unstable eigenvalue of $S$, and
$\mu_{2}$ and $\mu_{m}$ are defined in (5). Let $\Gamma$ be designed according
to

$\Gamma=\frac{2}{\mu_{2}+\mu_{m}}\frac{B^{T}\mathcal{P}S}{B^{T}\mathcal{P}B},$
(14)

where $\mathcal{P}>0$ is the solution to the following inequality

$\mathcal{P}-S^{T}\mathcal{P}S+\left(1-\zeta^{2}\right)\frac{S^{T}\mathcal{P}BB^{T}\mathcal{P}S}{B^{T}\mathcal{P}B}>0,$
(15)

with $\zeta$ satisfying
$\prod_{j}\left|\lambda_{j}^{u}(S)\right|<\zeta^{-1}\leq\frac{1+\mu_{2}/\mu_{m}}{1-\mu_{2}/\mu_{m}}.$
With the event-based controller (6), the network of agents reach
synchronization in the mean square sense..

##  4 Application to Event-based Distributed State Estimation

There are many practices where the synchronization of stochastic linear
systems is applicable. In this section, we will provide a non-trivial example
showing how the results developed in the previous sections can be applied to
solve the problem of distributed state estimation.

To this end, let us first introduce the formulation of distributed estimation
and consider the following discrete-time dynamical system:

$x(k+1)=Ax(k)+w(k),$ (16)

where $x(k)\in\mathbb{R}^{n}$ is the system state to be estimated,
$w(k)\sim\mathcal{N}(0,Q)$ is the independent and identically distributed
(i.i.d.) Gaussian disturbance with zero mean and covariance matrix $Q\geq 0$.
The initial state $x(0)$ also follows the Gaussian distribution which has zero
mean.

A sensor network monitors the system above. The measurement from each sensor
$i\in\\{1,2,...,m\\}$ is given by

$y_{i}(k)=C_{i}x(k)+v_{i}(k),$ (17)

where $y_{i}(k)\in\mathbb{R}$ is the measurement produced by sensor $i$,
$C_{i}\in\mathbb{R}^{1\times n}$, and $v_{i}(k)\in\mathbb{R}$ is the
measurement noise.

By collecting the measurements from all sensors, we have

$y(k)=Cx(k)+v(k),$ (18)

where

$\begin{split}y(k)\triangleq{\left[\begin{array}[]{c}y_{1}(k)\\\ \vdots\\\
y_{m}(k)\end{array}\right],}\;C\triangleq{\left[\begin{array}[]{c}C_{1}\\\
\vdots\\\
C_{m}\end{array}\right],}\;v(k)\triangleq{\left[\begin{array}[]{c}v_{1}(k)\\\
\vdots\\\ v_{m}(k)\end{array}\right]},\end{split}$ (19)

and $v(k)$ is zero-mean i.i.d. Gaussian noise with covariance $R\geq 0$ and is
independent of $w(k)$ and $x(0)$.

In this section, we make the following assumption on system observability,
which is the minimum requirement for solving the distributed estimation
problem:

######  Assumption 2 (Jointly observable)

The system is jointly observable, i.e., the pair $(A,C)$ is observable, while
$(A,C_{i})$ is not necessarily observable for each $i\in\\{1,\cdots,m\\}$.

###  4.1 Fundamental limit: Kalman filter

If there exists a single fusion center that has access to measurements from
all sensors, then the centralized Kalman filter offers the optimal solution
and provides a fundamental limit for all distributed estimation schemes. We
thus introduce Kalman filter before going to the distributed solution.

Let $P(k)$ be the error covariance of Kalman estimate at time $k$. Under
Assumption 2, the error covariance will converge to the steady state
exponentially fast (Anderson and Moore (2012)):

$\displaystyle P=\lim_{k\rightarrow\infty}P(k).$ (20)

Since a sensor network typically operates for a long period of time, we
consider the steady-state Kalman filter, which has the fixed gain

$\displaystyle K=PC^{T}\left(CPC^{T}+R\right)^{-1}.$ (21)

By using $K$, the optimal Kalman estimate is calculated recursively as

$\begin{split}\hat{x}(k+1)=(A-KCA)\hat{x}(k)+Ky(k+1).\end{split}$ (22)

###  4.2 An event-based distributed implementation of Kalman filter

It is clear that the optimal estimate (22) requires the information from all
sensors. However, in a distributed framework, each sensor is only capable of
communicating with immediate neighbors, which renders the centralized solution
impractical. Therefore, research attention has been devoted to approaching the
performance of Kalman estimate in a distributed manner.

In this work, we propose a new framework for achieving the distributed
estimation with an event-triggered communication strategy. Compared with
existing algorithms, the novelty of this framework is that we decouple the
procedure of local filters from that of the fusion process so that
communication among sensors occurs only in the fusion process and will not
affect the performance of local filters. Particularly, the update at each time
can be divided into two phases. In the first phase, a lossless local
decomposition of Kalman filter is proposed, based on which each sensor runs a
local filter solely using its own measurement. After that, the second phase
fuses the local knowledge by performing the event-based synchronization
algorithm from Section 2, to obtain the stable local estimators while
decreasing the transmission frequency.

In what follows, we will detail our solution for solving the distributed
estimation problem by respectively introducing these two phases.

####  4.2.1 Phase I: Implementation of local filters

In the first phase, a local decomposition of the Kalman filter is proposed,
where each sensor runs a local filter based on its own measurement without
communicating with others. It will be proved that the Kalman estimate (22) is
indeed a linear combination of the local filters. Therefore, by properly
fusing the local filters in the second phase, it is possible to approach the
performance of the Kalman filter.

Before going on, we make the following assumption:

###### Assumption 3

  1. 1.

all the eigenvalues of $A-KCA$ are distinct;

  2. 2.

$A-KCA$ and $A$ do not have any common eigenvalues.

###### Remark 1

In case that $A-KCA$ does not satisfy Assumption 3, we can always perturb $K$
to ensure the validity of these assumptions with little loss of the estimation
performance. This is feasible as the system is observable and thus the
eigenvalues of $A-KCA$ can be freely assigned.

Under Assumption 3, we can find a non singular matrix $V$, by which $A-KCA$ is
diagonalized as

$A-KCA=V\Lambda V^{-1},$ (23)

where $\Lambda\triangleq\operatorname*{diag}(\lambda_{1},\cdots,\lambda_{n})$
such that $\lambda_{1},\cdots,\lambda_{n}$ are the eigenvalues of $A-KCA$.

In this paper, we consider the following local filter, which is performed by
each sensor $i$ solely based on its local measurements:

$\begin{split}z_{i}(k)&=y_{i}(k+1)-\beta^{T}\hat{\xi}_{i}(k),\\\
\hat{\xi}_{i}(k+1)&=S\hat{\xi}_{i}(k)+\operatorname{\textbf{1}}_{n}z_{i}(k),\;\forall
i\in\\{1,\cdots,m\\},\end{split}$ (24)

where $\hat{\xi}_{i}(k)$ and $z_{i}(k)$ are respectively the state and output
of the local filter, $\beta\in\mathbb{R}^{n}$ and $G_{i}$ solve the following
equation:

$\displaystyle(G_{i}-\operatorname{\textbf{1}}_{n}C_{i})A=\Lambda
G_{i},\;\beta^{T}G_{i}=C_{i}A,$ (25)

and

$S\triangleq\Lambda+\operatorname{\textbf{1}}_{n}\beta^{T}.$ (26)

Before analyzing the performance of the local filter (24), some of its
properties are summarized below:

###### Lemma 2

Consider the local filter (24). There always exist $\beta$ and $G_{i}$ that
solve (25). Moreover, it holds for $k\geq 0$ that the covariance of $z_{i}(k)$
is bounded.

We shall next study the performance of local filter (24), namely, how it is
related to that of the Kalman filter. For simplicity, decompose the Kalman
gain as $K=[K_{1},\cdots,K_{m}]$, where $K_{j}$ denote the $j$th column of
$K$. In the next lemma, it is shown that the Kalman estimate (22) is indeed a
weighted sum of the local filters and thus can be perfectly recovered by the
collective information of them:

######  Lemma 3 ( Mo and Garone (2016))

Suppose Assumption 3 holds. The Kalman estimate (22) can be recovered by
$\hat{\xi}_{i}(k),i=1,...,m$, through a linear combination as below:

$\hat{x}(k)=\sum_{i=1}^{m}F_{i}\hat{\xi}_{i}(k),$ (27)

where222Notice $\operatorname*{diag}(V^{-1}K_{i})$ is a diagonal matrix, where
the $j$th diagonal entry of $\operatorname*{diag}(V^{-1}K_{i})$ equals the
$j$th entry of vector $V^{-1}K_{i}$.
$F_{i}=V\operatorname*{diag}(V^{-1}K_{i})$.

####  4.2.2 Phase II: Fusion of local filters via event-based synchronization
algorithm

By virtue of (27), it is clear that the optimal estimate needs to fuse
$\hat{\xi}_{j}(k)$ of all sensors. However, in the distributed formulation, a
single sensor can only access the local states within its neighborhood.
Therefore, in the second phase, each sensor $i$ runs a synchronization
algorithm, through which a stable local estimate is obtained by inferring
$\hat{\xi}_{j}(k)$ for all $j\in\mathcal{V}$. For the particular purpose of
decreasing the transmission frequency, we will adopt the event-based
synchronization algorithm proposed in Section 2.

To simplify notations, let $\eta_{i,j}(k)$ denote the inference from sensor
$i$ on $\hat{\xi}_{j}(k)$. By collecting its inference on all sensors, each
sensor $i$ contains a local state

$\eta_{i}(k)\triangleq{\left[\begin{array}[]{c}\eta_{i,1}(k)\\\ \vdots\\\
\eta_{i,m}(k)\end{array}\right]\in\mathbb{R}^{mn}.}$ (28)

Particularly, the local inference will be updated through the following
synchronization algorithm:

$\eta_{i}(k+1)=\tilde{S}\eta_{i}(k)+\tilde{L}_{i}z_{i}(k)+\tilde{B}\tilde{\Gamma}\sum_{j=1}^{m}a_{ij}(\hat{\eta}_{j}(k)-\hat{\eta}_{i}(k)),$
(29)

where $\tilde{S}\triangleq I_{m}\otimes S,\;\tilde{L}_{i}\triangleq
e_{i}\otimes\operatorname{\textbf{1}}_{n}$ with $e_{i}$ denoting the $i$th
canonical basis vector in $\mathbb{R}^{m}$, $\tilde{B}\triangleq I_{m}\otimes
B$ such that $(S,B)$ is controllable, $\tilde{\Gamma}=I_{m}\otimes\Gamma$, and
$\hat{\eta}_{i}(k)$ is the latest information broadcast by sensor $i$ and is
calculated by

$\hat{\eta}_{i}(k)=\tilde{S}^{(k-k_{s}^{i})}\eta_{i}(k_{s}^{i}).$ (30)

By collecting Phases I and II together, we summarize the update of any sensor
$i$ in Algorithm 1, where the initial values are set as $\hat{x}_{i}(0)=0$ and
$\eta_{i}(0)=0$. The information flow of Algorithm 1 is provided in Fig. 1,
which requires no fusion center and thus is achieved in a distributed manner.

1: Sorely using its own measurement, sensor $i$ computes $z_{i}(k)$ and
updates the state of the local filter by (24).
2: By fusing the most recently received information within its neighborhood,
sensor $i$ updates $\eta_{i}(k+1)$ according to the synchronization algorithm
(29)–(30) and (8)–(9).
3: Sensor $i$ updates the local estimate as

$\breve{x}_{i}(k+1)=mF\eta_{i}(k+1),$ (31)

where $F\triangleq\left[F_{1},F_{2},\cdots,F_{m}\right]$.
4: Sensor $i$ broadcasts the new state $\eta_{i}(k+1)$ to neighbors if and
only if the triggering function (8) exceeds $0$.

Algorithm 1 An event-based distributed estimation algorithm for sensor $i$ at
time $k>0$

$y_{i}(k)$$y_{j}(k)$Linear system$\eta_{i}(k)$Linear system$\eta_{j}(k)$Local
filter$\hat{\xi}_{i}(k)$Local
filter$\hat{\xi}_{j}(k)$$z_{i}(k)$$z_{j}(k)$$\hat{\eta}_{i}(k)$$\hat{\eta}_{j}(k)$Phase
IPhase II$\breve{x}_{i}(k)$$\breve{x}_{j}(k)$

Figure 1: The information flow of Algorithm 1.

###  4.3 Stability analysis of local estimators

As shown in Lemma 2, the covariance of $z_{i}(k)$ in (29) is bounded at any
time. As proved in Section 3, the synchronization algorithm (29) guarantees
that the local states $\eta_{i}(k)$’s achieve both the consistency and
consensus conditions. We shall, in the sequel of this part, show how these
conditions will help to achieve a stable local estimate at each sensor side.

###### Theorem 2

Suppose that the Mahler measure of $S$ meets condition (13), and $\Gamma$ is
designed based on (14) and (15). Algorithm 1 yields a stable estimator at each
sensor side. Specifically, it holds for any $k\geq 0$ that:

  1. 1.

The average of local estimates from all sensor equals the Kalman estimate.
That is,

$\frac{1}{m}\sum_{i=1}^{m}\breve{x}_{i}(k)=\hat{x}(k),\forall k\geq 0.$ (32)

  2. 2.

The error covariance of local estimate from each sensor $i$, i.e.,
$\operatorname*{cov}(\breve{x}_{i}(k)-x(k))$, is bounded.

The proof respectively hold by the consistency condition (3) and consensus
condition (4).

(1) Due to consistency condition (3), we obtain

$\sum_{i=1}^{m}\eta_{i}(k+1)=\tilde{S}\sum_{i=1}^{m}\eta_{i}(k)+\sum_{i=1}^{m}\tilde{L}_{i}z_{i}(k).$
(33)

Comparing it with (24), we conclude the following relation for any time $k$
and any $j\in\mathcal{V}$:

$\hat{\xi}_{j}(k)=\sum_{i=1}^{m}\eta_{i,j}(k).$ (34)

Thereby, the following equation is satisfied at any $k\geq 0$:

$\begin{split}\frac{1}{m}\sum_{i=1}^{m}\breve{x}_{i}(k)&=\sum_{i=1}^{m}F\eta_{i}(k)=\sum_{i=1}^{m}\sum_{j=1}^{m}F_{j}\eta_{i,j}(k)\\\
&=\sum_{j=1}^{m}F_{j}\Big{[}\sum_{i=1}^{m}\eta_{i,j}(k)\Big{]}=\sum_{j=1}^{m}F_{j}\hat{\xi}_{j}(k)=\hat{x}(k).\end{split}$
(35)

(2) Let us consider the local estimator of any sensor $i$. To begin with, by
the virtue of (4), we conclude that $\operatorname*{cov}(\delta_{i}(k))$ is
bounded at any time $k$, where

$\displaystyle\delta_{i}(k)=\eta_{i}(k)-\bar{\eta}(k).$

In order to prove the boundedness of error covariance at each sensor side, let
us denote

$\bar{e}_{i}(k)\triangleq\breve{x}_{i}(k)-\hat{x}(k),$ (36)

which is the bias from local estimate $\breve{x}_{i}(k)$ to the optimal Kalman
filter. Combining (27) and (34) yields

$\hat{x}(k)=F\sum_{i=1}^{m}\eta_{i}(k)=mF\bar{\eta}(k).$ (37)

One thus has

$\bar{e}_{i}(k)=mF(\eta_{i}(k)-\bar{\eta}(k))=mF\delta_{i}(k).$ (38)

Therefore, the estimation error of sensor $i$ is calculated as

$\begin{split}\breve{e}_{i}(k)&=\breve{x}_{i}(k)-x(k)=(\breve{x}_{i}(k)-\hat{x}(k))+(\hat{x}(k)-x(k))\\\
&=\bar{e}_{i}(k)+\hat{e}(k),\end{split}$

where $\hat{e}(k)$ is the estimation error of Kalman filter. Since Kalman
filter is optimal, $\bar{e}_{i}(k)$ is orthogonal to $\hat{e}(k)$. Therefore,
it follows that

$\displaystyle\operatorname*{cov}(\breve{e}_{i}(k))$
$\displaystyle=\operatorname*{cov}(\bar{e}_{i}(k))+\operatorname*{cov}(\hat{e}(k))$
(39) $\displaystyle=m^{2}F\operatorname*{cov}(\delta_{i}(k))F^{T}+P\leq
m^{2}F\Xi F^{T}+P,$

where the inequality holds by the consensus condition (4), and $P$ is steady-
state error convariance of Kalman filter as defined in (20). We therefore
complete the proof. $\square$

Therefore, by applying (29), the problem of distributed estimation is solved.

##  5 Numerical Example

In the numerical example, we consider the case where a network of four sensors
monitors the system with following parameters:

$\begin{split}&A=\begin{bmatrix}0.9&0\\\
0&1.1\end{bmatrix},\;C=\begin{bmatrix}1&0&1&1\\\
0&1&1&-1\end{bmatrix}^{T},\;Q=0.5I_{2},\;R=2I_{4}.\end{split}$ (40)

Suppose that the four sensors are connected as a ring, where each edge is
assigned with weight $1$. The Laplacian matrix is thus:

$\mathcal{L_{G}}=\begin{bmatrix}2&-1&0&-1\\\ -1&2&-1&0\\\ 0&-1&2&-1\\\
-1&0&-1&2\end{bmatrix}.$ (41)

Therefore, one has $\mu_{2}=2$, $\mu_{4}=4$. To guarantee the condition in
Theorem 2, let us choose $\zeta=0.5$. Based on these parameters, we calculate
$\Gamma=[0.80,-0.41].$

The initial states are set as $x(0)\sim\mathcal{N}(0,I)$ and
$\breve{x}_{i}(0)=0$ for each sensor $i$. By performing Algorithm 1, we
observe from Fig. 2 that the mean squared estimation error from each sensor is
stable during the operation. Moreover, the average communication rate over the
whole network is $66.5\%$ obtained by the $1000$-run Monte Carlo trials. On
the other hand, the performance loss is only $5.3\%$ compared to the full
transmission case.

$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$Estimate
error of
$x_{1}$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$TimeEstimate
error of $x_{2}$s1s2s3s4KF

Figure 2: Average mean square estimation error of system states in $1000$-run
Monte Carlo trials.

##  6 Conclusion

In this paper, the event-triggered synchronization problem has been addressed
for the discrete-time stochastic MASs which are subject to a general class of
noises. By using the stochastic Lyapunov stability theory, we propose an
event-based control protocol, which is proved to facilitate the mean-squared
synchronization among agents. In particular, through the local decomposition
of the Kalman filter, we transform the problem of distributed estimation to
the synchronization of stochastic linear systems, in which the proposed event-
based protocol is applied to yield stable local estimators while decreasing
the communication frequency.

## References

  * Anderson and Moore (2012) Anderson, B.D. and Moore, J.B. (2012).  _Optimal Filtering_.  Courier Corporation.
  * Dimarogonas et al. (2011) Dimarogonas, D.V., Frazzoli, E., and Johansson, K.H. (2011).  Distributed event-triggered control for multi-agent systems.  _IEEE Transactions on Automatic Control_ , 57(5), 1291–1297.
  * Ding et al. (2015) Ding, D., Wang, Z., Shen, B., and Wei, G. (2015).  Event-triggered consensus control for discrete-time stochastic multi-agent systems: The input-to-state stability in probability.  _Automatica_ , 62, 284–291.
  * Garcia et al. (2014) Garcia, E., Cao, Y., and Casbeer, D.W. (2014).  Decentralized event-triggered consensus with general linear dynamics.  _Automatica_ , 50(10), 2633–2640.
  * Gu et al. (2011) Gu, G., Marinovici, L., and Lewis, F.L. (2011).  Consensusability of discrete-time dynamic multiagent systems.  _IEEE Transactions on Automatic Control_ , 57(8), 2085–2089.
  * Kadowaki and Ishii (2014) Kadowaki, Y. and Ishii, H. (2014).  Event-based distributed clock synchronization for wireless sensor networks.  _IEEE Transactions on Automatic Control_ , 60(8), 2266–2271.
  * Ma and Zhang (2010) Ma, C.Q. and Zhang, J.F. (2010).  Necessary and sufficient conditions for consensusability of linear multi-agent systems.  _IEEE Transactions on Automatic Control_ , 55(5), 1263–1268.
  * Ma et al. (2016) Ma, L., Wang, Z., and Lam, H.K. (2016).  Event-triggered mean-square consensus control for time-varying stochastic multi-agent system with sensor saturations.  _IEEE Transactions on Automatic Control_ , 62(7), 3524–3531.
  * Mo and Garone (2016) Mo, Y. and Garone, E. (2016).  Secure dynamic state estimation via local estimators.  In _Proc. IEEE 55th Conference on Decision and Control (CDC)_ , 5073–5078.
  * Nowzari et al. (2019) Nowzari, C., Garcia, E., and Cortés, J. (2019).  Event-triggered communication and control of networked systems for multi-agent consensus.  _Automatica_ , 105, 1–27.
  * Qi et al. (2020) Qi, Y., Deng, H., Wang, J., and Tang, Y. (2020).  Passivity-based synchronization stability analysis for power-electronic-interfaced distributed generations.  _IEEE Transactions on Sustainable Energy_ , 12(2), 1141–1150.
  * Xu et al. (2019) Xu, L., Mo, Y., and Xie, L. (2019).  Distributed consensus over markovian packet loss channels.  _IEEE Transactions on Automatic Control_ , 65(1), 279–286.
  * Yan et al. (2021) Yan, J., Deng, C., and Wen, C. (2021).  Resilient output regulation in heterogeneous networked systems under Byzantine agents.  _Automatica_ , 133, 109872.
  * Yan et al. (2020) Yan, J., Li, X., Mo, Y., and Wen, C. (2020).  Resilient multi-dimensional consensus and optimization in adversarial environment.  _arXiv preprint arXiv:2001.00937_.
  * Yang et al. (2021) Yang, X., Yan, J., Mo, Y., and You, K. (2021).  A distributed implementation of steady-state kalman filter.  In _Proc. 40th Chinese Control Conference (CCC)_ , 5154–5159.
  * Yi et al. (2018) Yi, X., Liu, K., Dimarogonas, D.V., and Johansson, K.H. (2018).  Dynamic event-triggered and self-triggered control for multi-agent systems.  _IEEE Transactions on Automatic Control_ , 64(8), 3300–3307.
  * You and Xie (2011) You, K. and Xie, L. (2011).  Network topology and communication data rate for consensusability of discrete-time multi-agent systems.  _IEEE Transactions on Automatic Control_ , 56(10), 2262–2275.
  * Zhu et al. (2014) Zhu, W., Jiang, Z.P., and Feng, G. (2014).  Event-based consensus of multi-agent systems with general linear models.  _Automatica_ , 50(2), 552–558.


# On the suboptimality of the $p$-version discontinuous Galerkin methods for
first order hyperbolic problems

Zhaonan Dong and Lorenzo Mascotto  Inria, 2 rue Simone Iff, 75589 Paris,
France and CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée 2, France
<EMAIL_ADDRESS>Fakultät für Mathematik, Universität Wien, 1090 Vienna,
Austria<EMAIL_ADDRESS>
###### Abstract.

We address the issue of the suboptimality in the $p$-version discontinuous
Galerkin (dG) methods for first order hyperbolic problems. The convergence
rate is derived for the upwind dG scheme on tensor product meshes in any
dimension. The standard proof in seminal work [14] leads to suboptimal
convergence in terms of the polynomial degree by $3/2$ order for general
convection fields, with the exception of piecewise multi-linear convection
fields, which rather yield optimal convergence. Such suboptimality is not
observed numerically. Thus, it might be caused by a limitation of the
analysis, which we partially overcome: for a special class of convection
fields, we shall show that the dG method has a $p$-convergence rate suboptimal
by $1/2$ order only.

SUBMITTED AS A PROCEEDING ARTICLE

AMS subject classification: 65N12, 65N15, 65N30

Keywords: $hp$-finite element methods, discontinuous Galerkin methods,
hyperbolic problems

##  1\. Introduction

Discontinuous Galerkin (dG) finite element methods were introduced in the
early 1970s for the numerical solution of first-order hyperbolic problems [17]
and the weak imposition of inhomogeneous boundary conditions for elliptic
problems [16]. In the past several decades, dG methods have enjoyed
considerable success as a standard variational framework for the numerical
solution of many classes of problems involving partial differential equations
(PDEs); see, e.g., monographs [8, 9, 10] for reviews of some of the main
developments of dG methods. The interest in dG methods can be attributed to a
number of factors, including the great flexibility in dealing with
$hp$-adaptivity and general shaped elements [7, 5, 6], as well as in solving
convection-dominated PDEs; see, e.g., early works [3, 2] concerning hyperbolic
conservation laws and convection-diffusion problems.

Due to missing tools in the analysis, the convergence rate always contains
suboptimality in terms of the polynomial degree $p$. In [13], the first
optimal convergence rate of $hp$-dG methods is derived for linear convection
problems by using the SUPG stabilisation. However, the authors provide
numerical evidence that the $hp$-optimal convergence rate is achieved even
without such stabilisation. In seminal work [14], based on (back then) novel
optimal approximation results for the $L^{2}$-orthogonal projection, the
$hp$-optimal convergence rate is derived for dG methods applied to hyperbolic
problems, under the technical assumption that the convection field is
piecewise linear. Moreover, whenever the above assumption is violated, the
theoretical analysis in [14] leads to error bounds that are suboptimal in
terms of $p$ by $3/2$ order. Such suboptimality is yet not observed in the
numerical experiments. Over the last two decades, the above mentioned
technical assumption became standard in $hp$-dG methods for convection-
diffusion-reaction and hyperbolic problems; see, e.g., [12, 5, 6, 4]. It is
still an open question, whether the $p$-suboptimality for dG methods by $3/2$
order is true or not in general.

Our contribution represents a further step in shedding light on this issue.
Notably, we present the a priori error analysis for $hp$-dG methods applied to
pure hyperbolic problems employing a class of convection field, including
nonpolynomial cases. The new error is $h$-optimal and $p$-suboptimal by $1/2$
order only.

The rest of the paper is organized as follows: the continuous problem and its
dG discretization are addressed in Section 2; the classical analysis from [13]
is re-elaborated in Section 3, whereas the improved bounds under suitable
assumptions on the convection field are the topic of Section 4; we collect the
conclusions in Section 5.

Throughout, we employ a standard notation for Sobolev spaces [1].

##  2\. The continuous problem and its dG formulation

###  2.1. The continuous problem

Let $\Omega$ be a bounded polyhedral domain in $\mathbb{R}^{d}$,
$d\in\mathbb{N}$, with boundary $\Gamma$. We denote the unit outward normal
vector to $\Gamma$ at $x\in\Gamma$ by $\mathbf{n}_{\Gamma}(x)$ and introduce
the Fichera function $\mathbf{b}\cdot\mathbf{n}$ on $\Gamma$ to define

$\Gamma_{-}:=\\{x\in\Gamma:\mathbf{b}(x)\cdot\mathbf{n}(x)<0\\}\quad\text{and}\quad\Gamma_{+}:=\\{x\in\Gamma:\mathbf{b}(x)\cdot\mathbf{n}(x)\geq
0\\}.$ (1)

In the following, the sets $\Gamma_{-}$ and $\Gamma_{+}$ are referred to as
the _inflow_ and _outflow_ boundary, respectively, and clearly form a
nonoverlapping partition of $\Gamma$.

Given $\mathbf{b}\in[W^{1,\infty}(\Omega)]^{d}$, $c\in L^{\infty}(\Omega)$,
$f\in L^{2}(\Omega)$, and $g_{D}\in L^{2}(\Gamma_{-})$, we consider the
convection-reaction problem

$\displaystyle\mathbf{b}\cdot\nabla u+c\,u$
$\displaystyle=f\quad\quad\\!\\!\text{in }\Omega,$ (2) $\displaystyle u$
$\displaystyle=g_{D}\quad\text{on }\Gamma_{-}.$

Assuming the existence of a positive constant $c_{s}$ satisfying

$c_{0}^{2}:=c-\frac{1}{2}\nabla\cdot\mathbf{b}\geq c_{s}\qquad\text{for a.e.
}x\in\Omega,$ (3)

the well-posedness of problem (2) follows, e.g., as in [15].

###  2.2. The dG formulation on quadrilateral/hexahedral meshes

We are interested in discretizing solutions to (2) by means of a dG finite
element method. To the aim, consider sequences of meshes
$\\{\mathcal{T}_{n}\\}_{n}$ consisting of tensor product elements, which can
be defined through an affine mapping $\Phi_{K}$ on the reference
$d$-dimensional cube element $\widehat{K}:=(-1,1)^{d}$. For sake of
simplicity, we assume the all elements $K\in\mathcal{T}_{n}$ are shape
regular. We fix a uniform polynomial degree $p\in\mathbb{N}$ and denote the
space of tensor polynomials of degree $p$ over $\widehat{K}$ by
$\mathbb{Q}_{p}(\widehat{K})$. Next, we introduce the dG space

$V_{n}:=\\{v_{n}\in L^{2}(\Omega)\mid
v_{n}{}_{|K}\circ\Phi_{K}\in\mathbb{Q}_{p}(K)\;\forall
K\in\mathcal{T}_{n}\\}.$

Given an element $K\in\mathcal{T}_{n}$, we set its diameter and outward
pointing normal by $h_{K}$ and $\mathbf{n}_{K}$, respectively. Given boundary
$\Gamma^{K}$ of element $K$, we split it into the _inflow_ and _outflow_ parts
$\Gamma_{-}^{K}$ and $\Gamma_{+}^{K}$ defined as

$\Gamma_{-}^{K}:=\\{x\in\Gamma^{K}:\mathbf{b}(x)\cdot\mathbf{n}_{K}(x)<0\\}\quad\text{and}\quad\Gamma_{+}^{K}:=\\{x\in\Gamma^{K}:\mathbf{b}(x)\cdot\mathbf{n}_{K}(x)\geq
0\\}.$ (4)

Next, we define the classical _upwind jump_ operator. Given an internal face
$F$, let $K_{1}$ and $K_{2}$ be two elements in $\mathcal{T}_{n}$ sharing $F$.
Without loss of generality, we assume that $K_{1}$ is such that
$\mathbf{b}\cdot\mathbf{n}_{K}(x)<0$ for almost all $x$ in $F$. Then, we set

$\lfloor
v\rfloor_{F}:=(v{}_{|K_{1}}-v{}_{|K_{2}})_{|F}=v{}^{+}-v{}^{-}\quad\quad\forall
v\in H^{1}(\Omega,\mathcal{T}_{n}).$ (5)

In the rest of this work, when no confusion occurs, we shall write
$\mathbf{n}$ instead of $\mathbf{n}_{K}(x)$.

Following seminal work [14], we consider the upwind dG variational formulation
of (2). More precisely, introduce the dG bilinear form

$\begin{split}B_{n}(u_{n},v_{n})&:=\sum_{K\in\mathcal{T}_{n}}\Big{(}(\mathbf{b}\cdot\nabla
u_{n},v_{n})_{0,K}+(c\ u_{n},v_{n})_{0,K}\\\
&\quad\quad\quad\quad-((\mathbf{b}\cdot\mathbf{n})\lfloor
u_{n}\rfloor,v_{n}^{+})_{\Gamma_{-}^{K}\setminus\Gamma_{-}}-((\mathbf{b}\cdot\mathbf{n})u_{n}^{+},v_{n}^{+})_{\Gamma_{-}^{K}\cap\Gamma_{-}}\Big{)}\end{split}$
(6)

and the discrete right-hand side

$F_{n}(v_{n}):=\sum_{K\in\mathcal{T}_{n}}\Big{(}(f,v_{n})_{0,K}-((\mathbf{b}\cdot\mathbf{n})g_{D},v_{n}^{+})_{\Gamma_{-}^{K}\cap\Gamma_{-}}\Big{)}\quad\forall
v_{n}\in V_{n}.$ (7)

The dG method we consider reads

$\text{find }u_{n}\in V_{n}\quad\text{such that}\quad
B_{n}(u_{n},v_{n})=F_{n}(v_{n})\quad v_{n}\in V_{n}.$ (8)

Furthermore, we introduce the dG norm

$\begin{split}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}^{2}:=&\sum_{K\in\mathcal{T}_{n}}\Big{(}\|c_{0}v_{n}\|^{2}_{0,K}+\frac{1}{2}\|\sqrt{|\mathbf{b}\cdot\mathbf{n}|}v_{n}^{+}\|_{0,\Gamma_{-}^{K}\cap\Gamma}^{2}+\frac{1}{2}\|\sqrt{|\mathbf{b}\cdot\mathbf{n}|}v_{n}^{+}\|_{0,\Gamma_{+}^{K}\cap\Gamma}^{2}\\\
&\quad\quad\quad+\frac{1}{2}\|\sqrt{|\mathbf{b}\cdot\mathbf{n}|}\lfloor
v_{n}\rfloor\|^{2}_{\Gamma_{-}^{K}\setminus\Gamma_{-}}\Big{)}.\end{split}$ (9)

with $c_{0}$ defined as in (3).

It is easy to check that

$\begin{split}B_{n}(v_{n},v_{n})={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}^{2}\quad\quad\forall
v_{n}\in V_{n}.\end{split}$ (10)

The well-posedness of method (8) can be found, e.g., in [9, section 2.3].

##  3\. Standard analysis and suboptimality in terms of the polynomial degree

In this section, we recall the convergence analysis for method (8) from [14]
and where suboptimal estimates in terms of $p$ appear.

Preliminary, for all $K\in\mathcal{T}_{n}$, introduce the $L^{2}$ projector
$\Pi_{p}:L^{2}(K)\rightarrow\mathbb{Q}_{p}(K)$ through the affine mapping and
recall the standard $hp$-approximation estimates, see, e.g., [14, 11]: for any
function $u\in H^{\ell}(K)$ on the given element $K\in\mathcal{T}_{n}$, the
follow relations hold

$\begin{split}&\|u-\Pi_{p}u\|_{0,K}\lesssim\Big{(}\frac{h_{K}}{p+1}\Big{)}^{s}|u|_{s,K},\\\
&\|u-\Pi_{p}u\|_{0,\Gamma^{K}}\lesssim\Big{(}\frac{h_{K}}{p+1}\Big{)}^{s-\frac{1}{2}}|u|_{s,K},\end{split}$
(11)

with $s:=\min\\{p+1,\ell\\}$.

It is easy to check that method (8) is consistent, whence the following
Galerkin orthogonality follows:

$B_{n}(u-u_{n},v_{n})=0\quad\quad\forall v_{n}\in V_{n}.$ (12)

Then, we split error $u-u_{n}$ into $\eta+\xi$, where

$\eta:=u-\Pi_{p}u,\quad\quad\xi:=\Pi_{p}u-u_{n}.$ (13)

Using Galerkin orthogonality (12) and the properties of orthogonal projector
$\Pi_{p}$, we readily have the error equation

$0=B_{n}(u-u_{n},\xi)=B_{n}(\eta,\xi)+B_{n}(\xi,\xi)\quad\quad\Longrightarrow\quad\quad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\text{dG}}=-B_{n}(\eta,\xi).$
(14)

Since estimates on term $\eta$ are standard, error equation (14) allows us to
show a bound on term $\xi$, on which we now focus on.

We begin by computing the following error splitting:

$\begin{split}B_{n}(\eta,\xi)&=\sum_{K\in\mathcal{T}_{n}}\Big{(}(\mathbf{b}\cdot\nabla\eta,\xi)_{0,K}+(c\eta,\xi)_{0,K}\\\
&\quad\quad\quad-((\mathbf{b}\cdot\mathbf{n})\lfloor\eta\rfloor,\xi^{+})_{0,\Gamma_{-}^{K}\setminus\Gamma_{-}}-((\mathbf{b}\cdot\mathbf{n})\eta^{+},\xi^{+})_{0,\Gamma_{-}^{K}\cap\Gamma_{-}}\Big{)}\\\
&=\sum_{K\in\mathcal{T}_{n}}\Big{(}((c-\nabla\cdot\mathbf{b})\eta,\xi)_{0,K}-(\mathbf{b}\cdot\nabla\xi,\eta)_{0,K}\\\
&\quad\quad\quad+((\mathbf{b}\cdot\mathbf{n})\lfloor\xi\rfloor,\eta^{-})_{0,\Gamma_{-}^{K}\setminus\Gamma_{-}}+((\mathbf{b}\cdot\mathbf{n})\xi^{+},\eta^{+})_{0,\Gamma_{+}^{K}\cap\Gamma_{+}}\Big{)}=:T_{1}+T_{2}+T_{3}+T_{4}.\end{split}$
(15)

We show upper bounds for the four terms on the right-hand side of (15) and
anticipate that our analysis on term $T_{2}$ will lead to suboptimal bounds in
terms of $p$. Under further assumptions on vector $\mathbf{b}$, we shall
exhibit improved $p$-bounds in Section 4 below.

We begin with term $T_{1}$. Using that
$\mathbf{b}\in[W^{1,\infty}(\Omega)]^{2}$, $c\in L^{\infty}(\Omega)$, and
assumption (3), we obtain

$T_{1}\leq\sum_{K\in\mathcal{T}_{n}}\|c-\nabla\cdot\mathbf{b}\|_{\infty,K}\|\eta\|_{0,K}\|\xi\|_{0,K}\lesssim\|\eta\|_{0,\Omega}\|c_{0}\xi\|_{0,\Omega}\lesssim\|\eta\|_{0,\Omega}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}.$
(16)

As for terms $T_{3}$ and $T_{4}$, we have

$T_{3}+T_{4}\lesssim(\sum_{K\in\mathcal{T}_{n}}\|(\mathbf{b}\cdot\mathbf{n}_{K})^{\frac{1}{2}}\eta\|_{0,\Gamma^{K}}^{2})^{\frac{1}{2}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}\lesssim(\sum_{K\in\mathcal{T}_{n}}\|\eta\|_{0,\Gamma^{K}}^{2})^{\frac{1}{2}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}.$
(17)

As for term $T_{2}$, if we assume that $\mathbf{b}\cdot\nabla\xi\in V_{n}$ for
all $\xi\in V_{n}$, then $T_{2}=0$. Thus, inserting (16) and (17) in (15)
yields

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}\lesssim(\sum_{K\in\mathcal{T}_{n}}(\|\eta\|^{2}_{0,K}+\|\eta\|^{2}_{0,\Gamma^{K}}))^{\frac{1}{2}}\lesssim\Big{(}\frac{h_{K}}{p+1}\Big{)}^{s-\frac{1}{2}}|u|_{s,\Omega},$
(18)

Using a triangle inequality, and combining (11) with (18) leads to a
$p$-optimal error estimate.

Next, we focus on the case of nonzero $T_{2}$ to investigate the
$p$-suboptimality. Using the definition of $\eta$ in (13), and notably the
property of orthogonal projection $\Pi_{p}$, we can write

$T_{2}:=\sum_{K\in\mathcal{T}_{n}}\int_{K}(\mathbf{b}\cdot\nabla\xi)\eta=\sum_{K\in\mathcal{T}_{n}}\int_{K}(\mathbf{b}\cdot\nabla\xi-\mathbf{b}_{0}\cdot\nabla\xi)\eta,$

where $\mathbf{b}_{0}$ is the vector average over every $K$ of $\mathbf{b}$.
We deduce

$T_{2}\leq\sum_{K\in\mathcal{T}_{n}}\|\mathbf{b}-\mathbf{b}_{0}\|_{\infty,K}|\xi|_{1,K}\|\eta\|_{0,K}.$

On each element $K\in\mathcal{T}_{n}$, we have the following approximation
property and $hp$-polynomial inverse inequality:

$\|\mathbf{b}-\mathbf{b}_{0}\|_{\infty,K}\lesssim
h_{K}|\mathbf{b}|_{W^{1,\infty}(K)},\quad\quad|\xi|_{1,K}\lesssim\frac{p^{2}}{h_{K}}\|\xi\|_{0,K}.$

In the light of the two above bounds and (3), we have the following bound on
term $T_{2}$:

$T_{2}\lesssim
p^{2}\|\eta\|_{0,\Omega}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}.$
(19)

Inserting (16), (19), and (17) in (15), and using (11) yield

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}\lesssim(1+p^{2})\|\eta\|_{0,\Omega}+(\sum_{K\in\mathcal{T}_{n}}\|\eta\|^{2}_{0,\Gamma^{K}})^{\frac{1}{2}}\lesssim\left(\frac{h^{s}}{(p+1)^{s-2}}+\frac{h^{s-\frac{1}{2}}}{(p+1)^{s-\frac{1}{2}}}\right)|u|_{s,\Omega}.$
(20)

Using a triangle inequality, and combining (11) with (20) leads to the
following $p$-suboptimal error estimate:

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-u_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}\lesssim\left(\frac{h^{s-\frac{1}{2}}}{(p+1)^{s-2}}\right)|u|_{s,\Omega}.$
(21)

The above error bound is optimal in $h$ but suboptimal in terms of $p$ by
$3/2$ order, which is in accordance with [14, Remark $3.13$]. Notwithstanding,
such suboptimality is not observed in practice; see, e.g., [14, Numerical
Example $1$].

This motivates Section 4, where we shall exhibit improved estimates in terms
of $p$, under further assumptions on convection field $\mathbf{b}$.

##  4\. Improved bounds for special convection field

In this section, we show improved $p$-error estimates, under the following
assumption on convection field $\mathbf{b}$:

$\mathbf{b}\in[W^{2,\infty}(K)]^{d}\quad\forall
K\in\mathcal{T}_{n},\quad\quad\mathbf{b}=[b_{1}(x_{1}),b_{2}(x_{2}),\dots,b_{d}(x_{d})]^{T}.$
(22)

Since the $j$-th, $j=1,\dots,d$, component of $\mathbf{b}$ is assumed to be
single-valued in the $x_{j}$ variable, without loss of generality, we can
assume that element $K$ is the Cartesian product of intervals
$I_{j}:=(\alpha_{j},\beta_{j})$, $j=1,2,\dots,d$.

We can re-write term $T_{2}$ in (15) as

$T_{2}=\sum_{j=1}^{d}\int_{K}b_{j}(x_{j})\partial_{j}\xi\ \eta.$

Fix $j=1,\dots,d$. The $j$-th partial derivative of $\xi$ is a tensor
polynomial of degree $p-1$ along direction $x_{j}$ and $p$ along the others.
Define $\mathcal{I}_{j}b_{j}$ as the linear interpolant of $b_{j}$ at the end-
points of interval $I_{j}$:

$\big{(}b_{j}-\mathcal{I}_{j}b_{j}\Big{)}(\alpha_{j})=\big{(}b_{j}-\mathcal{I}_{j}b_{j}\Big{)}(\beta_{j})=0.$
(23)

Then, we clearly have that $\mathcal{I}_{j}b_{j}\
\partial_{i}\xi\in\mathbb{P}_{p}(I_{j})$ and $\mathcal{I}_{j}b_{j}\
\partial_{i}\xi\in V_{n}$. Consequently, the definition of orthogonal
projection $\Pi_{p}$ allows us to write

$\int_{K}b_{j}\partial_{j}\xi\
\eta=\int_{K}(b_{j}-\mathcal{I}_{j}b_{j})\partial_{j}\xi\ \eta.$

As $\alpha_{j}$ and $\beta_{j}$ denote the endpoints of interval $I_{j}$,
standard properties of one dimensional linear interpolation operators
guarantee the existence of $\widetilde{x}_{j}$ such that

$b_{j}(x_{j})-\mathcal{I}_{j}b_{j}(x_{j})=\frac{b_{j}^{(2)}(\widetilde{x}_{j})}{2}(x-a_{j})(x-b_{j})\quad\quad\forall
x_{j}\in I_{j}.$

Defining the standard quadratic bubble function on $I_{j}$

$w_{j}:=-(x-a_{j})(x-b_{j}),$

we can thus write

$\int_{K}(b_{j}-\mathcal{I}_{j}b_{j})\partial_{j}\xi\
\eta=\int_{K}\Big{(}\frac{b_{j}-\mathcal{I}_{j}b_{j}}{\sqrt{w_{j}}}\Big{)}\big{(}\sqrt{w_{j}}\partial_{j}\xi\big{)}\
\eta\leq\left\|\frac{b_{j}-\mathcal{I}_{j}b_{j}}{\sqrt{w_{j}}}\right\|_{\infty,K}\|\sqrt{w_{j}}\partial_{j}\xi\|_{0,K}\|\eta\|_{0,K}.$

Thanks to assumption (22) and an $hp$-polynomial inverse estimate involving
bubbles, see, e.g., [18, Lemma 3.42], we deduce

$\left\|\frac{b_{j}-\mathcal{I}_{j}b_{j}}{\sqrt{w_{j}}}\right\|_{\infty,K}\lesssim
h_{K}|b_{j}|_{W^{2,\infty}(K)},\quad\quad\|\sqrt{w_{j}}\partial_{j}\xi\|_{0,K}\lesssim\sqrt{p(p+1)}\|\xi\|_{0,K}.$

Collecting all the above estimates yields the following improved bound on term
$T_{2}$:

$T_{2}\lesssim
hp\|\eta\|_{0,\Omega}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}.$
(24)

Inserting (16), (24), and (17) in (15) gives the bound

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\xi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}\lesssim\frac{h^{s-\frac{1}{2}}}{(p+1)^{s-1}}|u|_{s,\Omega},$

which, combined with a triangle inequality and (11), eventually entails the
error estimate

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-u_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\text{dG}}\lesssim\frac{h^{s-\frac{1}{2}}}{(p+1)^{s-1}}|u|_{s,\Omega}.$

Compared with error estimate (21), the $p$-suboptimality improved by one
order.

##  5\. Conclusion

Employing a special class of convection fields, we derived an improved
$hp$-bound for dG methods discretising linear hyperbolic problems. The new
error bound is suboptimal by $1/2$ order in $p$ only, which improves the $3/2$
suboptimal order presented in [14]. Needless to say, the results in this work
do not provide full answers to the open questions in [14], notably on the
mismatch between the theoretical and numerical results on $p$-convergence.
However, they shed additional light on such issues and shows the possibility
of deriving sharper $hp$-error bounds for dG methods.

## References

  * [1] R. A. Adams and J. J. F. Fournier.  Sobolev Spaces, volume 140.  Academic Press, 2003.
  * [2] C. E. Baumann and J. T. Oden.  A discontinuous $hp$-finite element method for convection-diffusion problems.  Comput. Methods Appl. Mech. Engrg., 175(3-4):311–341, 1999.
  * [3] K. S. Bey and J. T. Oden.  $hp$-version discontinuous Galerkin methods for hyperbolic conservation laws.  Comput. Methods Appl. Mech. Engrg., 133(3-4):259–286, 1996.
  * [4] A. Cangiani, Z. Dong, and E. H. Georgoulis.  $hp$-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements.  arXiv:1906.01715, 2019.
  * [5] A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston.  $hp$-version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes.  ESAIM Math. Model. Numer. Anal., 50(3):699–725, 2016.
  * [6] A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston.  hp-Version Discontinuous Galerkin Methods on Polygonal and Polyhedral Meshes.  Springer, 2017.
  * [7] A. Cangiani, E. H Georgoulis, and P. Houston.  $hp$-version discontinuous Galerkin methods on polygonal and polyhedral meshes.  Math. Models Methods Appl. Sci., 24(10):2009–2041, 2014.
  * [8] B. Cockburn, G. E. Karniadakis, and C.-W. Shu.  The development of discontinuous Galerkin methods.  In Discontinuous Galerkin methods (Newport, RI, 1999), volume 11 of Lect. Notes Comput. Sci. Eng., pages 3–50. Springer, Berlin, 2000.
  * [9] D. A. Di Pietro and A. Ern.  Mathematical aspects of discontinuous Galerkin methods, volume 69.  Springer Science & Business Media, 2011.
  * [10] V. Dolejší and M. Feistauer.  Discontinuous Galerkin method, volume 48 of Springer Series in Computational Mathematics.  Springer, Cham, 2015.
  * [11] Z. Dong.  On the exponent of exponential convergence of $p$-version FEM spaces.  Adv. Comput. Math., 45(2):757–785, 2019.
  * [12] E. H. Georgoulis, E. Hall, and P. Houston.  Discontinuous Galerkin methods for advection-diffusion-reaction problems on anisotropically refined meshes.  SIAM J. Sci. Comput., 30(1):246–271, 2007/08.
  * [13] P. Houston, Ch. Schwab, and E. Süli.  Stabilized $hp$-finite element methods for first-order hyperbolic problems.  SIAM J. Numer. Anal., 37(5):1618–1643, 2000.
  * [14] P. Houston, Ch. Schwab, and E. Süli.  Discontinuous $hp$-finite element methods for advection-diffusion-reaction problems.  SIAM J. Numer. Anal., 39(6):2133–2163, 2002.
  * [15] P. Houston and E. Süli.  Stabilised $hp$-finite element approximation of partial differential equations with nonnegative characteristic form.  Computing, 66(2):99–119, 2001.
  * [16] J. Nitsche.  Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind.  Abh. Math. Sem. Univ. Hamburg, 36:9–15, 1971.
  * [17] W.H. Reed and T.R. Hill.  Triangular mesh methods for the neutron transport equation.  Technical Report LA-UR-73-479 Los Alamos Scientific Laboratory, 1973.
  * [18] R. Verfürth.  A posteriori error estimation techniques for finite element methods.  OUP Oxford, 2013.