Identification of non-causal systems with arbitrary switching modes
Yanxin Zhang Chengpu Yu Filippo Fabiani




                                              
                                              
                                              

 
                                                





Bit]Yanxin Zhangzhangyanxin@bit.edu.cn,   Bit]Chengpu Yuyuchengpu@bit.edu.cn, and MIT]Filippo Fabianifilippo.fabiani@imtlucca.it

[Bit]School of Automation, Beijing Institute of Technology, Beijing 100081, PR China  
[MIT]IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100, Lucca, Italy                                           
This work was supported by the National Natural Science Foundation of China (Grant No. 61991414, 62088101, 6193000461), Chongqing Natural Science Foundation CSTB
2023NSCQ-JQX0018, and Beijing Natural Science Foundation L221005. Corresponding author: Chengpu Yu.

                           
Switching systems; Non-causal systems; Expectation maximization; Kalman filter              
                             
                                          
                                          


§ ABSTRACT                           

	
	We consider the identification of non-causal systems with arbitrary switching modes (NCS-ASM), a class of models essential for describing typical power load management and department store inventory dynamics. The simultaneous identification of causal-and-anticausal subsystems, along with the presence of possibly random switching sequences, however, make the overall identification problem particularly challenging. To this end, we develop an expectation-maximization (EM) based system identification technique, where the E-step proposes a modified Kalman filter (KF) to estimate the states and switching sequences of causal-and-anticausal subsystems, while the M-step consists in a switching least-squares algorithm to estimate the parameters of individual subsystems. We establish the main convergence features of the proposed identification procedure, also providing bounds on the parameter estimation errors under mild conditions. Finally, the effectiveness of our identification method is validated through two numerical simulations.







§ INTRODUCTION

Non-causal switching dynamics denote a class of systems able to modeling several real-world scenarios, such as load management in power systems <cit.>, traffic signal control systems <cit.>, and robotic systems <cit.>. These systems encompass non-causality, meaning that their output signals depend not only on the current or past control actions, but also on future inputs. In addition, the systems exhibit switching characteristics, potentially transitioning among different operational states, and thus leading to variations in the system behavior. Understanding and managing the complexity of these systems is therefore crucial for enhancing efficiency, reliability, and adaptability, enabling them to better meet the demands of industrial production and operations. This essentially motivates the interest in modeling, analyzing, and controlling such type of systems.

In many parameter identification problems for dynamical systems, the input-output data are accompanied by temporal mode sequences. As the mode of the system changes over time, each data point corresponds to the mode at its corresponding time instant. In these circumstances, it is hence crucial to model the dynamics of different modes and estimate the transitions from one mode to another <cit.>. However, obtaining direct estimates of the dynamical system from input-output data is challenging and, in practice, the prior knowledge on the mode transitions is often unavailable. Therefore, estimating the switching behaviors poses a challenging, yet highly significant, problem addressed by several researchers. Available works indeed propose algorithms to estimate the individual system dynamics and the mode transition sequence based on observed behaviors <cit.>.



 §.§ Literature review

Several works consider the identification of switching models <cit.>. In <cit.>, a joint smoothing algorithm is proposed based on the expectation-maximization (EM) framework, where an E-step solution is introduced to effectively address issues related to the exponential complexity in the jump Markov linear model. In <cit.>, a numerically efficient, two-step estimation method was proposed, which iteratively updates the parameters and the switching sequence. The flexibility of this technique consists in its adaptability to different loss functions employed in the jump models, which significantly impact the overall shape and jumping behavior of the models. Furthermore, the identification of jump Box-Jenkins systems is investigated in <cit.>, where a maximum a-posteriori method is proposed to estimate the switching sequence of the model. Subsequently, the system parameters of the jump Box-Jenkins models are alternately estimated using the Gauss-Newton and the prediction error methods. In <cit.>, a switching least-squares algorithm for autonomous Markov jump linear systems is proposed. Here, the authors provided a formal proof characterizing the strong consistency of the underlying method, as well as established its convergence rate as 𝒪(√(log(T)/T)) a.s., where T is the time horizon. The aforementioned literature primarily focuses on linear systems and assumes that the switching behavior follows a Markov chain. However, these techniques are not applicable when the mode switching is random. Consequently, a number of approaches for identifying systems with random mode switching behavior have been proposed in the literature. For example, in <cit.> a kernel-based method is used to estimate the random switching system, which can solve both estimation and classification problems simultaneously. In <cit.>, a maximum-likelihood algorithm is presented for the identification of a random switching linear system, which combines the maximum-likelihood estimation criterion and the Kalman filtering technique to estimate the system modes in switched linear systems, resulting in a significant contribution to the estimation error stability of general switched linear systems.

All the studies mentioned above are focused on switching causal systems, where the system output is only related to current or past excitation. However, to the best of our knowledge, there is no literature on the system identification problem of non-causal systems with arbitrary switching modes (NCS-ASM). Nevertheless, NCS-ASM are widely present in real-world scenarios. For instance, in automatic control systems the time delay caused by sensors can lead to the switching non-causal characteristics of the system <cit.>. In financial markets, there exist interaction and feedback between investor decisions and market prices <cit.>. Although, there are some studies available on system identification for non-causal systems, such as the subspace <cit.> and the kernel methods <cit.>, these studies can only handle a single, non-causal system, rather than switching non-causal systems.



 §.§ Summary of contribution



In this paper, we focus on the identification of NCS-ASM. The proposed method is developed under the expectation-maximization (EM) framework, which can be divided into two main parts. Specifically, in the E-step we adopt a Bayesian rule to compute the posterior estimate of the switching sequence, along with a modified Kalman filter (KF) for estimating the state of the causal and anti-causal parts. In the M-step, instead, we propose a switching least-squares method to obtain the closed-form solution for the parameters and establish the convergence rate of the estimated parameters. Our main contributions can hence be summarized as follows:

	
  * To the best of our knowledge, this is the first work considering the idenfitication of NCS-ASM. In particular, in the E-step of the EM framework a modified KF is proposed to compute the posterior state estimates of the causal and anti-causal parts, which is crucial to deal with the identification of non-causal systems;
	
  * Compared with the system dynamics in <cit.>, where the switching behavior of the subsystems only occurs in the causal part, our methodology can cope with switching behaviors in both causal and anti-causal parts. Moreover, the switching sequences of the two directions is allowed to differ from each other. 




 §.§ Paper organization

The rest of the paper is organized as follows:

in Section <ref> we describe the considered system and formulate the related identification problem. In Section <ref>, instead, we discuss our EM method for the identification of the NCS-ASM, while in Section <ref> we provide its implementation details, as well as characterize the related convergence properties. Two simulation examples are finally discussed in Section <ref> to test the effectiveness of the proposed method numerically. 

The proofs of the technical results of the paper are all deferred to Appendix <ref>.


Notations:
ℤ and ℝ denote the set of integer and real numbers, respectively. 
Given a matrix X, ‖ X‖ and ‖ X‖_∞ denote respectively its spectral and infinity norms, λ_max(X) and λ_min(X) are respectively its maximum and minimum eigenvalues, and tr(X) denotes the trace. ℙ[·] and 𝔼[·] denote respectively a distribution probability and the related expected value. 𝕊^n is the space of n × n symmetric matrices and 𝕊_≻ 0^n (𝕊_≽ 0^n) is the cone of positive (semi-)definite matrices. Given two square matrices A, B of compatible dimension, A≽ B means that A-B is positive semidefinite. For a sequence {s_t}_t∈ℕ, s_T=𝒪(T) indicates that limsup_T→∞s_T/T∞, while s_T=o(T) that limsup_T→∞s_T/T=0. Finally, I identifies a standard identity matrix. 𝒩(μ,σ^2) denotes the normal distribution of a random variable with mean μ and standard deviation σ.
 



§ MATHEMATICAL FORMULATION

We now describe the system concerned in this paper, together with the main assumptions, and successively formalize the problem to be addressed.


 §.§ System model description

Consider the following discrete-time, non-causal system characterized by arbitrary switching modes:

	x_c(t)   =   A_c(s_c(t))x_c(t-1)+v_c(t),

	x_a(t)   =   A_a(s_a(t))x_a(t+1)+v_a(t),

	y(t)   =   C_c(s_c(t))x_c(t)+C_a(s_a(t))x_a(t)+v_m(t),


where t∈ℤ is the time instant, x_c(t)∈ℝ^n_x_c,x_a(t)∈ℝ^n_x_a are the causal and anti-causal state vectors, respectively, y(t)∈ℝ^n_y denotes the system output, while s_c(t)∈{1,2,…,m_c}≜Λ_c and s_a(t)∈{1,2,…,m_a}≜Λ_a are two discrete variables representing the possible switching modes. In addition, v_c(t)∈ℝ^n_x_c and v_a(t)∈ℝ^n_x_a are the system noise vectors, and v_m(t)∈ℝ^n_y is the measurement noise vector. Finally, A_c:Λ_c→ℝ^n_x_c× n_x_c and A_a:Λ_a→ℝ^n_x_a× n_x_a denote the matrix functions associated to the causal and anti-causal state dynamics, respectively, while C_c:Λ_c→ℝ^n_y× n_x_c and C_a:Λ_a→ℝ^n_y× n_x_a are those mapping the two state vectors to the measured output. 

Assume that the noise terms v_c(t), v_a(t) and v_m(t) are distributed according to a Gaussian distribution with zero mean and finite variance as follows:

    v_c(t)   ∼𝒩(0,Σ_c(s_c(t))),
    
    	v_a(t)   ∼𝒩(0,Σ_a(s_a(t))),
    
    	v_m(t)   ∼𝒩(0,Σ_m).









	The NCS-ASM (<ref>) is stable in the average sense, which means
	
    ∑_i=1^T‖ x_c(i)‖^2=𝒪(T),  ∑_i=1^T‖ x_a(i)‖^2=𝒪(T).

	where T is the sample size of the available dataset.



	Stability in the average sense is widely applied in linear systems <cit.>. Note that, compared to other commonly used notions, such as mean-square stability and almost sure stability, the assumption of stable in the average sense is weaker.


	The noise vectors v_c(t),v_a(t),v_m(t) are the martingale difference sequences with respect to an increasing sequence of σ-fields 𝒢_T generated by the history states, and satisfy the following conditions:
	
		
		
	
    liminf_N→∞1/T∑_i=1^Tv_c(i)v_c(i)^⊤=𝒞_1,
       liminf_N→∞1/T∑_i=1^Tv_a(i)v_a(i)^⊤=𝒞_2,
       liminf_N→∞1/T∑_i=1^Tv_m(i)v_m(i)^⊤=𝒞_3,

	where 𝒞_1∈𝕊_≻0^n_x_c, 𝒞_2∈𝕊_≻0^n_x_a, and 𝒞_3∈𝕊_≻0^n_y are also bounded.


	Standing Assumption <ref> denotes a common requirement for analyzing the convergence of system identification algorithms, enabling the noise process to exhibit non-stationary and heavy-tailed characteristics – see, e.g., <cit.>.


The NCS-ASM in (<ref>) thus consists of two state equations and one output equation. Specifically, the first state equation represents the dynamics of the causal state variables, while the second one the dynamics of the non-causal state variables. The system output is determined by both the causal and non-causal states. Furthermore, both the causal and non-causal parts of the system are composed of multiple subsystems, and their corresponding switching sequences are different. Given some T∈ℤ, which will denote the sample size of the available dataset, let the switching sequences of the causal and anti-causal parts being denoted by s_c≜{s_c(i)}_i=1^T and s_a≜{s_a(i)}_i=1^T, respectively. Each of them corresponds to a set of parameters, i.e., s_c(t)=i determines the model parameter θ^c_i≜{A_c(i),C_c(i),Σ_c(i)} that is active at the time instant t. In particular, the sequences s_c and s_a undergo random switches with certain (fixed) probabilities over time. Let us denote the underlying switching probabilities with ℙ[s_c=i]=π^c_i and ℙ[s_a=i]=π^a_i. We then have that ∑_i=1^m_cπ^c_i=1 and ∑_i=1^m_aπ^a_i=1. 

The complete set of model parameters that comprehensively describe the NCS-ASM can be conveniently encapsulated into a parameter object θ, defined as follows: 







    θ≜{{θ^c_i}_i=1^m_c,{θ^a_i}_i=1^m_a,{π^c_i}_i=1^m_c,{π^a_i}_i=1^m_a,Σ_m}.




 §.§ Problem statement

Our goal is hence to estimate the unknown model parameters θ characterizing the NCS-ASM (<ref>) with the known state dimension, number of causal system modes m_c and anti-causal system modes m_a, together with a collection of noisy output measurements y:

    y≜y_1:T={y(1),…,y(T)}.


It is worth remarking that the NCS-ASM consists of both causal and non-causal parts, and their switching sequences are different. Therefore, the problem we wish to addressing has two main challenges. First, it is difficult to obtain the parameters for each subsystem of the NCS-ASM, since the states of both parts are unknown, and the system output is determined by both causal and non-causal states. Identifying both parts simultaneously is hence challenging, especially due to the continuous switching pattern of different subsystems. Second, the system has two switching sequences, and the combination of subsystems varies at different time instants. For example, at time t the causal part may activate subsystem i (s_c(t)=i) while the non-causal part may activate the subsystem j (s_a(t)=j), thereby originating a total of m_a× m_c possible combinations. In addition, the switching  behavior of the subsystems is random and independent across different time instants:

    ℙ[s_c(t)| s_c(t-1),…,s_c(1)]=ℙ[s_c(t)],
       ℙ[s_a(t)| s_a(t+1),…,s_a(T)]=ℙ[s_a(t)],    t=1,…,T.


To deal with the identification problem of the NCS-ASM (<ref>), the EM framework is adopted, which is an iterative method that can yield an estimate of the parameters at each iteration <cit.>. Let us denote the parameter estimate at the k-th iteration of the underlying algorithm as θ^k. Then, the proposed method can be (qualitatively, for the moment) described by means of the following two steps:

	
  * In the E-step, we develop a modified KF to estimate the states of the causal and anti-causal parts. Furthermore, the Bayesian rule is used to obtain a posterior estimate of the switching sequence. Subsequently, the full-data likelihood function Q(θ,θ^k) can be calculated.
	
  * In the M-step, the likelihood function Q(θ,θ^k) is maximized with respect to the parameters θ. Then, the identification of the NCS-ASM is updated, yielding θ^k+1.


Next section will discuss in detail each step of the proposed technique for NCS-ASM identification.



§ THE EM METHOD FOR IDENTIFYING NCS-ASM

By making use of the dataset y, we aim at estimating the system parameters θ. To this end, a standard approach is to let coincide θ̂, i.e., our estimate of the true θ, with a maximizer of the likelihood function, namely:

    θ̂=θmax lnℙ_θ(y).

where we indicate with ℙ_θ(y) the probability density function of the output y given some sets of parameters θ. In the remainder, we will tacitly assume that the solution to the likelihood maximization problem (<ref>) is unique.

Given any collection of data y, note that the likelihood function lnℙ_θ(y), also called marginal density function of y, can be decomposed into the following form:












    lnℙ_θ(y)
       =lnℙ_θ[y(1)]+∑_i=2^T-1lnℙ_θ[y(i)]+lnℙ_θ[y(T)]
       =ln∫∫∑_s_c(1)∑_s_a(1)ℙ_θ[y(1)| x_c(1),x_a(1),s_c(1),s_a(1)]
       ℙ_θ[x_c(1),x_a(1),s_c(1),s_a(1)] dx_c(1)dx_a(1)
       +ln∫∫∑_s_c(T)∑_s_a(T)ℙ_θ[y(T)| x_c(T),x_a(T),s_c(T),s_a(T)]
       ℙ_θ[x_c(T),x_a(T),s_c(T),s_a(T)] dx_c(T)dx_a(T)
       +∑_i=2^T-1ln∫∫∑_s_c(i)∑_s_a(i)ℙ_θ[y(i)| x_c(i),x_a(i),s_c(i),s_a(i)]
       ℙ_θ[x_c(i),s_c(i)|y_1:i-1)p_θ(x_a(i),s_a(i)|y_i+1:T] dx_c(i)dx_a(i).


Let us denote the collection of state variables over T as x_c≜{x_c(i)}_i=1^T and x_a≜{x_a(i)}_i=1^T. Recall that in the NCS-ASM (<ref>) the state variables x_c and x_a are determined by the switching sequences s_c and s_a. Besides the potential nonconvexity of lnℙ_θ(y), which makes the direct maximization of lnℙ_θ(y) challenging (together with its high-dimensionality), from the decomposition above it is also clear that for calculating ℙ_θ(y) we need to sum over all possible values of s_a, s_c, thereby further complicating the solution of (<ref>). 


Another way to marginalize the latent variables (such as x_c, x_a, s_c, s_a) is by taking the expectation over these latter. Instead of maximizing the incomplete likelihood function lnℙ_θ(y), we can estimate the conditional density of the hidden variables given the observations y and an estimate of parameter θ̂. Then, parameter estimate θ̂ can be obtained by maximizing the complete likelihood function.

To stand out our technical contributions and contrast them with existing results, we will give a sample complexity analysis related to our EM-based technique for the identification of NCS-ASM. To achieve this, we have to further postulate the following:


	The following conditions hold true:
	
	
  * The switching sequences s_c, s_a, and the subsystem parameters θ^c, θ^a are all independent among them, i.e.,
	
    ℙ[s_c|θ^c]   =ℙ[s_c],  ℙ[θ^c|s_c]=ℙ[θ^c],
    ℙ[s_a|θ^a]   =ℙ[s_a],  ℙ[θ^a|s_a]=ℙ[θ^a].

	
  * The switching sequence follows a polynomial distribution, i.e.,
	
    ℙ[s_c=i]=π^c_i,    i=1,…,m_c,
    ℙ[s_a=i]=π^a_i,    i=1,…,m_a,

	with ∑_i=1^m_cπ^c_i=1, ∑_i=1^m_aπ^a_i=1.
	






Then, the full-data complete likelihood function can be expressed as follows:

    lnℙ_θ[y,x_c,s_c,x_a,s_a]=    lnℙ_θ[y]
       +lnℙ_θ[x_c,s_c,x_a,s_a|y].

This relation directly links ℙ_θ(y) and ℙ_θ[y,x_c,s_c,x_a,s_a], with the latter depending on the unknown states x_c, x_a and switching sequences s_c, s_a. The key step is then to approximate lnℙ_θ[y] by the above relation (<ref>), where x_c, s_c, x_a, and s_a can be approximated by their conditional expectations based on the observed data y. 
Therefore, at each iteration k of our EM-based algorithm, given the estimate θ^k the conditional expectation of lnℙ_θ^k[y,x_c,s_c,x_a,s_a] can be obtained based on the available data y as follows:

    𝔼_θ^k[lnℙ_θ[y,x_c,s_c,x_a,s_a]]
               =𝔼_θ^k[lnℙ_θ(y)]+𝔼_θ^k[lnℙ_θ[x_c,s_c,x_a,s_a|y]]
               =lnℙ_θ(y)+𝔼_θ^k[lnℙ_θ[x_c,s_c,x_a,s_a|y]].

Let us then define:

    Q(θ,θ^k)   =𝔼_θ^k[lnℙ_θ[y,x_c,s_c,x_a,s_a]],
    
    	V(θ,θ^k)   =𝔼_θ^k[lnℙ_θ[x_c,s_c,x_a,s_a|y]].

The EM approach iteratively estimates the parameters in the following two steps. First, we compute the expectation Q(θ,θ^k) based on θ^k obtained from the previous iteration. By the Bayesian rule and Markov properties, it can be inferred that:

    lnℙ_θ[y,x_c,s_c,x_a,s_a]
       =∑_i=1^Tlnℙ_θ[y(i)| x_c(i),x_a(i),s_c(i),s_a(i)]
       +lnℙ_θ[x_c(1),s_c(1)]+∑_i=2^Tlnℙ_θ[x_c(i)| x_c(i-1),s_c(i)]
       +lnℙ_θ[x_a(T),s_a(T)]+∑_i=1^T-1lnℙ_θ[x_a(i)| x_a(i+1),s_a(i)]
       =∑_i=1^T∑_j=1^m_c∑_l=1^m_alnℙ_θ[y(i)| x_c(i),x_a(i)]π_j^cπ_l^a
       +∑_j=1^m_clnℙ_θ[x_c(1)]π_j^c+∑_i=2^T∑_j=1^m_clnℙ_θ[x_c(i)| x_c(i-1)]π_j^c
       +∑_l=1^m_alnℙ_θ[x_a(T)]π_l^a+∑_i=1^T-1∑_l=1^m_alnℙ_θ[x_a(i)| x_a(i+1)]π_l^a.

In view of the white noise assumption characterizing the disturbance affecting both state variables and measured output, note that the distribution of the these variables, given the subsystem modes s_c(i)=j, s_a(i)=l, is Gaussian too and given as follows:

    ℙ_θ[y(i)|   x_c(i),x_a(i)]=| 2πΣ_m|^-1/2
       exp{(y(i)-μ_1(i))^⊤Σ_m^-1(y(i)-μ_1(i))},
    ℙ_θ[x_c(i)|   x_c(i-1)]=| 2πΣ_c(i)|^-1/2
       exp{(x_c(i)-μ_2(i))^⊤Σ_c^-1(i)(x_c(i)-μ_2(i))},
    ℙ_θ[x_a(i)|   x_a(i+1)]=| 2πΣ_a(i)|^-1/2
       exp{(x_a(i)-μ_3(i))^⊤Σ_a^-1(i)(x_a(i)-μ_3(i))},
 
where

    μ_1(i)   =y(i)-C_c(j)x_c(i)-C_a(l)x_a(i),
    μ_2(i)   =x_c(i)-A_c(j)x_c(i-1),
    μ_3(i)   =x_a(i)-A_a(l)x_a(i+1).

Let us indicate with w_ij^c the posterior probability of the switching sequence given that s_c(i)=j (w_il^a is defined similarly). Then, the objective function Q(θ,θ^k) assumes the following form:

    𝔼_θ^k[lnℙ_θ[y,x_c,s_c,x_a,s_a]|y]
    	=∑_i=1^5 Q_i(θ,θ^k),

 

		
		
		
		
		
		
		
		
		
where the closed form for the terms Q_i is in (<ref>). Note that the posterior densities w^c_ij and w^a_il, given the parameter θ^k and dataset y, can be computed as follows:

    w^c_ij= ℙ_θ^k[s_c(i)   =j|y]=ℙ_θ_i^k[y,s_c(i)=j]/ℙ_θ_i^k[y]
       =ℙ_θ_i^k[y| s_c(i)=j]π^c_j/∑_i=1^Tℙ_θ_i^k[y| s_c(i)=j]π^c_j,


    w^a_il= ℙ_θ^k[s_a(i)   =l|y]=ℙ_θ_i^k[y,s_a(i)=l]/ℙ_θ_i^k[y]
       =ℙ_θ_i^k[y| s_a(i)=l]π^a_l/∑_i=1^Tℙ_θ_i^k[y| s_a(i)=l]π^c_l.


Then, the second step is to maximize the Q(θ,θ^k) to obtain θ^k+1, formally defined as θ^k+1=max_θ Q(θ,θ^k).

Algorithm <ref> summarizes the two main steps of the proposed identification methodology for NCS-ASM. We characterize next the monotonic properties of the likelihood function in (<ref>) when the EM algorithm is iteratively applied to estimate the system parameters θ:

Given a dataset y, let {θ^k}_k∈ℤ be the sequence generated by Algorithm <ref>. Then, the likelihood function in (<ref>), evaluated along {θ^k}_k∈ℤ, is non-decreasing, thereby yielding lnℙ_θ^k+1[y]≥lnℙ_θ^k[y] for all k∈ℤ.






§ IMPLEMENTATION DETAILS OF THE EM ALGORITHM

We now delve into the details of the steps outlined in Algorithm <ref>, ultimately establishing our main technical result characterizing the sample complexity of the proposed identification technique for NCS-ASM.



 §.§ The E-step

This step requires the calculation of the objective function Q(θ,θ^k). Specifically, this shall be achieved on the basis of the parameter θ^k estimated in the previous iteration. Then, according to the expression of Q(θ,θ^k) in (<ref>), the expectations of states x_c,x_a and the switching sequences s_c,s_a given the data y are required. 

First, we calculate the posterior estimates of the switching sequences s_c and s_a by leveraging the Bayesian rule, namely ℙ_θ[s_c|y]=ℙ_θ[s_c,y]/ℙ_θ(y) and ℙ_θ[s_a|y]=ℙ_θ[s_a,y]/ℙ_θ(y). In addition, according to the formula of total probability one obtains:

    ℙ_θ(y)   =∑_j=1^m_cℙ_θ[y|s_c=j] ℙ[s_c=j]=∑_j=1^m_cℙ_θ[y|s_c=j]π_j^c,
    ℙ_θ(y)   =∑_l=1^m_aℙ_θ[y|s_a=l] ℙ[s_a=l]=∑_l=1^m_aℙ_θ[y|s_a=l]π_l^a.

 Then, the data point can be assigned to each subsystem at time i by solving the following optimization problem:

    ŝ_c(i)   =j∈{1,…,m_c}max ℙ_θ[y(i)| s_c(i)=j]π_j^c,
    ŝ_a(i)   =l∈{1,…,m_a}max ℙ_θ[y(i)| s_a(i)=l]π_l^a,

where maximizing ℙ_θ[y(i)| s_c(i)=j]π_j^c is equivalent to maximizing the posterior probability of ℙ_θ[s_c(i)=j| y(i)] which is commonly used for data classification.
After obtaining ŝ_c(i) and ŝ_a(i), for any (i,j)∈{1,…,T}×{1,…,m_c} (or (i,l)∈{1,…,T}×{1,…,m_a}), we obtained that 

    w_ij^c= {[           1 if ŝ_c(i)=j;           0        else ].,  
    	w_il^a= {[           1 if ŝ_a(i)=l;           0        else ].


Successively, we focus on the reconstruction of the state variables x_c and x_a, a task that is traditionally accomplished by means of a Kalman filter. The latter amounts to a recursive algorithm that estimates the system states by incorporating the information inferred from previous estimates and current measurements. 
Adapting the KF to our problem, however, requires few key modifications due to the following reason.
According to the system description in (<ref>), we note that the calculation of the posterior probability ℙ_θ[x_c(t)|y] is related to two equations. The first one requires to propagate the estimate from the previous time step, ℙ_θ[x_c(t-1)|y], to the current time step recursively, whereas the second one involves correcting the prior estimate with the measurement equation in a reverse manner (for the non-causal states happens exactly the opposite). Then, when correcting the prior prediction of the state variables x_c and x_a using the data y, we note that x_c and x_a are mutually correlated, thereby calling for a careful design of the KF as described below. To simplify notation we omit the dependency on the switching sequence, e.g., A_c=A_c(ŝ_c(t))):


First, we need to compute the prior state estimates of x_c and x_a, denoted as x̂_c^- and  x̂_a^-, through the first two relations in (<ref>), yielding x̂_c^-(t)=A_cx̂_c(t-1) and x̂_a^-(t)=A_ax̂_a(t+1). With this regard, note that the switching sequence for each step has already been calculated. Successively, the measurement equation in (<ref>) allows us to perform posterior corrections x̂_c and x̂_a on the underlying prior estimates x̂_c^- and  x̂_a^- as follows:

    x̂_c(t)   =x̂_c^-(t)+K_cC_c(C_c^-1(y(t)-C_ax̂_a^-(t))-C_cx̂_c^-(t))
       =x̂_c^-(t)+K_c(y(t)-C_ax̂_a^-(t)-C_cx̂_c^-(t)),
    x̂_a(t)   =x̂_a^-(t)+K_aC_a(C_a^-1(y(t)-C_ax̂_a^-(t))-C_cx̂_c^-(t))
       =x̂_a^-(t)+K_a(y(t)-C_ax̂_a^-(t)-C_cx̂_c^-(t)),

where K_c∈ℝ^n_x_c× n_y and K_a∈ℝ^n_x_a× n_y are the Kalman gains for the causal and anti-causal states, respectively, whose design is critical for the effectiveness of the KF.
Before delving into the derivation of K_c and K_a, let us first calculate the error covariance matrix for the prior state estimates based on the prior estimation errors e_c^-(t)=x_c(t)-x̂^-_c(t) and e_a^-(t)=x_a(t)-x̂^-_a(t), and the posterior estimation errors e_c(t)=x_c(t)-x̂_c(t) and e_a(t)=x_a(t)-x̂_a(t). Then, the prior error covariance and posterior estimation error covariance matrices can be defined as follows:

    P_c^-   ≜𝔼[(x_c(t)-x̂_c^-(t))(x_c(t)-x̂_c^-(t))^⊤],
    
    	P_a^-   ≜𝔼[(x_a(t)-x̂_a^-(t))(x_a(t)-x̂_a^-(t))^⊤],
    
    	P_c   ≜𝔼[(x_c(t)-x̂_c(t))(x_c(t)-x̂_c(t))^⊤],
    
    	P_a   ≜𝔼[(x_a(t)-x̂_a(t))(x_a(t)-x̂_a(t))^⊤].

Then, we have:

    P_c^-   =𝔼[(x_c(t)-x̂_c^-(t))(x_c(t)-x̂_c^-(t))^⊤]
       =𝔼[(A_c(x_c(t-1)-x̂_c(t-1))+v_c(t))]
       (A_c(x_c(t-1)-x̂_c(t-1))+v_c(t))^⊤]
       =A_cP_c(t-1)A_c^⊤+Σ_c,
    	P_a^-   =𝔼[(x_a(t)-x̂_a^-(t))(x_a(t)-x̂_a^-(t))^⊤]
       =𝔼[(A_a(x_a(t+1)-x̂_a(t+1))+v_a(t))
       (A_a(x_a(t+1)-x̂_a(t+1))+v_a(t))^⊤]
       =A_aP_a(t+1)A_a^⊤+Σ_a.

Then, the Kalman gains can be calculated so that the error covariance matrices for the posterior state estimates is minimized. The posterior estimation error can be rewritten as:

    e_c(t)   =x_c(t)-x̂_c(t)
       =x_c(t)-x̂_c^-(t)+K_c(y(t)-C_ax̂_a^-(t)-C_cx̂_c^-(t))
       =(I-K_cC_c)e_c^-(t)-K_cC_ae_a^-(t)-K_cv_m(t),
    
    e_a(t)   =x_a(t)-x̂_a(t)
       =x_a(t)-x̂_a^-(t)+K_a(y(t)-C_ax̂_a^-(t)-C_cx̂_c^-(t))
       =(I-K_aC_a)e_a^-(t)-K_aC_ce_c^-(t)-K_av_m(t),

while the error covariance matrices of the state estimates:

    P_c   =𝔼[((I-K_cC_c)e_c^-(t)-K_cC_ae_a^-(t)-K_cv_m(t))
       ((I-K_cC_c)e_c^-(t)-K_cC_ae_a^-(t)-K_cv_m(t))^⊤]
       =(I-K_cC_c)P_c^-(I-K_cC_c)^⊤+K_cC_aP_a^-C_a^⊤ K_c^⊤
       +K_cΣ_mK_c^⊤
       =P_c^–P_c^-C_c^⊤ K_c^⊤-K_cC_cP_c^-+K_cC_cP_c^-C_c^⊤ K_c^⊤
       +K_cC_aP_a^-C_a^⊤ K_c^⊤+K_cΣ_mK_c^⊤, 
    
    	P_a   =𝔼[((I-K_aC_a)e_a^-(t)-K_aC_ce_c^-(t)-K_av_m(t))
       ((I-K_aC_a)e_a^-(t)-K_aC_ce_c^-(t)-K_av_m(t))^⊤]
       =(I-K_aC_a)P_a^-(I-K_aC_a)^⊤+K_aC_cP_c^-C_c^⊤ K_a^⊤
       +K_aΣ_mK_a^⊤
       =P_a^–P_a^-C_a^⊤ K_a^⊤-K_aC_aP_a^-+K_aC_aP_a^-C_a^⊤ K_a^⊤
       +K_aC_cP_c^-C_c^⊤ K_a^⊤+K_aΣ_mK_a^⊤,

where the second equality in each derivation is established based on the independence of e_c^-(t), e_a^-(t) and v_m(t). Note that minimizing the variances of P_c and P_a is equivalent to minimizing their traces. Therefore, given the unconstrained nature of such trace minimization, the optimal Kalman gains K_c and K_a can be found as:

    d(tr(P_c))/dK_c   =0-2[C_cP_c^-]^⊤+2K_cC_cP_c^-C_c^⊤
       +2K_cC_aP_a^-C_a^⊤+2K_cΣ_m=0,
    ⇒ K_c   =(C_cP_c^-C_c^⊤+C_aP_a^-C_a^⊤+Σ_m)^-1(P_c^-C_c^⊤),
    d(tr(P_a))/dK_a   =0-2[C_aP_a^-]^⊤+2K_aC_aP_a^-C_a^⊤
       +2K_aC_cP_c^-C_c^⊤+2K_aΣ_m=0,
    ⇒ K_a   =(C_aP_a^-C_a^⊤+C_cP_c^-C_c^⊤+Σ_m)^-1(P_a^-C_a^⊤).


By substituting the Kalman gains above into (<ref>)–(<ref>), the updated error covariance matrices can be obtained as:

    P_c   =P_c^–P_c^-C_c^⊤ K_c^⊤-K_cC_cP_c^-+P_c^-C_c^⊤ K_c^⊤
       =(I-K_cC_c)P_c^-,
    
    	P_a   =P_a^–P_a^-C_a^⊤ K_a^⊤-K_aC_aP_a^-+P_a^-C_a^⊤ K_a^⊤
       =(I-K_aC_a)P_a^-.

By completing the steps of the modified KF, including the prediction, measurement update, and error covariance matrix update <cit.>, one can obtain all the posterior estimates of the state vectors x_c and x_a, which are optimal state estimates based on the available measurements and prior knowledge. In addition, to ensure the convergence of the proposed state estimation method, we establish the following properties of the state estimates:


	Let η_c(t)=x_c(t)-A(ŝ_c(t))x_c(t-1), η_a(t)=x_a(t)-A(ŝ_a(t))x_a(t+1), and η_m(t)=y(t)-C_c(ŝ_c(t))x_c(t)-C_a(ŝ_a(t))x_a(t). There exist α_1, α_2, α_3 > 0 so that η_c(t)^2≤α_1, η_a(t)^2≤α_2, and η_m(t)^2≤α_3, for all t∈ℤ.


Lemma <ref> states that the error of state estimation is bounded in the mean square sense, regardless of how the state trajectory evolves in time.
 
By making use of the steps described in this subsection one is able to find an expression for the objective function Q(θ,θ^k). We will discuss next how to actually maximize it with respect to its first argument.



 §.§ The M-step

The second step in Algorithm <ref> requires the maximization of Q(θ,θ^k) to update the parameters estimate θ^k:

    θ^k+1=θmax Q(θ,θ^k).

Let us first focus on the elements {{π^c_i}_i=1^m_c,{π^a_i}_i=1^m_a}, and recall that the objective function reads as:

    Q(θ,θ^k)
       =𝔼_θ^k[lnℙ_θ[y,x_c,s_c,x_a,s_a]|y],
       =∫∫∫∫lnℙ_θ[y,x_c,s_c,x_a,s_a] ℙ_θ^k[x_c|y] ℙ_θ^k[x_a|y]
       ℙ_θ^k[s_c|y] ℙ_θ^k[s_a|y] d(x_c)d(x_a)d(s_c)d(s_a)
       =∑_i=1^T∑_j=1^m_c∑_l=1^m_a∫∫ w^c_ijw^a_illnℙ_θ[y(i)| x_c(i),x_a(i)]
       ℙ_θ^k[x_c(i)|y] ℙ_θ^k[x_a(i)|y]π^c_jπ^a_l d(x_c(i))d(x_a(i)).

The (k+1)-th estimate of {{π^c_i}_i=1^m_c,{π^a_i}_i=1^m_a} can hence be obtained in closed-form by applying the first-order optimality conditions as follows:

    π^c_j   =π^c_jmax Q(θ,θ^k)=∑_i=1^Tw^c_ij/∑_i=1^T∑_j=1^m_cw^c_ij,
    π^a_l   =π^a_lmax Q(θ,θ^k)=∑_i=1^Tw^a_il/∑_i=1^T∑_l=1^m_aw^a_il.

Furthermore, also the expression for the parameters {{θ^c_i}_i=1^m_c,{θ^a_i}_i=1^m_a,Σ_m} can be computed in closed-form by using the switching least-squares approach as follows:

    A_c(j)=A_c(j)min ∑_i=1^Tw^c_ij‖ x_c(i)-A_c(j)x_c(i-1)‖^2,
       A_a(l)=A_a(l)min ∑_i=1^Tw^a_il‖ x_a(i)-A_a(l)x_a(i+1)‖^2,
       (C_c(j),C_a(l))=
       (C_c(j),C_a(l))min ∑_i=1^Tw^c_ijw^a_il‖ y(i)
    	-C_c(j)x_c(i)-C_a(l)x_a(i)‖^2.

Then, the covariance matrices related to the disturbances v_c, v_a, and v_m can also be estimated as:



    Σ_c(j)=
       ∑_i=1^T w^c_ij(x_c(i)-A_c(j)x_c(i-1))(x_c(i)-A_c(j)x_c(i-1))^⊤,
       Σ_a(l)=
       ∑_i=1^Tw^a_il(x_a(i)-A_a(l)x_a(i+1))(x_a(i)-A_a(l)x_a(i+1))^⊤,
       Σ_m=∑_i=1^T w^c_ijw^a_il(y(i)-C_c(j)x_c(i)-C_a(l)x_a(i))
       (y(i)-C_c(j)x_c(i)-C_a(l)x_a(i))^⊤.
































To show the convergence rate of the system matrices, we need the following definition of strong consistency of the parameter estimates, and few auxiliary results. Recall that θ̂ is the estimate of θ made by exploiting T samples. 

 The estimate θ̂ is strongly consistent if 
	lim_T→∞θ̂=θ.







	The following asymptotic relations hold true almost surely (a.s.):
	
    ‖∑_i=1^TA(s_c(i))x_c(i)v_c^⊤(i)+v_c(i)x_c^⊤(i)A(s_c(i))‖=o(T),
       ‖∑_i=1^TA(s_a(i))x_a(i)v_a^⊤(i)+v_a(i)x_a^⊤(i)A(s_a(i))‖=o(T).






















	The standard least-squares solution can be expressed as
	β̂_T=min_β‖ z-β^⊤ u‖^2,
	where u is the system inputs, z is the system output, β is the unknown parameters. Let U_T=∑_i=1^Tu(i)u(i)^⊤. If
	
		
  (C1) λ_min(U_T)→∞ a.s., and
		
  (C2) logλ_max(U_T)=o(λ_min(U_T)) a.s.,
	
	then the least-squares estimate β̂_T is strongly consistent with convergence rate
	
    ‖β̂_T-β‖_∞=𝒪(√(log(λ_max(U_T))/λ_min(U_T))) a. s.


Lemma <ref> indicates that the convergence rate of β̂_T depends on the covariance matrix of the system inputs.

We are now ready to establish the convergence rate for θ̂. Due to the possible different active subsystems at time t, it is convenient to define the following partition of the considered time interval {1, …, T} as 𝕋^c_j,T={t≤ T| s_c(t)=j} and 𝕋^a_l,T={t≤ T| s_a(t)=l}.

	Let W^c_j,T=∑_t∈𝕋^c_j,T x_c(t)x_c^⊤(t) and W^a_l,T=∑_t∈𝕋^a_l,Tx_a(t)x_a^⊤(t). Then, the estimate θ̂ generated by Algorithm <ref> is strongly consistent for any s_c∈Λ_c^T and s_a∈Λ_a^T, and the convergence rates are
	
    ‖Â_c(j)-A_c(j)‖_∞≤𝒪(√(log(λ_max(W^c_j,T))/λ_min(W^c_j,T))),
       ‖Â_a(l)-A_a(l)‖_∞≤𝒪(√(log(λ_max(W^a_l,T))/λ_min(W^a_l,T))),
       ‖Ĉ_c(j)-C_c(j)‖_∞≤𝒪(√(log(λ_max(W^c_j,T))/λ_min(W^c_j,T))),
       ‖Ĉ_a(l)-C_a(l)‖_∞≤𝒪(√(log(λ_max(W^a_l,T))/λ_min(W^a_l,T))).

	Furthermore, the convergence rate of Σ̂_c(j), Σ̂_a(l), Σ̂_m are
	
    ‖Σ̂_c(j)-Σ_c(j)‖_∞≤𝒪(log(T)/T),
       ‖Σ̂_a(l)-Σ_a(l)‖_∞≤𝒪(log(T)/T),
       ‖Σ̂_m-Σ_m‖_∞≤𝒪(log(T)/T).



	Theorem <ref> gives data-dependent upper bounds for the estimation errors of the parameter matrices. In order to have a data-independent characterization of the convergence rate for adaptive control or reinforcement learning purposes, in the proof of Theorem <ref>, specifically equation (<ref>), we provide with the corresponding convergence rate of the parameter estimate θ̂, which is equal to 𝒪(√(log(T)/T)).




§ NUMERICAL EXAMPLES

We now verify the effectiveness of the proposed methodology on two simulation examples. In both cases, we note that the true switching sequences s_c and s_a are only used to verify the accuracy of the estimated switching sequences, i.e., ŝ_c and ŝ_a. As performance index we make use of the mode match rate, defined as:

    L_mr= 1/T∑_i=1^Tι(s_c(i),ŝ_c(i)),

where ι(·,·) denotes the standard indicator function, i.e., ι(s_c(i),ŝ_c(i))=1 if s_c(i)=ŝ_c(i), 0 otherwise. 


 §.§ Example 1: Academic NCS-ASM

For illustrative purposes, we start by considering a simple non-causal system described in (<ref>) with m_c=m_a=2 modes and main parameters reported in Table <ref> (refer to the “True” columns). The dimensions of the outputs, causal states, and anti-causal states are n_y=1,n_x_c=n_x_a=2. The probabilities of all modes are π^c_1=0.7,π^c_2=0.3,π^a_1=π^a_2=0.5. The system is excited with white noise with zero mean and finite variance, and the data length is T=10^4.

The true and estimated parameters are reported in Table <ref>, which clearly shows that the parameter estimates are very close to their true values. In Fig. <ref> we report the partial estimation of the switching sequences s_c and s_a, where the mode match rates are 97.4% and 99.2%, respectively. Note that our method achieves an accurate parameter estimate, since each data point can be accurately assigned to the corresponding mode. 
To better validate the accuracy of the proposed algorithm in parameter estimation, Fig. <ref> illustrates the estimated states using the modified KF. The relative estimation errors, defined as δ_c=‖ x_c-x̂_c‖^2/‖ x_c‖^2 (δ_a has the same structure), are δ_c=3.74% and δ_a=3.14%, respectively.



 

For comparison purposes, we now implement the EM method for jump Markov linear systems proposed in <cit.>. The length of the data is set to T=10^4. The transition matrix in <cit.> is set to 𝒯=[[ 0.5 0.5; 0.5 0.5 ]], and the probability of the switching sequence in this paper is set to π^c_1=π^c_2=π^a_1=π^a_2=0.5. The subsystem match rates of the proposed method and <cit.> are compared at different noise levels by assuming Σ=Σ_c=Σ_a. The identification accuracy of the switching sequences are shown in Table <ref>.
 

To verify the robust performance of the proposed method against several noise levels, we run 100 Monte Carlo experiments under four different noise conditions, i.e., Σ∈{0.01I,0.1I,0.5I,I}. In Fig. <ref> we report the mean and the variance of the match rates in all the considered cases. We observe that the estimation accuracy of the switching sequence is not significantly affected by the noise variance, since even for high noise levels the estimation accuracy can still reach 98% due to the excellent performance of the modified KF.




 §.§ Example 2: The Department Store Inventory Price Index 

In this subsection we adopt “The Department Store Inventory Price Index"(DSIP) dataset from  The Bureau of Labor Statistics (BLS). These data come from inventory weighted price indices of goods carried by department stores.  

The department store inventory product prices can be influenced by various factors, such as supply and demand dynamics, market competition, seasonal variations and so on. These factors can collectively contribute to a complex relationship with the prices, thereby, suggesting that the inventory product prices may adhere to a mixed causal non-causal system. There are many causal relationships that affect the prices of department stores, such as supply and demand, cost of production, and market competition. These will directly lead to changes of future prices. There are also many non-causal factors that may not directly cause changes in the price of goods, such as economic conditions, consumer preferences, and trends. More importantly, seasonal changes can also lead to changes in its price index. Therefore, the system dynamics may change at different times to characterize DSIP. To sum up, a NCS-ASM  model (<ref>) is suitable for describing the DSIP.

In Fig. <ref> we show the true prices and the estimated prices with different number of subsystems. Specifically, we can infer that the larger the number of subsystems, the better the ability to describe the changes in the commodity price index. The estimation errors δ=‖y-ŷ‖/‖y‖ with different number of subsystems are shown in Table <ref>.
 


In conclusion, from Fig. <ref> and Table <ref> we note that switching systems with a larger number of modes can better describe the variations in DSIP, because seasonal changes can lead to different patterns in dynamic systems. The proposed identification method can therefore accurately capture the typical fluctuations of the DSIP.



§ CONCLUSION

We have proposed an expectation-maximization framework for identifying non-causal systems with arbitrary switching modes. In the E-step, we have embedded the reconstructed switching sequence into the modified Kalman filter so that the proposed algorithm can handle the joint state variable estimation for the causal and anti-causal parts. Furthermore, in M-step we have developed a switching least-squares algorithm  that can get the parameter estimates in closed-form. From a technical perspective, we have established the convergence of our identification methodology, also deriving an upper bound 𝒪(√(log(T)/T)) for the parameter errors. 

Note that the identification algorithm proposed in this paper can be adapted to the identification of switching linear descriptor systems with minor modifications, since a descriptor state-space model can be represented in the mixed causal and anti-causal form. When the subsystems are nonlinear, however, the identification task becomes more challenging, thus posing greater difficulties. This aspect will be further investigated in our future work. In addition, addressing the joint identification of structured subsystems and piecewise constant switching sequences is an interesting future research direction.





















99     





Tan2023 K. Tan, W. J. Parquette, & M. Tao. (2023). A predictive algorithm for maximum power point tracking in solar photovoltaic systems through load management. Solar Energy, 265, 112127.
Liao2024 S. Liao, Y. Wu, K. Ma, & Y. Niu, (2024). Ant Colony Optimization With Look-Ahead Mechanism for Dynamic Traffic Signal Control of IoV Systems. IEEE Internet of Things Journal, 11(1), 366-377.
Carloni2007 R. Carloni, R. G. Sanfelice, A. R. Teel, & C. Melchiorri. (2007). A hybrid control strategy for robust contact detection and force regulation. In Proc. American Control Conf., New York City, USA, 1461-1466.
Schlegl2003 T. Schlegl, M. Buss, & G. Schmidt. (2003). A hybrid systems approach toward modeling and dynamical simulation of dextrous manipulation. IEEE/ASME Trans. on Mechatronics, 8(3), 352-361.




Chan2008 Chan, A. B., & Vasconcelos, N. (2008). Modeling, clustering, and segmenting video
with mixtures of dynamic textures. IEEE Transactions on Pattern Analysis and
Machine Intelligence, 30(5), 909–926.

Ferrari2003 Ferrari-Trecate, G., Muselli, M., Liberati, D., & Morari, M. (2003). A clustering
technique for the identification of piecewise affine systems. Automatica, 39(2),
205–217.

Garulli2012 Garulli, Andrea, Paoletti, Simone, & Vicino, Antonio (2012). A survey on switched
and piecewise affine system identification. In 16th IFAC symposium on system
identification, Brussels, Belgium (pp. 344–355).

Bianchi2021 Bianchi, Federico, Breschi, Valentina, Piga, Dario, & Piroddi, Luigi (2021). Model
structure selection for switched NARX system identification: A randomized
approach. Automatica, 125, Article 109415.

Mark2022 Mark P. Balenzuela, Adrian G. Wills, Christopher Renton, & Brett Ninness. (2022). Parameter estimation for Jump Markov Linear Systems. Automatica, 135 109949.

Bemporad2018 Alberto Bemporad, Valentina Breschi, Dario Piga, & Stephen P. Boyd. (2018). Fitting jump models. Automatica, 96, 11-21.

Piga2020 Dario Piga, Valentina Breschi, & Alberto Bemporad. (2020). Estimation of jump Box–Jenkins models. Automatica, 120 109126.

Sayedana2024 Borna Sayedana, Mohammad Afshari, Peter E. Caines, & Aditya Mahajan. (2024). Strong Consistency and Rate of Convergence of Switched Least Squares System Identification for Autonomous Markov Jump Linear Systems. IEEE transactions on Automatic Control, 1-8.

Anna2018 Anna Scampicchio, Alberto Giaretta, & Gianluigi Pillonetto. (2018). Nonlinear Hybrid Systems Identification using Kernel-Based Techniques. In IFAC-PapersOnline, 51(15), 269-274.

Angelo2010 Angelo Alessandri, Marco Baglietto, & Giorgio Battistelli. (2010). A maximum-likelihood Kalman filter for switching discrete-time linear systems. Automatica, 46, 1870-1876.

Tsurumoto2022 K. Tsurumoto, W. Ohnishi, T. Koseki, N. Strijbosch, & T. Oomen. (2022). A non-causal approach for suppressing the estimation delay of state observer. 022 American Control Conference (ACC), 3356-3356.

El2023 El Ammari Anis, & Terzi, Chokri. (2023). Causal Nexus Between Ownership Structure, Dividend Policy and Financial Performance: A Bootstrap Panel Granger non-causality Analysis. Journal of African Business, 24(4), 562-579.


Verhaegen1996 Verhaegen, M. (1996). A subspace model identification solution to the identification of mixed causal, anti-causal LTI systems. SIAM Journal on Matrix Analysis and Applications, 17(2), 332–347.

Fang2024 X. Fang, & T. Chen. (2024). On kernel design for regularized non-causal system identification. Automatica, 159, 111335.

Blanken2020 Blanken, L., & Oomen, T. (2020). Kernel-based identification of non-causal systems with application to inverse model control. Automatica, 114.

sta_as1 T. E. Duncan, & B. Pasik-Duncan. (1990). Adaptive control of continuoustime linear stochastic systems. Math. Control signals systems, 3(1), 45–60.
sta_as2 M. K. S. Faradonbeh, A. Tewari, & G. Michailidis. (2020). On adaptive linear–quadratic regulators. Automatica, 117, 108982.

Dempster1977 Dempster, Arthur P., Laird, Nan M., & Rubin, Donald B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B. Statistical Methodology, 1–38.


















Lai1982 T. L. Lai, & C. Z. Wei. (1982). Least squares estimates in stochastic regression
models with applications to identification and control of dynamic
systems. Ann. Statist., 10(1), 154–166.

PE2018 P. E. Caines. (2018). Linear stochastic systems. SIAM.
Chen1986 H. F. Chen, & L. Guo. (1986). Convergence rate of least-squares identification and adaptive control for stochastic systems. Int J Control, 44(5), 1459–1476.

Stout1974 W. F. Stout. (1974). Almost Sure Convergence. Academic Press.

Gibson2005 Stuart Gibson, & Brett Ninness. (2005). Robust maximum-likelihood estimation of multivariable dynamic
systems. Automatica, 41, 1667-1682.

Kalman Kalman RE. (1960). A new approach to linear fltering and prediction problems. J Basic Eng, 82(1), 35–45.




























§ TECHNICAL PROOFS

Proof of Lemma <ref>: The log likelihood difference between the θ and θ^k can be expressed as 

    lnℙ_θ(y)-lnℙ_θ^k[y]=Q(θ,θ^k)-   Q(θ^k,θ^k)
       +V(θ,θ^k)-V(θ^k,θ^k),

where the difference V(θ,θ^k)-V(θ^k,θ^k) coincides with the Kullback–Leibler distance that possess an important property, i.e., being non-negative. Therefore, the maximization of Q(θ,θ^k) can yield an increase in the log-likelihood function lnℙ_θ(y), namely

    Q(θ,θ^k+1)≥ Q(θ,θ^k)⇒lnℙ_θ^k+1[y]≥lnℙ_θ^k[y],

thus concluding the proof.


Proof of Lemma <ref>: Only the boundedness of η_c(t) will be proven in detail, since that of η_a(t) and η_m(t) can be derived in a similar way.

First, we note that x_c(t-1) can be equivalently expressed  as follows:

    x_c(t-1)=φ_1(s_c)x_c(1)+φ_2(s_c)v_c(1:t-1),

where φ_1(s_c) and φ_2(s_c) are shown in A.1, A.2. v_c(1:t-1)≜[v_c(1),⋯,v_c(t-1)]. Both matrices are uniquely determined by the switching sequence s_c and system matrices A_c. Then, one obtains that:
TempEqCnt 
 
 


    η_c(t)   =x_c(t)-A(ŝ_c(t))x_c(t-1)
       =(A_c(s_c(t))-A(ŝ_c(t))) x_c(t-1)+v_c(t)
       =(A_c(s_c(t))-A(ŝ_c(t)))φ_1(s_c)x_c(1)
       +φ_3(s_c)v_c(1:t),

where φ_3(s_c)=[φ_2(s_c),1]. Passing to the (squared) norm in the expression above we note that, in view of the fact that the noise v_c has a bounded covariance, the last term is bounded too. For what concerns the first term, instead, we have:

    ‖ [A_c(s_c(t))-A(ŝ_c(t))]φ_1(s)x_c(1)‖^2≤λ_1‖ x_c(1)‖^2,

where 

    λ_1≜λ_max(φ^⊤_1(s)   (A_c(s_c(t))-A(ŝ_c(t)))^⊤
       (A_c(s_c(t))-A(ŝ_c(t)))φ_1(s)),

which concludes the proof.


Proof of Theorem <ref>: In the interest of space, we establish the convergence rate for Â_a(l) only, since the other bounds on the system matrices can be derived similarly.

In view of Lemma <ref>, sufficient conditions for establishing the convergence rate of Â_a(l) are (C1) λ_min(W^a_l,T)→∞, a.s., and (C2) logλ_max(W^a_l,T)=o(λ_min(W^a_l,N)), a.s.. We therefore have to show that these two conditions are verified in our case. Then, for what concerns (C1), one has:

    x_a(t)x_a(t)^⊤   =(Â_a(l)x_a(t+1)+v_a(t))
       (Â_a(l)x_a(t+1)+v_a(t))^⊤
       =Â_a(l)x_a(t+1)x_a^⊤(t+1)Â^⊤_a(l)
       +2v_a(t)x^⊤_a(t+1)Â^⊤_a(l)+v_a(t)v_a^⊤(t).

Since Â_a(l)x_a(t+1)x_a^⊤(t+1)Â^⊤_a(l) is a positive semidefinite matrix, by relying on Lemma <ref> we can infer that

    W^a_l,T   =∑_t∈𝕋^a_l,Tx_a(t)x^⊤_a(t)
       ≽∑_t∈𝕋^a_l,Tv_a(t)v_a^⊤(t)+x_a(T)x^⊤_a(T)
       +∑_t∈𝕋^a_l,T(Â_a(l)x_a(t+1)v^⊤_a(t)+v_a(t)x^⊤_a(t+1)Â^⊤_a(l))
       ≽∑_t∈𝕋^a_l,Tv_a(t)v_a^⊤(t)+o(T).

Then, we readily obtain:

    lim_|𝕋^a_l,T|→∞inf   ∑_t∈𝕋^a_l,Tx_a(t)x^⊤_a(t)/|𝕋^a_l,T|
       ≽lim_|𝕋^a_l,T|→∞inf∑_t∈𝕋^a_l,Tv_a(t)v_a^⊤(t)/|𝕋^a_l,T|≻0.

Therefore, we can conclude that λ_min(W^a_l,T)→∞ a.s..

To prove (C2) we note that:

    λ_max(∑_t∈𝕋^a_l,Tx_a(t)x^⊤_a(t))   ≤tr(∑_t∈𝕋^a_l,Tx_a(t)x^⊤_a(t))
       ≤∑_i=1^T‖ x_a(t)‖^2=𝒪(N),
 
where the last equality follows in view of the stability, in average sense, of the NCS-ASM in (<ref>).
Then, one can readily obtain that

    lim_T→∞log(λ_max(W^a_l,T))/λ_min(W^a_l,T)≤lim_T→∞log(T)/|𝕋_i,T|=log(T)/𝒪(T)=0.


We are now able to establish the convergence rate for the covariance matrices. Specifically, we will give the detailed proof for Σ̂_c(j) only, since the remaining ones follow similarly.

From the NCA-ASM in (<ref>), the true covariance matrix for v_c can be expressed as:

    Σ_c(j)=1/|𝕋^c_j,T|∑_t∈𝕋^c_j,T(x_c(t)-   A_c(j)x_c(t-1))
       (x_c(t)-A_c(j)x_c(t-1))^⊤.

Then, the estimation error can take the following form: 

    Σ̂_c(j)-Σ_c(j)=   1/|𝕋^c_j,T|∑_t∈𝕋^c_j,T((A_c(j)-Â_c(j))x_c(t-1))
       ((A_c(j)-Â_c(j))x_c(t-1))^⊤.

Therefore, the convergence rate for Σ̂_c(j) reads as:

    ‖Σ̂_c(j)-Σ_c(j)‖_∞≤   ∑_t∈𝕋^c_j,Tx_c(t-1)x_c(t-1)^⊤/|𝕋^c_j,T|
       ‖ (A_c(j)-Â_c(j))(A_c(j)-Â_c(j))^⊤‖_∞
       ≤𝒪(log(T)/T),

which completes the proof.


