1
COIN-GP: Cooperative Online Learning in Networked Distributed Systems with Partial Measurements via Gaussian Process Regression
arXiv:2609.20598v1 [cs.LG] 17 Sep 2026
Zewen Yang1 , Member, IEEE, Xiaobing Dai*2 , Zhenxiao Yin3 , Hang Zhao3 , Zhijun Li4 , Fellow, IEEE, C.C. Chan5 , Life Fellow, IEEE
Abstract—In this paper, we tackle the problem of jointly estimating the system states and partially unknown dynamics within distributed sensor-equipped networks, particularly in scenarios where only partial state observations are available. To address this issue, we propose an observer-based dynamic cooperative learning framework incorporating online distributed Gaussian Process (GP) regression, which enables accurate estimation despite incomplete in measurements and deficient GP models. In addition, a novel data collection strategy is introduced, with theoretical conditions ensuring feasible data acquisition. Moreover, we also derive an error upper bound encompassing state estimation and model estimation, leveraging the deterministic error bounds of GPs. Empirical simulations demonstrate the superiority of our approach compared to existing distributed GP-based methods. Index Terms—Cooperative learning, online learning, networked systems, Gaussian process regression, system estimation.
I. I NTRODUCTION OOPERATIVE learning in distributed dynamical systems has gained significant attention due to its broad applications in areas such as sensor networks, multi-robot systems, and autonomous underwater vehicles [1]–[3]. The primary goal of cooperative learning is to enable a group of nodes or agents to collaboratively infer the unknown pattern by sharing information and leveraging neighbors’ learned results. In this study, we investigate the problem of handling partial measurements in distributed learning systems while ensuring performance guarantees. The major challenge lies in the joint estimation of the system states and the unknown function. Although state estimation techniques such as the Luenberger observer, high-gain observer, and sliding mode observer have been widely employed in feedback systems, they face significant limitations for discrete-time systems with unknown nonlinear dynamics components. In particular, the Luenberger observer, despite its efficacy in reconstructing states from partial measurements, relies on a fully known
C
*Corresponding author: Xiaobing Dai <[email protected]> 1 Zewen Yang is with the Chair of Robotics and Systems Intelligence (RSI), Munich Institute of Robotics and Machine Intelligence (MIRMI), Technical University of Munich (TUM), 80992 Munich, Germany. 2 Xiaobing Dai is with the Technical University of Munich, 80333 Munich, Germany. 3 Zhenxiao Yin and Hang Zhao are with The Hong Kong University of Science and Technology (Guangzhou). 4 Zhijun Li is with the School of Mechanical Engineering, Translational Research Center, Tongji University, Shanghai 201804, China, affiliated with Shanghai Yangzhi Rehabilitation Hospital and also with Shanghai Key Laboratory of Wearable Robotics and Human-Machine Interaction, and also with Department of Automation, University of Science and Technology of China, Hefei 230026, China. 5 C.C. Chan is with The Hong Kong Polytechnic University.
system model, rendering it impractical when the nonlinear component is unknown, even if the linear part of the dynamics is known [4]. Similarly, high-gain observers, primarily designed for continuous-time systems, encounter stability issues and exacerbated peaking phenomena when discretized, due to finite sampling rates [5]. Sliding mode observers, while robust in continuous-time applications, struggle with the implementation of infinitely fast switching in discrete-time systems, leading to chattering and reduced accuracy [6]. To address these challenges, we employ Gaussian process (GP) regression, a machine learning technique particularly well-suited for unknown function modeling in dynamical systems [7]. As a Bayesian learning method, GP regression is widely recognized for its rigorous prediction error quantification [8] and is frequently used in model identification. Its capability for online learning offers a distinct advantage, enabling real-time model updates based on new measurements and thereby facilitating adaptation to changing environments. In addition to the problems posed by the partial measurements and unknown nonlinear dynamics components, we consider a distributed sensor network comprising multiple subsystems. These subsystems are equipped with sensors that measure only partial states of the overall system, adding another layer of complexity to the estimation process. Furthermore, we consider the scenario where only a limited number of subsystems have access to sufficient training data for learning. Our proposed approach for cooperative online learning in networked distributed systems under partial measurements via GP regression (COIN-GP) enables the exchange of predictions among interconnected subsystems, which not only enhances the efficiency of the learning process but also improves the overall state estimation performance, particularly for subsystems with limited or poor datasets. Through dynamic cooperative learning, subsystems can leverage information from their neighbors, significantly enhancing prediction accuracy and system robustness. A. Related Work In addressing the joint estimation problem within nonlinear systems, [9] introduced an adaptive state observer capable of both state estimation and model identification. Building upon this framework, [10] further enhanced it by integrating GP models. These observers rely on high-gain mechanisms, which present limitations when applied to discretetime nonlinear systems as mentioned above. Nevertheless, the existing methodologies primarily focus on single systems, thereby presenting challenges when extending their application to distributed systems, such as multi-agent systems (MASs)
2
or multi-robot systems, due to the complexities inherent in interconnections between subsystems. In MASs, distributed observers are essential for estimating unknown dynamics or states. Adaptive distributed observers have been proposed for leader-following consensus problems [11]–[13], including scenarios with uncertain leaders in a linear form [14] and reduced-order design [15]. For formation tracking, distributed extended state observers have been used to model the nonlinear unknown functions of agents [16], but these proposed methods are limited to second-order strictfeedback dynamics. Even with significant progress in the development of observers for MASs, current methods primarily address either state estimation or unknown function modeling separately, lacking strategies for joint estimation problems in MASs. To fill this gap, recent work has explored the use of neural networks for learning unknown functions in MASs [17]. For instance, learning-based distributed observers have been proposed to estimate both the states and parameters of the leader node [18], [19]. However, these approaches often rely on the assumption that the unknown function is in a vector form of a linear combination of known features. GP-based methods provide deterministic error bounds and naturally handle agents with limited or no training data through uncertaintyweighted cooperation, whereas NN-based approaches can offer greater scalability and capacity at the cost of interpretability and formal guarantees. Therefore, we consider GP regression for handling unknown dynamics that do not presuppose a predefined model structure. In order to adapt to the distributed systems, distributed Gaussian processes (DGPs), where a central entity aggregates predictions from distributed experts using divided sub-datasets [20]. Advanced methods for DGPs in MASs include optimizing aggregation weights through variations in the product of experts (PoE) method [21] and applying the dynamical average consensus algorithm [22] or maximal information coefficient [23] to expedite convergence in joint predictions. However, these approaches encounter challenges when applied to scenarios involving online learning. Specifically, the offline collected dataset, which does not update dynamically during the operation with new observations, limits their adaptability in environments that are dynamic. Recent research [24]–[26] underscores the critical importance of integrating streaming data to enhance the adaptability of DGPs. However, these methodologies require complete state observation in order to work. While current approaches are wellsuited for continuous-time systems, their direct applicability to discrete-time systems is constrained due to different stability conditions. More recently, cooperative GP-based learning has been advanced by elective learning [27], [28], control under switching topologies [29], and event-triggered online learning [26], while distributed impulsive observers handle aperiodically sampled outputs [30], interconnected systems with input and output disturbances [31], and uncertaintyaware learning is strengthened in [32], [33]. Nevertheless, the learning-based approaches assume full state measurements, whereas the estimation-oriented schemes do not simultaneously learn the unknown dynamics. Therefore, to the best of our knowledge, there exists no approach designed explicitly to
resolve the problems associated with unknown dynamics and incomplete measurements in online distributed learning with GPs in discrete-time systems while also allowing subsystems without inference capabilities. B. Contribution and Structure The summary of the contributions of the paper is listed below. • Data acquisition strategy: For constructing the training dataset with limited measurements, the data collection strategy is proposed along with sufficient conditions of data collectability. • Distributed learning-based estimation: The COIN-GP method is presented not only to facilitate full state estimation with incomplete observations but also to enable each individual system to infer the unknown function. Additionally, the proposed dynamic cooperative learning approach allows the distributed system to operate without the need to establish a GP model for every subsystem. • Learning and control performance analysis: The joint error upper bound of the state observation and model prediction is derived by employing the deterministic error bound inherited from GP regression, which ensures guaranteed performance for safety-critical applications. The paper is organized as follows: In Section II, we lay out the preliminaries and define the problem setting. Section III then elaborates on the data collection strategy and the COIN-GP approach, along with a comprehensive performance analysis. In Section IV, we demonstrate the effectiveness of our proposed method by conducting comparative analyses with existing distributed GP-based aggregation frameworks through simulation experiments. Finally, Section V presents the conclusions drawn from this study. II. P RELIMINARIES AND P ROBLEM F ORMULATION A. Notation and Graph Theory We denote R by the set of real numbers and let R>0 represents positive reals (0, ∞) and R≥0 non-negative reals [0, ∞). Natural numbers are symbolized by N and N0 with zero. I[a,b] denotes the integer interval from a to b. The vector 1N ∈ RN is defined as 1N = [1, · · · , 1]T and the N × N identity matrix as IN . Define the function λi (·) returns the i-th eigenvalue of a square matrix, while λ(·) returns the minimum eigenvalue, and λ̄(·) returns the maximum eigenvalue of the matrix. Additionally, we represent a diagonal matrix as diag(a1 , . . . , aN ) and concatenated block matrices diagonally as blkdiag(A1 , . . . , AN ). The operator ⊗ denotes the Kronecker product. The nodes (systems/agents/sensors) are connected through a communication network, which is described as an undirected graph G = {V, E} with the vertex set V = {1, 2, . . . , N } with N ∈ N. The edge set E ∈ V × V depicts the connection between each node, where the pair (i, j) ∈ E indicates the existence of the communication channel from node j to node i for ∀i, j ∈ V, vice versa. The adjacent matrix A is defined based on the topology G, where ai,j = 1 if (j, i) ∈ E and ai,j = 0 otherwise. Moreover, we define the set Ni is the neighbor set of node i, i.e., Ni = {j ∈ V|(j, i) ∈ E}.
3
B. Problem Formulation In this paper, we investigate the estimation of a class of discrete-time Lur’e systems [34]. These systems, characterized by a feedback structure involving a linear dynamical component and a nonlinear element, are governed by the dynamics x0 (k + 1) = Ax0 (k) + Bf (x0 (k)),
∀k ∈ N0 ,
(1)
n
where x0 (k) ∈ R denotes the system state at time step k, with n ∈ N representing the state dimension. The matrices A ∈ Rn×n and B ∈ Rn×m are known reflecting the system structure, with m ∈ N being the dimension of the nonlinear output. However, the unknown nonlinear function f (·) : Rn → Rm embodies the latent trends governing the evolution of the system, which needs to be inferred. Remark 1: The system described in (1) can also be interpreted as a discrete-time control-affine system with matched nonlinearities, expressed as x0 (k + 1) = Ax0 (k) + B (u(x0 (k)) + f (x0 (k))), where u(·) : Rn → Rm represents a state-dependent control input. Employing a feedback controller u = Kx0 (k), the system (1) is retrieved from x0 (k + 1) = Ãx0 (k) + Bf (x0 (k)) with à = A + K. Furthermore, in the absence of linear dynamics, i.e., when A = 0, the system simplifies to a discrete-time nonlinear autonomous form: x(k + 1) = Bf (x(k)). Therefore, this versatility enables the framework in (1) to model a diverse array of dynamical systems, ranging from linear controlled systems to fully nonlinear autonomous ones, enhancing its applicability across various domains. In order to obtain the behaviors of the partially unknown system (1) over time, specifically concerning the states x0 and the dynamics f (·), we assume a collective of N ∈ R>0 networked sensors is deployed. This arrangement gives rise to a multi-agent or multi-robot configuration, where each agent or robot is equipped with its own sensor. Considering the partial observation of x0 , each subsystem’s measurement is given by yi (k) = Ci x0 (k) + vi (k),
∀i ∈ V,
(2)
where yi ∈ Rpi with pi ∈ I[1,n] is considered as the system output of the i-th sensor with Ci ∈ Rpi ×n as known, such that the pair (A, Ci ) is observable. The measurement noise denoted by vi is bounded by v̄i ∈ R>0 , i.e., ∥vi (k)∥ ≤ v̄i . Remark 2: Notice that the dimension of the output states yi in (2) is not strictly equal to the dimension of the state x0 , which indicates the incompleteness of the observations obtained through the sensor measurements. Given the challenging circumstance posed by the partial observations, the principal aim is to empower the multi-sensor distributed system with the capability to cooperatively infer both the system states x0 (k) and the unknown dynamics f (·). Let the state estimation of the system states x0 for sensor i ∈ V denote as xi and function estimation f (·) as fˆi (·), respectively. Then the desired performance defined as practical convergence is formulated as follows. Definition 1: The networked sensor system is said to achieve jointly practical convergence for the observation of x0 and the inference of f (·), if ∥x(k)−1N ⊗x0 (k)∥ ≤ ϵx ∧∥fˆ(x(k))−1N ⊗f (x0 (k))∥ ≤ ϵf
for k > 0k̄ with small constants ϵx , ϵf ∈ R>0 and k̄ ∈ N0 , T where x = [xT1 ,· · ·, xTN ]T , fˆ(x) = [fˆ1T (x1 ),· · ·, fˆN (xN )]T are the concatenated estimation for x0 and f (x0 ), respectively. For ease of notation, we assume that f (·) is scalar, i.e., m = 1. Consequently, the matrix B is simplified to a vector b ∈ Rn . Note that the derived results in the following sections can be extended to the multi-dimensional case using multioutput machine learning techniques, by assuming each output is independent [29], [35]. C. Gaussian Process Regression In order to infer the unknown function f (·), Gaussian process regression is employed, inducing a Gaussian distribution over function characterized by a kernel function κ(·, ·) : Rn × Rn → R≥0 and a mean function m(·) : Rn → R≥0 . Specifically, the representation is formalized as f (·) ∼ GP(m(·), κ(·, ·)). The mean function embodies prior knowledge of the unknown function and is commonly set to zero in the absence of observable trends of the system in advance. Moreover, the kernel function delineates the covariance between two samples and adheres to the following assumption. Assumption 1: The unknown function f belongs to the reproduced kernel Hilbert space (RKHS) Hκ corresponding to a given stationary and differentiable Lipschitz-continuous kernel κ(·, ·) with Lipschitz constant Lκ ∈ R>0 induced by Euclidean norm. Moreover, the kernel satisfies κ(ξ, ξ) = σf2 and σf ∈ R+ for ξ ∈ Rn . The induced inner product defined on RKHS p denotes ⟨·, ·⟩κ , and the RKHS norm defined as ∥f ∥2κ = ⟨f, f ⟩κ is bounded by Γ ∈ R+ , i.e., ∥f ∥κ ≤ Γ. The RKHS norm ∥ · ∥κ indicates the smoothness of the function f (·), such that the existence of the upper bound Γ is equivalent to the bounded derivative of f (·) requiring f (·) to be Lipschitz. Note that the Lipschitz system dynamics in (1) is a common prerequisite in continuous-time systems in reality, guaranteeing the dynamics (1) has a unique solution [36]. For instance, common nonlinear dynamics in robotics and control, such as friction models, aerodynamic drag, and smooth manipulator gravitational/Coriolis terms, are Lipschitz continuous over compact operating regions, ensuring that Assumption 1 is satisfied. The upper bound Γ of the RKHS norm can be approximated using the method in [37]. More recently, data-driven methods for over-estimating the RKHS norm bound have been developed in [38], which leverage neural network approximators and collected data. Therefore, Assumption 1 does not impose restrictive conditions and is broadly satisfied in real-world dynamical systems. As modeling f using GPs requires the capability to reconstruct the system state xi and the value of f (xi ) using the output yi from each sensor i, we introduce the data set Di contains Mi ∈ N samples {(ξ (ι) , φ(ι) )} for ι = 1, · · · , Mi described in the following assumption. Assumption 2: A scalar state φ(ι) is an observation of f (·) measured at the state ξ (ι) ∈ Rn for ι ∈ I[1,M ] , i.e., φ(ι) = f (ξ(ι) )+w(ι) . And w(ι) is the related observation noise bounded by a constant w̄ ∈ R>0 , i.e., |w(ι) | ≤ w̄ for ∀ι ∈ N0 . Assumption 2 admits the presence of noise in the training samples, particularly φ(ι) . This mirrors the real-world scenario,
4
where φ(ι) is influenced by environmental disturbances and uncertainties. The upper bound is derivable from the sensor measurements in (2) by considering the bounded nature of measurement noise, which is elucidated later in Theorem 1. Notably, extracting the information from the measurement yi is inaccessible for some sensors, even though (A, Ci ) is observable due to the unknown function f . In order to address this issue, we present a data collection strategy for building the training data set for GPs as discussed in Section III-A. Before that, we first present an estimated function for f (·) using the data set D with complete states. GP predicts the value of f (·) and point ξ as a Gaussian distribution with posterior mean µ(·) and variance σ 2 (·) as µ(ξ) = kDT (ξ)(KD + w̄2 IM )−1 φ, 2
σ (ξ) = κ(ξ, ξ) − kDT (ξ)(KD + w̄2 IM )−1 kD (ξ),
collection strategy and an analysis of noise propagation. We present the property of observability for the sensors as follows. Lemma 2: Consider a system (1) with states x0 ∈ Rn , which is measured by sensor i following (2) with the output dimension pi . If (A, Ci ) is observable, then there exists a set of matrices {Hi,d }d∈I[0,n−1] such that Ci An =
yi (k + n) = Ci An x0 (k) + vi (k + n) + Ci
Within this section, we initially present the data collection strategy in Section III-A. Subsequently, we propose the COINGP algorithm in Section III-B, followed by a rigorous analysis of its learning performance in Section III-C, ensuring the prediction guarantee for safety-critical applications.
Owing to the challenge induced by the combination of incomplete measurements yi of states x0 from sensor i ∈ V and the unknown dynamics f (·), the observability requirement of (A, Ci ) is insufficient, resulting in the inadequacy of usable training data pairs. However, local observability of (A, Ci ) of every agent is not required when agents do not maintain individual GP models and instead rely on predictions received from their neighbors through the cooperative learning approach described in Section III-B. Moreover, while the work [41] addresses only a specific case, we need to consider arbitrarily structured matrices A, b, and Ci for constructing the datasets. To solve this problem, a condition for sufficient data acquisition is introduced, along with a corresponding data
An−d−1 bf (x0 (k + d)).
(7)
Combining the outcomes in Lemma 2, Eq. (7) is further reformulated as n−1 X Hi,d Ci Ad x0 (k) + vi (k + n) yi (k + n) = d=0
+ Ci
n−1 X
An−d−1 bf (x0 (k + d)),
(8)
d=0
which is equivalent to yi (k + n) =
n−1 X
Hi,d yi (k + d) − vi (k + d) + vi (k + n)
d=0 n−1 X
−
Hi,d Ci
d=1 n−1 X
d−1 X
Ad−l−1 bf (x0 (k + l))
l=0
An−d−1 bf (x0 (k + d)).
(9)
d=0
Due to the fact that the double summation1 n−1 d−1 X X Hi,d Ci Ad−l−1 bf (x0 (k + l)) l=0 d=1 n−1 X n−1 X
Hi,l Ci Al−d−1 bf (x0 (k + d)),
=
d=0
A. Data Collection Strategy
n−1 X d=0
+ Ci III. C OOPERATIVE L EARNING WITH PARTIAL M EASUREMENTS
(6)
Proof: See the supplementary material. Considering the system dynamics defined in (1) along with the measurement function in (2), the output yi (k + n) for all k ∈ N is expressed as
(4)
where β = Γ2 − φT (KD + w̄2 IM )−1 φ + M . Lemma 1 establishes a theoretical error bound in Gaussian processes, a methodology widely applied in safety-critical contexts, as exemplified in [39]. While the coefficient β scales with the dataset size, given by β ≤ Γ2 + M , the work in [40] assures a concurrent reduction in the posterior variance σ 2 (·).
Hi,d Ci Ad .
d=0
(3)
where kD (ξ) = [κ(ξ, ξ (1) ), · · · , κ(ξ, ξ(M ) )]T ∈ RM , φ = [φ1 , · · · , φM ]T ∈ RM and KD = [κ(ξ (i) , ξ(j) )]i,j=1,··· ,M ∈ RM ×M . Therefore, the posterior mean µ(·) is adopted to estimate f (·), and the posterior variance σ 2 (·) quantifies the prediction error shown in the following lemma. Lemma 1 ( [37]): Assume the unknown function f (·) satisfying Assumption 1 and predicted by GP regression using data set D with M training samples following Assumption 2. Then, the prediction error is bounded as p ∀ξ ∈ Rn , (5) |µ(ξ) − f (ξ)| ≤ η(ξ) := βσ(ξ),
n−1 X
(10)
l=d+1
then yi (k + n) is written as yi (k + n) =
n−1 X
Hi,d yi (k + d)−vi (k + d) +vi (k + n)
d=0 n−1 X
n−1 X
d=0
l=d+1
+
Ci An−d−1−
Hi,l Ci Al−d−1 bf (x0 (k+d)). (11)
For facilitating the construction of the training dataset and notional simplicity, we define an auxiliary state ỹi (k) from measurements as n−1 X ỹi (k) = yi (k + n) − Hi,d yi (k + d) (12) d=0 1 The proof of Eq. (10) is provided in the supplementary material.
5
and an auxiliary vector from matrices as ρi,d = Ci An−d−1 −
n−1 X
Hi,l Ci Al−d−1 b.
(13)
l=d+1
Based on the defined auxiliary variables (12) and (13), then one has n−1 X ỹi (k) = ρi,d f (x0 (k + d)) − Hi,d vi (k + d) d=0
+ vi (k + n).
(14)
To assess the feasibility of constructing a training dataset conducive to the inference of the unknown function f and the estimated state vector x0 from measurements yi , we introduce the concept of collectability for a system. This concept indicates the sensor’s capability in acquiring a training dataset Di associated with the sensor i ∈ V that satisfies Assumption 2, which is described in the following definition. Definition 2: The collectability of the training dataset of a system i ∈ V is affirmed if and only if its corresponding dataset Di can be established by obtaining measurements yi and satisfies Assumption 2. Considering Definition 2, the sufficient condition of the collectability of a system is given as follows. Lemma 3: Consider the system (1) measured by sensor i ∈ V following (2). If (A, Ci ) is observable and there exists one d∗i ∈ I[0,n−1] such that ( ̸= 0pi ×1 , if d = d∗i ρi,d , d ∈ I0,n−1 , (15) = 0pi ×1 , otherwise then the data set Di satisfying Assumption 2 is obtainable and constructible. Proof: See the supplementary material. Specifically, with ti ∈ Rpi satisfying tTi ρi,d∗i ̸= 0 and Ti ∈ Rn×npi rendering Ti Oi invertible, the training pair is constructed from the measurements as −1 T φi (k) = tTi ρi,d∗i ti ỹi (k − d∗i ) = f (x0 (k)) + wi,1 (k), (16) −1 ξi (k) = Ti Oi Ti Yi (k) − Gi φi (k) = x0 (k) + wi,2 (k), (17) where Yi (k) = [yiT (k), · · · , yiT (k + n − 1)]T , φi (k) = [φi (k), · · · , φi (k + n − 2)]T , Gi is the block lower-triangular Toeplitz matrix whose (r, c)-th block is Ci Ar−c−1 b for r > c and 0pi ×1 otherwise, and the transformed noises wi,1 (k), wi,2 (k) collect the measurement noise contributions. It is worth mentioning that the sensor measurements and the training datasets of individual systems are not transmitted among the interconnected subsystems, thereby alleviating communication overhead. Moreover, despite the inability of all systems to assemble the data set Di , cooperative learning mechanisms are leveraged to mitigate this limitation. Through cooperative learning, neighboring systems contribute to compensating for the unavailability of data by sharing predictions derived from GP models trained on non-empty datasets, which is elaborated in Section III-B.
Remark 3: Condition (15) indicates the auxiliary variable ỹ(k) can be written as the function of f (x(k + d∗i )), and is irrelevant to the unknown function f (x(·)) and system state x at other time steps, which means the value of f (x(k)) can be recovered using only ỹ(k − d∗i ) as shown in (16). Note that this condition depends on the selection of Hi,d , which may have multiple choices when pi > 1. For sensors whose valid Hi,d for collectability check is hard to obtain, we assume they are not collectable and apply the cooperative learning strategy to enhance their prediction performance. The geometric picture is that observability ensures the observation space “fully captures” the state space dynamics, making the redundancy encoding via {Hi,d } possible. For each row j ∈ I[1,pi ] of Ci An , solve the linear system: gj C i A n =
pi n−1 XX
d hj,i d,q gq Ci A .
d=0 q=1 n Find coefficients {hj,i d,q } such that the j-th row of Ci A is expressed as a linear combination of rows from Oi . Stack these into matrices Hi,d ∈ Rpi ×pi as defined in the proof.
Since the data collection condition in Lemma 3 is a sufficient condition, to find such Hi,d with d = 0, · · · , n − 1 satisfying (15) is nontrivial. Nevertheless, the search is not exhaustive. For a fixed candidate index d∗i , both (6) and the vanishing conditions ρi,d = 0, d ̸= d∗i , are linear in the entries of {Hi,d }, so collectability reduces to a linear feasibility problem with np2i unknowns and pi (2n − 1) equations over the n candidates d∗i , which can be solved before deployment. Since this problem is overdetermined for pi = 1 and underdetermined for pi ≥ 2, multi-output sensors generically admit a valid {Hi,d }, whereas scalar-output sensors are collectable only for structured Ci . In this light, we discuss some special cases. • For single output measurement with pi = 1, the matrices Hi,d degraded to scalar are unique, since Oi ∈ Rn×n . Then, the evaluation of (15) is fast for unique Hi,d , but no freedom is left to enforce (15). • For a discrete-time canonical system in control and observation with sampling time T ∈ R+ as In−1 0(n−1)×1 0 T In−1 , A= + (n−1)×1 01×(n−1) 0 0 01×(n−1) b = [01×(n−1) , 1]T , Ci = [1, 01×(n−1) ], it has Ci Ad b = 0 for d = 0, · · · , n − 2 and Ci An−1 b = 1 satisfying (15) without considering Hi,d . Remark 4: Considering (12) and (16), obtaining φi (k) for computing ξi (k) via (17) requires access to measurements yi (k), · · · , yi (k +2n−1). Consequently, under the worst-case scenario, acquiring the training data pair {ξi (k), φi (k)} demands a sequential time interval of k + 2n − 1 steps. However, when Ci is an identity matrix, it is straightforward to derive that ξi (k) = yi (k) and φi (k) = bT (y(k+1)−Ay(k))/(bT b). Considering the process of derivation of Lemma 3, for a system i ∈ V satisfying (15), the approach for data collection is summarized in Algorithm 1. Remark 5: Algorithm 1 allows online learning for each agent satisfying (15) by adding new samples into the training
6
Algorithm 1 Data acquisition on subsystem i
B. Observer-based Cooperative Learning
1: Condition: (A, Ci ) is observable satisfying (15); 2: Initialization: Di = ∅, Yi = ∅, φi = ∅, k = 0; 3: while Terminate condition for collection is not satisfied
do Collect y(k); if |Yi | < n then Yi ← [Yi , yi (k)] else Yi ← [Yi (pi + 1 : npi ); yi (k)]; Calculate φi (k − n + d∗i ) from (16); 10: φi ← [φi , φi (k − n + d∗i )] 11: if |φi | ≥ n − 1 then 12: Calculate ξi (k − n + 1) from (17); 13: Di ← {Di , {ξi (k − n + 1), φi (k − n + d∗i )}}; 14: end if 15: end if 16: k ←k+1 17: end while 4: 5: 6: 7: 8: 9:
data set, which continuously improves the observation and prediction performance of x0 and f (·) during operation as detailed in Section III-B. Since one valid training pair requires measurements spanning up to 2n − 1 steps, the newest sample in Di lags the current state by up to 2n − 1 steps; as f (·) is time-invariant, such delayed pairs remain valid samples in the sense of Assumption 2, and the delay only postpones the arrival of samples without altering their distribution along the trajectory of x0 . Consequently, the delay affects only the posterior variance σi (xi (k)) in the error bound of Section III-C, which grows at most linearly with the distance travelled by x0 within 2n − 1 steps, and hence with n for fast-varying systems, when x0 enters a region not yet covered by data, while the influence of the delay vanishes in revisited regions, since the GP exploits the entire data set rather than the most recent sample. The ultimate boundedness of the observation and prediction errors in Theorem 2 is unaffected, as it holds for any Di satisfying Assumption 2, including an empty one. Given the collectible dataset, we now investigated the noise bound of the measurements for f in the following theorem. Theorem 1: Considering the dynamical system (1) and the sensors i ∈ V equipped with measurement systems described by (2), the collection of the dataset Di satisfies the conditions in Lemma 3. The dataset Di is collected through Algorithm 1 for GP prediction of f (·) under Assumption 1, then Assumption 2 holds for Di with the noise wi , which is bounded by p |wi (k)| ≤ w̄i = w̄i,1 + Lf w̄i,2 , ∀k ∈ N, (18) √ where Lf = 2Lκ Γ, and w̄i,1 = tTi ρi,d∗i
−1
n−1 X ∥ti ∥ 1 + ∥Hi,d ∥ v̄i ,
(19)
d=0
w̄i,2 = ∥ Ti Oi
−1
√ √ Ti ∥ ∥Gi ∥ n−1w̄i,1 + nv̄i .
(20)
with ti and Ti chosen according to (16) and (17), respectively. Proof: See the supplementary material.
To address the potential absence of some training datasets for GP models, we introduce the COIN-GP approach for the estimated function of the system i ∈ V, expressed as fˆi (k + 1) = γ1
N X
ãi,j (k)(fˆi (k) − fˆj (k))
(21)
j=1
+ γ2 ϖi (k)(fˆi (k) − µi (xi (k))) + µi (xi (k + 1)), where γ1 , γ2 ∈ R are estimation gains to be designed later, fˆi (k) represents the estimation of f (xi (k)), ∀k ∈ N0 for simplicity. xi (k + 1) = Axi (k) + bfˆi (k) + Li (Ci xi (k)−yi (k)). (22) Remark 6: The communication graph is assumed undirected because the proposed framework relies on bidirectional information exchange between neighboring agents. This requirement can be relaxed by utilizing a directed graph in which each agent only broadcasts its own predictions, though this typically leads to degraded performance. Specifically, Li ∈ Rn×pi denotes the observer gain matrix selected such that A+Li Ci is Schur, that is, all eigenvalues lie strictly within the unit circle in the complex plane. The design of Li is performed locally at each node, so that the observer gain design is implemented in a fully distributed manner. If (A, Ci ) is observable, such a matrix Li always exists and can be obtained using Ackermann’s formula with given desired eigenvalues. The local prediction from the GP model for system i is denoted by µi (·), leveraging the individual dataset Di if available. The gains γ1 , γ2 ∈ R in (21) are designed constants to balance the weights of the cooperative and individual learning. The consensus weights ãi,j and ϖi are designed as ãi,j (k) =
ai,j σi2 (xi (k)) , σj2 (xj (k))+ w̄i2
ϖi (k) =
σi2 (xi (k)) , 2 σf,i
(23)
where ai,j is the entry of the adjacent matrix A at the i-th row and j-th column. Therefore, the computation outlined in (21) solely relies on neighboring information. Consequently, a distributed computation framework is realized, wherein the scalability of the proposed methodology is attainable, rendering it well-suited for large-scale sensor networks. Moreover, (23) ensure bounded consensus weights, i.e., ãi,j (k) ∈ 2 {0, σf,i /w̄i2 } and ϖi (k) ∈ (0, 1]. In the early exploration phase, all posterior variances remain close to the priors, so (23) reduces to uniform-weight cooperation. As data accumulate asymmetrically, a data-rich agent i assigns a small weight ãi,j to a data-poor neighbor j while the latter keeps a large weight ãj,i , so that on a connected graph the predictions of agents without data follow those of data-rich agents, directly or through intermediate agents. Remark 7: The evaluation of the proposed estimator in (21) on agent i ∈ V only requires the information exchange for the estimated prediction fˆj (k) and posterior variance σj2 (xj (k)) from j ∈ Ni . Without sharing the training samples or the local prediction µj (xj (k)), the proposed method in (21) protects the data privacy of the systems. Moreover, the proposed method
7
saves the computation of µj (xi ) required for the conventional aggregation strategy used in MAS as [21], [25], [26]. Throughout this work, the inter-agent communication is assumed synchronous and delay-free: at every sampling instant, each agent exchanges its current prediction and posterior variance with its neighbors before the updates are evaluated. For imperfect timing, COIN-GP can be integrated with delaycompensation strategies [42], asynchronous distributed GP updates [43], or dynamic average consensus [22], which is left for future work.
collected data set Di satisfies Assumption 2. Use the observerbased cooperative learning strategy proposed in (21) satisfying Assumption1 and choose proper γ1 ∈ R<0 , γ2 , α1 ∈ R>0 and symmetric positive matrix Q ∈ R2×2 such that λ∗f =
max
|λi (γ1 L(k)+γ2 B(k))| < 1,
C. Performance Analysis
λ∗s =
max
To analyze the estimation and identification performance for system states x0 and dynamics f (·), we define the their associated prediction and observation error for f (·) and x0 as
Q=−
efi (k) = fˆi (k) − f (x0 (k)),
esi (k) = xi (k) − x0 (k), (24)
respectively. Combining the system dynamics (1) and output observation (2) with estimated function (21), the dynamics of prediction and observation error are reformulated as efi (k + 1) = γ1
N X
ãi,j (k)(efi (k) − efj (k))
(25)
j=1
+ γ2 ϖi,i (k)(efi (k) + µ̃i (k)) − µ̃i (k + 1), es (k + 1) = Axi (k) + bfˆi (k) + Li (Ci xi (k) − yi (k)) i
− Ax0 (k) − bf (x0 (k))
(26)
= (A + Li Ci )esi (k) + befi (k + 1) − Li vi (k), with µ̃i (k) = f (x0 (k)) − µi (xi (k)). Moreover, we define the concatenated prediction and observation error as ef (k) = [ef1 (k), · · · , efN (k)]T and es (k) = [es1 (k)T , · · · , esN (k)T ]T respectively, whose dynamics are written as ef (k + 1) = γ1 L(k) + γ2 B(k) ef (k) + γ2 B(k)µ̃(k) − µ̃(k + 1), s
s
(27) f
e (k + 1) = Θe (k) + (IN ⊗ b)e (k) − Υv(k),
(28)
with the concatenated prediction error denotes by µ̃(k) = [µ̃1 (k), · · · , µ̃N (k)]T , and the concatenated measurement T noise written as v(k) = [v1T(k),· · ·, vN (k)]T , where 2 Σ(k) = diag(σ12 (x1 (k)), · · · , σN (xN (k))), 2 2 Σ̄(k) = diag(σ12 (x1 (k)) + w̄12 , · · · , σN (xN (k)) + w̄N ),
L(k) = D(k) − Σ(k)AΣ̄−1 (k), Υ = blkdiag(L1 , · · ·, LN ), D(k) = diag(Σ(k)AΣ̄−1 (k)1N ), −2 −2 B(k) = diag(σf,1 , · · · , σf,N )Σ(k),
Θ = blkdiag(A + L1 C1 , · · · , A + LN CN ). Given the error dynamics (27) and (28) of prediction and observation of the distributed system, the joint learning performance for the distributed system is shown in the following theorem. Theorem 2: Consider a system (1) measured by a distributed system (2) and bounded measurement noise v̄i for i ∈ V. Each subsystem is observable and under the condition (15) with the training samples using Algorithm 1, such that the
(29) ζ0 = λ(Q)−2α1 max(λ∗f , λ∗s ) max γ22 + 1, ∥Υ∥2 ∥P ∥ 2 (30) − α12 max γ22 + 1, ∥Υ∥2 ∥P ∥ > 0, i=1···N,k∈N
with |λi (A + Lj Cj )|, T ∗ 0 λf 0 P + P, λ∗s ∥b∥ λ∗s
j=1···N,i=1···n
λ∗f ∥b∥
(31) (32)
where the matrix P is the solution of discrete-time Lyapunov equation ΦT1 P Φ1 − P = −Q. Then there exists k̄ ∈ N such that the observation and prediction error are ultimately bounded by |efi (k)| ≤ ∥ef (k)∥ ≤ ē,
∥esi (k)∥ ≤ ∥es (k)∥ ≤ ē,
(33)
respectively, for ∀k > k̄, with the ultimate error bound ē formulated as q q (34) ē = λ̄(P )/λ(P )ζ0−1 ζ1 + ζ12 + ζ0 ζ2 α2 , where the coefficients ζ1 and ζ2 are defined by ζ1 = τ (λ∗ + α1 τ )∥P ∥ and ζ2 = τ 2 ∥P ∥. The scalar α2 is defined as √ α2 = α1−1 L2f N + ∥η̄∥ + ∥η̄ ∗ ∥ + ∥v̄∥, (35) where ∥v̄∥ = ∥[v̄1 , · · · , v̄N ]T ∥, the prediction accuracy related terms η̄ and η̄ ∗ are defined as η̄ = [η̄1 , · · · , η̄N ]T ∗ T ] with η̄i = supk∈N0 ηi (xi (k)) and and η̄ ∗ = [η̄1∗ , · · · , η̄N ∗ 2 −1 3 η̄i = (βσf,i ) η̄i , respectively. Moreover, the constants λ∗ = max(λ∗f , λ∗s ) and τ = max(γ22 + 1, ∥Υ∥2 ) denote the maximal singular values of Φ1 and Φ2 respectively, in which λ∗s = maxj=1,··· ,N maxi=1,··· ,n |λi (A + Lj Cj )| ∈ [0, 1) due to the block triangular structure of Λ and the designed Li . Proof: See the supplementary material. Remark 8: Notably, the stability guarantees in Theorem 2 rely on the global spectral properties of the graph, particularly the eigenvalues of γ1 L(k) + γ2 B(k). While exact verification requires global topology information, practitioners can derive a conservative sufficient condition for ad-hoc networks using the Gershgorin Circle Theorem and Weyl’s inequality. Since λmax (L) ≤ 2N and λmax (B) ≤ 1, choosing parameters 1 such as |γ1 | < 4N and γ2 < 12 easily satisfies the bound λmax (γ1 L + γ2 B) < 1 without requiring full graph knowledge, at the price of a conservatism that grows with N : the gain is scaled by the worst case 4N , whereas the spectrum of a sparse graph grows only with the node degree, so the resulting design converges more slowly and settles at a larger ultimate error than a globally tuned one. Nevertheless, extending the current framework to a fully distributed, topology-agnostic adaptive tuning mechanism, where agents do not even require knowledge of the network upper bounds (e.g., N ), remains a highly relevant direction for future work.
8
TABLE I: Sensor configurations
fˆ2 (x2 ) σ2 (x2 )
1 y1
y2
fˆ3 (x3 ) fˆ2 (x2 ) σ3 (x3 )
fˆ4 (x4 ) fˆ (x ) 1
2
fˆ1 (x1 ) σ1 (x1 )
1
2
σ4 (x4 ) σ1 (x1 )1
σ2 (x2 )
fˆ3 (x3 ) σ3 (x3 )
4
3
fˆ4 (x4 ) σ4 (x4 )
y4
y1
y3
y2
Properties Sensor 1 Sensor 2 Sensor 3 Sensor 4 Ci in (2) [1,0] [1,0; 2,1] [1,1; 1,2] [0,1; 3,2] Hi,0 in (6) [0] − 0.5[1,0; 1,0] − Hi,1 in (6) [1] − 0.5[0,1; 0,1] − t in (16) 1 − [2; 1] − T in (17) I2 − [I2 , 02×2 ] − λ1,2 (A+Li Ci ) [0.4,0.5] [0.4,−0.5] [0.5,−0.3] [0.4,−0.3]
3
4 y3
y4
0
Fig. 1: Communication topology among sensors/agents. Teal nodes represent the GP model instances, and magenta nodes represent the actual sensing nodes that transmit measurements from system node 0. The black numbers on the nodes indicate they lack individual training datasets, i.e., the node 2 and 4, while the white numbers imply the existence of datasets on these nodes, i.e., node 1 and 3. Theorem 2 shows that the observation and prediction using the proposed observer-based cooperative learning in (21) will converge to a small area around the zero, achieving practical convergence defined in Definition 1. The additional condition requiring ζ0 > 0 is practically not restrictive, which can be easily achieved by choosing sufficiently small α1 . However, smaller α1 induces larger ultimate bound ē by considering the inverse relationship between α1 and α2 . Unlike the setting in continuous-time without measurement noise such as in [44], increasing the control gains by choosing eigenvalues of γ1 L(k) + γ2 B(k) or A + Li Ci for i ∈ V closer to 0 may not tighter the error bound ē. This is because the induced high gains increase the sensitivity of the multi-sensor system to the model uncertainty and measurement noise reflected by a larger singular value of Φ2 . IV. S IMULATIONS In Section IV-A, we describe the simulation setting. The proposed approach is then compared with prevailing distributed GP frameworks and an adaptive neural-network baseline through Monte Carlo tests in Sections IV-B and IV-C, and its robustness under unreliable communication is examined in Section IV-D. A. Simulation Setting We consider a system represented by (1) with state dimension n = 2.2 In this case, the systems state is denoted T as x0 = [x 0,1 , x0,2 ] , and the system matrices are defined 1 1 as A = and b = [0 1]T . The state trajectory, 0 0 which the distributed system aims to estimate, is given by x0,1 (k) = arctan(a1 k) sin(a2 k), where a1 and a2 ∈ R chosen as 0.01 and 0.05, respectively. The overall COIN-GP framework is illustrated in Fig. 1. The simulation is conducted over a maximum time step of k = 500. In order to estimate the system states x0 and unknown dynamics f (·), the distributed system consisting of N = 4 2 The full simulation configuration is given in the supplementary material.
sensors with an undirected topology shown in Fig. 1 is employed. Each subsystem follows (2) with noise variance v̄i = 0.001, ∀i ∈ V. The specific configurations for each sensor are detailed in Table I. Note that Li is chosen according to λ1,2 (A + Li Ci ) using robust pole placement [45]. For agents 2 and 4, data collection is deliberately deactivated, so that their data sets remain empty and the matrices Hi,d , ti and Ti are left unspecified in Table I. This setting serves to demonstrate the efficacy of the proposed cooperative learning approach, which does not require all agents to collect data or participate in modeling the unknown function. For agent 1 and 3, the training data set D1 and D3 are collected online using Algorithm 1 satisfying Assumption 2. Additionally, each Gaussian process model established is initialized at k = 0 with an empty data set, i.e., Di = ∅, and prior knowledge regarding the kernel is chosen as automatic relevance determination 2 squared exponential kernel with κi (ξ, ξ ′ ) = σf,i exp − P n 1 ′ 2 2 T d=1 |ξd − ξd | /ld , where ξ = [ξ1 , ξ2 ] . The hyperparam2 eters are set to σf,i = 0.017, l1 = 1.8 and l2 = 0.025, which are obtained through hyperparameter optimization using the non-empty dataset D1 and D3 with k = 500. B. Performance comparison To demonstrate the superiority of the proposed observerbased online cooperative learning (21) with the devised state estimation method (22), we conduct a comparative simulation for analyzing the observation and prediction performance against the prevailing distributed GP frameworks. It’s crucial to note that due to the existing methods rely on complete observations necessitating full measurement of system states x0 and are unsuitable for output-based control scenarios. In order to ensure a fair comparison, we employ the same observer (22), across all compared approaches. The details of our experimental setup are elucidated below. The compared methods are the proposed cooperative learning (21) withP γ2 = 1 and the statedependent gain γ1,i (k) = −γ2 ϖi (k)/ j∈Ni ãi,j (k), under which the consensus term in (21) reduces to the weighted average of the neighbors’ estimates with the weights (23); ˆ local learning [46] Pwith fi (xi (k)) = µi (xi (k)); MoE [47] with fˆi (xi (k)) = j∈{i,Ni } ai,j µj (xj (k)); and PoE [21] with P −2 j∈{i,Ni } ai,j σj (xj (k))µj (xj (k)) fˆi (xi (k)) = . P −2 j∈{i,Ni } ai,j σj (xj (k)) Moreover, with ai,i = 1, the shorthands µj = µj (xj (k)) and σj2 = σj2(xj (k)), and the log-variance weight ρj (k) = 1 2 2 2 ln σf,j /σj , we evaluate GPoE [20] with P −2 j∈{i,Ni } ai,j ρj (k)σj µj ˆ fi (xi (k)) = P −2 , j∈{i,Ni } ai,j ρj (k)σj
9
−2 j∈{i,Ni } ai,j σj µj −2 , fˆi (xi (k)) = P P ai,j σj−2 + 1− ai,j σf,i j∈{i,Ni } j∈{i,Ni }
P
GP Error
BCM [48] with
∥µ(x)−1N ⊗f (x0 )∥ 10
10−4 100
and RBCM [20] with
Finally, to assess the trade-off between capacity and formal guarantees, an adaptive radial-basis-function neural network baseline is evaluated, replacing each agent’s GP by an RBFNN with linear output weights, fˆiNN (x) = ŵiT ϕi (x), n 1X (xd − cim,d )2 , m = 1, . . . , Mϕ , [ϕi (x)]m = exp − 2 s2d d=1
with Mϕ = 100 Gaussian features whose centers cim are assigned randomly over the operating region, since no prior data is available to place them, mirroring the a-priori fixed GP kernel hyperparameters, and whose widths sd equal the GP lengthscales. Upon each reconstructed pair (ξi , φi ), the weights are updated online by normalized least-mean-squares with σ-modification, ϕi (ξi ) φi − ŵiT ϕi (ξi ) ŵi ← ŵi + η − ησm ŵi , ε + ∥ϕi (ξi )∥2 with η = 0.5, σm = 10−4 , ε = 10−8 , and a replay buffer sharing the GP data budget, under the identical observer, datacollection pipeline, topology, and noise realizations. RBFNNLocal is non-cooperative, while RBFNN-Coop applies the consensus law (21) with µi (·) replaced by the NN prediction; since the NN provides no posterior variance, the weights (23) degenerate to uniform weights (ãi,j = ai,j , ϖi = 1, gains retuned to γ1 = −0.05, γ2 = 0.2). All NN hyperparameters and consensus gains are grid-searched in favor of the baseline. All sensors are initialized at k = 0 with same state estimation for each methods, specifically x1 (0) = [0.6294, 0.0406]T , x2 (0) = [−0.7460, 0.0413]T , x3 (0) = [0.2647, −0.0402]T and x4 (0) = [−0.4430, 0.0047]T . On the individual-agent level, the benefit of cooperation concentrates on the agents without data collectability: over the 100 Monte Carlo runs of Section IV-C, the per-agent mean prediction error of agents 2 and 4 drops from 0.067 (Local) to 0.031 (COIN-GP), while for the data-collectable agents COINGP performs on par with or better than local learning; MoE instead weighs empty and trained models equally, deteriorating the accuracy of the reliable agents. Fig. 2 shows the prediction error from GP, which is bounded by its theoretical upper bound derived in Lemma 1. Moreover, the RKHS norm bound Γ = 5 in computing η(·) is estimated by using a data-driven method provided in [37]. C. Monte Carlo Test To show the generalizability of our method, the Monte Carlo tests are conducted 100 times for each method. The
200
300
400
k
Fig. 2: Prediction errors of local GP models and their theoretical bounds. ∥fˆ(x)−1N f (x0 )∥ ∥x−1N ⊗x0 ∥
−2 j∈{i,Ni } ai,j ρj (k)σj µj −2 . P P ai,j ρj (k)σj−2 + 1− ai,j ρj (k) σf,i j∈{i,Ni } j∈{i,Ni }
P
fˆi (xi (k)) =
Error bound
−2
100
10−1
10−1
10−2
100
200
300
400
k COIN-GP BCM
PoE RBCM
MoE RBFNN-Local
Local RBFNN-Coop
GPoE
Fig. 3: State observation (top) and prediction (bottom) errors over 100 Monte Carlo runs (M = 20, p = 0.2). initial observer states xi (0) for i ∈ V are sampled randomly from a uniform distribution in [−1, 1] × [−0.05, 0.05], and the system dynamics is also randomized with a1 following a uniform distribution between 0.01 to 0.05 and a2 satisfying a uniform distribution between 0.05 to 0.1. To reflect realistic long-term networked deployments, the main comparison is conducted with the streaming implementation, i.e., slidingwindow datasets with budget M = 20 ensuring a bounded per-step complexity, under mild communication unreliability (per-step link failure probability p = 0.2). In addition to the GP-based baselines, the adaptive RBFNN baselines introduced in the experimental setup are included in the same figures and tables. As illustrated in Fig. 3, where solid lines denote the mean and shaded areas the ±1 standard deviation (shown for COINGP, PoE, MoE, and Local; the bands of the newly added baselines are omitted for legibility), COIN-GP achieves the lowest observation and prediction errors among all nine methods, including the strongest uncertainty-weighted aggregators (GPoE, BCM) and the RBFNN baselines, while RBCM diverges under link failures in this streaming setting and leaves the displayed range. The same ranking holds for every link failure probability (Section IV-D) and for every data budget M ∈ {20, 100, 500}. As reported in the last column of Table II, one full COIN-GP time step takes 0.33 ms on a single CPU core, on par with purely local learning (0.34 ms) and only marginally above the static and RBFNN baselines (0.28– 0.31 ms). For a clear comparison, we have numerically listed several
10
observation error
Methods COIN-GP PoE MoE Local GPoE BCM RBCM RBFNN-Coop RBFNN-Local
Observation error Mean Median RMSE 0.144 0.129 0.170 0.197 0.193 0.214 0.248 0.257 0.278 0.237 0.233 0.257 0.187 0.185 0.199 0.187 0.185 0.199 div. 0.928 div. 0.226 0.218 0.248 0.223 0.210 0.240
Prediction error Mean Median RMSE 0.075 0.066 0.089 0.099 0.097 0.110 0.115 0.120 0.130 0.125 0.127 0.138 0.096 0.094 0.105 0.096 0.094 0.105 div. 0.103 div. 0.123 0.119 0.137 0.124 0.119 0.134
Time [ms] 0.33 0.29 0.29 0.34 0.30 0.30 0.31 0.28 0.28
0.14
prediction error
TABLE II: Observation and prediction errors and wall-clock time per step using different methods (k ≥ 100, 100 runs, M = 20, p = 0.2, single CPU core)
0.25
0.12
0.2
0.1 0.08
0.15 0
0.2
0.4
0.6
0.8 0
0.4
0.6
0.06 0.8
p
p COIN-GP GPoE
0.2
PoE BCM
MoE RBCM
Local
Fig. 4: Mean observation (left) and prediction (right) errors (k ≥ 100, 20 Monte Carlo runs, M = 20) versus the link failure probability p. evaluation metrics for both observation error and prediction error in Table II. To mitigate the impact of initial errors, all values are calculated using data starting from k = 100. COINGP attains the best mean, median, and RMSE for both the observation and the prediction error, while RBCM diverges under link failures in the streaming setting (denoted div.). Therefore, the effectiveness and superiority of the proposed method are clearly demonstrated. Regarding the neural-network baseline, three observations substantiate the capacity-versus-guarantee trade-off: the RBFNN matches the local GP almost exactly, confirming comparable capacity; uniform-weight cooperation provides the NN with at best marginal benefit, since without calibrated uncertainty the consensus incorporates the predictions of agents with empty datasets, whereas the variance-weighted fusion of COIN-GP reduces the error of local learning by 40%; and the GP prediction is accompanied by the deterministic bound of Lemma 1 at every time step (cf. Fig. 2), whereas no comparable bound exists for the RBFNN error. D. Unreliable Communication To assess the robustness of the cooperative mechanism under degraded connectivity, every communication link fails independently at each time step with probability p ∈ {0, 0.2, 0.4, 0.6, 0.8}; a failed link removes the neighbor from the consensus term of COIN-GP and from the aggregation of the static baselines at that step, using the streaming setting of the main comparison (M = 20). As shown in Fig. 4, COIN-GP is the most accurate method for every p, since its recursion retains the last consensus estimate whenever a neighbor is unavailable, whereas the
memoryless aggregators approach the non-cooperative baseline as p grows. The dashed BCM curve is close to GPoE in Fig. 4, with mean errors differing by less than 2 × 10−4 for every p. The prior-variance experts of the non-collecting agents are exactly cancelled by the BCM prior-correction term and receive zero weight ρj in GPoE, so both aggregators fuse the same informative experts and differ only in their weighting. RBCM diverges for every p > 0 (mean errors above 103 ), as its prior-correction denominator vanishes whenever the effective neighborhood contains only near-prior experts, a situation that link failures readily create under tight budgets, confirming that the robustness of COIN-GP stems from its uncertainty-aware recursive structure. V. C ONCLUSION AND D ISCUSSION This paper introduced COIN-GP, an observer-based cooperative learning framework that combines distributed state estimation with Gaussian process regression for systems with known linear dynamics, an unknown nonlinear component, and partial state measurements. The framework consists of three elements, namely a data collection strategy together with sufficient conditions under which a sensor is able to reconstruct training data from its own incomplete measurements, an uncertainty-weighted consensus law that lets sensors without collectable data benefit from the models of their neighbors, and a joint analysis of the coupled estimation and learning loop that yields a single error bound covering both the state estimate and the model prediction. The central conceptual message is that the posterior variance of the Gaussian process is not merely an accuracy indicator but the quantity that renders cooperation meaningful, since it tells each sensor how much of a neighbor’s prediction to trust. It is precisely this calibrated weighting, which is absent both in aggregation schemes without a probabilistic model and in deterministic function approximators, that makes the interconnection of observer and learner analyzable and keeps the cooperative estimate reliable when parts of the network hold no data at all. Several limitations delineate the scope of the present results and suggest directions for future research. First, for longterm real-time deployment, the computational complexity of the Gaussian process must be prevented from growing unboundedly. A sliding-window budget already bounds the perstep cost uniformly over time at a moderate accuracy loss, and more elaborate streaming techniques, such as logGP [49] and SkyGP [50], can be seamlessly embedded into the datacollecting agents. Second, the system matrices (A, B) are assumed to be known. Bounded perturbations (∆A, ∆B) split into a matched part in the range of B, which is absorbed into f via f˜(x) = f (x) + B † (∆Ax + ∆Bf (x)) and hence learned online, and an unmatched remainder, which acts as a bounded disturbance and enlarges the error bound ē in proportion to its magnitude. A quantitative robustness analysis and the extension to fully unknown (A, B) are left for future work. Third, the collectability condition (15) is sufficient but not necessary, and it is verified per sensor for a fixed structure of Hi,d . Characterizing collectability by necessary
11
and sufficient conditions, and exploiting the network structure so that sensors jointly satisfy a collectability condition that none of them meets individually, would considerably widen the class of admissible sensor configurations. Beyond these points, the framework is developed for estimation, whereas the joint error bound is exactly the ingredient required by safety-critical control; embedding COIN-GP into distributed predictive or safety-filter architectures is therefore a natural next step. R EFERENCES [1] Y. Huang, P. Li, S. Ma, S. Yan, M. Tan, J. Yu, and Z. Wu, “VisualInertial-Acoustic Sensor Fusion for Accurate Autonomous Localization of Underwater Vehicles,” IEEE Transactions on Cybernetics, vol. 55, no. 2, pp. 880–896, 2025. [2] T. Ma, F. Hu, X. Su, C. Shen, and X. Ma, “Adaptive Neural Cooperative Control of Multirobot Systems With Input Quantization,” IEEE Transactions on Cybernetics, vol. 54, no. 9, pp. 5518–5528, 2024. [3] Z. Yan, Z. Yang, X. Pan, J. Zhou, and D. Wu, “Virtual leader based path tracking control for Multi-UUV considering sampled-data delays and packet losses,” Ocean Engineering, vol. 216, p. 108065, 2020. [4] M. Zeitz, “The extended Luenberger observer for nonlinear systems,” Systems & Control Letters, vol. 9, no. 2, pp. 149–156, 1987. [5] H. K. Khalil and L. Praly, “High-gain observers in nonlinear feedback control,” International Journal of Robust and Nonlinear Control, vol. 24, no. 6, pp. 993–1015, 2014. [6] S. K. Spurgeon, “Sliding mode observers: a survey,” International Journal of Systems Science, vol. 39, no. 8, pp. 751–764, 2008. [7] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning, Cambridge, Mass: MIT Press, 2006. [8] E. T. Maddalena, P. Scharnhorst, and C. N. Jones, “Deterministic error bounds for kernel-based learning techniques under bounded noise,” Automatica, vol. 134, p. 109896, 2021. [9] M. Bin and L. Marconi, “Model Identification and Adaptive State Observation for a Class of Nonlinear Systems,” IEEE Transactions on Automatic Control, vol. 66, no. 12, pp. 5621–5636, 2021. [10] M. Buisson-Fenet, V. Morgenthaler, S. Trimpe, and F. di Meglio, “Joint state and dynamics estimation with high-gain observers and Gaussian process models,” in 2021 American Control Conference (ACC), pp. 4027–4032, 2021. [11] H. Cai and J. Huang, “The Leader-Following Consensus for Multiple Uncertain Euler-Lagrange Systems With an Adaptive Distributed Observer,” IEEE Transactions on Automatic Control, vol. 61, no. 10, pp. 3152–3157, 2016. [12] W. Yu, Y. Li, G. Wen, X. Yu, and J. Cao, “Observer Design for Tracking Consensus in Second-Order Multi-Agent Systems: Fractional Order Less Than Two,” IEEE Transactions on Automatic Control, vol. 62, no. 2, pp. 894–900, 2017. [13] H. Wang, G. Wen, W. Yu, and X. Yu, “Designing Event-Triggered Observers for Distributed Tracking Consensus of Higher-Order Multiagent Systems,” IEEE Transactions on Cybernetics, vol. 52, no. 5, pp. 3302– 3313, 2022. [14] S. Wang and X. Meng, “Adaptive Consensus and Parameter Estimation of Multiagent Systems With an Uncertain Leader,” IEEE Transactions on Automatic Control, vol. 66, no. 9, pp. 4393–4400, 2021. [15] X. Wang, G.-P. Jiang, H. Su, and Z. Zeng, “Consensus-Based Distributed Reduced-Order Observer Design for LTI Systems,” IEEE Transactions on Cybernetics, vol. 52, no. 7, pp. 6331–6341, 2022. [16] J. Lv, Y. Kao, C. Wang, and F. Wang, “Fully predefined-time distributed observer design for second-order strict-feedback nonlinear systems,” International Journal of Robust and Nonlinear Control, vol. 33, no. 18, pp. 11553–11567, 2023. [17] A. Tahoun and M. Arafa, “Adaptive leader–follower control for nonlinear uncertain multi-agent systems with an uncertain leader and unknown tracking paths,” ISA Transactions, vol. 131, pp. 61–72, 2022. [18] S. Wang, X. Meng, H. Zhang, and F. L. Lewis, “Learning nonlinear dynamics in synchronization of knowledge-based leader-following networks,” Automatica, vol. 166, p. 111695, 2024. [19] C. Wang, J. Lv, Y. kao, and Y. Jiang, “Neural network based adaptive finite-time distributed estimation for an uncertain leader,” Information Sciences, vol. 655, p. 119894, 2024.
[20] M. Deisenroth and J. W. Ng, “Distributed Gaussian Processes,” in Proceedings of the 32nd International Conference on Machine Learning (F. Bach and D. Blei, eds.), vol. 37 of Proceedings of Machine Learning Research, (Lille, France), pp. 1481–1490, PMLR, 07–09 Jul 2015. [21] Z. Yang, S. Sosnowski, Q. Liu, J. Jiao, A. Lederer, and S. Hirche, “Distributed Learning Consensus Control for Unknown Nonlinear MultiAgent Systems based on Gaussian Processes,” in 2021 60th IEEE Conference on Decision and Control (CDC), pp. 4406–4411, 2021. [22] A. Lederer, Z. Yang, J. Jiao, and S. Hirche, “Cooperative Control of Uncertain Multiagent Systems via Distributed Gaussian Processes,” IEEE Transactions on Automatic Control, vol. 68, no. 5, pp. 3091–3098, 2023. [23] Z. Ge, “Structure Learning of Deep Gaussian and Non-Gaussian Information Fusion Framework for Automated Predictive Data Analytics,” IEEE Transactions on Cybernetics, pp. 1–12, 2025. [24] T. N. Hoang, Q. M. Hoang, K. H. Low, and J. How, “Collective Online Learning of Gaussian Processes in Massive Multi-Agent Systems,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 33, pp. 7850–7857, Jul. 2019. [25] X. Dai, Z. Yang, M. Xu, S. Zhang, F. Liu, G. Hattab, and S. Hirche, “Decentralized event-triggered online learning for safe consensus control of multi-agent systems with Gaussian process regression,” European Journal of Control, vol. 80, p. 101058, 2024. 2024 European Control Conference Special Issue. [26] X. Dai, Z. Yang, S. Zhang, D.-H. Zhai, Y. Xia, and S. Hirche, “Cooperative Online Learning for Multiagent System Control via Gaussian Processes With Event-Triggered Mechanism,” IEEE Transactions on Neural Networks and Learning Systems, vol. 36, no. 7, pp. 13304–13318, 2025. [27] Z. Yang, X. Dai, A. Dubey, Y. Wu, X. Chen, E. O. Awotoro, and S. Hirche, “Whom to Trust? Selective Online Learning in Multi-Agent Systems with Prior-Aware Gaussian Process Regression,” Autonomous Agents and Multi-Agent Systems, vol. 40, no. 2, p. 38, 2026. [28] Z. Yang, X. Dai, J. Cheng, Y. Huang, and P. Shi, “Quality or Quantity? Error-Informed Selective Online Learning with Gaussian Processes in Multi-Agent Systems,” IEEE/CAA Journal of Automatica Sinica, vol. 13, no. 6, pp. 1325–1338, 2026. [29] Z. Yang, S. Dong, A. Lederer, X. Dai, S. Chen, S. Sosnowski, G. Hattab, and S. Hirche, “Cooperative Learning with Gaussian Processes for EulerLagrange Systems Tracking Control Under Switching Topologies,” in 2024 American Control Conference (ACC), pp. 560–567, 2024. [30] D. Zhou, W.-H. Chen, and X. Lu, “Leader-Following Consensus of Linear Multiagent Systems With Aperiodically Sampled Outputs: A Distributed Impulsive-Observer-Based Approach,” IEEE Transactions on Cybernetics, vol. 55, no. 1, pp. 161–171, 2025. [31] Y. Mu, H. Zhang, Y. Yan, and X. Xie, “Distributed observer-based robust fault estimation design for discrete-time interconnected systems with disturbances,” IEEE Transactions on Cybernetics, vol. 53, no. 10, pp. 6737–6747, 2023. [32] Y. Yang, H. Modares, K. G. Vamvoudakis, and F. L. Lewis, “Cooperative Finitely Excited Learning for Dynamical Games,” IEEE Transactions on Cybernetics, vol. 54, no. 2, pp. 797–810, 2024. [33] T. Liu, J. Lu, Z. Yan, and G. Zhang, “Robust Gaussian Process Regression With Input Uncertainty: A PAC-Bayes Perspective,” IEEE Transactions on Cybernetics, vol. 54, no. 2, pp. 962–973, 2024. [34] A. I. Lur’e and V. N. Postnikov, “On the theory of stability of control systems,” Applied mathematics and mechanics, vol. 8, no. 3, pp. 246– 248, 1944. [35] D. Xu, Y. Shi, I. W. Tsang, Y.-S. Ong, C. Gong, and X. Shen, “Survey on Multi-Output Learning,” IEEE Transactions on Neural Networks and Learning Systems, vol. 31, no. 7, pp. 2409–2429, 2020. [36] H. K. Khalil, Nonlinear Systems. Prentice-Hall, 2002. [37] K. Hashimoto, A. Saoud, M. Kishida, T. Ushio, and D. V. Dimarogonas, “Learning-based symbolic abstractions for nonlinear control systems,” Automatica, vol. 146, p. 110646, 2022. [38] A. Tokmak, K. G. Krishnan, T. B. Schön, and D. Baumann, “Safe exploration in reproducing kernel Hilbert spaces,” in Proceedings of The 28th International Conference on Artificial Intelligence and Statistics (Y. Li, S. Mandt, S. Agrawal, and E. Khan, eds.), vol. 258 of Proceedings of Machine Learning Research, pp. 784–792, PMLR, 03–05 May 2025. [39] S. Zhang, D.-H. Zhai, Y. Xiong, J. Lin, and Y. Xia, “Safety-critical control for robotic systems with uncertain model via control barrier function,” International Journal of Robust and Nonlinear Control, vol. 33, no. 6, pp. 3661–3676, 2023. [40] C. K. Williams and C. E. Rasmussen, Gaussian Processes for Machine Learning, vol. 2. MIT press Cambridge, MA, 2006.
12
[41] Z. Yang, X. Dai, W. Yang, B. Ilgen, A. Anzel, and G. Hattab, “Kernelbased Learning for Safe Control of Discrete-Time Unknown Systems under Incomplete Observations,” in 2024 43rd Chinese Control Conference (CCC), pp. 6627–6632, 2024. [42] M. Doostmohammadian and T. Charalambous, “Distributed target tracking based on localization with linear time-difference-of-arrival measurements: A delay-tolerant networked estimation approach,” Systems & Control Letters, vol. 196, p. 106009, 2025. [43] Z. Yang, X. Dai, and S. Hirche, “Asynchronous Distributed Gaussian Process Regression,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 39, pp. 22065–22073, Apr. 2025. [44] J. Umlauft and S. Hirche, “Feedback Linearization Based on Gaussian Processes With Event-Triggered Online Learning,” IEEE Transactions on Automatic Control, vol. 65, no. 10, pp. 4154–4169, 2020. [45] J. Kautsky, N. K. Nichols, and P. Van Dooren, “Robust Pole Assignment in Linear State Feedback,” International Journal of control, vol. 41, no. 5, pp. 1129–1155, 1985. [46] D. Nguyen-Tuong, J. Peters, and M. Seeger, “Local Gaussian Process Regression for Real Time Online Model Learning,” Advances in neural information processing systems, vol. 21, 2008. [47] V. Tresp, “Mixtures of Gaussian processes,” Advances in neural information processing systems, vol. 13, 2000. [48] V. Tresp, “A Bayesian committee machine,” Neural computation, vol. 12, no. 11, pp. 2719–2741, 2000. [49] A. Lederer, A. J. O. Conejo, K. A. Maier, W. Xiao, J. Umlauft, and S. Hirche, “Gaussian Process-Based Real-Time Learning for Safety Critical Applications,” in Proceedings of the 38th International Conference on Machine Learning, vol. 139 of Proceedings of Machine Learning Research, pp. 6055–6064, PMLR, 18–24 Jul 2021. [50] Z. Yang, D. Zhang, X. Dai, F. Yu, C. Zhang, B. Huang, H. Sadeghian, and S. Haddadin, “Streaming Generated Gaussian Process Experts for Online Learning and Control,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 40, pp. 27719–27727, Mar. 2026.