# Interaction-Aware Personalized Vehicle Trajectory Prediction Using Temporal
Graph Neural Networks

Amr Abdelraouf, Rohit Gupta, Kyungtae Han A. Abdelraouf, R. Gupta, and K. Han
are with InfoTech Labs, Toyota Motor North America R&D, 465 Bernardo Avenue,
Mountain View, CA 94043 (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>kyungtae.han@toyota.com).

###### Abstract

Accurate prediction of vehicle trajectories is vital for advanced driver
assistance systems and autonomous vehicles. Existing methods mainly rely on
generic trajectory predictions derived from large datasets, overlooking the
personalized driving patterns of individual drivers. To address this gap, we
propose an approach for interaction-aware personalized vehicle trajectory
prediction that incorporates temporal graph neural networks. Our method
utilizes Graph Convolution Networks (GCN) and Long Short-Term Memory (LSTM) to
model the spatio-temporal interactions between target vehicles and their
surrounding traffic. To personalize the predictions, we establish a pipeline
that leverages transfer learning: the model is initially pre-trained on a
large-scale trajectory dataset and then fine-tuned for each driver using their
specific driving data. We employ human-in-the-loop simulation to collect
personalized naturalistic driving trajectories and corresponding surrounding
vehicle trajectories. Experimental results demonstrate the superior
performance of our personalized GCN-LSTM model, particularly for longer
prediction horizons, compared to its generic counterpart. Moreover, the
personalized model outperforms individual models created without pre-training,
emphasizing the significance of pre-training on a large dataset to avoid
overfitting. By incorporating personalization, our approach enhances
trajectory prediction accuracy.

##  I Introduction

The advent of vehicle-generated big data has sparked considerable interest in
data-driven personalized advanced driver assistance systems (ADAS) [1]. By
leveraging the wealth of personalized driving patterns and insights expressed
by drivers, personalization has the potential to enhance the performance of
ADAS systems greatly. This, in turn, leads to improved driving experiences,
increased driver acceptance, and greater utilization of ADAS functionalities.
In recent years, many personalized ADAS applications have been proposed,
including Adaptive Cruise Control [2], Forward Collision Warning [3], Lane
Keeping Assistance[4], among many others. Furthermore, personalization extends
to battery electric vehicle (BEV) ADAS applications. Personalized range
estimation, in particular, holds significant promise in mitigating range
anxiety, as it enables more accurate predictions of available driving range
for individual drivers’ given their unique driving styles [5].

Vehicle trajectory prediction is a foundational component in many ADAS
applications and autonomous vehicle (AV) systems [6]. It plays a significant
role in various safety-critical ADAS applications, such as collision warning
[7], automated braking [8, 9], and lane change prediction [10]. The accuracy
of trajectory prediction is crucial for enhancing the effectiveness of these
safety-critical applications [11], enabling them to provide timely warnings
while minimizing false positives. Moreover, in the context of autonomous
vehicles, trajectory prediction plays a vital role in ensuring safe motion
planning and overall system safety [12]. Anticipating and responding to
dynamic traffic situations helps to reduce the risk of accidents and instills
trust in AV systems. Additionally, trajectory prediction can significantly
contribute to Vehicle-to-Vehicle (V2V) communication, particularly in intent-
sharing [13] and negotiation applications [14]. By facilitating vehicular
communication of their intentions to others, trajectory prediction for intent-
sharing promotes cooperative behavior in mixed traffic, improves traffic flow,
and enhances overall safety.

Recent advancements in trajectory prediction have emphasized the significance
of vehicle interactions [12]. Interaction-aware trajectory prediction
considers the spatio-temporal relationship between the target vehicle and its
surrounding neighbors. In the recent past, some of the most effective
approaches for short-term trajectory prediction combined graph neural networks
(GNNs) to model the spatial domain and recurrent neural networks (RNNs) to
capture the temporal domain for vehicle interactions [15, 16].

Despite the substantial body of research dedicated to vehicle trajectory
prediction, a notable research gap persists in personalized vehicle trajectory
prediction. Specifically, there is limited exploration of interaction-aware
personalized trajectory prediction methods [17]. A major challenge stems from
the lack of individual driver trajectory data which captures the driver’s
trajectory in addition to the driver’s surrounding vehicle trajectories for an
extended period of time.

This paper presents a novel approach for generating personalized interaction-
aware vehicle trajectory predictions. Our method builds upon state-of-the-art
techniques and combines graph convolution networks (GCN) with Long Short-Term
Memory (LSTM) models to effectively capture the dynamics of vehicle
interactions. To obtain individual driver trajectories which can support
interaction-aware modeling, we utilized human-in-the-loop simulation using the
CARLA driving simulator [18]. The proposed personalization framework leverages
transfer learning. Initially, a base model was pre-trained on a large-scale
vehicle trajectory dataset named CitySim [19]. Subsequently, the model was
fine-tuned for each driver using their individual driving data, resulting in a
unique model for each driver tailored to their specific driving
characteristics. Experiments indicated that our proposed personalized approach
generates more accurate future trajectories compared to generic trajectory
prediction. The proposed method can be used to personalize various ADAS
applications, such as collision warning, enhancing their accuracy while
reducing their false positive rate. The improvements can contribute to a
better driving experience by improving safety, reliability, and driver trust.

In summary, the contribution of this paper is three-fold:

  1. 1.

We introduced a new method for predicting future vehicle trajectories that
incorporates personalized driving characteristics. Our approach employed
transfer learning to first learn generic interaction-aware driving patterns,
which were then fine-tuned using driver-specific trajectories.

  2. 2.

To obtain driver-specific trajectories which can support interaction-aware
trajectory prediction, we conducted a human-in-the-loop simulation experiment
that captured naturalistic driving behavior and surrounding vehicle
trajectories.

  3. 3.

We demonstrated the effectiveness of our approach through experiments and
analysis, which showed that our interaction-aware personalized trajectory
prediction outperforms generic trajectory prediction methods.

The remainder of this paper is structured as follows: Section II provides an
overview of related work, Section III describes the personalization scheme and
model architecture, Section IV describes the personalized driving data and
large-scale trajectory dataset used for experimentation, Section V presents
the experimental results, and Section VI concludes the paper and discusses
future work.

##  II Related Work

###  II-A Vehicle Trajectory Prediction

In recent years, there have been several comprehensive surveys focusing on
vehicle trajectory prediction [20, 21, 12]. Notably, Huang et al. [12]
conducted a recent survey that examined various trajectory prediction
methodologies. These approaches were categorized into physics-based methods,
classic machine learning methods, reinforcement learning-based methods, and
deep learning-based methods. Physics-based models are often regarded as the
simplest approach to trajectory prediction. While these methods may be
suitable for short prediction horizons, their accuracy deteriorates as the
output horizon increases. Physics-based methods encompass a range of models,
from basic ones like constant velocity and constant acceleration to more
sophisticated techniques like Kalman Filtering (KF) methods [22]. Classic
machine learning methods utilize a data-driven approach for trajectory
prediction. These algorithms incorporate various techniques such as Gaussian
Process (GP) [23, 24], Support Vector Machine (SVM) [25, 26], and Hidden
Markov Models (HMMs) [27, 28]. Reinforcement learning approaches leverage
expert demonstrations to learn optimal driving policies by maximizing the
expected reward. Example algorithms include Inverse Reinforcement Learning
(IRL) [29] and Generative Adversarial Imitation Learning (GAIL) [30].

In the recent past, deep learning-based approaches have demonstrated state-of-
the-art performance for vehicle trajectory prediction, particularly given the
increased availability of large trajectory datasets. The superior performance
can be attributed to the innate effectiveness of deep learning at capturing
the intricate patterns present in large datasets [21, 31]. Previous research
has focused on utilizing sequential deep neural networks to address the
temporal nature of the vehicle trajectory prediction problem. Methods such as
Recurrent Neural Networks (RNNs) [32, 33] and sequence-to-sequence
Transformers [34] have been employed for this purpose.

Vehicle interactions play a significant role in vehicle trajectory prediction.
To effectively incorporate these interactions within the deep learning
modeling framework, prior research has frequently integrated a sequential
neural network module with a spatial counterpart to create an end-to-end
spatio-temporal trajectory prediction network. For instance, Convolutional
Neural Networks (CNNs) have often been utilized to model the spatial context
for the trajectory prediction problem [35, 36, 37]. Deo et al. [35] employed
CNNs and LSTM to create a convolutional social pooling algorithm. The authors
utilized a CNN-LSTM architecture to capture the trajectories of neighboring
vehicles, extracting embeddings that were then utilized to generate multi-
modal trajectories based on vehicle maneuvers. In recent publications, Graph
Neural Networks (GNNs) have gained popularity for modeling the spatial domain
of vehicle interactions [15, 16, 38, 39]. GNNs are well-suited for this task
as they effectively capture the relationships among vehicles in the
neighborhood. Each vehicle is represented as a graph node, and the
interactions between vehicles are represented by edges in the graph.

Figure 1: Proposed trajectory prediction network architecture

###  II-B Personalized Trajectory Prediction

A few research efforts proposed methods for personalized trajectory
prediction. However, due to lack of long-term interaction-aware personal
trajectory data, many of these researches proposed to first group driving
behavior into distinct clusters [40, 41]. For example, Xing et al. [40]
introduced an approach for trajectory prediction that starts by clustering the
driving style of the target vehicle into one of three categories:
conservative, moderate, or aggressive using Gaussian Mixture Models (GMM). The
method utilized a pre-trained LSTM network and subsequently fine-tuned a
regression head tailored to each of the three driving styles.

Other approaches aimed to learn individual driver profiles by constructing
personalized reward functions using Inverse Reinforcement Learning (IRL) [42,
43]. However, the approach in [42] was developed without leveraging long-term
personalized trajectories and was tested on a few seconds of individual
driving periods. Moreover, the approach presented in [43] did not consider the
interactions between the target vehicle and its neighboring vehicles.

This paper aimed to address the literature gap in personalized trajectory
prediction by first collecting extended personalized vehicle trajectory. To
facilitate interaction-aware trajectory prediction, the dataset captured
target vehicles’ surrounding traffic trajectory data as well. Next, an
interaction-aware approach which utilized temporal graph neural networks was
introduced. Transfer learning techniques were incorporated to create
personalized models for each driver. Finally, the effectiveness of the
approach was evaluated using the collected individual driving data.

##  III Methodology

The proposed model is illustrated in Figure 1. It presents an architecture
designed to predict the future trajectory of the target vehicle by processing
a traffic graph centered around it. The model follows an encoder-decoder
structure, with a GCN-LSTM network serving as the encoder and an LSTM network
serving as the decoder.

###  III-A Traffic Graph

The input graph is represented as an undirected graph $G=(V,E)$, where the set
of nodes $v_{i}\in V$ of size $N$ consists of the target vehicle and the
vehicles in its neighborhood. The set of edges $(v_{i},v_{j})\in E$ represents
the connections between vehicles in close proximity, satisfying two
conditions: they are within 1 lane of each other and their distance is less
than or equal to a threshold $\tau$.

Each vehicle $v_{i}$ is associated with a node feature vector
$X_{i}\in\mathbb{R}^{T_{\text{in}}+1\times F}$, containing its historical
trajectory for the past $T_{\text{in}}$ time steps in addition to the current
time step. Each element within a time step contains the vehicle’s
$x$-coordinate, $y$-coordinate, and speed. The historical trajectory $X_{i}$
of a vehicle $v_{i}$ is defined as:

$X_{i}=\\{C^{t-T_{in}}_{i},\ldots,C^{t-1}_{i},C^{t}_{i}\\}$ (1)

where each element $C^{t}_{i}$ is defined as:

$C^{t}_{i}=\\{x^{t}_{i},y^{t}_{i},s^{t}_{i}\\}$ (2)

###  III-B Interaction Modeling

The proposed encoder model combined an LSTM network for temporal feature
embedding and a GCN network for capturing the spatial interaction between
neighboring vehicles. The LSTM network was specifically designed to capture
sequential dependencies, enabling it to extract valuable features from each
vehicle’s historical trajectory data. By sharing weights across LSTMs, the
model ensured consistent representations. The LSTM generated 1-dimensional
embeddings $X^{embed}$, which served as node feature inputs for the GCN
layers.

Graph Convolution Networks, initially introduced in [44], were specifically
designed as an alternative to CNNs for effective modeling of graph-based
structures. GCNs require two inputs: an adjacency matrix and a set of
features. In our implementation, the adjacency matrix $A\in\mathbb{R}^{N\times
N}$ represents the connectivity between nodes in the graph, where $N$
corresponds to the number of nodes. $A_{i,j}=A_{j,i}=1$ if there exists an
edge $(v_{i},v_{j})$, and 0 otherwise. The set of features
$X\in\mathbb{R}^{N\times D}$ represents the input feature vector for each node
in the graph, where $D$ denotes the number of features. Notably, in our
implementation, $X$ is constructed from the output of the LSTM embeddings
$X^{embed}$ and $D$ is equal to the size of the embedding output from the
encoder LSTM.

Building upon the methodology described in [16], our approach incorporates a
2-layer GCN module, which performs the following operation:

$\text{GCN}(X^{embed},A)=\hat{D}^{-\frac{1}{2}}\hat{A}\hat{D}^{-\frac{1}{2}}H^{(0)}W^{(1)}$
(3)

Here, $H^{(0)}$ is computed as:

$H^{(0)}=\text{RELU}(\hat{D}^{-\frac{1}{2}}\hat{A}\hat{D}^{-\frac{1}{2}}X^{embed}W^{(0)})$
(4)

In the equations above, a node-wise self-connection was established by adding
the input adjacency matrix $A$ to the identity matrix $I$, denoted as
$\hat{A}=A+I$. This self-connection allows vehicle nodes to incorporate their
own historical trajectories during GCN propagation. Furthermore, the weights
of the adjacency matrix were normalized by the node degree using the diagonal
matrix $\hat{D}_{i,i}=\sum_{j}^{N}A_{i,j}$. Finally, $W^{(0)}$ and $W^{(1)}$
represented the trainable parameters of the first and second GCN layers,
respectively.

###  III-C Future Trajectory Generation

The decoder layer of the proposed model is responsible for generating the
future trajectory of the target vehicle. It begins by extracting the GCN
embedding state of the target node and utilizing it as input for an LSTM
module. The decoder LSTM captures the temporal dependencies between the output
time steps. Subsequently, the model generates the output
$Y_{i}\in\mathbb{R}^{T_{\text{out}}\times 2}$, which contains the $x$ and $y$
coordinates of the future trajectory for each time step in the output horizon
$T_{\text{out}}$.

The future trajectory $Y_{i}$ for a target vehicle $v_{i}$ is defined as
follows:

$Y_{i}=\\{\\{x^{t+1}_{i},y^{t+1}_{i}\\},\ldots,\\{x^{t+T_{\text{out}}}_{i},y^{t+T_{\text{out}}}_{i}\\}\\}$
(5)

###  III-D Model Training

The proposed model was trained using backpropagation, with the loss function
$\mathcal{L}$ defined as:

$\mathcal{L}(Y,\hat{Y})=\frac{1}{T_{\text{out}}\cdot
M}\sum_{i=1}^{M}\sum_{t=1}^{T_{\text{out}}}\left|x_{i}^{t}-\hat{x}_{i}^{t}\right|+\left|y_{i}^{t}-\hat{y}_{i}^{t}\right|$
(6)

In this equation, $x_{i}^{t}$ and $y_{i}^{t}$ represent the true
$x$-coordinate and $y$-coordinate of example $i$ at time step $t$,
respectively. The predicted $x$-coordinate and $y$-coordinate at time step $t$
are denoted as $\hat{x}_{i}^{t}$ and $\hat{y}_{i}^{t}$, respectively. $M$ is
the number of examples used for training and $T_{out}$ is the output
prediction horizon.

###  III-E Personalization

Figure 2: Transfer learning-based trajectory prediction personalization

The proposed personalization scheme, depicted in Figure 2, involves two main
steps. Initially, a base model is constructed by training the proposed
architecture on a large-scale vehicle trajectory dataset. Subsequently,
personalized models are created for each driver using transfer learning. This
is achieved by loading the weights of the pre-trained base model and fine-
tuning it using the driver’s own trajectories. During fine-tuning, the encoder
network weights remain frozen to preserve the learned spatio-temporal
interaction embedding representations, while the decoder weights are updated.
This approach allows for personalized trajectory prediction by leveraging the
base model’s learned features and adapting them to individual driver
characteristics.

##  IV Data Description

TABLE I: Trajectory prediction RMSE (m) |  |  |  |  |
---|---|---|---|---|---
Prediction Horizon (s) |  CV |  Seq2seq LSTM |  Generic GCN-LSTM |  Individual GCN-LSTM |  Personalized GCN-LSTM
1 |  1.427 |  1.168 |  0.834 |  1.063 |  0.763
2 |  2.360 |  1.546 |  1.344 |  1.533 |  1.191
3 |  3.573 |  2.741 |  1.820 |  1.965 |  1.529
4 |  5.439 |  3.275 |  2.547 |  2.670 |  2.022
5 |  6.109 |  5.885 |  5.351 |  4.486 |  3.024

###  IV-A Personalized Trajectories

To overcome the limited availability of long-term naturalistic driving data
that can support interaction-aware trajectory prediction, we adopted a human-
in-the-loop driver simulation approach. The CARLA driving simulator [18]
served as our platform, providing pre-implemented vehicle models, control and
physics, road geometry, and traffic management modules. In our hardware setup,
we utilized a G29 Logitech steering wheel, floor pedals, and a driving seat.

Our experimentation focused on highway scenarios. Therefore, we utilized Town
04 in CARLA, which offers a continuous 3-lane highway configured in the shape
of a figure-eight. We extracted the trajectories from the non-curved segments
of the highway for use in our experiment as illustrated in Figure 3.

Figure 3: CARLA Town 04 map annotated with the portions of the road used for
data collection

We gathered data from a total of five distinct drivers for our study. Each
driver received instructions to adhere to traffic laws and drive in a typical
manner as they would on public roads. To capture a comprehensive range of
driving behaviors, each participant completed four driving rounds, with each
round lasting 10 minutes. Thus, we obtained a total of 40 minutes of
personalized naturalistic driving data per driver. To capture the range of
driving behaviors under diverse traffic conditions, we systematically varied
the number of vehicles in the background across the rounds. Specifically, we
set the number of vehicles in the background traffic with three different
settings: 100, 200, and 300 vehicles. This variation enabled us to capture how
different drivers adapt and respond to different traffic densities in our
personalized dataset.

###  IV-B Large-Scale Trajectories

The CitySim dataset [19], specifically the Freeway B scenario, was used for
base model pre-training. This dataset includes bird’s-eye view vehicle
trajectories obtained from drone videos. Multiple drones were flown
concurrently to capture continuous trajectories along an extended highway
segment. The Freeway B scenario is a long continuous stretch of a non-curved
highway. The dataset covered a duration of 35 minutes and included trajectory
data from a total of 7,307 vehicles in various traffic conditions. The average
trajectory length was 45 seconds, with a standard deviation of 24 seconds. The
dataset provides a diverse range of traffic conditions, allowing the base
model to learn from a wide range of scenarios encountered by numerous drivers
on the freeway.

##  V Experiments and Results

###  V-A Experimental Setup

In this experiment, 10 minutes from each driver’s personalized trajectories
were used for evaluation. To generate future trajectories, 5 seconds of
historical trajectories were utilized to predict the subsequent 5 seconds of
future trajectories for the target vehicle. The time step interval employed
was 0.5 seconds. Thus, the proposed model parameters $T_{in}=T_{out}=10$.
Furthermore, in the input traffic graph, the distance threshold for edge
creation $\tau$ was set to 100 feet.

###  V-B Evaluation Metrics

Similar to previous trajectory prediction works, we adopted Root Mean Square
Error (RMSE) for evaluation. RMSE measures the difference between the ground
truth and the predicted value at time step $t$ using the following equation:

$\mathrm{RMSE}^{t}=\sqrt{\frac{1}{M}\sum_{i=1}^{M}(x_{i}^{t}-\hat{x}_{i}^{t})^{2}+(y_{i}^{t}-\hat{y}_{i}^{t})^{2}}$
(7)

In this equation, $x_{i}^{t}$ and $y_{i}^{t}$ represent the true coordinates
at time step $t$, and $\hat{x}_{i}^{t}$ and $\hat{y}_{i}^{t}$ represent the
predicted coordinates at time step $t$. $M$ is the number of samples. We
compare the RMSE values at output horizons of 1 second to 5 seconds using
1-second intervals. The RMSE values were reported in meters.

###  V-C Baseline Models

Several baseline models were used to compare the performance of the proposed
model. The results compare between the following models:

  * •

Constant Velocity (CV): A constant velocity Kalman filter was employed to
predict the future trajectory. It assumes constant velocity motion.

  * •

Sequence-to-Sequence (Seq2seq) LSTM: This approach utilizes an LSTM encoder
and LSTM decoder network. It does not incorporate a graph module for spatial
interaction modeling.

  * •

Generic GCN-LSTM: The base GCN-LSTM model was employed without any fine-tuning
for individual drivers.

  * •

Individual GCN-LSTM: Separate GCN-LSTM models were trained for each individual
driver using their respective data. No pre-training was conducted for these
models.

###  V-D Results

Table I compares the results of the proposed model against the baseline
models. It contains output RMSEs for different prediction horizons. The
results indicate that, while the output of the CV model is relatively
comparable for shorter prediction horizons, it rapidly deteriorates as the
prediction horizon increases. The generic GCN-LSTM outperforms the CV model
and the Seq2seq LSTM model, highlighting the effectiveness of the spatio-
temporal interaction-aware modeling approach.

The results in Table I indicate that the personalized GCN-LSTM surpasses its
generic counterpart, emphasizing the importance of personalization in vehicle
trajectory prediction. Furthermore, the corresponding graph presented in
Figure 4 illustrates the RMSE reduction percentage achieved by the
personalized GCN-LSTM model compared to the generic model. The results
indicate that the further out the prediction horizon, the higher the error
reduction margin. Notably, longer prediction horizons are more challenging to
forecast given the higher level of uncertainty. Nevertheless, the results
indicate that the personalized model improvement is more pronounced for longer
horizons.

In order to assess the significance of the pre-training phase for the
personalized GCN-LSTM model, its performance was compared to individual GCN-
LSTM models trained exclusively on each driver’s data. The findings presented
in Table I demonstrate that the personalized GCN-LSTM model surpasses the
individual model results. This observation implies that the individual model
may encounter overfitting issues when solely trained on the driver’s data. By
leveraging a diverse range of driving scenarios during pre-training, the
personalized model showcases superior generalization capabilities compared to
the individual model.

Figure 4: Percentage RMSE reduction of the personalized GCN-LSTM model
compared to the generic GCN-LSTM model for different prediction horizons

###  V-E Impact of Personalization

Figure 5: Impact of driving duration used for personalization on trajectory
prediction RMSE

To further investigate the impact of personalization on vehicle trajectory
prediction, the number of driving minutes used for fine-tuning was varied. As
mentioned earlier, 40 minutes of personalized naturalistic driving data were
collected for each driver, with 10 minutes allocated for testing.
Subsequently, the number of driving minutes for model fine-tuning was varied
from 5 to 30 minutes, and the RMSE was computed for each variation.

Figure 5 illustrates the relationship between the driving duration used for
fine-tuning and the trajectory prediction RMSE for different horizons. The
graph demonstrates a clear trend where an increase in the number of fine-
tuning minutes leads to a decrease in RMSE. Furthermore, in alignment with the
results depicted in Figure 4, it can be noticed that personalization has a
more pronounced effect on longer prediction horizons. The graph trends suggest
that a longer duration of personalized driving data positively impacts the
trajectory prediction accuracy, particularly for longer prediction horizons.

##  VI Conclusion and Future Work

In this paper, we proposed a personalized interaction-aware trajectory
prediction method that leverages GCNs and LSTMs to model the spatio-temporal
vehicle interactions between the target vehicle and its neighboring vehicles.
By employing transfer learning, we personalized the proposed model by pre-
training it on a large-scale vehicle trajectory dataset and subsequently fine-
tuning it for each unique driver, thus creating a personalized model that
captures the unique patterns of individual drivers. To overcome the lack of
availability of extended interaction-based personalized trajectories, we
employed human-in-the-loop simulation to collect interaction-based
personalized trajectories.

Through experiments, we demonstrated that the personalized GCN-LSTM model
outperforms its generic counterpart, particularly for longer prediction
horizons. This highlights the effectiveness of our personalized approach in
capturing the driver-specific characteristics and improving trajectory
predictions. Furthermore, our results indicated that the personalized model
surpasses individual models created without pre-training, emphasizing the
significance of pre-training on a large dataset to mitigate overfitting issues
and enhance prediction accuracy.

The implications of personalized trajectory prediction are profound for
various ADAS systems such as collision warning and lane keeping assist. By
tailoring the predictions to individual patterns, the accuracy is improved,
and the false positive rate is reduced. Consequently, these enhancements lead
to greater safety and reliability within the systems, while instilling driver
trust. Furthermore, personalized trajectory prediction can affect vehicle-to-
vehicle (V2V) communication applications, specifically intent-sharing and road
resource negotiations between vehicles. This is particularly relevant in mixed
traffic scenarios, where both connected automated vehicles (CAVs) and manually
driven connected vehicles interact.

Future work should focus on further validating the significance of
personalization in interaction-aware trajectory prediction by collecting a
larger dataset of personalized trajectories. In addition, although simulation
data offers a convenient way to demonstrate the effectiveness of
personalization, it is crucial to address the domain shift between simulation
and real driving. Integrating real-world data into the training process will
enhance the model’s ability to generalize and accurately predict trajectories
in practical driving scenarios.

## Acknowledgment

The contents of this work only reflect the views of the authors, who are
responsible for the facts and the accuracy of the data presented herein. The
contents do not necessarily reflect the official views of Toyota Motor North
America.

The authors would like to sincerely thank Yongkang Liu and Ahmadreza
Moradipari for their participation in the human-in-the-loop experiments and
Hazem Abdelkawy for his support.

## References

  * [1] M. Hasenjäger, M. Heckmann, and H. Wersing, “A survey of personalization for advanced driver assistance systems,” _IEEE Transactions on Intelligent Vehicles_ , vol. 5, no. 2, pp. 335–344, 2019.
  * [2] Z. Zhao, Z. Wang, K. Han, R. Gupta, P. Tiwari, G. Wu, and M. J. Barth, “Personalized car following for autonomous driving with inverse reinforcement learning,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 2891–2897.
  * [3] J. Wang, C. Yu, S. E. Li, and L. Wang, “A forward collision warning algorithm with adaptation to driver behaviors,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 17, no. 4, pp. 1157–1167, 2015.
  * [4] W. Wang, D. Zhao, W. Han, and J. Xi, “A learning-based approach for lane departure warning systems with a personalized driver model,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 10, pp. 9145–9157, 2018.
  * [5] H. Shen, X. Zhou, H. Ahn, M. Lamantia, P. Chen, and J. Wang, “Personalized velocity and energy prediction for electric vehicles with road features in consideration,” _IEEE Transactions on Transportation Electrification_ , 2023\.
  * [6] A. Abdelraouf, M. Abdel-Aty, Z. Wang, and O. Zheng, “Trajectory prediction for vehicle conflict identification at intersections using sequence-to-sequence recurrent neural networks,” _arXiv preprint arXiv:2210.08009_ , 2022.
  * [7] N. Lyu, J. Wen, Z. Duan, and C. Wu, “Vehicle trajectory prediction and cut-in collision warning model in a connected vehicle environment,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 2, pp. 966–981, 2020.
  * [8] L. Zhang, J. Yuan, K. Yuan, J. Hong, H. Ding, and H. Chen, “Automated braking decision and control for pedestrian collision avoidance based on risk assessment,” _IEEE Intelligent Transportation Systems Magazine_ , vol. 14, no. 3, pp. 41–61, 2021.
  * [9] Z. Wang, O. Zheng, L. Li, M. Abdel-Aty, C. Cruz-Neira, and Z. Islam, “Towards next generation of pedestrian and connected vehicle in-the-loop research: A digital twin co-simulation framework,” _IEEE Transactions on Intelligent Vehicles_ , 2023.
  * [10] X. Liao, Z. Wang, X. Zhao, Z. Zhao, K. Han, P. Tiwari, M. J. Barth, and G. Wu, “Online prediction of lane change with a hierarchical learning-based approach,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 948–954.
  * [11] N. Mahmoud, M. Abdel-Aty, and A. Abdelraouf, “The impact of target speed on pedestrian, bike, and speeding crash frequencies,” _Accident Analysis & Prevention_, vol. 192, p. 107263, 2023.
  * [12] Y. Huang, J. Du, Z. Yang, Z. Zhou, L. Zhang, and H. Chen, “A survey on trajectory-prediction methods for autonomous driving,” _IEEE Transactions on Intelligent Vehicles_ , vol. 7, no. 3, pp. 652–674, 2022.
  * [13] H. M. Wang, S. S. Avedisov, T. G. Molnár, A. H. Sakr, O. Altintas, and G. Orosz, “Conflict analysis for cooperative maneuvering with status and intent sharing via v2x communication,” _IEEE Transactions on Intelligent Vehicles_ , vol. 8, no. 2, pp. 1105–1118, 2022.
  * [14] A. Correa, S. S. Avedisov, M. Sepulcre, A. H. Sakr, R. Molina-Masegosa, O. Altintas, and J. Gozalvez, “On the impact of v2x-based maneuver coordination on the traffic,” in _2021 IEEE 93rd Vehicular Technology Conference (VTC2021-Spring)_. IEEE, 2021, pp. 1–5.
  * [15] X. Li, X. Ying, and M. C. Chuah, “Grip++: Enhanced graph-based interaction-aware trajectory prediction for autonomous driving,” _arXiv preprint arXiv:1907.07792_ , 2019.
  * [16] Z. Zhao, H. Fang, Z. Jin, and Q. Qiu, “Gisnet: Graph-based information sharing network for vehicle trajectory prediction,” in _2020 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2020, pp. 1–7.
  * [17] Y. Ma, W. Ye, X. Cao, A. Abdelraouf, K. Han, R. Gupta, and Z. Wang, “Cemformer: Learning to predict driver intentions from in-cabin and external cameras via spatial-temporal transformers,” _arXiv preprint arXiv:2305.07840_ , 2023.
  * [18] A. Dosovitskiy, G. Ros, F. Codevilla, A. Lopez, and V. Koltun, “Carla: An open urban driving simulator,” in _Conference on robot learning_. PMLR, 2017, pp. 1–16.
  * [19] O. Zheng, M. Abdel-Aty, L. Yue, A. Abdelraouf, Z. Wang, and N. Mahmoud, “Citysim: A drone-based vehicle trajectory dataset for safety oriented research and digital twins,” _arXiv preprint arXiv:2208.11036_ , 2022.
  * [20] S. Lefèvre, D. Vasquez, and C. Laugier, “A survey on motion prediction and risk assessment for intelligent vehicles,” _ROBOMECH journal_ , vol. 1, no. 1, pp. 1–14, 2014.
  * [21] S. Mozaffari, O. Y. Al-Jarrah, M. Dianati, P. Jennings, and A. Mouzakitis, “Deep learning-based vehicle behavior prediction for autonomous driving applications: A review,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 1, pp. 33–47, 2020.
  * [22] S. Ammoun and F. Nashashibi, “Real time trajectory prediction for collision risk estimation between vehicles,” in _2009 IEEE 5th International Conference on Intelligent Computer Communication and Processing_. IEEE, 2009, pp. 417–422.
  * [23] Q. Tran and J. Firl, “Online maneuver recognition and multimodal trajectory prediction for intersection assistance using non-parametric regression,” in _2014 ieee intelligent vehicles symposium proceedings_. IEEE, 2014, pp. 918–923.
  * [24] D. Vasquez and T. Fraichard, “Motion prediction for moving objects: a statistical approach,” in _IEEE International Conference on Robotics and Automation, 2004. Proceedings. ICRA’04. 2004_ , vol. 4. IEEE, 2004, pp. 3931–3936.
  * [25] G. S. Aoude, B. D. Luders, K. K. Lee, D. S. Levine, and J. P. How, “Threat assessment design for driver assistance system at intersections,” in _13th international ieee conference on intelligent transportation systems_. IEEE, 2010, pp. 1855–1862.
  * [26] P. Kumar, M. Perrollaz, S. Lefevre, and C. Laugier, “Learning-based approach for online lane change intention prediction,” in _2013 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2013, pp. 797–802.
  * [27] N. Deo, A. Rangesh, and M. M. Trivedi, “How would surround vehicles move? a unified framework for maneuver classification and motion prediction,” _IEEE Transactions on Intelligent Vehicles_ , vol. 3, no. 2, pp. 129–140, 2018.
  * [28] S. Zhang, Y. Zhi, R. He, and J. Li, “Research on traffic vehicle behavior prediction method based on game theory and hmm,” _IEEE Access_ , vol. 8, pp. 30 210–30 222, 2020.
  * [29] D. Xu, Z. Ding, X. He, H. Zhao, M. Moze, F. Aioun, and F. Guillemard, “Learning from naturalistic driving data for human-like autonomous highway driving,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 12, pp. 7341–7354, 2020.
  * [30] A. Kuefler, J. Morton, T. Wheeler, and M. Kochenderfer, “Imitating driver behavior with generative adversarial networks,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2017, pp. 204–211.
  * [31] A. Abdelraouf, M. Abdel-Aty, and J. Yuan, “Utilizing attention-based multi-encoder-decoder neural networks for freeway traffic speed prediction,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 8, pp. 11 960–11 969, 2021.
  * [32] L. Xin, P. Wang, C.-Y. Chan, J. Chen, S. E. Li, and B. Cheng, “Intention-aware long horizon trajectory prediction of surrounding vehicles using dual lstm networks,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2018, pp. 1441–1446.
  * [33] S. H. Park, B. Kim, C. M. Kang, C. C. Chung, and J. W. Choi, “Sequence-to-sequence prediction of vehicle trajectory via lstm encoder-decoder architecture,” in _2018 IEEE intelligent vehicles symposium (IV)_. IEEE, 2018, pp. 1672–1678.
  * [34] F. Giuliari, I. Hasan, M. Cristani, and F. Galasso, “Transformer networks for trajectory forecasting,” in _2020 25th international conference on pattern recognition (ICPR)_. IEEE, 2021, pp. 10 335–10 342.
  * [35] N. Deo and M. M. Trivedi, “Convolutional social pooling for vehicle trajectory prediction,” in _Proceedings of the IEEE conference on computer vision and pattern recognition workshops_ , 2018, pp. 1468–1476.
  * [36] R. Chandra, U. Bhattacharya, A. Bera, and D. Manocha, “Traphic: Trajectory prediction in dense and heterogeneous traffic using weighted interactions,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 8483–8492.
  * [37] G. Xie, A. Shangguan, R. Fei, W. Ji, W. Ma, and X. Hei, “Motion trajectory prediction based on a cnn-lstm sequential model,” _Science China Information Sciences_ , vol. 63, pp. 1–21, 2020.
  * [38] W. Zeng, M. Liang, R. Liao, and R. Urtasun, “Lanercnn: Distributed representations for graph-centric motion forecasting,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 532–539.
  * [39] J. Su, P. A. Beling, R. Guo, and K. Han, “Graph convolution networks for probabilistic modeling of driving acceleration,” in _2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2020, pp. 1–8.
  * [40] Y. Xing, C. Lv, and D. Cao, “Personalized vehicle trajectory prediction based on joint time-series modeling for connected vehicles,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 2, pp. 1341–1352, 2019.
  * [41] D. Yi, J. Su, C. Liu, and W.-H. Chen, “Trajectory clustering aided personalized driver intention prediction for intelligent vehicles,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 6, pp. 3693–3702, 2018.
  * [42] Z. Huang, J. Wu, and C. Lv, “Driving behavior modeling using naturalistic human driving data with inverse reinforcement learning,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 8, pp. 10 239–10 251, 2021.
  * [43] X. Liao, X. Zhao, Z. Wang, Z. Zhao, K. Han, R. Gupta, M. J. Barth, and G. Wu, “Driver digital twin for online prediction of personalized lane change behavior,” _IEEE Internet of Things Journal_ , 2023.
  * [44] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _arXiv preprint arXiv:1609.02907_ , 2016.


# Computation of open-loop inputs for uniformly ensemble controllable systems

Michael<EMAIL_ADDRESS>[ University
of Würzburg  Institute for Mathematics
University of Würzburg
Emil-Fischer Straße 40
97074 Würzburg, Germany

E-mail:

###### Abstract

This paper presents computational methods for families of linear systems
depending on a parameter. Such a family is called ensemble controllable if for
any family of parameter-dependent target states and any neighborhood of it
there is a parameter-independent input steering the origin into the
neighborhood. Assuming that a family of systems is ensemble controllable we
present methods to construct suitable open-loop input functions. Our approach
to solve this infinite-dimensional task is based on a combination of methods
from the theory of linear integral equations and finite-dimensional control
theory.

45A05,

93B05,

93B40,

Linear integral equation of first kind,

Moment discretization,

Convergence rates,

Reproducing kernel Hilbert spaces,

Controllability,

Input computation,

Consensus methods,

###### keywords:

[class=AMS]

###### keywords:

This work is accepted for publication in _Mathematical Control and Related
Fields_ and was supported by the German Research Foundation (DFG)
grant SCHO 1780/1-1.
Received 30. September 2019; revised 30. August 2020.

##  1 Introduction

Recently, the topic of ensemble control has become an emerging field in
mathematical systems and control theory. An essential characteristic is the
task of controlling a large, potentially infinite, number of states, or
systems at once, using a single open-loop input function or a single feedback
controller. Although the notion ensemble control is relatively new, such
problems have been studied in the 1970s and 1980s under different names. For
instance, in [34] the blending problem refers to the task of stabilizing a
family of scalar transfer functions by a common dynamic feedback compensator.
Also, in a series of papers Ghosh and co-workers consider the problem of
simultaneous controllability and stabilizability for families of linear
systems, cf. [15, 16, 17, 18].

In principle, ensemble control embraces several facets. A current nice
exposition for the applicability of ensemble control is given in [6, Section
2.4]. For instance, it contains problems where states or observations of a
system are subject to uncertainties and only a probability distribution on the
state space or output snapshots are available. In this case, the ensemble
control problem is to morph a given initial probability density function of
the states into a desired one by transporting it along a linear system. This
yields controllability problems for the Liouville equation and the Fokker-
Plank equation. We refer to, e.g. [6, 7, 10, 13, 14] for recent works on this
topic and note that this list is not intended to be complete. Also, in [39,
41] the problem of reconstructing an initial probability density function from
output snapshots is addressed.

Moreover, ensemble control also deals with families of systems, usually
defined by parameter-dependent systems. Here, the goal is to control the
entire family by finitely many open-loop inputs or feedback controllers. That
is, given families of initial states and desired terminal states the task is
to steer the initial family towards the desired family in some finite time by
a parameter-independent input function. Problems of his kind arise for
instance in robotics and quantum control, cf. [3, 26, 37]. Of course,
characterizations of ensemble controllability depend on the structure of the
parameter-dependent system. For necessary and sufficient conditions for
ensemble controllability for linear systems, bilinear systems and nonlinear
systems we refer to [11, 20, 24, 25, 31, 40] and [4, 8, 9] and [1],
respectively.

In this paper we consider ensembles defined by families of parameter-dependent
linear systems. Assuming that it is ensemble controllable and given a family
of target states, we tackle the question how to derive suitable open-loop
input functions. The rest of the paper is organized as follows. In Section 2
we outline the basic problem studied in this paper and provide the necessary
notations and definitions concerning ensemble controllability. In Section 3 we
extend a known discretization method to obtain approximate solutions to
vector-valued linear integral equations of first kind. To get such a
collocation method, we provide the reproducing kernel structure for Hilbert
spaces of vector-valued functions. Moreover, we derive a convergence rate
under the weaker regularity assumption that the reproducing kernel is
Lipschitz continuous. Note that the standard result by Nashed and Wahba [28,
Theorem 3.1] requires that the reproducing kernel satisfies a
differentiability condition. In Section 4 we present new methods that provide
algorithms for the computation of suitable open-loop input functions.

##  2 Problem statement and computational method

In this paper we consider families of parameter-dependent linear control
systems

$\begin{split}\tfrac{\partial x}{\partial
t}(t,\theta)&=A(\theta)x(t,\theta)+B(\theta)u(t)\\\
x(0,\theta)&=x_{0}(\theta)\in\mathbb{R}^{n},\end{split}$ (1)

where the matrices $A(\theta)\in\mathbb{R}^{n\times n}$ and
$B(\theta)\in\mathbb{R}^{n\times m}$ depend continuously on a parameter
$\theta$ which is varying over a nonempty compact interval
$\mathbf{P}\subset\mathbb{R}$. Throughout the paper we will use the short
notation $(A,B)\in C_{n,n}(\mathbf{P})\times C_{n,m}(\mathbf{P})$ for $A\in
C(\mathbf{P},\mathbb{R}^{n\times n})$ and $B\in
C(\mathbf{P},\mathbb{R}^{n\times m})$. Similarily we shortly write
$L^{q}_{m}([0,T])$ instead of $L^{q}([0,T],\mathbb{R}^{m})$. For $u\in
L^{1}_{m}([0,T])$ let $\varphi\big{(}t,u,x_{0}(\theta),\theta\big{)}$ denote
the solution to (1). To express that (1) defines an infinite-dimensional
system whose states are functions mapping the compact interval $\mathbf{P}$ to
$\mathbb{R}^{n}$ we write
$\varphi(t,u,x_{0})(\theta):=\varphi\big{(}t,u,x_{0}(\theta),\theta\big{)}$.

The precise definition of ensemble controllability is as follows. A pair
$(A,B)\in C_{n,n}(\mathbf{P})\times C_{n,m}(\mathbf{P})$ is called _uniformly
ensemble controllable_ if for every $\varepsilon>0$ and for every $f\in
C_{n}(\mathbf{P})$ there are $T>0$ and $u\in L^{1}_{m}([0,T])$ such that

$\|\varphi(T,u,x_{0})-f\|_{\infty}=\sup_{\theta\in\mathbf{P}}\|\varphi(T,u,x_{0}(\theta),\theta)-f(\theta)\|_{\mathbb{R}^{n}}<\varepsilon.$

As shown in [36, Thm. 3.1.1, Remark 3.1.2] in continuous-time the final time
$T>0$ can be chosen arbitrary. Thus, from now we fix $T>0$. We note that for
discrete-time systems the final time $T>0$ depends on the target state $f$ and
the desired accuracy $\varepsilon>0$.

Most of the literature mentioned in the introduction considers the problem of
deriving convenient necessary and sufficient conditions for uniform ensemble
controllability in terms of $(A,B)$. Almost all results are existence results,
i.e. the existence of an appropriate input function is guaranteed but explicit
expressions are not available in general. Here, we tackle the next step,
namely supposing a pair $(A,B)$ is uniformly ensemble controllable, this paper
considers the problem of determining suitable input functions.

In order to make use of Hilbert space properties we are looking for input
functions in $L^{2}_{m}([0,T])$. If the pair $(A,B)\in
C_{n,n}(\mathbf{P})\times C_{n,m}(\mathbf{P})$ is uniformly ensemble
controllable, the range of the compact reachability integral operator

$\mathcal{R}\colon L^{2}_{m}([0,T])\to
L^{2}_{n}(\mathbf{P}),\quad\mathcal{R}u(\theta)=\int_{0}^{T}e^{A(\theta)(T-s)}B(\theta)\,u(s)\operatorname{d}\\!s$
(2)

is dense in $C_{n}(\mathbf{P})$, cf. [11]. We consider how to determine an
appropriate input function $u\in L^{2}_{m}([0,T])$ such that for given $f\in
C_{n}(\mathbf{P})$ and $\varepsilon>0$ one has

$\|\mathcal{R}u-f\|_{C_{n}(\mathbf{P})}<\varepsilon.$ (3)

In other words, the problem of computing an appropriate open-loop input
function is equivalent to finding an approximate solution to the linear
operator equation

$\mathcal{R}u=f.$ (4)

This is a classical ill-posed inverse problem. Well-known necessary and
sufficient conditions such that (4) has a solution is given by Picard’s
criterion, cf. [12, 22, 23]. It also provides an explicit representation of
the solution. This, however, has the drawback that it depends on the singular
system of the compact reachability operator $\mathcal{R}$. The results in [25]
are based on Picard’s criterion. Other related approaches can be found in [35,
42].

In this paper we follow a different approach, namely a discretization method.
Based on the assumption that $(A,B)\in C_{n,n}(\mathbf{P})\times
C_{n,m}(\mathbf{P})$ is uniformly ensemble controllable, we have

$\displaystyle\overline{\operatorname{Im}\mathcal{R}}=C_{n}(\mathbf{P})\qquad\text{
and }\qquad(\operatorname{Im}\mathcal{R})^{\bot}=\\{0\\}.$

The uniform ensemble reachability assumption yields that for every $f\in
C_{n}(\mathbf{P})$ and every $\varepsilon>0$ there is a function
$g\in\operatorname{Im}\operatorname{\mathcal{R}}$ such that
$\|f-g\|_{\infty}<\varepsilon$. Thus, without loss of generality we treat the
case that $f$ is an element of the range of $\mathcal{R}$. Despite the
regularization techniques and projection methods available in the literature
above, we follow a different approach called moment collocation, which
inherits well-known techniques from control theory, e.g. the controllability
Gramian. To obtain computational methods we will incorporate ideas originating
from decentralized methods to solve linear finite-dimensional equations, as
presented, e.g. in [2, 27, 32]. Finally, we note that our subsequent methods
do not depend on the singular system of the reachability operator
$\mathcal{R}$, which is used in [25].

##  3 Moment collocation

In this section we present the moment collocation method to derive approximate
solutions to the linear equation (4). The basic idea is simply to take
$N\in\mathbb{N}$ different parameters out of the compact parameter interval
$\mathbf{P}=[\theta^{-},\theta^{+}]$, i.e. we pick

$\theta^{-}=\theta_{1}<\theta_{2}<\cdots<\theta_{N}=\theta^{+}.$

The $\theta_{i}$’s are sometimes called _moments_. For the pair $(A,B)\in
C_{n,n}(\mathbf{P})\times C_{n,m}(\mathbf{P})$ evaluated at the moments we use
the short notation

$(A_{k},B_{k}):=(A(\theta_{k}),B(\theta_{k})),\quad k=1,...,N.$

Then, we define the interpolation operator $\mathcal{I}_{N}\colon
L^{2}_{m}([0,T])\to\mathbb{R}^{Nn}$,

$\mathcal{I}_{N}u=\begin{pmatrix}\int_{0}^{T}e^{A_{1}(T-s)}B_{1}\,u(s)\operatorname{d}\\!s\\\
\vdots\\\
\int_{0}^{T}e^{A_{N}(T-s)}B_{N}\,u(s)\operatorname{d}\\!s\end{pmatrix}.$

A justification for the term interpolation operator can be found e.g. in [23,
Section 11.2]. As in the proof of [11, Theorem 3 (c)], the application of the
restriction property to $\\{\theta_{1},...,\theta_{N}\\}$ yieds that the
parallel connection

$\left(\begin{pmatrix}A_{1}&&\\\ &\ddots&\\\
&&A_{N}\end{pmatrix},\begin{pmatrix}B_{1}\\\ \vdots\\\
B_{N}\end{pmatrix}\right)$ (5)

of the finite-dimensional pairs $(A_{k},B_{k})$, $k=1,...,N$ is controllable.
Thus, for each $f\in C_{n}(\mathbf{P})$ and

$F_{N}=\begin{pmatrix}f(\theta_{1})\\\ \vdots\\\ f(\theta_{N})\end{pmatrix}$

the linear equation

$\mathcal{I}_{N}u=F_{N}$ (6)

has a unique minimal $L^{2}$-norm solution, denoted by $u_{\parallel N}$. It
is well-known, cf. [33, Section 3.5, Theorem 5] that the solution can be
written as

$u_{\parallel
N}(s)=\begin{pmatrix}B_{1}^{\mathsf{T}}&\cdots&B_{N}^{\mathsf{T}}\end{pmatrix}\begin{pmatrix}e^{A_{1}^{\mathsf{T}}(T-s)}&&\\\
&\ddots&\\\
&&{}^{A_{N}^{\mathsf{T}}(T-s)}\end{pmatrix}W^{-1}_{\parallel}F_{N},$

where $W_{\parallel}$ denotes the controllability Gramian of the parallel
connetion (5) given by

$\begin{split}\int_{0}^{T}\begin{pmatrix}e^{A_{1}(T-s)}B_{1}&\cdots&e^{A_{N}(T-s)}B_{N}\end{pmatrix}\begin{pmatrix}B_{1}^{\mathsf{T}}e^{A_{1}^{\mathsf{T}}(T-s)}&\cdots&B_{N}^{\mathsf{T}}e^{A_{N}^{\mathsf{T}}(T-s)}\end{pmatrix}\operatorname{d}\\!s\\\
=\left(\int_{0}^{T}e^{A_{l}(T-s)}B_{l}B_{k}^{\mathsf{T}}e^{A_{k}^{\mathsf{T}}(T-s)}\operatorname{d}\\!s\right)_{l,k=1,...,N}.\end{split}$
(7)

In the following we are interested in, whether for given $f\in
C_{n}(\mathbf{P})$ and $\varepsilon>0$, the input $u_{\parallel N}$ can be
used to obtain a suitable input satisfying (3) by choosing suitable moment
parameters $\\{\theta_{1},...,\theta_{N}\\}$. To this end, we define

$\Delta_{N}:=\sup_{\theta\in\mathbf{P}}\left(\inf_{i=1,...,N}|\theta_{i}-\theta|\right).$

Note, if the moments are chosen equidistantly, i.e. for $\delta>0$ and
$N\in\mathbb{N}$ the moments

$\theta^{-}=\theta_{1}\leq\cdots<\theta_{N}=\theta_{+}\quad\text{ are such
that }\quad|\theta_{i+1}-\theta_{i}|=\delta,$

it holds

$\Delta_{N}\leq\frac{\delta}{2}.$

In order to carry out the convergence analysis we need some preparations. In
doing so, we introduce Hilbert spaces with reproducing kernels. Let
$X_{n}(\mathbf{P})$ denote a Hilbert space of functions from $\mathbf{P}$ to
$\mathbb{R}^{n}$. The following definition is taken from [30] and is in line
with the definition in [28] which treats the scalar case.

######  Definition 3.1 (Reproducing kernel).

A reproducing kernel for the Hilbert space $X_{n}(\mathbf{P})$ is a family
$Q(\theta,\tilde{\theta})$ of transformations from $\mathbb{R}^{n}$ to
$\mathbb{R}^{n}$ defined for each $\theta,\tilde{\theta}\in\mathbf{P}$ and
having the properties:

  1. 1.

for every $\tilde{\theta}\in\mathbf{P}$ and $v\in\mathbb{R}^{n}$ it holds
$Q(\cdot,\tilde{\theta})v\in X_{n}(\mathbf{P})$.

  2. 2.

for every $\tilde{\theta}\in\mathbf{P},f\in X_{n}(\mathbf{P})$ and
$v\in\mathbb{R}^{n}$ it holds

$\langle v,f(\tilde{\theta})\rangle_{\mathbb{R}^{n}}=\langle
Q(\cdot,\tilde{\theta})v,f\rangle_{X_{n}(\mathbf{P})}.$

To get a convergence result, we equip the range $\operatorname{Im}\mathcal{R}$
with an inner product. Let
$\mathcal{R}^{\sharp}\colon\operatorname{Im}\mathcal{R}\to L^{2}_{m}([0,T])$
denote the Moore-Penrose inverse of $\mathcal{R}$, that is for
$f\in\operatorname{Im}\mathcal{R}$ the element $\mathcal{R}^{\sharp}f\in
L^{2}_{m}([0,T])$ is the unique solution of $\mathcal{R}u=f$ with minimal
$L^{2}_{m}([0,T])$-norm, cf. [12]. Then, an inner product on
$\operatorname{Im}\mathcal{R}$ is given by

$\langle
f,g\rangle_{\operatorname{Im}\mathcal{R}}:=\langle\mathcal{R}^{\sharp}f,\mathcal{R}^{\sharp}g\rangle_{L^{2}_{m}([0,T])}=\int_{0}^{T}\langle\mathcal{R}^{\sharp}f(s),\mathcal{R}^{\sharp}g(s)\rangle_{\mathbb{R}^{m}}\operatorname{d}\\!s\qquad
f,g\in\operatorname{Im}\mathcal{R}.$

######  Lemma 3.1.

$(\operatorname{Im}\mathcal{R},\langle\cdot,\cdot\rangle_{\operatorname{Im}\mathcal{R}})$
is a reproducing kernel Hilbert space.

Proof. That $\operatorname{Im}\mathcal{R}$ defines a Hilbert space is shown in
[28, Prop. 2.1]. To show the second claim, we verify that

$Q(\theta,\tilde{\theta})=\int_{0}^{T}e^{A(\theta)(T-s)}B(\theta)B(\tilde{\theta})^{\mathsf{T}}\,e^{A(\tilde{\theta})^{\mathsf{T}}(T-s)}\operatorname{d}\\!s$

defines a reproducing kernel.

To conclude the first property, let $v\in\mathbb{R}^{n}$ and
$\tilde{\theta}\in\mathbf{P}$. Then, by definition we have

$Q(\theta,\tilde{\theta})v=\int_{0}^{T}e^{A(\theta)(T-s)}B(\theta)B(\tilde{\theta})^{\mathsf{T}}\,e^{A(\tilde{\theta})^{\mathsf{T}}(T-s)}v\operatorname{d}\\!s.$

Since $s\mapsto
B(\tilde{\theta})^{\mathsf{T}}\,e^{A(\tilde{\theta})^{\mathsf{T}}(T-s)}v$ is
an $L^{2}_{m}([0,T])$-function, the first property follows.

To conclude the second property take $f\in\operatorname{Im}\mathcal{R}$ and
let $u_{f}:=\mathcal{R}^{\sharp}f\in L^{2}_{m}([0,T])$. Then, for every
$v\in\mathbb{R}^{n}$ and $\tilde{\theta}\in\mathbf{P}$ one has on one hand

$\langle
v,f(\tilde{\theta})\rangle_{\mathbb{R}^{n}}=\int_{0}^{T}v^{\mathsf{T}}\,e^{A(\tilde{\theta})(T-s)}B(\tilde{\theta})u_{f}(s)\operatorname{d}\\!s$

and on the other hand

$\begin{split}\langle
Q(\cdot,\tilde{\theta})v,f\rangle_{\operatorname{Im}\mathcal{R}}&=\langle\mathcal{R}^{\sharp}Q(\cdot,\tilde{\theta})v,\mathcal{R}^{\sharp}f\rangle_{L^{2}_{m}([0,T])}\\\
&=\int_{0}^{T}\langle
B(\tilde{\theta})^{\mathsf{T}}\,e^{A(\tilde{\theta})^{\mathsf{T}}(T-s)}v,u_{f}(s)\rangle_{\mathbb{R}^{n}}\operatorname{d}\\!s\\\
&=\int_{0}^{T}v^{\mathsf{T}}\,e^{A(\tilde{\theta})(T-s)}B(\tilde{\theta})u_{f}(s)\operatorname{d}\\!s.\end{split}$

The assertion follows by comparing the right hand sides. $\Box$

Next, we present the essential structure given by the reproducing kernel. Let
$\operatorname{e}_{i}$ denote the $i$th standard basis vector in
$\mathbb{R}^{n}$ and

$\displaystyle
U:=\overline{\operatorname{span}\\{u_{i,\theta}=B(\theta)^{\mathsf{T}}e^{A(\theta)^{\mathsf{T}}(T-\cdot)}\operatorname{e}_{i}\,\,|\,\,i=1,...,n\,,\theta\in\mathbf{P}\\}}\subset
L_{m}^{2}([0,T])$

and note that

$\displaystyle\operatorname{Im}\mathcal{R}=\operatorname{span}\\{Q(\cdot,\tilde{\theta})\operatorname{e}_{i}\,|\,i=1,...,n,\,\tilde{\theta}\in\mathbf{P}\\}.$

Then, by the reproducing kernel structure we get

$\begin{split}\langle\mathcal{R}u_{i,\tilde{\theta}},\mathcal{R}u_{j,\hat{\theta}}\rangle_{\operatorname{Im}\mathcal{R}}&=\langle
Q(\cdot,\tilde{\theta})\operatorname{e}_{i},Q(\cdot,\hat{\theta})\operatorname{e}_{j}\rangle_{\operatorname{Im}\mathcal{R}}=\langle\operatorname{e}_{i},Q(\tilde{\theta},\hat{\theta})\operatorname{e}_{j}\rangle_{\mathbb{R}^{n}}\\\
&=\int_{0}^{T}\langle
u_{i,\tilde{\theta}}(s),u_{j,\hat{\theta}}(s)\rangle_{\mathbb{R}^{n}}\operatorname{d}\\!s=\langle
u_{i,\tilde{\theta}},u_{j,\hat{\theta}}\rangle_{L_{m}^{2}([0,T])}.\end{split}$

That is, $\mathcal{R}\colon U\to\operatorname{Im}\mathcal{R}$ is an isometric
isomorphism and using $\mathcal{R}\mathcal{R}^{\sharp}f=f$ for all
$f\in\operatorname{Im}\mathcal{R}$, we have

$\begin{split}\|\mathcal{R}u-f\|_{\operatorname{Im}\mathcal{R}}^{2}&=\langle\mathcal{R}u,\mathcal{R}u\rangle_{\operatorname{Im}\mathcal{R}}-2\langle\mathcal{R}u,\mathcal{R}\mathcal{R}^{\sharp}f\rangle_{\operatorname{Im}\mathcal{R}}+\langle
f,f\rangle_{\operatorname{Im}\mathcal{R}}\\\ &=\langle
u,u\rangle_{L^{2}_{m}([0,T])}-2\langle u,\
\mathcal{R}^{\sharp}f\rangle_{L^{2}_{m}([0,T])}+\langle\mathcal{R}^{\sharp}f,\mathcal{R}^{\sharp}f\rangle_{L^{2}_{m}([0,T])}\\\
&=\|u-\mathcal{R}^{\sharp}f\|_{L^{2}_{m}([0,T])}^{2}\end{split}$ (8)

for all $u\in U$.

The significance here is that the convergence analysis of the moment
collocation method in the Hilbert space $L^{2}_{m}([0,T])$ can be traced back
to a convergence analysis in
$(\operatorname{Im}\mathcal{R},\langle\cdot,\cdot\rangle_{\operatorname{Im}\mathcal{R}})$.
Next, we construct the solution $u_{\parallel N}$ by following a classical
approach in the theory of integral equations, cf. [23, 28]. To do so, we take
finite-dimensional subspaces of $\operatorname{Im}\mathcal{R}$ by choosing for
$i\in\\{1,...,n\\}$ and $k\in\\{1,...,N\\}$ the $i$th column of the matrix
$B_{k}^{\mathsf{T}}e^{A^{\mathsf{T}}_{k}(T-s)}$ as an element of
$L^{2}_{m}([0,T])$. That is, we define

$\displaystyle v_{i,k}(\cdot):=Q(\cdot,\theta_{k})\operatorname{e}_{i},\quad
u_{i,k}(\cdot):=B(\theta_{k})^{\mathsf{T}}e^{A(\theta_{k})^{\mathsf{T}}(T-\cdot)}\operatorname{e}_{i}$

and

$\displaystyle V_{N}$
$\displaystyle:=\operatorname{span}\\{v_{1,k},...,v_{n,k}\,\,|\,\,k=1,...,N\\}\subset\operatorname{Im}\mathcal{R}$
$\displaystyle U_{N}$
$\displaystyle:=\operatorname{span}\\{u_{1,k},...,u_{n,k}\,\,|\,\,k=1,...,N\\}\subset
L_{m}^{2}([0,T]).$

We note explicitly that, by definition it holds

$\displaystyle\mathcal{R}u_{i,k}=v_{i,k}\qquad\text{ and }\qquad
u_{i,k}=\mathcal{R}^{\sharp}v_{i,k}$ (9)

for all $i=1,...,n$ and $k=1,...,N$. Without loss of generality, we assume
that the functions $v_{i,k}$ are linearly independent. Consequently, $U_{N}$
and $V_{N}$ are $nN$-dimensional subspaces of $L^{2}_{m}([0,T])$ and
$\operatorname{Im}\mathcal{R}$, respectively.

Now, we return to the solution $u_{\parallel N}$ of (6) and shall show that it
lies in $U_{N}$. To see this, let $u\in U_{N}$, i.e. there are real numbers
$\alpha_{1,1},...,\alpha_{n,N}$ such that

$u=\sum_{k=1}^{N}\sum_{i=1}^{n}\alpha_{i,k}u_{i,k}.$ (10)

Let $\mathcal{R}_{k}\colon L^{2}_{m}([0,T])\to\mathbb{R}^{n}$,

$\mathcal{R}_{k}u=\int_{0}^{T}e^{A_{k}(T-s)}B_{k}\,u(s)\operatorname{d}\\!s.$

Then, the linear equation (6) can also be written as follows

$\begin{pmatrix}\mathcal{R}_{1}u_{1,1}&\cdots&\mathcal{R}_{1}u_{n,N}\\\
\vdots&&\vdots\\\
\mathcal{R}_{N}u_{1,1}&\cdots&\mathcal{R}_{N}u_{n,N}\end{pmatrix}\begin{pmatrix}\alpha_{1,1}\\\
\vdots\\\ \alpha_{n,N}\end{pmatrix}=\begin{pmatrix}f(\theta_{1})\\\ \vdots\\\
f(\theta_{N})\end{pmatrix}.$

In terms of the reproducing kernel $Q$, the matrix

$\begin{pmatrix}\mathcal{R}_{1}u_{1,1}&\cdots&\mathcal{R}_{1}u_{n,N}\\\
\vdots&&\vdots\\\
\mathcal{R}_{N}u_{1,1}&\cdots&\mathcal{R}_{N}u_{n,N}\end{pmatrix}$

can also be written in block form

$Q:=\begin{pmatrix}Q(\theta_{1},\theta_{1})&\cdots&Q(\theta_{1},\theta_{N})\\\
\vdots&&\vdots\\\
Q(\theta_{N},\theta_{1})&\cdots&Q(\theta_{N},\theta_{N})\end{pmatrix}.$ (11)

Recall that,

$Q(\theta_{l},\theta_{k})=\int_{0}^{T}e^{A_{l}(T-s)}B_{l}B_{k}^{\mathsf{T}}\,e^{A_{k}^{\mathsf{T}}(T-s)}\operatorname{d}\\!s$

and thus $Q$ is just the controllability Gramian (7), i.e. $Q$ is positiv
definite. Hence, (6) is equally written as $Q\alpha=F_{N}$ and the real
numbers $\alpha_{1,1},...,\alpha_{n,N}$ are given by

$\begin{pmatrix}\alpha_{1,1}\\\ \vdots\\\
\alpha_{n,N}\end{pmatrix}=Q^{-1}\begin{pmatrix}f(\theta_{1})\\\ \vdots\\\
f(\theta_{N})\end{pmatrix}.$

Plugging this into (10) we get that

$u_{\parallel N}=\sum_{k=1}^{N}\sum_{i=1}^{n}\alpha_{i,k}\,u_{i,k}$

is another explicit representation of the minimal $L^{2}_{m}([0,T])$-norm
solution to (6). The next result states conditions under which the solution
$u_{\parallel N}$ also approximates the minimal $L^{2}$-norm solution to
$\mathcal{R}u=f$. Before we state and prove this, we need some notation. For a
matrix $M=(M_{ij})\in\mathbb{R}^{n\times n}$ we will use the following matrix
norm:

$\|M\|_{n,\infty}:=n\max_{1\leq i,j\leq n}|M_{ij}|$

and recall that it is consistent and submultiplicative (cf. [21, Chapter 5]),
i.e. it holds

$\|Mx\|\leq\|M\|_{n,\infty}\|x\|\quad\text{ and
}\quad\|M\,N\|_{n,\infty}=\|M\|_{n,\infty}\,\|N\|_{n,\infty}$

for all $M,N\in\mathbb{R}^{n\times n}$ and for all $x\in\mathbb{R}^{n}$.
Moreover, we say that the reproducing kernel $Q$ is Lipschitz continuous on
$\mathbf{P}\times\mathbf{P}$ if there is an $L_{Q}>0$ such that

$\|Q(\eta_{1},\tau_{1})-Q(\eta_{2},\tau_{2})\|_{n,\infty}\leq
L_{Q}\left(\,|\eta_{1}-\eta_{2}|+|\tau_{1}-\tau_{2}|\,\right)$

for all
$(\eta_{1},\tau_{1}),(\eta_{2},\tau_{2})\in\mathbf{P}\times\mathbf{P}$.
Moreover, let

$\Delta_{\max}:=\max_{k=1,...,N-1}|\theta_{k+1}-\theta_{k}|.$

Then, we have the following refinement of the convergence result obtained in
[28, Theorem 3.1] and [38].

######  Theorem 3.1.

Let $(A,B)\in C_{n,n}(\mathbf{P})\times C_{n,m}(\mathbf{P})$ and
$f\in\operatorname{Im}\mathcal{R}$.

  1. (a)

Then, it holds

$\displaystyle\lim_{\Delta_{N}\to 0}\|u_{\parallel
N}-\mathcal{R}^{\sharp}f\|_{L^{2}([0,T],\mathbb{R}^{m})}=0.$

  2. (b)

If $Q$ is Lipschitz contiunuous on $\mathbf{P}\times\mathbf{P}$ and there is
an $g\in L^{2}_{n}(\mathbf{P})$ such that
$\mathcal{R}^{\sharp}f=\mathcal{R}^{*}g$ (i.e.
$f\in\mathcal{R}\mathcal{R}^{*}(L^{2}_{n}(\mathbf{P})$) then, we have

$\displaystyle\|u_{\parallel
N}-\mathcal{R}^{\sharp}f\|_{L^{2}_{m}([0,T])}=\sqrt{\Delta_{\max}\,3\,n\left(4\,M_{Q}e^{-\tfrac{2}{\Delta_{\max}}}+L_{Q}\,\sqrt{\Delta_{\max}}\right)}\,\|g\|_{L_{n}^{2}(\mathbf{P})},$

where
$M_{Q}=\max_{i,j=1,...,n}\max_{(\tau,\eta)\in\mathbf{P}\times\mathbf{P}}|Q(\tau,\eta)_{ij}|$.

###### Proof.

(a): In a first step we show that $u_{\parallel N}$ satisfies

$\mathcal{R}u_{\parallel N}=P_{V_{N}}f,$ (12)

where $P_{V_{N}}$ denotes the orthogonal projection in
$\operatorname{Im}\mathcal{R}$ onto $V_{N}$. To see this, we shall show that

$\langle\mathcal{R}u_{\parallel
N}-f,v\rangle_{\operatorname{Im}\mathcal{R}}=0\qquad\forall\,\,v\in V_{N}.$

Recall that, by construction it holds $u_{\parallel
N}=I_{N}^{\sharp}F_{N}=P_{U_{N}}\mathcal{R}^{\sharp}f$, i.e.

$\langle u_{\parallel
N}-\mathcal{R}^{\sharp}f,u\rangle_{L^{2}_{m}([0,T])}=0\qquad\forall\,\,u\in
U_{N}.$

Let $v\in V_{N}$. Since $\mathcal{R}\colon U_{N}\to V_{N}$ is isometric, there
is an $u\in U_{N}$ so that $\mathcal{R}u=v$ and since
$\mathcal{R}\mathcal{R}^{\sharp}f=f$ we have

$\displaystyle\langle\mathcal{R}u_{\parallel
N}-f,v\rangle_{\operatorname{Im}\mathcal{R}}$
$\displaystyle=\langle\mathcal{R}u_{\parallel
N}-f,\mathcal{R}u\rangle_{\operatorname{Im}\mathcal{R}}=\langle\mathcal{R}(u_{\parallel
N}-\mathcal{R}^{\sharp}f),\mathcal{R}u\rangle_{\operatorname{Im}\mathcal{R}}$
$\displaystyle=\langle u_{\parallel
N}-\mathcal{R}^{\sharp}f,u\rangle_{L^{2}_{m}([0,T])}=0.$

Let

$\displaystyle
V:=\operatorname{span}\left\\{v_{1,\tilde{\theta}},...,v_{n,\tilde{\theta}}\,\,|\,\,\tilde{\theta}\in\bigcup_{N\in\mathbb{N}}\\{\theta_{1},...,\theta_{N}\\}\right\\}.$

Then, as $Q$ is continuous and $\lim_{N\to\infty}\Delta_{N}=0$, it follows
that $V$ is dense in $\operatorname{Im}\mathcal{R}$, cf. [23, Proof of Theorem
17.13]. Thus, using (8) and (12) we have

$\displaystyle\lim_{\Delta_{N}\to 0}\|u_{\parallel
N}-\mathcal{R}^{\sharp}f\|_{L^{2}_{m}([0,T])}$
$\displaystyle=\lim_{\Delta_{N}\to 0}\|\mathcal{R}u_{\parallel
N}-f\|_{\operatorname{Im}\mathcal{R}}$ $\displaystyle=\lim_{\Delta_{N}\to
0}\|\mathcal{P}_{V_{N}}f-f\|_{\operatorname{Im}\mathcal{R}}=0.$

(b): The assumption $f\in\mathcal{R}\mathcal{R}^{*}(L^{2}_{n}(\mathbf{P}))$
implies that there is a function $g=\begin{pmatrix}g_{1}\\\ \vdots\\\
g_{N}\end{pmatrix}\in L^{2}_{n}(\mathbf{P})$ such that

$f(\theta)=\int_{\mathbf{P}}Q(\theta,\eta)g(\eta)\operatorname{d}\\!\eta.$
(13)

Since for any $\tilde{f}\in V_{N}$ it holds

$\displaystyle\|f-\mathcal{P}_{V_{N}}f\|_{\operatorname{Im}\mathcal{R}}\leq\|f-\tilde{f}\|_{\operatorname{Im}\mathcal{R}},$

the verification of (b) is simply to find appropriate functions
$p_{k}\colon\mathbf{P}\to\mathbb{R}$, $k=1,...,N$ such that the function

$\displaystyle\tilde{f}(\cdot):=\sum_{k=1}^{N}Q(\cdot,\theta_{k})c_{k},\qquad
c_{k}=\begin{pmatrix}\int_{\mathbf{P}}p_{k}(\tau)g_{1}(\tau)\operatorname{d}\\!\tau\\\
\vdots\\\
\int_{\mathbf{P}}p_{k}(\tau)g_{n}(\tau)\operatorname{d}\\!\tau\end{pmatrix}$

is in $V_{N}$ and $\|f-\tilde{f}\|_{\operatorname{Im}\mathcal{R}}$ can be
bounded from above.

It holds,

$\displaystyle\|f-\tilde{f}\|^{2}_{\operatorname{Im}\mathcal{R}}=\langle
f,f\rangle_{\operatorname{Im}\mathcal{R}}-2\langle
f,\tilde{f}\rangle_{\operatorname{Im}\mathcal{R}}+\langle\tilde{f},\tilde{f}\rangle_{\operatorname{Im}\mathcal{R}}.$

Next, we consider the terms individually. By assumption it holds
$\mathcal{R}^{\sharp}f=\mathcal{R}^{*}g$ and one has

$\displaystyle\langle f,f\rangle_{\operatorname{Im}\mathcal{R}}$
$\displaystyle=\langle\mathcal{R}^{\sharp}f,\mathcal{R}^{\sharp}f\rangle_{L^{2}_{m}([0,T])}=\langle\mathcal{R}^{*}g,\mathcal{R}^{*}g\rangle_{L^{2}_{m}([0,T])}$
$\displaystyle=\int_{0}^{T}\left(\int_{\mathbf{P}}B(\tau)^{\mathsf{T}}e^{A(\tau)^{\mathsf{T}}(T-s)}g(\eta)\operatorname{d}\\!\tau\right)^{\mathsf{T}}\left(\int_{\mathbf{P}}B(\eta)^{\mathsf{T}}e^{A(\eta)^{\mathsf{T}}(T-s)}g(\eta)\operatorname{d}\\!\eta\right)\operatorname{d}\\!s$
$\displaystyle=\int_{\mathbf{P}}\int_{\mathbf{P}}g(\tau)^{\mathsf{T}}\left(\int_{0}^{T}e^{A(\tau)(T-s)}B(\tau)B(\eta)^{\mathsf{T}}e^{A(\eta)^{\mathsf{T}}(T-s)}\operatorname{d}\\!s\right)g(\eta)\operatorname{d}\\!\tau\operatorname{d}\\!\eta$
$\displaystyle=\int_{\mathbf{P}}\int_{\mathbf{P}}g(\tau)^{\mathsf{T}}Q(\tau,\eta)g(\eta)\operatorname{d}\\!\tau\operatorname{d}\\!\eta.$

Furthermore, exploiting the reproducing kernel structure and (13) we get

$\displaystyle\langle f,\tilde{f}\rangle_{\operatorname{Im}\mathcal{R}}$
$\displaystyle=\sum_{k=1}^{N}\langle
f,Q(\cdot,\theta_{k})c_{k}\rangle_{\operatorname{Im}\mathcal{R}}=\sum_{k=1}^{N}\langle
f(\theta_{k}),c_{k}\rangle_{\mathbb{R}^{n}}=\sum_{k=1}^{N}c_{k}^{\mathsf{T}}f(\theta_{k})$
$\displaystyle=\sum_{k=1}^{N}\int_{\mathbf{P}}\int_{\mathbf{P}}p_{k}(\tau)g(\tau)^{\mathsf{T}}Q(\theta_{k},\eta)g(\eta)\operatorname{d}\\!\tau\operatorname{d}\\!\eta$

and

$\displaystyle\langle\tilde{f},\tilde{f}\rangle_{\operatorname{Im}\mathcal{R}}$
$\displaystyle=\sum_{k=1}^{N}\sum_{l=1}^{N}\left\langle
Q(\cdot,\theta_{k})c_{k},Q(\cdot,\theta_{l})c_{l}\right\rangle_{\operatorname{Im}\mathcal{R}}=\sum_{k=1}^{N}\sum_{l=1}^{N}\langle
c_{k},Q(\theta_{k},\theta_{l})c_{l}\rangle_{\mathbb{R}^{n}}$
$\displaystyle=\int_{\mathbf{P}}\int_{\mathbf{P}}\sum_{k=1}^{N}\sum_{l=1}^{N}g(\tau)^{\mathsf{T}}p_{k}(\tau)Q(\theta_{k},\theta_{l})p_{l}(\eta)g(\eta)\operatorname{d}\\!\tau\operatorname{d}\\!\eta.$

This yields

$\|f-\tilde{f}\|^{2}_{\operatorname{Im}\mathcal{R}}=\\\
\int_{\mathbf{P}}\int_{\mathbf{P}}g(\tau)^{\mathsf{T}}\left\\{Q(\tau,\eta)-\sum_{k=1}^{N}\Big{(}2\,p_{k}(\tau)Q(\theta_{k},\eta)-\sum_{l=1}^{N}p_{k}(\tau)Q(\theta_{k},\theta_{l})p_{l}(\eta)\Big{)}\right\\}g(\eta)\operatorname{d}\\!\tau\operatorname{d}\\!\eta.$

In the next step we construct linear interpolations of the matrices
$Q(\eta,\tau)$. To this end, we partition the interval $\mathbf{P}$ according
to the moments

$\displaystyle\mathbf{P}=[\theta_{1},\theta_{2}]\cup[\theta_{2},\theta_{3}]\cup\cdots\cup[\theta_{N-1},\theta_{N}]=:\mathbf{P}_{1}\cup\cdots\cup\mathbf{P}_{N-1}.$

Then, on each subinterval $\mathbf{P}_{k}$ we define the polynomials $l_{k,1}$
and $l_{k,2}$ as follows

$\displaystyle\begin{split}l_{k,1}(\theta)=\frac{\theta-\theta_{k+1}}{\theta_{k}-\theta_{k+1}}\qquad\text{
and }\qquad
l_{k,2}(\theta)=\frac{\theta-\theta_{k}}{\theta_{k+1}-\theta_{k}}.\end{split}$
(14)

Thus, on $\mathbf{P}_{k}$ and for fixed $\eta\in\mathbf{P}$ the linear
interpolation of the matrix-valued function $\tau\mapsto Q(\tau,\eta)$ is
given by the matrix-valued Lagrange polynomial

$\displaystyle\tau\mapsto
B_{1,Q}(\tau,\eta):=l_{k,1}(\tau)\,Q(\theta_{k},\eta)+l_{k,2}(\tau)\,Q(\theta_{k+1},\eta).$

The notation $B_{1,Q}$ is intended to express that the interpolation
polynomial is equal to the first Bernstein polynomial. Hence, for fixed
$\eta>0$ and $\tau\in\mathbf{P}_{k}$, as shown in the proof of Theorem 1 in
[19, Eq. (4)] we get for each entry $Q(\tau,\eta)_{ij}$ the following estimate

$\displaystyle|Q(\tau,\eta)_{ij}-B_{1,Q}(\tau,\eta)_{ij}|\leq
4\,e^{-\tfrac{2}{\theta_{k+1}-\theta_{k}}}\,\max_{\tau\in\mathbf{P}_{k}}|Q(\tau,\eta)_{ij}|\,+L_{Q}\sqrt{\theta_{k+1}-\theta_{k}},$

where $L_{Q}>0$ denotes the Lipschitz constant of $Q$. In order to apply the
previous estimation, we define for $k=1,...,N$ the functions

$\displaystyle\begin{split}p_{k}(\theta)=\begin{cases}l_{k,1}(\theta)&\quad\theta\in\mathbf{P}_{k}\\\
0&\quad\theta\not\in\mathbf{P}_{k}\end{cases}\quad\text{ and }\quad
p_{k+1}(\theta)=\begin{cases}l_{k,2}(\theta)&\quad\theta\in\mathbf{P}_{k}\\\
0&\quad\theta\not\in\mathbf{P}_{k}.\end{cases}\end{split}$ (15)

Then, for $\tau\in\mathbf{P}_{k}$ and $\eta\in\mathbf{P}_{l}$ it holds

$\displaystyle Q(\tau,\eta)$
$\displaystyle-\sum_{k=1}^{N}\left(2\,p_{k}(\tau)Q(\theta_{k},\eta)-\sum_{l=1}^{N}p_{k}(\tau)Q(\theta_{k},\theta_{l})p_{l}(\eta)\right)=$
$\displaystyle Q(\tau,\eta)$ $\displaystyle-
l_{k,1}(\tau)\,Q(\theta_{k},\eta)-l_{k,2}(\tau)\,Q(\theta_{k+1},\eta)$
$\displaystyle-
l_{k,1}(\tau)\Big{(}Q(\theta_{k},\eta)-l_{l,1}(\eta)Q(\theta_{k},\theta_{l})-l_{l,2}(\eta)Q(\theta_{k},\theta_{l+1})\Big{)}$
$\displaystyle-
l_{k,2}(\tau)\Big{(}Q(\theta_{k+1},\eta)-l_{l,1}(\eta)Q(\theta_{k+1},\theta_{l})-l_{l,2}(\eta)Q(\theta_{k+1},\theta_{l+1})\Big{)}.$

Defining

$M_{Q}=\max_{i,j=1,...,n}\max_{(\tau,\eta)\in\mathbf{P}\times\mathbf{P}}|Q(\tau,\eta)_{ij}|$

one has for each entry the estimate

$\displaystyle\left|Q(\tau,\eta)_{ij}-\sum_{k=1}^{N}\left(2\,p_{k}(\tau)Q(\theta_{k},\eta)_{ij}-\sum_{l=1}^{N}p_{k}(\tau)Q(\theta_{k},\theta_{l})_{ij}p_{l}(\eta)\right)\right|$
$\displaystyle\leq|Q(\tau,\eta)_{ij}-l_{k,1}(\tau)\,Q(\theta_{k},\eta)_{ij}-l_{k,2}(\tau)\,Q(\theta_{k+1},\eta)_{ij}|$
$\displaystyle\quad+|l_{k,1}(\tau)|\,|Q(\theta_{k},\eta)_{ij}-l_{l,1}(\eta)Q(\theta_{k},\theta_{l})_{ij}-l_{l,2}(\eta)Q(\theta_{k},\theta_{l+1})_{ij}|$
$\displaystyle\quad+|l_{k,2}(\tau)|\,|Q(\theta_{k+1},\eta)_{ij}-l_{l,1}(\eta)Q(\theta_{k+1},\theta_{l})_{ij}-l_{l,2}(\eta)Q(\theta_{k+1},\theta_{l+1})_{ij}|$
$\displaystyle\leq
12\,e^{-\tfrac{2}{\Delta_{\max}}}\,M_{Q}\,+3\,L_{Q}\,\sqrt{\Delta_{\max}}.$

That is, using
$C=3\,n\left(4\,M_{Q}e^{-\tfrac{2}{\Delta_{\max}}}+L_{Q}\,\sqrt{\Delta_{\max}}\right)$
for a moment, we have

$\displaystyle\Bigg{\|}Q(\tau,\eta)$
$\displaystyle-\sum_{k=1}^{N}\left(2\,p_{k}(\tau)Q(\theta_{k},\eta)-\sum_{l=1}^{N}p_{k}(\tau)Q(\theta_{k},\theta_{l})p_{l}(\eta)\right)\Bigg{\|}_{n,\infty}\leq
C.$

Then, using the Cauchy-Schwarz inequality we get

$\displaystyle\|f-\tilde{f}\|^{2}_{\operatorname{Im}\mathcal{R}}$
$\displaystyle\leq\sum_{\nu=1}^{N}\int_{\mathbf{P}_{\nu}}\int_{\mathbf{P}_{\nu}}\|g(\tau)\|_{\mathbb{R}^{n}}\,C\,\|g(\eta)\|_{\mathbb{R}^{n}}\operatorname{d}\\!\tau\operatorname{d}\\!\eta$
$\displaystyle\leq\sum_{\nu=1}^{N}\,C\left(\int_{\mathbf{P}_{\nu}}\,\|g(\tau)\|_{\mathbb{R}^{n}}^{2}\,\operatorname{d}\\!\tau\right)^{\frac{1}{2}}\left(\int_{\mathbf{P}_{\nu}}\,1\,\operatorname{d}\\!\eta\right)\left(\int_{\mathbf{P}_{\nu}}\,\|g(\eta)\|_{\mathbb{R}^{n}}^{2}\,\operatorname{d}\\!\eta\right)^{\frac{1}{2}}$
$\displaystyle\leq
C\,\Delta_{\max}\,\sum_{\nu=1}^{N}\left(\int_{\mathbf{P}_{\nu}}\|g(\tau)\|_{\mathbb{R}^{n}}^{2}\operatorname{d}\\!\tau\right)\leq
C\,\Delta_{\max}\|g\|_{L_{n}^{2}(\mathbf{P})}^{2}$

and finally

$\displaystyle\|f-\tilde{f}\|^{2}_{\operatorname{Im}\mathcal{R}}\leq\Delta_{\max}\,3\,n\left(4\,M_{Q}e^{-\tfrac{2}{\Delta_{\max}}}+L_{Q}\,\sqrt{\Delta_{\max}}\right)\,\|g\|^{2}_{L_{n}^{2}(\mathbf{P})}.$

This shows the assertion. ∎

Note that, since we neither used the particular kernel
$K(\theta,s)=e^{A(\theta)(T-s)}B(\theta)$ of the reachability operator nor
that the operator is compact, Theorem 3.1 is valid for general integral
equations of first kind with a continuous kernel $K(\theta,s)$ and a Lipschitz
continuous reproducing kernel

$Q(\theta,\tilde{\theta})=\int_{0}^{T}K(\theta,s)K(\tilde{\theta},s)^{\mathsf{T}}\operatorname{d}\\!s.$

We will now formulate a different version of the result, which is more
convenient for the purpose of this paper. Since for
$f\in\operatorname{Im}\mathcal{R}$ it holds
$\mathcal{R}\mathcal{R}^{\sharp}f=f$ one has

$\|\mathcal{R}u_{N}-f\|_{\infty}\leq\|\mathcal{R}\|\,\|u_{N}-\mathcal{R}^{\sharp}f\|_{L^{2}_{m}([0,T])},$

where $\|\mathcal{R}\|$ denotes the operator norm of the reachability operator
$\mathcal{R}$ defined in (2). Then, for given $\varepsilon>0$ and appropriate
$f$ Theorem 3.1 yields the subsequent lower bound on the number of moment
parameters that are required such that $u_{\parallel N}$ is a suitable input.

######  Theorem 3.2.

Suppose that $(A,B)\in C^{1}_{n,n}(\mathbf{P})\times C^{1}_{n,m}(\mathbf{P})$
is uniformly ensemble controllable. Let $\varepsilon>0$ and
$f\in\operatorname{Im}\mathcal{R}$ such that there is a function $g\in
L_{n}^{2}(\mathbf{P})$ with $\mathcal{R}^{\sharp}f=\mathcal{R}^{*}g$. Then,
the input $u_{\parallel N}$ defined in (3) satisfies

$\|\mathcal{R}u_{\parallel N}-f\|_{\infty}\leq\varepsilon$

if the moments $\theta_{1}<\cdots<\theta_{N}$ are chosen such that

$\Delta_{\max}<\tfrac{\varepsilon^{2}}{\|\mathcal{R}\|^{2}\,\|g\|_{L_{n}^{2}(\mathbf{P})}^{2}\,3\,n\,\big{(}4\,M_{Q}e^{-\tfrac{2}{\theta_{N}-\theta_{1}}}+L_{Q}\sqrt{|\theta_{N}-\theta_{1}|}\big{)}}.$

##  4 Consensus methods

In this section, we present additional methods to obtain solutions to the
linear equation

$\mathcal{I}_{N}u=F_{N}.$

Recall that uniform ensemble controllability implies that each finite-
dimensional moment system $(A_{i},B_{i}),i=1,...,N$ is controllable, cf. [20].
Thus, for each $f_{i}:=f(\theta_{i})$ there is an open-loop input $u_{i}\in
L^{2}_{m}([0,T])$ steering the origin to $f_{i}$. Using the controllability
Gramian

$W_{i}=\int_{0}^{T}e^{A_{i}\,(T-s)}B_{i}B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}\,(T-s)}\,\operatorname{d}\\!s$

the $L^{2}$-norm minimal input function $u_{i}$ is given by

$u_{i}=B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}\,(T-\cdot)}\,W_{i}^{-1}f_{i}.$
(16)

We will now use the moment input functions $u_{i}$, $i=1,...,N$ to derive a
global input function which is an approximate solution to (4). We note that
the subsequent methods do not require the inversion of the matrix $Q$ defined
in (11).

###  4.1 Projected gradient method

In order to define our method, we consider for $k=1,...,N$ the operators
$\mathcal{R}_{k}\colon L^{2}_{m}([0,T])\to\mathbb{R}^{n}$,

$\mathcal{R}_{k}u=\int_{0}^{T}e^{A_{k}(T-s)}B_{k}\,u(s)\operatorname{d}\\!s$

and their adjoints $\mathcal{R}_{k}^{*}\colon\mathbb{R}^{n}\to
L^{2}_{m}([0,T])$,

$\mathcal{R}_{k}^{*}v=B_{k}^{\mathsf{T}}e^{A^{\mathsf{T}}_{k}(T-s)}v.$

For $i=1,...,N$ we consider the system of differential equations

$\displaystyle\dot{u}_{i}(t)=\frac{1}{N}\sum_{k=1}^{N}\big{(}u_{k}(t)-u_{i}(t))-\frac{1}{N}\sum_{k=1}^{N}\mathcal{R}_{k}^{*}\mathcal{R}_{k}u_{k}(t)-\mathcal{R}_{k}^{*}f_{k}$
(17)

defined on the Hilbert space $L^{2}_{m}([0,T])$. As initial values $u_{i}(0)$
we pick the known moment control input functions $u_{i}$ defined in (16). The
first part of the right-hand side is a consensus term and the second part is a
gradient term so that $u_{i}(t)$ converges to a solution of
$\mathcal{R}_{i}u=f_{i}$. Before stating and proving the convergence result,
we first recall that a sequence $(u_{n})_{n\in\mathbb{N}}$ in
$L^{2}_{m}([0,T])$ is said to converge weakly to $u$ if for all $v\in
L^{2}_{m}([0,T])$ it holds

$\lim_{n\to\infty}\langle u_{n},v\rangle_{L^{2}_{m}([0,T])}=\langle
u,v\rangle_{L^{2}_{m}([0,T])}.$

We shortly write $u_{n}\rightharpoonup u$ if $(u_{n})_{n\in\mathbb{N}}$
converges weakly to $u$. The convergence result is as follows.

######  Theorem 4.1.

For each $i=1,...,N$ the solution of the differential equation (17) exists for
all $t\geq 0$ and converges weakly to a global solution, i.e. for all
$i=1,...,N$ one has

$u_{i}(t)\rightharpoonup u\qquad\text{ as }\quad t\to\infty$

for some $u\in L^{2}_{m}([0,T])$. In particular, we have

$\mathcal{I}_{N}u=\begin{pmatrix}\mathcal{R}_{1}u\\\ \vdots\\\
\mathcal{R}_{N}u\end{pmatrix}=F_{N}.$

Proof. First we define the consensus subspace

$C=\\{u=(u_{1},...,u_{N})\,|\,u_{1}=\cdots=u_{N}\\},$

which is nonempty, closed and convex. Furthermore, let $V\colon
L^{2}_{m}([0,T])^{N}\to[0,\infty)$,

$V(u)=\sum_{k=1}^{N}\frac{1}{2}\|\mathcal{R}_{k}u_{k}-f_{k}\|_{\mathbb{R}^{n}}^{2}.$

The functional $V$ is Frechét differentiable with

$\nabla
V(u)=\begin{pmatrix}\mathcal{R}_{1}^{*}\mathcal{R}_{1}u_{1}-\mathcal{R}_{1}^{*}f_{1}\\\
\vdots\\\
\mathcal{R}_{N}^{*}\mathcal{R}_{N}u_{N}-\mathcal{R}_{N}^{*}f_{N}\end{pmatrix}.$

Next, we note that the orthogonal projection $P_{C}\colon
L^{2}_{m}([0,T])^{N}\to C$ is given by

$P_{C}\,u=\begin{pmatrix}\tfrac{1}{N}(u_{1}+\cdots+u_{N})\\\ \vdots\\\
\tfrac{1}{N}(u_{1}+\cdots+u_{N})\end{pmatrix}.$

Using this notation and stacking the equations (17) for $i=1,...,N$ we get

$\dot{u}(t)=-u(t)+P_{C}\left(u(t)-\nabla V\big{(}u(t)\big{)}\right).$ (18)

The assertion then follows from [5, Theorem 3.2]. $\Box$

In a next step we modify the system of differential equations (17) to get
strong convergence of the corresponding trajectories. To this end, let
$\eta:[0,\infty)\to[0,\infty)$ be a nonincreasing continuously differentiable
function with bounded derivative satisfying

$\lim_{t\to\infty}\eta(t)=0,\qquad\int_{0}^{\infty}\eta(t)\operatorname{d}\\!t=\infty.$
(19)

Then, for $i=1,...,N$ we consider the modified ordinary differential equations

$\displaystyle\dot{u}_{i}(t)=\frac{1}{N}\sum_{k=1}^{N}\Big{(}\big{(}1-\eta(t)\big{)}u_{k}(t)-u_{i}(t)\Big{)}-\frac{1}{N}\sum_{k=1}^{N}\mathcal{R}_{k}^{*}\mathcal{R}_{k}u_{k}(t)-\mathcal{R}_{k}^{*}f_{k}.$
(20)

As initial conditions we pick one of the moment input functions given in (16),
i.e

$u_{i}(0)=u_{0}\in\left\\{u_{1}(s)=B_{1}^{\mathsf{T}}e^{A_{1}^{\mathsf{T}}\,(T-s)}\,W_{1}^{-1}f_{1},...,u_{N}(s)=B_{N}^{\mathsf{T}}e^{A_{N}^{\mathsf{T}}\,(T-s)}\,W_{N}^{-1}f_{N}\right\\}.$

Collocating for $i=1,...,N$, the equations (20) can be written as

$\dot{u}(t)=-u(t)+P_{C}\left(u(t)-\nabla
V\big{(}u(t)\big{)}-\eta(t)\,u(t)\right).$ (21)

Using [5, Theorem 5.1] we conclude the following strong convergence result.

######  Theorem 4.2.

For each $i=1,...,N$ the solution of the differential equation (20) exists for
all $t\geq 0$ and converges strongly to a global solution, i.e. for all
$i=1,...,N$ one has

$\lim_{t\to\infty}u_{i}(t)=u,$

for some $u\in L^{2}_{m}([0,T])$. In particular, we have
$\mathcal{I}_{N}u=F_{N}$.

###  4.2 Averaging Method

Another method to solve the infinite-dimensional linear integral equation

$\mathcal{I}_{N}u=F_{N}$

using the known moment input functions $u_{i}\in M_{i}$ is inspired by [32],
where the authors present distributed algorithms for solving finite-
dimensional linear equations. Let

$M_{i}:=\\{u\in L^{2}_{m}([0,T])\,\,|\,\mathcal{R}_{i}u=f_{i}\,\\}$

denote the set of all solutions to

$\mathcal{R}_{i}u=f_{i}.$ (22)

Following the ideas for distributed solving of linear finite-dimensional
equations, cf. [27, 32], we consider the orthogonal projection
$P_{M_{i}}\colon L^{2}_{m}([0,T])\to M_{i}$. In order to compute $P_{M_{i}}$,
we note that $M_{i}$ is an affine subspace and so $P_{M_{i}}$ is given by

$\displaystyle g\mapsto P_{M_{i}}g$
$\displaystyle=u_{i}+P_{\operatorname{ker}\mathcal{R}_{i}}(g-u_{i})$
$\displaystyle=u_{i}+(I-P_{(\operatorname{ker}\mathcal{R}_{i})^{\perp}})(g-u_{i}),$

with $u_{i}\in M_{i}$ defined in (16). Since the range of $\mathcal{R}_{i}$ is
finite-dimensional, it is also closed. Then, by standard functional analysis
one has

$(\operatorname{ker}\mathcal{R}_{i})^{\perp}=\operatorname{Im}\mathcal{R}_{i}^{*}.$

Consequently, it remains to determine the orthogonal projection onto
$\operatorname{Im}\mathcal{R}_{i}^{*}$. This is derived as follows. Let $h\in
L_{m}^{2}([0,T])$ and $y=P_{\operatorname{Im}\mathcal{R}_{i}^{*}}h$. Then,
there is a $v_{y}\in\mathbb{R}^{n}$ such that

$y(s)=B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}(T-s)}v_{y},\quad\forall\,s\in[0,T].$

Furthermore, it holds

$\langle
h-y,g\rangle_{L^{2}_{m}([0,T])}=0\qquad\forall\,g\in\operatorname{Im}\mathcal{R}_{i}^{*},$

or equivalently

$\langle
h-B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}(T-\cdot)}v_{y},B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}(T-\cdot)}v\rangle_{L^{2}_{m}([0,T])}=0\qquad\forall\,v\in\mathbb{R}^{n}.$

It follows that

$v_{y}=W_{i}^{-1}\int_{0}^{T}e^{A_{i}(T-s)}B_{i}h(s)\operatorname{d}\\!s$

and therefore, we obtain

$\displaystyle P_{(\operatorname{ker}\mathcal{R}_{i})^{\perp}}h(s)$
$\displaystyle=B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}(T-s)}W_{i}^{-1}\int_{0}^{T}e^{A_{i}(T-\tau)}B_{i}h(\tau)\operatorname{d}\\!\tau$
$\displaystyle=\mathcal{R}_{i}^{*}(W_{i}^{-1}\mathcal{R}_{i}h).$

As a result, the orthogonal projection onto $M_{i}$ is given by

$P_{M_{i}}g=g-\mathcal{R}_{i}^{*}(W_{i}^{-1}\mathcal{R}_{i}g)+\mathcal{R}^{*}_{i}(W_{i}^{-1}f_{i}).$
(23)

A simple calculation shows that the orthogonal projection onto $M_{i}$ is
affine linear, i.e.

$P_{M_{i}}(g+h)=P_{M_{i}}g+P_{M_{i}}h-P_{M_{i}}(0)$

and the following identity holds

$P_{M_{i}}\left(\frac{g_{1}+\cdots+g_{N}}{N}\right)=\frac{P_{M_{i}}g_{1}+\cdots+P_{M_{i}}g_{N}}{N}.$

Now, we are in the position to present our method. In order to obtain a global
solution we define the projection-consensus system of differential equations

$\dot{u}_{i}(t)=\sum_{k=1}^{N}\left(P_{M_{i}}u_{k}(t)-u_{i}(t)\right).$ (24)

As initial conditions we chose again the known moment solutions
$u_{i}(0)=u_{i}\in M_{i}$ defined in (16). The following result shows that the
mean-value of the moment solutions converges and its limit is a global
solution. More precisely, we have the following

######  Theorem 4.3.

The solutions to (24) satisfy

$\lim_{t\to\infty}\frac{1}{N}\sum_{i=1}^{N}u_{i}(t)=u,$

and $u$ is a solution to $\mathcal{I}_{N}u=F_{N}$.

Proof. We begin with denoting the mean-value by $z(t)$, i.e.

$z(t):=\frac{1}{N}\sum_{i=1}^{N}u_{i}(t).$

The evolution of the mean-value is given by the following ordinary
differential equation

$\displaystyle\dot{z}(t)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\dot{u}_{i}(t)=\frac{1}{N}\sum_{i=1}^{N}\left(\sum_{k=1}^{N}P_{M_{i}}u_{k}(t)-u_{i}(t)\right)$
$\displaystyle=\sum_{i=1}^{N}\left(\frac{1}{N}\sum_{k=1}^{N}P_{M_{i}}u_{k}(t)-\frac{1}{N}\sum_{l=1}^{N}u_{l}(t)\right)$
$\displaystyle=\sum_{i=1}^{N}\left(P_{M_{i}}\left(\frac{1}{N}\sum_{k=1}^{N}u_{k}(t)\right)-\frac{1}{N}\sum_{l=1}^{N}u_{l}(t)\right)$
$\displaystyle=\sum_{i=1}^{N}\left(P_{M_{i}}z(t)-z(t)\right).$

Spelling out the orthogonal projection (23) the latter differential equation
also takes the form

$\dot{z}(t)=-\sum_{i=1}^{N}\mathcal{R}_{i}^{*}\left(W_{i}^{-1}\mathcal{R}_{i}z(t)\right)+\sum_{i=1}^{N}\mathcal{R}_{i}^{*}\left(W_{i}^{-1}f_{i}\right).$
(25)

Since $W_{i}^{-1}$ is symmetric and positiv definite, by [21, Theorem 7.2.6]
there is a symmetric and positiv definite matrix $S_{i}$ such that $S_{i}\cdot
S_{i}=W_{i}^{-1}$. Further, we consider the modified interpolation operator
$\widetilde{\mathcal{I}}_{N}\colon L^{2}_{m}([0,T])\to\mathbb{R}^{nN}$,

$\widetilde{\mathcal{I}}_{N}u=\begin{pmatrix}S_{1}\int_{0}^{T}e^{A_{1}(T-s)}B_{1}\,u(s)\operatorname{d}\\!s\\\
\vdots\\\
S_{N}\int_{0}^{T}e^{A_{N}(T-s)}B_{N}\,u(s)\operatorname{d}\\!s\end{pmatrix}$

and its adjoint $\widetilde{\mathcal{I}}_{N}^{*}\colon\mathbb{R}^{nN}\to
L^{2}_{m}([0,T])$,

$v=\begin{pmatrix}v_{1}\\\ \vdots\\\
v_{N}\end{pmatrix}\mapsto\widetilde{\mathcal{I}}_{N}^{*}v(s)=\sum_{i=1}^{N}B_{i}^{\mathsf{T}}e^{A_{i}^{\mathsf{T}}(T-s)}S_{i}v_{i}.$

Using
$\widetilde{F}_{n}:=\begin{pmatrix}S_{1}f_{1}&\cdots&S_{N}f_{N}\end{pmatrix}^{\mathsf{T}}$,
the linear differential equation (25) also reads as

$\dot{z}(t)=-\widetilde{\mathcal{I}}_{N}^{*}\widetilde{\mathcal{I}}_{N}z(t)+\widetilde{\mathcal{I}}_{N}^{*}\widetilde{F}_{N}.$

We note that for $V\colon L^{2}_{m}([0,T])\to[0,\infty)$,

$V(z):=\frac{1}{2}\|\widetilde{\mathcal{I}}_{N}z-\widetilde{F}_{N}\|_{\mathbb{R}^{nN}}^{2}$

one has

$\dot{z}(t)=-\nabla
V(z)=-\widetilde{\mathcal{I}}_{N}^{*}\widetilde{\mathcal{I}}_{N}z+\widetilde{\mathcal{I}}_{N}^{*}\widetilde{F}_{N}.$

Thus, applying [29, Theorem 3.1] we conclude that $\lim_{t\to\infty}z(t)=u$
exists and satisfies $\widetilde{\mathcal{I}}_{N}u=\widetilde{F}_{N}$, or
equivalently

$\int_{0}^{T}e^{A_{i}(T-s)}B_{i}\,u(s)\operatorname{d}\\!s=f_{i}\qquad\text{
for all }i=1,...,N.$

This shows the assertion. $\Box$

## References

  * [1] A. Agrachev, Y. Baryshnikov and A. Sarychev,  Ensemble controllability by Lie algebraic methods,  _ESAIM Control Optim. Calc. Var._ , 22 (2016), 921–938.
  * [2] B. Anderson, S. Mou, A. S. Morse and U. Helmke,  Decentralized gradient algorithm for solution of a linear equation,  _Numerical Algebra, Control and Optimization_ , 6 (2017), 319–328.
  * [3] A. Becker and T. Bretl,  Approximate steering of a unicycle under bounded model perturbation using ensemble control,  _IEEE Transactions on Robotics_ , 28 (2012), 580–591.
  * [4] M. Belhadj, J. Salomon and G. Turinici,  Ensemble controllability and discrimination of perturbed bilinear control systems on connected, simple, compact Lie groups,  _Eur. J. Control_ , 22 (2015), 23–29.
  * [5] J. Bolte,  Continuous gradient projection method in Hilbert spaces,  _J. Optim. Theory Appl._ , 119 (2003), 235–259.
  * [6] R. Brockett,  Notes on the control of the Liouville equation,  in _Control of Partial Differential Equations_ (eds. F. Alabau-Boussouira, R. Brockett, O. Glass, J. LeRousseau and E. Zuazua), vol. 2048 of Lecture Notes in Mathematics,  Springer, Heidelberg, 2012,  101–129.
  * [7] Y. Chen, T. T. Georgiou and M. Pavon,  Optimal transport over a linear dynamical system,  _IEEE Transactions on Automatic Control_ , 62 (2017), 2137–2152.
  * [8] F. C. Chittaro and J.-P. Gauthier,  Asymptotic ensemble stabilizability of the Bloch equation,  _Systems Control Lett._ , 113 (2018), 36–44.
  * [9] G. Dirr,  Ensemble controllability of bilinear systems,  _Oberwolfach Reports_ , vol. 9, no. 1, pp. 674–676, 2012.
  * [10] G. Dirr, U. Helmke and M. Schönlein,  Controlling mean and variance in ensembles of linear systems,  _IFAC-PapersOnLine_ , 49 (2016), 1018–1023.
  * [11] G. Dirr and M. Schönlein,  Uniform and $L^{q}$-ensemble ensemble reachability of parameter-dependent linear systems,  _Journal of Differential Equations_ , accepted for publication.
  * [12] H. W. Engl, M. Hanke and A. Neubauer,  _Regularization of Inverse Problems_ ,  Kluwer Academic Publishers Group, Dordrecht, 1996.
  * [13] A. Fleig and L. Grüne,  Estimates on the minimal stabilizing horizon length in model predictive control for the Fokker-Planck equation,  _IFAC-PapersOnLine_ , 49 (2016), 260–265.
  * [14] P. A. Fuhrmann and U. Helmke,  _The Mathematics of Networks of Linear Systems_ ,  Springer International Publishing, Switzerland, 2015.
  * [15] B. K. Ghosh,  Some new results on the simultaneous stabilizability of a family of single input, single output systems,  _Systems & Control Letters_, 6 (1985), 39–45,
  * [16] B. K. Ghosh,  An approach to simultaneous system design. I: Semialgebraic geometric methods.,  _SIAM J. Control Optim._ , 24 (1986), 480–496.
  * [17] B. K. Ghosh,  An approach to simultaneous system design. II: Nonswitching gain and dynamic feedback compensation by algebraic geometric methods,  _SIAM J. Control Optim._ , 26 (1988), 919–963.
  * [18] B. Ghosh,  Transcendental and interpolation methods in simultaneous stabilization and simultaneous partial pole placement problems,  _SIAM J. Control Optim._ , 24 (1986), 1091–1109.
  * [19] H. Gzyl and J. L. Palacios,  The Weierstrass approximation theorem and large deviations,  _Am. Math. Mon._ , 104 (1997), 650–653.
  * [20] U. Helmke and M. Schönlein,  Uniform ensemble controllability for one-parameter families of time-invariant linear systems,  _Systems & Control Letters_, 71 (2014), 69–77.
  * [21] R. A. Horn and C. R. Johnson,  _Matrix Analysis_ ,  Cambridge University Press, 1990.
  * [22] A. Kirsch,  _An Introduction to the Mathematical Theory of Inverse Problems_ ,  Springer Science & Business Media, 2011.
  * [23] R. Kress,  _Linear Integral Equations_ ,  3rd edition,  Springer, New York, 2014,
  * [24] J. S. Li and J. Qi,  Ensemble control of time-invariant linear systems with linear parameter variation,  _IEEE Transactions on Automatic Control_ , 61 (2016), 2808–2820.
  * [25] J.-S. Li,  Ensemble control of finite-dimensional time-varying linear systems,  _IEEE Transactions on Automatic Control_ , 56 (2011), 345–357.
  * [26] J.-S. Li and N. Khaneja,  Control of inhomogeneous quantum ensembles,  _Physical review A_ , 73 (2006), 030302.
  * [27] S. Mou and A. Morse,  A fixed-neighbor, distributed algorithm for solving a linear algebraic equation,  in _Proc. IEEE European Control Conference (ECC), 2013_ ,  pp. 2269–2273.
  * [28] M. Z. Nashed and G. Wahba,  Convergence rates of approximate least squares solutions of linear integral and operator equations of the first kind,  _Math. Comp._ , 28 (1974), 69–80.
  * [29] J. W. Neuberger,  _Sobolev Gradients and Differential Equations_ ,  2nd edition,  Springer-Verlag, Berlin, 2010.
  * [30] G. Pedrick,  _Theory of Reproducing Kernels for Hilbert Spaces of Vector Valued Functions_ ,  PhD. Thesis, University of Kansas, 1958.
  * [31] M. Schönlein and U. Helmke,  Controllability of ensembles of linear dynamical systems,  _Mathematics and Computers in Simulation_ , 125 (2016), 3–14.
  * [32] G. Shi, B. D. Anderson and U. Helmke,  Network flows that solve linear equations,  _IEEE Transactions on Automatic Control_ , 62 (2017), 2659–2674.
  * [33] E.D. Sontag,  _Mathematical Control Theory. Deterministic Finite Dimensional Systems_ ,  2nd Ed.,  Springer, New York, 1998.
  * [34] A. Tannenbaum,  _Invariance and System Theory: Algebraic and Geometric Aspects._ ,  Lecture Notes in Mathematics. 845. Berlin-Heidelberg-New York: Springer-Verlag, 1981.
  * [35] L. Tie, W. Zhang, S. Zeng, and J.-S. Li, “Explicit input signal design for stable linear ensemble systems,” _IFAC-PapersOnLine_ , 50 (2017), 3051 – 3056.
  * [36] R. Triggiani,  Controllability and observability in Banach space with bounded operators,  _SIAM Journal on Control_ , 13 (1975), 462–491.
  * [37] G. Turinici, V. Ramakhrishna, B. Li and H. Rabitz,  Optimal discrimination of multiple quantum systems: Controllability analysis,  _Journal of Physics A: Mathematical and General_ , 37 (2004), 273–282.
  * [38] G. Wahba,  Convergence rates of certain approximate solutions to Fredholm integral equations of the first kind,  _J. Approximation Theory_ , 7 (1973), 167–185.
  * [39] S. Zeng, H. Ishii and F. Allgöwer,  Sampled observability and state estimation of discrete ensembles,  _IEEE Trans. Autom. Contr._ , 62 (2017), 2406–2418.
  * [40] S. Zeng and F. Allgöwer,  A moment-based approach to ensemble controllability of linear systems,  _Systems & Control Letters_, 98 (2016), 49 – 56.
  * [41] S. Zeng, S. Waldherr, C. Ebenbauer and F. Allgöwer,  Ensemble observability of linear systems,  _IEEE Transactions on Automatic Control_ , 61 (2016), 1452–1465.
  * [42] S. Zeng, W. Zhang, and J.-S. Li,  On the computation of control inputs for linear ensembles,  Proc. _2018 Annual American Control Conference (ACC)_ , (2018), pp. 6101–6107.