Design of Optimal Intervention Based on a Generative Structural Causal Model

Abstract

In the industrial sector, malfunctions of equipment that occur during the production and operation process typically necessitate human intervention to restore normal functionality. However, the question that follows is how to design and optimize the intervention measures based on the modeling of actual intervention scenarios, thereby effectively resolving the faults. In order to address the aforementioned issue, we propose an improved heuristic method based on a causal generative model for the design of optimal intervention, aiming to determine the best intervention measure by analyzing the causal effects among variables. We first construct a dual-layer mapping model grounded in the causal relationships among interrelated variables to generate counterfactual data and assess the effectiveness of intervention measures. Subsequently, given the developed fault intervention scenarios, an adaptive large neighborhood search (ALNS) algorithm employing the expected improvement strategy is utilized to optimize the interventions. This method provides guidance for decision-making during equipment operation and maintenance, and the effectiveness of the proposed model and search strategy have been validated through tests on the synthetic datasets and satellite attitude control system dataset.

1. Introduction

Intervention is typically defined as a deliberate action or strategy aimed at changing the state of one or several variables, thereby affecting the overall behavior or output of the system. In fields such as medicine and social science, intervention often takes the form of controlled experiments, such as administering a specific dosage of a drug to patients in order to study its therapeutic effects. When measuring the impact of teaching methods on students, the specific content of the intervention might involve choosing between traditional and digital teaching approaches. As the theory of experimental design has evolved, the forms of experiment used to evaluate the effectiveness of interventions have become increasingly diverse, including methods such as randomized controlled trials (RCT) [1], among others.

As a term with rich meanings, intervention is often abstractly expressed as an operator $do(X=x)$ in statistics and causal reasoning, as proposed by Judea Pearl. This notation reflects a general meaning of intervention, which is to force one or more variables in the model set as a specific value or distribution [2]. The introduction of the do-calculus simplifies and essentializes the original intervention measure, greatly broadening the boundary of intervention research, especially when the intervention action is difficult to perform in reality due to ethical and legal restrictions, or higher implementation costs are required, the do operator converts the actual impossibility into theoretical possibility, allowing researchers to evaluate potential results through theoretical models without doing actual interventions. This counterfactual analysis method has been widely used in causal reasoning.

We all know that in the real world, there are complex connections between different entities. Broadly speaking, these connections can be categorized into two types: correlation and causality. Research on interventions and their impacts typically requires a thorough understanding of real causality. If we simply observe a certain correlation between the variables and rashly intervene, we may find that the measure does not have the desired effect due to the false relationship caused by the existence of underlying factors and thus results in a waste of resources [3].

For instance, suppose there are three variables that conform to the relationships illustrated in Figure 1, where X is the common cause of Y and Z, but X may not be observed due to incomplete prior knowledge. It can be found that since Y and Z are simultaneously affected by X, changes in X will invariably result in correlated fluctuations in the output data of Y and Z. If we do not know the existence of potential factor X in advance, we might be misled by this apparent correlation and make incorrect intervention decisions.

In order to accurately describe the real relationship between variables, Pearl introduced the concept of a structural causal model (SCM). SCMs are mathematical models based on the relationships between potential and observed variables through a series of structural equations [4]. These models provide a systematic approach to formally articulate the causal relationships among variables and have become increasingly significant in the realm of causal reasoning. In particular, when it is challenging to assess the effectiveness of an intervention by practical means, constructing a SCM is crucial. Through a SCM, the do operator can link the observed distribution to the intervention distribution, thereby estimating the nominal post-intervention value of a variable [5].

Furthermore, a key question for decision-makers is how to select a good intervention. Lee et al. [6] pointed out that the design of optimal intervention can be regarded as optimizing an expected strategy that contains two dimensions: observation and action. Simply put, it involves making the best choice based on the observed information when there are alternative interventions. Current methods that combine causal models with decision frameworks for the design of intervention primarily include two categories: Bayesian optimization and the causal bandit algorithm.

Bayesian optimization [7,8], as a global optimization method, employs a prior probability model and an acquisition function to optimize computationally expensive and non-analytical objective functions. With the increasing use of causal Bayesian networks (CBN) [9,10], incorporating causal information into the Bayesian optimization (BO) can better guide the search process for interventions. Causal Bayesian optimization (CBO) formulates a series of interventions by maximizing a target function. This method leverages causal relationships between variables to enhance the reliability of decision-making, and it effectively handles uncertainties through Bayesian techniques, delivering more robust optimization results. Meanwhile, it can adapt strategies in dynamic environments to accommodate new data and changes. Due to these advantages, causal Bayesian optimization has garnered considerable attention. Aglietti et al. [11] introduced a causal Bayesian optimization algorithm to infer optimal decision strategies using causal information available, aiming to address a new class of problems called causal global optimization (CGO). Further, due to the potential constraints in actual intervention scenarios, methods like constrained causal Bayesian optimization (cCBO) [12] have been proposed based on the theory of CBO. Considering that the choice of intervention measures is affected by the uncertainty of the causal structure, Branchini et al. [13] extended the causal Bayesian optimization framework and proposed a causal entropy optimization algorithm (CEO) that incorporates uncertainty information into the decision-making process. Zhang et al. [14] developed a causal active learning framework and designed a closed-form acquisition function to facilitate the rapid optimization of interventions. Similarly, Li et al. [15] designed the acquisition function in Bayesian optimization to make the causal reasoning process more robust to noise interference. Scott et al. [16] proposed a model-based CBO (MCBO) method that models the uncertainties of SCMs. Ren et al. [17] introduced a CBO method based on learning the distribution of exogenous variables to address the common issue of neglecting exogenous variables in existing methods. However, when using the CBO method, it is usually assumed that the causal graph of variables is known, while the structural equations are unknown, and probabilistic surrogate models, often linear Gaussian models, are employed as substitutes. Simple linear Gaussian models cannot model complex causal relationships, which can lead to unreasonable intervention designs. Additionally, CBO involves sequential decision-making, which has a high computational complexity that may exceed the data sampling time, rendering the algorithm impractical for real-world applications. In contrast, heuristic optimization offers rapid convergence and higher computational efficiency, effectively addressing these issues.

The design of intervention measures is essentially a decision optimization problem. Therefore, it was conceptualized as a stochastic multi-armed bandit problem with side information, where a good intervention is defined as one that maximizes the expected reward [18,19]. In the context of the causal bandit problem, Kroon et al. [20] introduced an estimator based on the property of separating sets into the discrete causal bandit algorithm, allowing them to learn better interventions without requiring the incorporation of prior causal knowledge. Bilodeau et al. [21] considered the case where the observed variables might be d-separators after intervention and proposed an adaptive optimal regret algorithm, extending its adaptability to conditions such as the front-door criterion. Moreover, methods like dynamic programming [22] and reinforcement learning [23] have also been employed to find superior interventions. The causal bandit algorithm utilizes causal reasoning to understand the impact of actions, thereby making more effective decisions based on causal relationships rather than mere correlations. It can also tailor actions to specific user environments by considering individual characteristics and previous interactions to improve outcomes. By focusing on causal effects, this type of algorithm can reduce the exploration needed to identify optimal actions. However, the effectiveness of such methods depends on the validity of the causal models used. Incorrect assumptions can lead to misleading results, and these methods can be computationally intensive, especially when dealing with large datasets and complex models, which may limit their scalability.

Although a large number of methods have been developed to study the optimization design of interventions based on causal models, there are still several shortcomings in the current approaches:

• Model-based approaches often employ simplistic causal modeling such as the Gaussian process (GP) [9,11,13] and linear Gaussian models [15,18]. However, in practical problems, target variables may exhibit non-Gaussian distributions, making traditional Gaussian-based models insufficient for capturing complex nonlinear relationships and non-stationary noise distributions among variables.
• The construction of intervention scenarios remains incomplete. For instance, causal bandit algorithms for multi-armed bandit problems typically focus on discrete or limited interventions. But, in reality, any value within the feasible range is possible, not just limited discrete ones. Furthermore, current researchers seldom consider optimization decisions under constraints [13], especially when we do not want other unmanipulated variables to deviate from our expectations after intervening on the target variable. This consideration cannot be ignored in practical applications.
• Heuristic algorithms are increasingly employed for intervention optimization [24,25]. However, in fault intervention scenarios, due to the large solution space of feasible interventions, applying optimization algorithms directly to such problems often results in becoming stuck in local optima, making it difficult to find better solutions quickly. Therefore, improvement in the optimization process of the algorithm is necessary.

Based on the discussion above, we propose an improved heuristic algorithm based on a causal generative model with the aim of designing optimal interventions and conducting corresponding experiments to demonstrate the effectiveness of this method. The main contributions of this paper are summarized as follows:

• A dual-layer mapping generative SCM is developed to characterize the causal structure among fault-related variables, along with the updating process of the counterfactual SCM following an intervention.
• The fault intervention problem is modeled based on real-world scenarios and formally described as a constrained optimization problem aimed at maximizing intervention utility.
• An ALNS algorithm based on the expected improvement strategy (EI-ALNS) is proposed to address the optimization problem of fault intervention, and experiments demonstrate the superiority of this heuristic method in optimizing intervention measures.

The value of the method for intervention design in this paper lies in its ability to provide guidance for decision-making during equipment operation and maintenance. Specifically, through real-time monitoring data, operators can assess whether the equipment is likely to fail in the near future reflected in changes in relevant parameters and proactively gain insights into potential failure information using predictive analytics and other methods. Based on this, operators can evaluate the effects of specific interventions (such as adjusting certain parameters within the system) in advance using causal generative models and optimize these interventions with improved heuristic algorithms, enabling them to implement effective measures when a failure occurs.

The remainder of this paper is organized as follows: Section 2 presents some background konwledge. Section 3 provides a detailed description of the implementation process of the method. Section 4 validates the proposed model and search strategy on both synthetic data and satellite attitude control system data. Section 5 summarizes the advantages and disadvantages.

2. Background

2.1. Generative Structural Causal Model

Let ${\left\{{\mathit{x}}_t\right\}}_{t=1}^M$ denote time-series data where ${\mathit{x}}_t={\left\{x_t^i\right\}}_{i=1}^n\in {ℝ}^n$, and $n$ represents the number of variables. Suppose the causal relationships among variables are known as a priori, thus allowing us to describe the generation mechanism of these variables through SCMs. According to Strobl [26], a complete SCM $\mathcal{M}=〈{\mathit{V}}_{\mathit{a}},\mathcal{F},P,G〉$ typically includes the following components:

1.

Set of variables ${\mathit{V}}_{\mathit{a}}=\{\mathit{V},\mathit{U}\}$: When considering causal relationships among system variables, it is generally believed that there are two types of variables—endogenous variables $\mathit{V}$ and exogenous variables $\mathit{U}$. The endogenous variables $\mathit{V}={\left\{V_i\right\}}_{i=1}^n$ constitute the main body of the SCM and are directly influenced by other variables within the model, while the exogenous variables $\mathit{U}={\left\{U_i\right\}}_{i=1}^n$ indicate the sources that are external to the causal relationships, typically representing random factors or dependencies not considered in the model.

2.

Structural equations $\mathcal{F}={\left\{{\mathcal{F}}_i\right\}}_{i=1}^n$: Set of functions that characterize the causal relationships among the endogenous variables:

$$V_i={\mathcal{F}}_i\left(Pa(V_i),U_i\right)$$

(1)

Here, ${\mathcal{F}}_i$ is the structural function corresponding to the endogenous variable $V_i$, and $Pa(V_i)$ represents the parent variables of $V_i$, i.e., the set of variables that directly influence $V_i$.

3.

Joint probability density of exogenous variables $P(\mathit{U})$: Determines the randomness of the entire model and impacts the accuracy of analysis and prediction when using the model.

4.

Causal graph $G=({\mathit{V}}_a,\mathit{E})$: Visually describes the causal relationships between variables, which satisfies the constraints of a directed acyclic graph (DAG). $\mathit{E}$ indicates the directed edges in $G$.

In this paper, SCM is treated as a generative model in order to answer the following question: ”What would the data look like after the intervention?”. The role of the generative SCM primarily encompasses the following two aspects. On one hand, it can capture complex relationships within the causal structure and simulate the impacts of intervened variables on other variables in the system. On the other hand, it provides a formalized counterfactual perspective, which allows the model to be trained to generate data under hypothetical scenarios, thereby facilitating the analysis and evaluation of hypothetical actions.

It should be noted that the generative SCM described in this paper is a surrogate model. It requires training with observational data to approximately simulate the data generation process for each variable. The causal structure $G$ is a priori, predetermined by expert information, which is a critical assumption in constructing the model.

2.2. Post-Nonlinear Causal Model

The design of the generative model in this paper, which incorporates a dual-layer mapping structure, draws formal inspiration from the nonlinear acyclic causal model proposed by Zhang et al. [27], which is also known as post-nonlinear causal (PNL) model. PNL is one type of SCM, but it has a more complex functional form. Mathematically, the representation of the PNL model for each variable can typically be articulated as follows:

$$x^{(i)}={\mathcal{F}}_{i,2}\left({\mathcal{F}}_{i,1}\left(P{a}_G^{(i)}\right)+{\epsilon }^{(i)}\right),i=1,2,\dots ,n$$

(2)

where ${\mathcal{F}}_{i,1}$ denotes the nonlinear effect of the causal variables $P{a}_G^{(i)}$, ${\epsilon }^{(i)}$ is the independent disturbance, and ${\mathcal{F}}_{i,2}$ denotes post-nonlinear distortion in variable $x^{(i)}$.

Specifically, the two nonlinear processing layers of the PNL model can be viewed as a “progressive” relationship, where the data generation process for each variable can be divided into the following two stages:

1.

The first stage

The primary goal of the first layer is to extract preliminary nonlinear features from the input causal variables. In this layer, we construct a nonlinear transformation of its parents $P{a}_G^{(i)}$, denoted by ${\mathcal{F}}_{i,1}\left(P{a}_G^{(i)}\right)$, plus some disturbance ${\epsilon }^{(i)}$, thereby obtaining the intermediate outputs $u^{(i)}={\mathcal{F}}_{i,1}\left(P{a}_G^{(i)}\right)+{\epsilon }^{(i)}$ generated by this layer.

2.

The second stage

Nonlinear distortion is applied to the output of the first stage $u^{(i)}$, which is used as input for the transformation of the second layer $x^{(i)}={\mathcal{F}}_{i,2}\left(u^{(i)}\right)$, thereby producing the final output $x^{(i)}$.

Compared to the traditional linear SCMs [14] and single-layer nonlinear SCMs [28], the expressiveness and superiority of the PNL model are primarily reflected in the following aspects:

• Modeling capability of nonlinear relationships

Traditional linear causal models $x^{(i)}={\displaystyle \sum B_{ik}{x}_k+{\epsilon }^{(i)}}$ can only describe the direct linear relationships between variables, while single-layer nonlinear causal models $x^{(i)}={\mathcal{F}}_i\left(P{a}_G^{(i)}\right)+{\epsilon }^{(i)}$ often only depict some relatively simple dependency structures. The PNL model, through the introduction of two nonlinear functions ${\mathcal{F}}_{i,1}$ and ${\mathcal{F}}_{i,2}$, is able to more finely decompose complex nonlinear causal relationships, thereby capturing the nonlinear effect of the causes [29].

• Flexible control of error structures

It is known that in linear and single-layer nonlinear models, the error terms are additive noise. However, in the PNL model, the error term is encompassed by the posterior nonlinear function, allowing the model to handle non-stationary noise distributions and external disturbances in a more complex manner [30], depicting the noise effect in the transmission process from $P{a}_G^{(i)}$ to $x^{(i)}$, which follows an unconventional distribution.

• Adaptability and robustness of the model

The nonlinear component can be customized to adapt to specific data structures, which are often regarded as nonlinear distortion caused by the sensor or measurement. This makes the model more adaptable and robust when facing varied and heterogeneous real-world data.

2.3. Counterfactual Updating in SCM

When we intervene on variables in endogenous variables $\mathit{V}$, the original SCM may no longer be able to effectively describe the new actual situation. This is because the intervention operation forces the variable to take on a specific value, causing the intervened variable not to be influenced by its parent variables, which is reflected in the DAG by deleting incoming edges of the corresponding variable. The causal structure is revised in the “counterfactual world”, leading to a change in the data generation mechanism, thereby simulating counterfactual occurrences that do not exist in the real world. Formally, suppose we perform the intervention $do(\mathit{P}{\mathit{a}}_G^{(i)}={\mathit{X}}_{iv})$, replacing $\mathit{P}{\mathit{a}}_G^{(i)}$ with intervention value ${\mathit{X}}_{iv}$. The original SCM $\mathcal{M}$ then transforms into the counterfactual SCM $\tilde{\mathcal{M}}=〈{\mathit{V}}_{\mathit{a}},{\mathcal{F}}^{\prime },P^{\prime },G^{\prime }〉$, where ${\mathcal{F}}^{\prime }$ represents the counterfactual distribution that needs to be relearned, $P^{\prime }$ is the probability density of posterior noise, $G^{\prime }$ is the post-intervention DAG.

2.4. Intervention Optimization

In practical operations, when intervening at a fault, we typically consider how to balance the cost and benefit of the intervention measures, because the application of intervention measures will inevitably lead to the consumption of resources. Specifically, if the intervention intensity is too low, the control over the fault may not achieve the desired effect. On the other hand, excessive intervention could lead to unnecessary waste of manpower and material resources. Therefore, a good intervention measure should have a high utility value, i.e., achieving better intervention results at a lower cost. To achieve the optimal intervention intensity, we formulate the intervention optimization as a classic constrained optimization problem, and then we can develop a heuristic algorithm to solve it.

3. Proposed Method

The method for designing optimal intervention proposed in this paper primarily consists of the following two components, as illustrated in Figure 2:

• Construction of generative model

Before devising intervention measures, the first consideration is how to evaluate the impact of intervention measures on the fault. To solve this, we introduce a generative SCM based on a dual-layer mapping structure called graph-based warped multi-task Gaussian process (G-WMGP), which includes two sub-models: G-MTGP and normalizing flows. We also detail the updating process through which the G-WMGP model is adapted into the counterfactual G-WMGP model in the context of an intervention, thereby providing a theoretical foundation for the generation of post-intervention data.

2.

Formulation of intervention measures

To transform the intervention decision into a constrained optimization problem, we define the objective function and the specific form of the intervention measures, where the intervention measures align with the “3W” intuitive thinking of decision-makers in real-world scenarios: the target of the intervention (“which”), the timing and duration of the intervention (“when”), and the specific values of the intervention (“what”). Given this, we consider employing a heuristic optimization algorithm to identify high-quality intervention measures, where an intervention measure is defined as an individual solution. The objective function is calculated based on the counterfactual data, which are generated by the generative SCM, and the intervention measure generated in each iteration will act as input to “urge” the update of the SCM. Moreover, we propose a search strategy aimed at addressing efficiency (i.e., how to quickly find an optimal solution) and accuracy (i.e., how to enhance the quality of the solution) issues during the algorithm’s execution.

3.1. Data Generation Based on SCM

3.1.1. Generative Model Construction

In order to facilitate the construction of the model, we assume that the causal relationships among variables satisfy the condition of causal stationarity [31], which posits that the causal structure $G$ does not change over time. Further, considering the causal links between any variables exhibit a relatively complex nonlinear pattern, we employ the structure of a post-nonlinear causal model [32] to simulate the data generation process of the variables of interest, i.e., for each variable $V_i$:

$${\mathit{x}}^{(i)}=g_i\left(f_i(\mathit{P}{\mathit{a}}_G^{(i)})+{\epsilon }^{(i)}\right), i=1,2,\dots ,n$$

(3)

Here, both $g_i$ and $f_i$ represent nonlinear mapping relationships, ${\epsilon }^{(i)}$ represents mutually independent noise that illustrate the potential influence of the exogenous variable $U_i$ on the system variable $V_i$, where ${\epsilon }^{(i)}\sim \mathcal{N}(0,{\sigma }_i^2)$. The parent variables of $V_i$ in the causal structure $G$ are denoted as $P{a}_G^{(i)}$, where $P{a}_G^{(i)}\subseteq \mathit{V}$, and the bold form $\mathit{P}{\mathit{a}}_G^{(i)}$ represents the values of the corresponding variables. Intuitively, serving as the input to the generative model, $P{a}_G^{(i)}$ directly impacts the generation data of $V_i$ [30]. Therefore, we hope to construct a transformation ${\mathcal{F}}_i=g_i(f_i(•))$ that can effectively characterize the data generation mechanism when $V_i$ is influenced by $P{a}_G^{(i)}$.

To achieve this, we introduce a generative SCM called G-WMGP, which is based on the PNL dual-layer mapping architecture for designing the specific functional form of transformation ${\mathcal{F}}_i\in \mathcal{F}$. The principle involves mapping the original input through a war-ping transformation to a new input space, where a graph-based multi-task Gaussian process (G-MTGP) is established. The inner mapping $f_i$ uses the G-MTGP model, which is an extension of GP and denoted as $f_i\sim \mathcal{G}\mathcal{P}(\mu ,{\mathbf{K}}_{\theta ,i})$, where the inputs and outputs are based on the many-to-one relationships among variables in the causal graph $G$. In the multi-task learning framework, the correlations among tasks are utilized to improve the learning performance; thus, the overall dataset can be viewed as a collection containing $S$ tasks, with each task $r_j\in \left\{r_j,j\in P{a}_G^{(i)}\right\}$ comprising $M$ samples. The covariance matrix ${\mathbf{K}}_{\theta ,i}$ encompasses the covariance between different tasks and among different data points within the same task, characterized by a coregionalization kernel in this paper as follows:

$${\mathrm{k}}_{\theta ,i}{\left(\mathit{P}{\mathit{a}}_G^{(i)}(m),\mathit{P}{\mathit{a}}_G^{(i)}(m^{\prime })\right)}_{j,j^{\prime }}={\sigma }^2\left(1+{\displaystyle \frac{\sqrt{3}}{l}}\alpha \right)\exp \left(-{\displaystyle \frac{\sqrt{3}}{l}}\alpha \right)\cdot Q(r_j,r_{j^{\prime }})$$

(4)

$$Q(r_j,r_{j^{\prime }})={\left(w{w}^T+diag(u)\right)}_{j,j^{\prime }}$$

(5)

Here, $m,m^{\prime }=1,2,\dots ,M$. $\alpha ={‖\mathit{P}{\mathit{a}}_G^{(i)}(m)-\mathit{P}{\mathit{a}}_G^{(i)}(m^{\prime })‖}_2$ denotes the Euclidean distance between data points at $m$ and $m^{\prime }$, where $\mathit{P}{\mathit{a}}_G^{(i)}(m),\mathit{P}{\mathit{a}}_G^{(i)}(m^{\prime })\in {ℝ}^n$. ${\sigma }^2$ and $l$ represent scale factors, and $\theta$ summarizes parameters from the multi-task Gaussian process. The positive definite matrix $Q(r_j,r_{j^{\prime }})$ describes the correlation between tasks $r_j$ and $r_{j^{\prime }}$ and corresponds to the variance of task $r_j$ while $j=j^{\prime }$.

The outer mapping $g_i$ employs a generative architecture of normalizing flows, aiming at constructing flexible probability distributions on continuous random variables. The main idea is to express ${\mathit{x}}^{(i)}$ as the transformation $g_i$ of a real vector ${\mathit{u}}^{(i)}$ sampled from $p_u$:

$${\mathit{x}}^{(i)}=g_i\left({\mathit{u}}^{(i)}\right)$$

(6)

where $p_{\mathrm{u}}$ is regarded as the base distribution of the flow-based model, which equals $\mathcal{N}(\mu ,{\mathbf{K}}_{\theta ,i}+{\sigma }_i^2)$ due to the construction of the G-MTGP model. To enhance the learning capability of the flow model, $g_i$ is typically composed of $n^{\prime }$ sub-flows, i.e., $g_i=T_1\circ T_2\cdots \circ T_{n^{\prime }}$, while each of them can be a different flow type. When building the flow-based model, it is necessary to ensure that $g_i$ is invertible and both $g_i$ and its inverse $g_i^{-1}$ are differentiable. According to the property of the inverse, $g_i^{-1}$ are given by the following:

$$g_i^{-1}=T_{n^{\prime }}^{-1}\circ T_{n^{\prime }-1}^{-1}\cdots \circ T_1^{-1}$$

(7)

Under these conditions, the density of ${\mathit{x}}^{(i)}$ can be obtained through the change in variables:

$$p_{\mathrm{x}}({\mathit{x}}^{(i)})=p_{\mathrm{u}}({\mathit{u}}^{(i)}){\left|\det {\mathcal{J}}_{g_i}({\mathit{u}}^{(i)})\right|}^{-1}$$

(8)

where $\left|\det {\mathcal{J}}_{g_i}({\mathit{u}}^{(i)})\right|$ denotes the absolute Jacobian determinant of $g_i$, which quantifies the relative change in a small neighborhood around ${\mathit{u}}^{(i)}$ due to the transformation of $g_i$. Because ${\mathit{u}}^{(i)}=g_i^{-1}({\mathit{x}}^{(i)})$, according to Equation (6), $p_{\mathrm{x}}({\mathit{x}}^{(i)})$ can be written as follows:

$$p_{\mathrm{x}}({\mathit{x}}^{(i)})=p_{\mathrm{u}}(g_i^{-1}({\mathit{x}}^{(i)}))\left|\det {\mathcal{J}}_{g_i^{-1}}({\mathit{x}}^{(i)})\right|$$

(9)

Similarly, when containing composable sub-flows, the Jacobian determinant of $g_i$ is calculated as follows:

$$\det {\mathcal{J}}_{g_i}({\mathit{u}}^{(i)})=\det {\mathcal{J}}_{T_1}(T_2({\mathit{u}}^{(i)}))\cdot \det {\mathcal{J}}_{T_2}(T_3({\mathit{u}}^{(i)}))\cdots \det {\mathcal{J}}_{T_{n^{\prime }-1}}(T_{n^{\prime }}({\mathit{u}}^{(i)}))$$

(10)

The loss function used for model training can be characterized by the divergence between the flow model $g_i({\mathit{x}}^{(i)};\phi )$ and the target distribution $p_{\mathrm{x}}^{\ast }({\mathit{x}}^{(i)})$, specifically employing the forward Kullback–Leibler divergence and approximated using the Monte Carlo method [33]:

$$\begin{array}{l}\mathcal{L}(\theta ,\phi )=-{\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M\frac{1}{2}}\left({\mathit{y}}^T{({\mathbf{K}}_{\theta ,i}+{\sigma }_i^2\mathbb{I})}^{-1}\mathit{y}+\log \left|{\mathbf{K}}_{\theta ,i}+{\sigma }_i^2\mathbb{I}\right|+MS\log (2\pi )\right)\\ +\log \left|\det {\mathcal{J}}_{g_i^{-1}}\left(x_m^{(i)};\phi \right)\right|+\mathrm{const}\end{array}$$

(11)

where $\mathit{y}:=g_i^{-1}\left({\mathit{x}}_m^{(i)};\phi \right)-\mu$, $\phi$ defines the parameters of the flow model. $\left|\det {\mathcal{J}}_{g_i^{-1}}\right|$ represents the Jacobian determinant of the inverse transformation $g_i^{-1}$, calculated as follows:

$$\left|\det {\mathcal{J}}_{g_i^{-1}}({\mathit{x}}_m^{(i)};\phi )\right|={\displaystyle \prod _{m=1}^M\left|\frac{\partial g_i^{-1}({\mathit{x}}_m^{(i)};\phi )}{\partial x_i}\right|}$$

(12)

Through Equation (5), parameters $\{\theta ,\phi \}$ can be learned with the observational data ${\mathit{X}}_F^{(i)}={\left\{{\mathit{x}}^{(i)},\mathit{P}{\mathit{a}}_G^{(i)}\right\}}_F$.

By constructing a multi-layered nested architecture, the G-WMGP model not only more accurately captures the latent structures and nonlinear dependencies in the data, but also provides high-quality data generation and explicit modeling of complex distributions, while maintaining the flexibility and interpretability of the model.

3.1.2. Counterfactual SCM Calculation

According to the three-step method proposed by Pearl, the post-intervention counterfactual data ${\mathit{x}}_{CF}^{(i)}$ can be obtained through the following process [28,33]:

1.

Abduction: First, calculate the posterior noise distribution for the transformation ${\mathcal{F}}_i$ based on ${\mathit{X}}_F^{(i)}$, which conforms to a Gaussian distribution according to Weilbach et al. [29], as follows:

$${\epsilon }_F^{(i)}\left|{\mathit{X}}_F^{(i)}\right.\sim \mathcal{N}\left({\widehat{\mu }}_{\epsilon ,i},{\widehat{\sigma }}_{\epsilon ,i}^2\right)$$

(13)

where

$${\widehat{\mu }}_{\epsilon ,i}={\sigma }_i^2{\left({{\mathbf{K}}^{\prime }}_{\theta ,i}+{\sigma }_i^2\mathbb{I}\right)}^{-1}\left(g_i^{-1}({\mathit{X}}_F^{(i)})-\mu \right)$$

(14)

$${\widehat{\sigma }}_{\epsilon ,i}^2={\sigma }_i^2{{\mathbf{K}}^{\prime }}_{\theta ,i}{\left({{\mathbf{K}}^{\prime }}_{\theta ,i}+{\sigma }_i^2\mathbb{I}\right)}^{-1}$$

(15)

Here, ${{\mathbf{K}}^{\prime }}_{\theta ,i}$ is calculated with ${\left\{\mathit{P}{\mathit{a}}_G^{(i)}\right\}}_F$.

2.

Action: Then, update the structural equation ${\mathcal{F}}_i$ with the intervention measure $do(\mathit{P}{\mathit{a}}_G^{(i)}={\mathit{X}}_{iv})$, i.e., replacing $\mathit{P}{\mathit{a}}_G^{(i)}$ with value ${\mathit{X}}_{iv}$:

$${\mathit{x}}_{CF}^{(i)}\left(do(\mathit{P}{\mathit{a}}_G^{(i)}={\mathit{X}}_{iv})\right)=g_i\left(f_i({\mathit{X}}_{iv})+{\epsilon }_F^{(i)}\right)$$

(16)

3.

Prediction: Finally, calculate the predictive posterior distribution of $Z_i=f_i({\mathit{X}}_{iv})+{\epsilon }_F^{(i)}$ at ${\mathit{X}}_{iv}$:

$$Z_i\left|{\mathit{X}}_F^{(i)}\right.\sim \mathcal{N}\left({\tilde{\mu }}_i+{\widehat{\mu }}_{\epsilon ,i},{\tilde{\sigma }}_i^2+{\widehat{\sigma }}_{\epsilon ,i}^2\right)$$

(17)

where ${\tilde{\mu }}_i$ and ${\tilde{\sigma }}_i$ are given by

$${\tilde{\mu }}_i=\mu +{\tilde{\mathbb{K}}}_i{\left({{\mathbf{K}}^{\prime }}_{\theta ,i}+{\sigma }_i^2\mathbb{I}\right)}^{-1}\left({\mathit{x}}_F^{(i)}-\mu \right)$$

(18)

$${\tilde{\sigma }}_i^2={\tilde{\mathrm{k}}}_i-{\tilde{\mathbb{K}}}_i{\left({{\mathbf{K}}^{\prime }}_{\theta ,i}+{\sigma }_i^2\mathbb{I}\right)}^{-1}{\tilde{\mathbb{K}}}_i^{\top }$$

(19)

${\tilde{\mathrm{k}}}_i$ is the covariance matrix among ${\left\{{\mathit{X}}_{iv}(m)\right\}}_{m=1}^M$, while ${\tilde{\mathbb{K}}}_i$ is calculated as another between ${\mathit{X}}_{iv}$ and ${\left\{\mathit{P}{\mathit{a}}_G^{(i)}\right\}}_F$. Then, we can obtain the post-intervention counterfactual values through samples from the posterior distribution of $Z_i$:

$${\mathit{x}}_{CF}^{(i)}=g_i\left({\mathit{z}}_{CF}^{(i)}\right)$$

(20)

By constructing SCM $\tilde{\mathcal{M}}$, it is possible to derive a counterfactual mapping ${\widetilde{{\mathcal{F}}^{\prime }}}_i$ from the causes of $V_i$ to $V_i$ if an intervention measure is implemented on $P{a}_G^{(i)}$.

3.2. Design of Optimal Intervention Based on EI-ALNS

In this section, we primarily discuss the process of seeking better solutions for intervention measures based on the optimization algorithm and the improved search strategy, which aims at enhancing the optimization capability of the algorithm. We first elaborate on the fault intervention issue in Section 3.2.1, and construct the scenario based on real-world fault conditions and intervention needs, while also defining the optimization objective and constraints. In Section 3.2.2, we detail the principle of the expected improvement search strategy which is abbreviated as EI, and then provide the implementation process of the EI-ALNS algorithm in the context of intervention optimization.

3.2.1. Problem Modeling

Assume that the faulty variable is denoted as $V_k$, and the interval during which the fault occurs is represented as $\mathcal{S}$, with a length of $M$. Our goal is to implement interventions to eliminate the abnormal state of $V_k$ as much as possible. To achieve this, an intervention measure $\mathcal{R}=({\mathit{V}}_{iv},t_{iv},{\mathcal{L}}_{iv},{\mathit{X}}_{iv})$ should include the following elements:

• Which: Intervened variable ${\mathit{V}}_{iv}$, a discrete set. Given the causal structure, we can clearly identify which variables directly or indirectly influence $V_k$. Considering the constraints of cost and operability in real-world scenarios, and decision-makers may hope to eradicate the fault at its source rather than merely addressing the symptoms, it is common practice to intervene on a subset ${\mathit{V}}_{iv}$ of $P{a}_{\mathrm{G}}^{(k)}$.
• When: Intervention start time $t_{iv}\in \mathcal{S}$ and step length ${\mathcal{L}}_{iv}$, which lies within the range $[1,M-{\mathcal{L}}_{iv}]$, discrete values. These two constitute the time interval of intervention, indicating when the decision-makers start the intervention and how long the intervention measures last.
• How: Specific values for intervention ${\mathit{X}}_{iv}\in {ℝ}^{{\mathit{V}}_{iv}\times {\mathcal{L}}_{i\nu }}$, continuous values. The intervention values ${\mathit{X}}_{iv}$, which should be within the normal range of the variable, reflect the intensity of the intervention through their difference from the original data.

In order to transform the intervention decision problem into a classic constrained optimization problem, we first define the following variables of interest within $G$:

• Intervention variables: Variables about to be intervened denoted as $P{a}_{\mathrm{G}}^{(k)}$.
• Target variable: The faulty variable $V_k$. Ensuring that the faulty variable returns to normal functioning as much as possible is one of the core purposes of designing intervention measures.
• Constraint variables: For certain critical variables ${\mathit{V}}_c\subseteq \mathit{V}$ in the system, it is necessary to avoid causing new faults when intervening in order to control the adverse effects of the intervention.

Intervention values should be continuously derived within the normal range of the variables, as mentioned above. In this paper, $x_{out,t}^{(i)}$ is defined as a deviation from the expected interval $Q(i,t)$:

$$x_{out,t}^{(i)}\notin Q(i,t):=\left({\mu }_{tar,t}^{(i)}-{\delta }_{i,t},{\mu }_{tar,t}^{(i)}+{\delta }_{i,t}\right)$$

(21)

where $Q(i,t)$ represents the normal interval of $V_i$ at time $t$, and ${\mu }_{tar}^{(i)}$ denotes the target values under non-intervention conditions. ${\delta }_{i,t}$ represents the permissible fluctuation range for the variable in factual scenarios, which can be determined based on expert knowledge or approximated using the statistical interquartile range method. Thus, we can calculate the abnormal state of $V_k$ at time $t$:

$${\mathcal{C}}_{state,t}^{(k)}=\left|{\mu }_t^{(k)}\left(do(\mathit{P}{\mathit{a}}_G^{(k)}={\mathit{X}}_{iv})\right)-{\mu }_{tar,t}^{(k)}\right|-{\delta }_{k,t}, if {\mu }_t^{(k)}\left(do(\mathit{P}{\mathit{a}}_G^{(k)}={\mathit{X}}_{iv})\right)\notin Q(k,t)$$

(22)

The post-intervention expected values of $V_k$, denoted as ${\mu }_t^{(k)}\left(do(\mathit{P}{\mathit{a}}_G^{(k)}={\mathit{X}}_{iv})\right)$, are derived from the aforementioned counterfactual update process. Given these definitions, the relative abnormal value can be expressed as follows:

$${\mathcal{C}}_{state}^{(k)}={\displaystyle \frac{1}{\left|{\overline{\mathit{x}}}^{(k)}\right|}}{\displaystyle \sum _{t\in \mathcal{S}}{\mathcal{C}}_{state,t}^{(k)}}$$

(23)

where ${\overline{\mathit{x}}}^{(k)}$ represents the mean of outlier data within $\mathcal{S}$. Taking practical context into account, we not only aim to mitigate the abnormal state of variables through intervention measures, but also seek to ensure that the cost of intervention is not excessively high. The cost of each intervened variable $V_b$ can be represented by the following:

$${\mathcal{C}}_{cost}^{(b)}={\displaystyle \sum _{t\in \mathcal{S}}\left|\frac{1}{x_t^{(b)}}\left(x_{iv,t}^{(b)}-x_t^{(b)}\right)\right|}, V_b\in {\mathit{V}}_{iv}$$

(24)

And the total cost of the intervention is the sum, as follows:

$${\mathcal{C}}_{cost}={\displaystyle \sum _{V_b\in V_{iv}}{\mathcal{C}}_{cost}^{(b)}}$$

(25)

Taking relative changes aims to eliminate differences in dimensions among variables. In summary, the objective function and the constraints of this optimization problem can be formalized as follows:

$$F(\mathcal{R})={\mathcal{C}}_{cost}+{\mathcal{C}}_{state}^{(k)}$$

(26)

$$\mathrm{s}.\mathrm{t}.\left\{\begin{array}{l}{\mu }_t^{(i)}\left(do(\mathit{P}{\mathit{a}}_G^{(k)}={\mathit{X}}_{iv})\right)\in Q(i,t), for V_i\in {\mathit{V}}_c\\ {\mathit{x}}_{iv,t}^{(b)}\in Q(b,t), for V_b\in {\mathit{V}}_{iv}\end{array}\right.$$

(27)

The optimal intervention measure is thus determined by solving the following:

$${\mathcal{R}}^{\ast }=\arg \underset{s.t.}{\min }F(\mathcal{R})$$

(28)

3.2.2. Expected Improvement Search Strategy

The ALNS algorithm is a powerful heuristic method that relies on its flexible neighborhood search mechanism and adaptive strategy adjustments, aiming at solving large-scale and complex constrained optimization problems [34]. The basic idea of the ALNS algorithm is to start from the current solution and dynamically search for new solutions in the solution space through a process of destruction and reconstruction, combined with adaptive mechanisms and acceptance criteria.

Although the ALNS algorithm itself has a strong global optimization capability, directly applying it to the optimization of intervention measures reveals challenges. Specifically, when there are many intervention variables and large intervention intervals, due to the solution form $\mathcal{R}=({\mathit{V}}_{iv},t_{iv},{\mathcal{L}}_{iv},{\mathit{X}}_{iv})$, which includes continuous optimization of intervention values, the solution space becomes very large, making it more difficult to find high-quality solutions without extensive search time.

Considering the aforementioned issues and the difficulty of obtaining superior solutions using the basic form of the ALNS algorithm, we aim to refine the optimization process by combining the ideas of global and local neighborhood searches, thereby avoiding the pitfalls of local optima. Since the complex relationships among variables in the scenario have already been formalized using a SCM, we consider how to integrate information provided by the causal model to assist in guiding the search process of the algorithm. We have found that the cCBO approach improves the intervention optimization process by incorporating information from the causal surrogate model, which has been proven effective both theoretically and practically [12]. This idea has greatly inspired us to improve the algorithm.

We propose a search strategy with the core principle that during the iterative optimization of intervention values, the direction and step size of the neighborhood search are adaptively adjusted based on the effectiveness of the current intervention, thereby enhancing the algorithm’s optimization capability. As a classic optimization search method, the gradient descent method enables the solution to be updated towards the current negative gradient direction. Considering that the objective function (20) is smooth and differentiable, we can perform local optimization on the intervention values ${\mathit{X}}_{iv}$ based on gradient information. Formally, for the current intervention measure $\mathcal{R}=({\mathit{V}}_{iv},t_{iv},{\mathcal{L}}_{iv},{\mathit{X}}_{iv})$, we derive a new one ${\mathcal{R}}^{\prime }=({\mathit{V}}_{iv},t_{iv},{\mathcal{L}}_{iv},{{\mathit{X}}^{\prime }}_{iv})$ based on the following improvement rule:

$$\mathrm{\Delta }{\mathit{X}}_{iv,t}={{\mathit{X}}^{\prime }}_{iv,t}-{\mathit{X}}_{iv,t}=-{\displaystyle \sum _{n=1}^{n_0}\beta {\nabla }_{x}F{({\mathcal{R}}^{(n)})}_t-{\mathcal{Z}}_t, t\in (t_{iv},t_{iv}+{\mathcal{L}}_{iv})}$$

(29)

where $n_0$ represents the number of iterations for executing gradient descent, and $\beta$ is the learning rate. ${\mathcal{Z}}_t$ is the step size control term, which serves to adjust the step size in the search process during the iterations of the algorithm. Increase the search step size when the current intervention $\mathcal{R}$ significantly deviates from the theoretically optimal intervention ${\mathcal{R}}^{\ast }$. Conversely, the step size is gradually reduced to achieve a more refined search. For the intervened variable $V_q\in {\mathit{V}}_{iv}$ at time $t$, we have the following:

$${\mathcal{Z}}_t^{(q)}=\alpha \left|x_{iv,t}^{(q)}-x_t^{(q)}\right|$$

(30)

Here, $\alpha$ is the control coefficient, which is adaptively adjusted during the optimization process. Given the current solution $\mathcal{R}$, we can calculate ${\mathcal{C}}_{state}^{(k)}$ for $V_k$ and ${\mathcal{C}}_{cost}^{(q)}$ for $V_q$ while taking their ratio as the value of $\alpha$:

$$\alpha ={\mathcal{C}}_{state}^{(k)}\times {\mathcal{C}}_{cost}^{(q)}$$

(31)

This ratio intuitively represents the effectiveness of the intervention measure in fault clearance, which incorporates the intervention information from the generative SCM. If the current intervention measure $\mathcal{R}$ is of low utility, meaning the intervention input is excessive but the effect is poor, and continuing in the current search direction may be more likely to find an inferior solution. At this point, the coefficient $\alpha$ will increase and change the search direction based on the evaluative information on $\mathcal{R}$ from the generative SCM, thereby helping to escape the local area of inferior solutions. Conversely, when $\mathcal{R}$ has a high utility, it is appropriate to decrease the step size $\mathrm{\Delta }{\mathit{X}}_{iv,t}$ controlled by $\alpha$ conducting more detailed searches in its neighborhood. During the search process, $\alpha$ can adaptively adjust without the need to set values manually. This strategy enhances the ability of the ALNS algorithm to escape local optima through a “feedback regulation” mechanism, thereby increasing the opportunity for finding superior solutions (Algorithm 1).

Algorithm 1: Expected Improvement Search

The process of the ALNS algorithm based on the expected improvement strategy is shown in Figure 3. First, set the initial parameters of the algorithm and generate an intervention ${\mathcal{R}}_0$. Then, one of the operators, which is defined as the destruction operator and reconstruction operator, is selected to iteratively update ${\mathcal{R}}_0$ through the weights of operators, which results in ${\mathcal{R}}_n$ for $n=1,2,\dots$. After the n-th intervention measure ${\mathcal{R}}_n$ is obtained, the following steps are executed:

• Update $SCM \mathcal{M}$ to the counterfactual $SCM \tilde{\mathcal{M}}$ through performing $do(\mathit{P}{\mathit{a}}_G^{(k)}={\mathit{X}}_{iv})$ with observed facta ${\mathit{X}}_F^{(k)}$.
• Calculate the expected step size $\mathrm{\Delta }{\mathit{X}}_{iv}={\left\{\mathrm{\Delta }{\mathit{X}}_{iv,t}\right\}}_{t\in \mathcal{S}}$ and update ${{\mathit{X}}^{\prime }}_{iv}={\mathit{X}}_{iv}+\mathrm{\Delta }{\mathit{X}}_{iv}$ where the objective function $F(\mathcal{R})$ contains the gradient information of the counterfactual generation transformation ${\widetilde{{\mathcal{F}}^{\prime }}}_k$ at ${\mathit{X}}_{iv}$.
• Replace ${\mathit{X}}_{iv}$ with ${{\mathit{X}}^{\prime }}_{iv}$ and obtain new intervention measure ${{\mathcal{R}}^{\prime }}_n$.

Following this, we then compare ${{\mathcal{R}}^{\prime }}_n$ with the intervention from the previous iteration ${{\mathcal{R}}^{\prime }}_{n-1}$ and decide whether to accept the inferior one based on the Metropolis criterion, thereby updating the current intervention and weights of operators for the next iteration.

4. Experimental Results

4.1. Models Comparison

Generating accurate post-intervention counterfactual data is one of the key objectives of generative models based on SCM. To verify the capability of the G-WMGP model to capture the factual distribution, we compare several GP-based baselines described below:

• GP-SCM [28]: A SCM based on a single-task Gaussian process. For each variable $V_i$, a GP-SCM is constructed as follows:

  $${\mathit{x}}^{(i)}=f_i(\mathit{P}{\mathit{a}}_G^{(i)})+{\epsilon }^{(i)}$$

  (32)

  where $f_i\sim \mathcal{G}\mathcal{P}(0,{\mathbf{k}}_{\theta ,i})$. The covariance matrix ${\mathbf{k}}_{\theta ,i}={[{\mathrm{k}}_{\theta ,i}]}_{M\times M}$ is characterized using the RBF kernel.
• G-MTGP [12]: Graph-based multi-task Gaussian process. The relationship between $V_i$ and $P{a}_G^{(i)}$ is similar to Equation (26), but the kernel function is selected as a coregionalization kernel, which is formulated in Equations (3) and (4), where the target GP is considered as a linear combination of independent GPs.
• G-WGP [29]: Graph-based warped Gaussian process. A dual-layer mapping architecture based on a warped transformation, where the transformation is modeled by a normalizing flow similar to the G-WMGP model proposed in this paper, except that the base distribution is represented by a standard GP-SCM.

The reason why we chose these three SCMs as baseline models is as follows: GP-SCM is one of the most classic nonlinear SCMs, the G-MTGP model is the core generative model of the cCBO approach, and the flow-based G-WGP model is also constructed as the generative SCM in solving the algorithmic recourse problem. The G-MTGP model and G-WGP model have been demonstrated to have sufficiently superior performance in the relevant literature through comparisons with other intervention models.

To compare the ability to uncover the underlying distribution and capture nonlinear effects between variables, we use Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and Signal-to-Noise Ratio (SNR) as evaluation metrics. Here, SNR is calculated by the following:

$$SNR={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}10\log \left(\frac{{\sigma }_{s,i}^2}{{\sigma }_{n,i}^2}\right)}$$

(33)

where ${\sigma }_{s,i}^2$ represents the variance of the predicted values for variable $V_i$, ${\sigma }_{n,i}^2$ represents the variance of the prediction errors for variable $V_i$. The purpose of taking the average is to consider the different dimensions of variables. A higher SNR indicates that the model is more effective at extracting useful information from a noisy environment, meaning it is more robust under noisy conditions.

4.1.1. Synthetic Datasets

We simulate some data using the structured synthetic data generation process proposed by Lawrence et al. [35] to test the G-WMGP model, as well as the aforementioned comparative methods. The advantage of using synthetic data lies in the ability to conduct experiments with accurate causal prior information, thereby avoiding the adverse effects that errors in causal knowledge, which are non-critical factors, might have on the experiments.

Before generating the data, it is necessary to specify the causal structure of the variables and initialize the SCM, including the functional dependencies between variables and the noise distribution. First, a nonlinear SCM with additive noise is constructed, where the functional dependencies and noise distribution are randomly selected from a specified collection of distributions, with the assumption that no time-lag effects and no influence from latent factors are considered. Given this SCM, we vary the number of variables $n$ to generate synthetic datasets labeled as Synthetic 1, Synthetic 2, and Synthetic 3. Additionally, to verify whether the length of the fault interval $\mathcal{S}$ has an impact on the performance of these surrogate models, we set different values of $M$, denoted as Reg 1, Reg 2, and Reg 3, respectively. In this paper, $M,n\in [20,30,40]$. The experimental results are presented in

Table 1. Comparison with baseline models on synthetic datasets.

Baseline Models			GP-SCM	G-MTGP	G-WGP	G-WMGP
Dataset	Region	Metrics
Synthetic 1	Reg 1	MAE	0.1941	0.1975	0.1958	0.1912
		RMSE	0.2340	0.2390	0.2365	0.2299
		SNR	6.5542	7.2749	7.4981	7.9175
	Reg 2	MAE	0.1911	0.1836	0.1867	0.1829
		RMSE	0.2307	0.2235	0.2281	0.2212
		SNR	6.5309	7.6280	7.4668	7.8326
	Reg 3	MAE	0.2025	0.2018	0.2002	0.1978
		RMSE	0.2472	0.2487	0.2463	0.2403
		SNR	6.4119	6.7236	7.2403	7.5264
Synthetic 2	Reg 1	MAE	0.2055	0.1952	0.1927	0.1894
		RMSE	0.2498	0.2359	0.2318	0.2302
		SNR	8.3081	10.0011	10.0951	10.2510
	Reg 2	MAE	0.2068	0.1968	0.1952	0.1917
		RMSE	0.2512	0.2385	0.2355	0.2341
		SNR	9.0110	9.7833	9.9564	10.3954
	Reg 3	MAE	0.2068	0.2011	0.2017	0.1994
		RMSE	0.2510	0.2478	0.2490	0.2438
		SNR	9.2907	10.0112	10.0269	10.3497
Synthetic 3	Reg 1	MAE	0.2113	0.1974	0.1976	0.1935
		RMSE	0.2597	0.2401	0.2392	0.2336
		SNR	5.1458	6.7793	6.9453	7.1746
	Reg 2	MAE	0.1959	0.1854	0.1870	0.1857
		RMSE	0.2375	0.2276	0.2280	0.2250
		SNR	6.9043	8.0212	8.0459	8.1312
	Reg 3	MAE	0.2012	0.1947	0.1928	0.1928
		RMSE	0.2463	0.2358	0.2325	0.2332
		SNR	6.7114	7.7421	7.9201	7.9916

, as follows:

Experimental comparisons with several baselines reveal that the G-WMGP model has the smallest MAE and RMSE on Synthetic 1 and Synthetic 2, showing a better performance in learning the data distribution for datasets with fewer variables. As for Synthetic 3, it is observed that the G-WMGP model still outperforms traditional single-mapping GP-SCM, while only showing a distinct advantage over the G-MTGP model and G-WGP model when the value of $M$ is relatively small. However, it is noteworthy that in almost all cases, the SNR values of the G-WMGP model are significantly higher than those of other models, demonstrating the strong robustness of the G-WMGP model in noisy environments. Considering the performance across all these cases, we can see that GP-SCM is relatively the worst, both in terms of prediction and robustness to noise. The G-WGP model exhibits better capability of prediction than the G-MTGP model, while its main advantage lies in its stronger robustness compared to the G-MTGP model. The G-WMGP model proposed in this paper undoubtedly has the best performance in both learning power and noise robustness.

4.1.2. Satellite Attitude Control System Dataset

Considering that synthetic datasets might not adequately reflect the complexity and variability of real-world data, we further verified the effectiveness of the proposed model on a satellite attitude control system dataset. This dataset comprises 68,713 samples, including 19 telemetry variables, with a sampling frequency of 5 min (

Table 2. Information of satellite attitude control system dataset.

Attribute	Value
Number of variables	19
Sampling frequency	5 min
Data length	68,713

). It also contains some fault information, and the intervals corresponding to faulty variables have been labeled, with a label value of 1 indicating the occurrence of a fault and 0 otherwise. The causal relationships among these 19 variables are already known and displayed in Figure 4.

Similarly, we compare the MAE, RMSE, and SNR values of different models on a real dataset, with the results as follows:

From

Table 3. Comparison with baseline models on a real dataset.

Baseline Models		GP-SCM	G-MTGP	G-WGP	G-WMGP
Region	Metrics
Reg 1	MAE	0.8308	0.8023	0.7984	0.7701
	RMSE	1.2641	1.1871	1.1872	1.1235
	SNR	−6.9671	−3.5833	−3.3389	−3.0929
Reg 2	MAE	1.1328	1.0106	1.0302	0.9751
	RMSE	2.9878	2.9163	2.9069	2.9104
	SNR	−1.9956	−1.4988	−1.5842	−1.3890
Reg 3	MAE	0.9988	1.0358	1.0053	0.9602
	RMSE	1.4118	1.4231	1.4594	1.3906
	SNR	0.2838	3.4263	0.9867	0.0655

, we can observe negative SNR values on Reg 1 and Reg 2, indicating that the noise interference in this dataset is more intense compared to the previous synthetic datasets. In such conditions, the G-WMGP model remains reliable in terms of noise robustness, particularly on Reg 3 where its SNR value is only 0.0655. Except for the RMSE value on Reg 2, which is higher than that of the G-WGP model, the other metrics demonstrate the dominant performance of the G-WMGP model on the real dataset, further proving the effectiveness of the proposed model.

4.2. Algorithms Comparison

In this section, we still employ the satellite attitude control system dataset for the comparison of intervention optimization algorithms. We extract the fault information from the dataset, which includes the faulty variables, the start time of the fault, the duration of the fault, and the fault values of the variables. These 19 variables are assigned labels from A to T, and based on the varying lengths of the fault intervals, we categorize the fault scenarios, selecting specific ones as the basis of the intervention experiment. The detailed fault information is as follows (

Table 4. Fault scenarios.

Scenario ID	Faulty Variable	Fault Start Time	Fault Duration
1	I	57,496	25
2	H	25,528	31
3	S	7265	48

):

We first compare the performance of different heuristic algorithms on the intervention optimization problem, including the genetic algorithm (GA), simulated annealing algorithm (SA), and particle swarm optimization algorithm (PSO). We run each algorithm ten times, providing the mean and standard error of the convergence values from these runs as follows:

As shown in

Table 5. Comparison of ALNS, GA, SA, and PSO.

Scenario ID	Objective Value
	ALNS	GA	SA	PSO
1	3.290 ± 2.0 × 10−3	3.298 ± 1.7 × 10−3	3.387 ± 3.3 × 10−3	3.314 ± 2.5 × 10−3
2	4.256 ± 2.1 × 10−3	4.267 ± 2.4 × 10−3	4.273 ± 2.9 × 10−3	4.264 ± 1.9 × 10−3
3	7.774 ± 3.4 × 10−3	7.792 ± 2.8 × 10−3	7.799 ± 3.8 × 10−3	7.790 ± 2.4 × 10−3

, although the performance of the ALNS algorithm may not be as stable as other algorithms, particularly evident in scenario 3, we can still find that the ALNS algorithm has the smallest convergence mean in all three scenarios, demonstrating its superiority in solving this intervention optimization problem.

To further validate the effectiveness of the improved algorithm, we select two related baselines—random minimization baseline and greedy search baseline—and their specific meanings are as follows:

• Random minimization baseline (RM): Given the current intervention $\mathcal{R}$, a random walk is performed $N$ times in a local region. Both the direction and step size are randomly generated. The improved intervention ${\mathcal{R}}^{\prime }$ is chosen such that

  $${\mathcal{R}}^{\prime }=\underset{r\in N}{\arg \min }F({\mathcal{R}}_r)$$

  (34)

  where ${\mathcal{R}}_r\in {\left\{{\mathcal{R}}_i\right\}}_{i=1}^N$.
• Greedy search baseline (GS): The improvement of $\mathcal{R}$ is purely based on the gradient information of the objective function. This involves taking $n_1$ steps in the direction of the gradient descent:

  $$\mathrm{\Delta }{\mathit{X}}_{iv,t}=-{\displaystyle \sum _{m=1}^{n_1}\beta {\nabla }_{x}F{({\mathcal{R}}^{(m)})}_t}$$

  (35)

As the gradient descent algorithm approaches a local optimum, it tends to gradually converge. Therefore, for the optimization strategies that utilize gradient information to improve $\mathcal{R}$, i.e., EI and GS, it is crucial to determine the appropriate number of gradient descent steps $n_0$ and $n_1$, which ensures that significant improvements can be achieved with a relatively low time cost. To illustrate the impact that $n_0$ and $n_1$ have on the performance of the algorithm, we conduct experiments by altering these parameters and observing the changes in the algorithm’s runtime and the values of the objective function.

Taking scenario 1 as an example, as shown in Figure 5. Figure 5a and Figure 5b, respectively, illustrate the performance of the I-ALNS algorithm and GS-ALNS algorithm under different values of $n_0$ and $n_1$. The red bar chart represents the runtime of the algorithm, while the blue curve indicates the objective value calculated by Equation (20) at the maximum number of iterations. It can be observed from Figure 5a that, when $n_0=20$, the optimal objective value that the algorithm can find begins to converge and continues to increase the value of $n_0$, which has a negligible effect on the improvement of the solution. Similarly, in Figure 5b, the blue curve shows an inflection point at $n_1=60$. Therefore, considering both the benefit of finding a better solution and the computational cost of the algorithm, $n_0=20$ and $n_1=60$ are deemed appropriate. For other fault scenarios, the values of $n_0$ and $n_1$ can also be determined using the inflection point method.

To compare the capability of ALNS algorithms based on different search strategies to find high-quality solutions, we apply them to the three scenarios described in Table 3 and compare them with the ALNS algorithm that has no improvement strategy, in which the Adam optimizer with a learning rate of $\beta ={10}^{-2}$ is used under the EI and GS strategy, the number of random walks is given $N=50$ under the RM strategy, and the maximum number of iterations is 200. The objective values of different algorithms are plotted as follows:

As shown in Figure 6, the algorithm based on the RM strategy, due to its inherent uncertainty, may produce inferior solutions at the initial stage. However, when the number of random walks is sufficiently large, it can converge to a comparatively good result. Experimental data indicate that the two methods using gradient information are able to find better solutions than the random baseline and converge to a stable state more quickly, and the EI-ALNS algorithm, represented by the blue curve, consistently outperforms other baselines and achieves the smallest convergence value.

A diagram of the optimal intervention measure when the EI-ALNS algorithm converges on scenario 1 is illustrated in Figure 7. For the causal graph on the left, the yellow node represents the faulty variable, while red nodes denote the causes of the fault, i.e., from which the variables to be intervened are selected. By using the EI-ALNS algorithm, we obtain the optimal intervention measure ${\mathcal{R}}^{\ast }=({\mathit{V}}_{iv}^{\ast },t_{iv}^{\ast },{\mathcal{L}}_{iv}^{\ast },{\mathit{X}}_{iv}^{\ast })$, where the intervened variables ${\mathit{V}}_{iv}^{\ast }$ include variable F, variable H, and variable N, indicated in red. The middle graphs display the intervention values ${\mathit{X}}_{iv}^{\ast }$ for the corresponding variables, which are represented by blue curves, where the red area indicates the intervention interval determined by start time $t_{iv}^{\ast }$ and step length ${\mathcal{L}}_{iv}^{\ast }$. The rightmost graph shows the post-intervention desired values of the faulty variable and is displayed as a red curve. It can be observed that the abnormal condition of variable I within $\mathcal{S}$ is significantly reduced, indicating a notable intervention effect.

Furthermore, considering that the generative model might be influenced by external noise, which could lead to deviations in the experiment results when evaluating the effectiveness of intervention measures, we run each algorithm 10 times across the example scenarios to minimize the interference of extraneous factors and obtain more representative results. The convergence values of these 10 runs for each algorithm are presented in a box plot, as depicted in Figure 8. We can observe that the optimization capability of the method based on the RM strategy exhibits significant instability, with noticeable fluctuations in convergence results, and a higher likelihood of extreme performance. In contrast, the results of the gradient-based methods show less variation with more consistent algorithmic performance, making them more applicable in real-world scenarios. Additionally, we compare the average runtime spent by different algorithms during these 10 runs, as illustrated in Figure 9. It can be seen that the original ALNS algorithm runs the fastest but produces lower-quality solutions. The RM-ALNS and GS-ALNS algorithms improve solution quality at the cost of higher time consumption, while the EI-ALNS algorithm strikes a balance between time cost and solution quality. From this, we can conclude that the EI strategy has a greater advantage over the other two strategies.

5. Discussion

This paper proposes a heuristic method based on a generative SCM and validates the effectiveness of the method on both synthetic data and satellite attitude control system data. The proposed G-WMGP model, constructed as a dual-layer mapping model based on causal coupling knowledge, effectively captures the complex nonlinear relationships when the fault occurs. The dual-layer structure also exhibits a strong capability to capture uncertainties in the data. Thus, we embed the SCM into the ALNS algorithm for optimizing interventions, enabling the generation of the post-intervention desired values. Additionally, we provide the computational process of the counterfactual SCM to assess the effects of performing intervention measures more accurately. In the search process for the optimal intervention, we introduce a search strategy with the purpose of enhancing the optimization capability. To validate the superiority of the proposed method, in Section 4.1, we demonstrate that the G-WMGP model has a strong learning ability and robustness to noise by comparing both synthetic and real datasets with baseline models. In Section 4.2, through comparison with various heuristic algorithms and different search strategies, we argue that the superior performance of the improved ALNS algorithm lies in its ability to strike a balance between time cost and solution quality. Compared to other search strategies, it not only finds better interventions but also requires less search time.

However, when constructing the generative model, this paper considers assumptions of causal stationarity and the absence of latent confounding factors, which may partially limit the applicability of the model as follows:

• The causal stationarity assumption implies that causal relationships in the system remain unchanged over time. Due to the need for counterfactual generation, we only consider scenarios where causal structure changes when the variables in the system are disturbed by external forces, such as interventions. However, in real-world background, causal relationships may not be static. When real-world systems are dynamic and have feedback effects, static SCMs may not be able to capture the complex and dynamic interactions between variables in the system.
• No latent confounding factors assumption requires that all potential confounding variables are identified and controlled, assuming complete causal knowledge of the system. However, even if we believe that all factors related to failures have been observed and controlled, there may still be unidentified confounding variables in reality. Missing confounding factors could potentially mislead us in identifying the true causes of failures, resulting in biases when estimating causal effects, as shown in Figure 1.

The ALNS algorithm is effective in solving the intervention optimization issue, but it still has certain limitations in the following aspects:

• Scalability: When the problem scale significantly increases, the performance of the ALNS algorithm may decline. This is because the algorithm needs to find effective solutions in a larger search space, and each neighborhood operation requires more computing resources and time, which may lead to an inefficient search process. Meanwhile, the ALNS algorithm typically relies on heuristic rules to guide the searching process, which may not be sufficient to capture all important dependencies and interactions in the face of complex causal structures, thus affecting the convergence of the algorithm and the quality of the final solution.
• Computational complexity: when dealing with dynamically changing datasets, the ALNS algorithm may need to frequently update its search strategies or adjust parameters, which could affect the computational efficiency and the feasibility of real-time applications. Additionally, if there are complex causal relationships among variables in the dataset, it may be necessary to design specific neighborhood structures or customized improvement strategies, thereby increasing the complexity of the algorithm.

Future work could focus on developing more reasonable causal models or capturing causal relationships that are dynamically changing. For intervention optimization algorithms, it is important to consider whether the computational efficiency of the algorithms can be improved, for instance, by addressing complex causal structures through the design of neighborhood structures. Additionally, investigating whether these algorithms can be deployed for real-time optimization is also worth exploring.