Failure-averse Active Learning for Physics-constrained Systems

Active learning is a subfield of machine learning that is devised for design and modeling of systems with highly expensive sampling costs. Industrial and engineering systems are generally subject to physics constraints that may induce fatal failures when they are violated, while such constraints are frequently underestimated in active learning. In this paper, we develop a novel active learning method that avoids failures considering implicit physics constraints that govern the system. The proposed approach is driven by two tasks: the safe variance reduction explores the safe region to reduce the variance of the target model, and the safe region expansion aims to extend the explorable region exploiting the probabilistic model of constraints. The global acquisition function is devised to judiciously optimize acquisition functions of two tasks, and its theoretical properties are provided. The proposed method is applied to the composite fuselage assembly process with consideration of material failure using the Tsai-wu criterion, and it is able to achieve zero-failure without the knowledge of explicit failure regions.


Failure-averse Active Learning for
Physics-constrained Systems Cheolhei Lee, Xing Wang, Jianguo Wu, and Xiaowei Yue, Senior Member, IEEE Abstract-Active learning is a subfield of machine learning that is devised for design and modeling of systems with highly expensive sampling costs. Industrial and engineering systems are generally subject to physics constraints that may induce fatal failures when they are violated, while such constraints are frequently underestimated in active learning. In this paper, we develop a novel active learning method that avoids failures considering implicit physics constraints that govern the system. The proposed approach is driven by two tasks: the safe variance reduction explores the safe region to reduce the variance of the target model, and the safe region expansion aims to extend the explorable region exploiting the probabilistic model of constraints. The global acquisition function is devised to judiciously optimize acquisition functions of two tasks, and its theoretical properties are provided. The proposed method is applied to the composite fuselage assembly process with consideration of material failure using the Tsai-wu criterion, and it is able to achieve zero-failure without the knowledge of explicit failure regions.
Index Terms-Active Learning, Physics-integrated Machine Learning, Composite Structures Assembly

I. INTRODUCTION
A CTIVE learning is a subfield of machine learning that maximizes information acquisition to reduce the labelling cost in supervised learning [1]. Contrary to passive learning such as factorial, maximum entropy, Latin Hypercube design (LHD) [2], active learning optimizes the acquisition function that quantifies the potential importance of unlabelled data, and interactively queries the most informative design point to the oracle. Generally, acquisition functions refer to the up-to-date model or labelled data, while there are various strategies that can be adopted according to the preference in information criteria and characteristics of systems. The engineering domain is one of beneficiaries of active learning due to the high complexity of systems and the expensive cost of evaluation [3]- [5].
However, active learning in engineering applications has been mostly utilized without considering coexisting or inherent constraints that may have different processes thereof. It is very Manuscript received 00 00, 2021; revised 00 00, 2021. crucial to consider such constraints in engineering systems, since most of them are subject to physics constraints that may induce fatal and irreversible failures. For example, design of the automatic shape control system in composite aircraft manufacturing need to consider potential material failures such as crack, buckling, and delamination caused by intolerable inputs [6]. The inverse partial differential equation (PDE) problem is another example that is used to calibrate parameters in physics models based on observations. It usually involves physics constraints comprised of different PDEs, and the constraints must be satisfied in order to build first principled models. In both cases, the application of active learning without considering physics constraints may induce fatal failures, so constraints must be considered in order to avoid failures in the application of active learning in engineering problems.
To consider physics constraints in active learning, it is straightforward to define a safe region where design points satisfy safe conditions (with high probability at least), and conduct active learning within the safe region. However, it is not always possible in practice, since physics constraints cannot be explicitly attained due to the complex nature of the system. Examples include the crash damage analysis of commercial vehicles composed of different materials, and physicochemical interactions in corrosion of alloys. In these cases, fundamental physics laws and equations cannot be directly applied or are insufficient to accommodate the complexities of mechanisms. Physics-based numerical methods such as finite difference methods and finite element methods (FEM) [7] are well-established and convincing to analyze large classes of complex structures including failures, while they are too computationally demanding to identify the entire safe region. Moreover, their deterministic solutions are vulnerable to various uncertainty sources such as material properties, geometries, and loads.
In order to circumvent aforementioned limitations of physics-based approaches, machine learning models have been widely used in the engineering domain due to their flexibility, inexpensive prediction, and capability of uncertainty quantification (UQ). Physics information can be highly advantageous for machine learning in several aspects such as generalization and physical consistency [8]. Especially, Gaussian processes (GPs) have shown remarkable performance in stochastic analysis of structural reliability that aims to evaluate the probability of system failure [9]. A common reliability analysis approach employs the GP surrogate model of performance function associated with the system failure, and uses the acquisition function (e.g., [10], [11]) that leads to sampling near the boundary of safe and failure regions. The boundary is called the limit-state, which is the margin of acceptable structural design. However, the reliability analysis is mainly interested in the response surface associated with failure, and it is time-consuming due to the requirement of a large number of samples to estimate the underlying distribution at the limit-state. Hence, it can be data-inefficient to implement reliability analysis prior to the estimation of safe region. So the development of flexible active learning that takes account of target and failure processes simultaneously is promising for systems with implicit physics constraints.
The principle of active learning is exploitation of knowledge from observations, and exploration by tackling the knowledge such as the design point with maximum entropy or the most disagreeable point in the set of hypotheses. However, if unknown constraints exist in the design space, active learning can be very challenging, since the most informative design point may be located in the unknown failure region. Conversely, if active learning is too conservative to avoid failures, the resulted model will be vulnerable to unexplored safe regions. Consequently, active learning for physics-constrained systems should simultaneously take into account the following objectives: 1) maximizing the information acquisition for the target model; 2) expanding explorable safe region by focusing on the constraint function, and they are must be achieved safely. Definitely, two objectives are at odds since they are associated with different functions, so the active learning strategy must be judiciously controlled. In this paper, we propose an active learning methodology for systems that are constrained by unknown failure processes. Our proposed method mainly focuses on the target process, while it simultaneously considers the failure process based on the GP model informed by physics-based models. The GP is adopted to incorporate uncertainty in the prediction of failure model and reduce the computational cost of the physics-based models. Two confidence levels are considered for safe exploration and safe region expansion, and a new acquisition function is developed to integrate two objectives of which the policy can be flexibly adjusted by the decision maker's preference. Our contributions in this paper are as follows.
1) We develop the safe active learning strategy under the regime of implicit physics constraints in the system. 2) A new acquisition function is proposed to enable flexible exploration and exploitation in the constrained design space by referring to two heterogeneous processes. 3) Theoretical properties of the proposed method are provided to guarantee the performance. This paper is organized as follows. In section II, we review literature related to active learning applications in which physics involved and machine learning with unknown constraints. In section III, we elucidate our active learning strategy considering implicit constraints to avoid failures. Section IV illustrates that how the proposed strategy works under the regime of unknown constraints with the simulation study. The real-world application to predictive modeling of composite fuselage deformation considering structural failure is presented in section V. Lastly, a summary of this paper is provided in section VI.

II. LITERATURE REVIEW
In this section, we discuss related literature dichotomizing into (i) active learning for engineering systems; and (ii) machine learning with unknown constraints. While the literature of two topics may not be related, the integration of two approaches can be the cornerstone of our approach. In the application of active learning to engineering systems, we focus on that how physics in the engineering system influences active learning. In machine learning with unknown constraints, we do not restrict the constraints therein being related to physics, yet focus on the ways of considering unknown constraints in machine learning. Since there are not many references whose main contribution lies in active learning, we extend our attention to more general machine learning approaches including the Bayesian optimization, and the inverse design in structural reliability analysis.

A. Active Learning for Engineering Systems
In the engineering domain, active learning is vastly utilized for surrogate modeling of expensive-to-evaluate systems, and the limit-state estimation in the structural reliability analysis. For the surrogate modeling, active learning aims to sample design points that minimizes the generalization error of the model. One of the most applied area is computational physics that often involves surrogate models of systems of which physics-based models are costly or absent. For response surface method (RSM), Alaeddini et al. [12] proposed an active learning strategy adopting the variance reduction of Laplacian regularized parameters. Since RSM is often too restrictive to approximate complex pheomena, deep neural networks (DNNs) and GPs are substantially used for the surrogate modeling of PDE-based models. For DNNs, Costabal et al. [13] proposed a physics-informed neural network for cardiac activation mapping. The neural network model is guided by a loss function that involves the physics equation, one of the most widely used approaches to infuse physics into machine learning. Thereby the resulted model can be physically consistent. In their active learning, they chose the design point with the largest uncertainty quantified by the network model. Lye et al. [14] proposed active learning for surrogate modeling of PDE solutions that chooses the next query point minimizing the cost function with the sequentially updated DNN model. In order to ensure the feasibility, they confined the explorable settings with the feasible region known a priori. Pestourie et al. [15] employed the ensemble neural network model for photonic-device model, and used the quantified uncertainty of network model for active learning. However, their strategy does not consider physically infeasible settings, and they may induce failures if there exist failure conditions in the design space.
For GP models, Yang et al. [16], [17] proposed a physicsinformed GP for the stochastic PDE simulator. The GP model is informed by replicated observations of the PDE simulator, and the predictive variance is referred for active learning. However, their strategy is unconstrained, so the next query point may not satisfy the preliminary conditions. Chen et al. [18] developed the GP that incorporates linear and nonlinear PDE information. They involve the active learning strategy to determine the PDE points for their GP model by employing the integrated mean-squared error (IMSE) criterion [2]. However, the criterion only refers to physics information that does not constrain the learning of the target process, so their approach also may lead to infeasible design points in the constrained system. Yue et al. [4] proposed the variance-based and the Fisher information criteria for GP considering uncertainty, and applied to the modeling of composite fuselage deformation. However, they also did not involve physics constraints during the learning process.
The GP is also frequently employed in structural reliability analysis to model the system reliability, and acquisition functions thereof give more weights on sampling from the vicinity of the limit-state. Echard et al. [11] proposed an acquisition function to estimate the limit-state with GP models, and suggested the framework that encompasses the limitstate estimation and the Monte-Carlo simulation for density estimation. Bichon et al. [10] proposed an acquisition function called expected feasibility function that quantifies the closeness of design points to the limit-state considering quantified uncertainty of the GP. More explicitly, the vicinity of the limit-state was defined with the predictive uncertainty of GP for every candidate design point. Bect et al. [19] proposed the stepwise uncertainty reduction approach that employs the IMSE criterion associated with the estimated safe region. Wang et al. [20] proposed the maximum confidence enhancement method whose acquisition function aggregates the distance to the limit-state, input density, and the predictive uncertainty by multiplying them. The components are the same to [10], while the formulation of acquisition function is different. Sadoughi et al. [21] proposed the dynamically adjustable acquisition function using smooth weighting the closeness unlike to the aforementioned approaches. However, these surrogate-based structural reliability analysis and their variations are free from failures, since they observe systems using simulators. Therefore, their approaches cannot be used in the presence of constraints.

B. Machine Learning with Unknown Constraints
Constraints are critical in every optimization problem that must be satisfied for every solution, although they may be unknown a priori in practice. It is important to consider constraints in machine learning as well, since it is one way to incorporate background domain knowledge and thereby make the model more consistent with reality. There are mainly two distinct manners in considering unknown constraints in machine learning. One is that any evaluation at an infeasible design point is never allowed, and the other is that constraints can be disregarded in the trajectory of solution. We refer the first type as safe exploration since any infeasible design is directly related to the system failure which is critical for continuing the system analysis. The second type is usually ad-dressed in Bayesian optimization, which is a sequential design strategy for the global optimization of black-box functions.
For safe exploration, Schreiter et al. [22] proposed active learning for GPs whose motivation is closest to ours. They used the nuisance functions of GP classifiers for constraints, and used the lower confidence interval to ensure safety during the maximum entropy based active learning. However, their approach only focuses on the entropy of target function, so it could be problematic if the constraint models are lack of accuracy due to the target process-focused samples. Furthermore, the nuisance function may distort the numerical information from the constraint observations that are usually related to closeness to the boundary of the safe region. Sui et al. [23] proposed a Bayesian optimization algorithm with unknown constraints. Their principle is to expand and conquer. That is, they expand the safe region in the first phase, and then the ultimate learning objective is achieved within the disclosed safe region. However, their safe expansion can be subject to slow convergence to the ground-truth unless the Lipschitz constants of unknown constraint functions are known beforehand. Considering the limitation of the safe expansion, Turchetta et al. [24] suggested the safe exploration for interactive machine learning that can be adopted for active learning and Bayesian optimization. Their principle is that the safe region expansion is conducted when the next query position lies in the uncertain safe region. However, their method mainly considers dynamic systems that involve decisions and state transitions, so their approach is unsuitable to ours.
In Bayesian optimization, the GP is mostly considered for the approximation of black-box function, and one of the most widely used acquisition function is the expected improvement (EI) function. The function is straightforwardly formulated upon tractable Gaussian distribution of GP models: it chooses the point with the greatest expectation of improvement in the optimization solution considering uncertainty. Schonlau et al. [25] multiplied the probability of constraint satisfaction (also called expected feasibility) to the EI function. Gelbart et al. [26] employs the latent GP models of unknown constraints that may not be observable simultaneously with the target process. They also used the EI function weighted by expected feasibility induced by the latent GP model. Gramacy et al. [27] assumed infeasible designs also can be informative to the solution of problem, so they proposed the integrated expected conditional improvement that is also weighted by expected feasibility. Later, the augmented Lagrangian approach coupled with the EI function was proposed in [28]. The Lagrangian penalty is optimized in the outer loop of the algorithm so that the solution can be led to the feasible region. Picheny [29] suggested to choose the design point that induces the largest expected volume reduction, quantified by integrating the production of expected feasibility and probability of the improvement function. Hernández et al. [30] proposed the predictive entropy search with constraints that refers to the expected entropy reduction at the minimum associated with observations from the objective function and constraints. It automatically focuses one of objective and constraints by merging them into the integrated entropy. Basudhar et al. [31] used the probabilistic support vector classifier to discriminate the safe region. Then, the probability of constraint satisfaction is used to weight the EI function as well. However, the aforementioned approaches for unknown constraints do not enforce sampling from the safe region, so they are incapable of safe exploration.

III. METHODOLOGY
In this section, we propose our failure-averse active learning method for physics-constrained systems, the overview of which is shown in Fig. 1. The system of interest contains the target function, and there are governing constraints. The predictive modeling of target function is initialized with a passive design (e.g., space-filling), and the constrained model is introduced to inform the active learning module comprised of two objectives: the safe variance reduction and the safe region expansion. Then, both models are updated sequentially with data queried by the active learning strategy and labelled by the oracle. We begin with specifying our problem whose constraints can be evaluated along with the target process. Then, the safe exploration incorporating constraints is described in detail with its theoretical properties.

A. Problem Statement and Gaussian Process Priors
Consider a system defined over a compact and connected design space Ω ⊆ R D . The system includes the target function f : Ω → R, we want to predict, and the constraint function h : Ω → R, related to the system failure and assumed to be independent of f . Let ξ ∈ R be a tolerable failure threshold associated with h, which should be defined conservatively considering intrinsic uncertainty of the process. For any design point x = [ x 1 · · · x D ] ∈ Ω, the system failure occurs when and safe otherwise. For example, f can be the dimensional deformation of a solid structure given the force vector x, and h can be the resulted von Mises stress in the structure. Both functions can be observed or evaluated with the costly real experiment and relevant physics-based models as We assume that ν f and ν h are assumed to be known, while they can be also estimated with GP priors on f and h to be described later with nugget effects. The safe region Ω S , the subset of design space comprised of non-failure design settings, is unknown and difficult to obtain due to the prohibitive cost of the evaluation. We refer the complementary of safe region as the failure region such that Ω F = Ω \ Ω S .
In our case, due to the high cost of evaluation of h and f , we prefer to sample the most informative set of design points that minimizes the generalization error associated with the target function. Suppose we have n samples from the system, , andf n andĥ n are our predictors of f and h trained with D n , respectively. Then, the expected risk off n using L p loss is where L p (y, y ) = |y − y | p is the loss function, and λ(x) is a probability measure defined over Ω. In our problem, any violation of (1) may incur prohibitive cost of failure in the system, so any x ∈ Ω F will not be considered for the system.
We assume that there exists a reproducing kernel Hilbert space (RKHS) for each of f and h, and they are bounded therein. It allows us to model both functions with GPs with corresponding kernels such that k f : Ω 2 → R and k h : Ω 2 → R [32]. In this paper, we consider the automatic relevance determination using radial basis function (RBF) kernel for x, x ∈ Ω as where κ f is the nonnegative scale hyperparameter, M f is the diagonal matrix of nonnegative length hyperparameters θ f = [ θ f,1 , . . . , θ f,D ] , and δ is the Kronecker delta function for the nugget effect. By defining k h in the same manner, we can write f and h as where µ f and µ h are mean functions assumed to be zero without loss of generality. Instead of employing a discriminative function for estimating the safe region, the GP regressor is more suitable for physics constraints since the output of h is numerically informative. More explicitly, as h(x) is closer to the failure threshold, we may notice that x is closer to the safe boundary. Moreover, discriminative functions require observations from both safe and failure regions, while regressors are not subject to such imbalance or absence of one class. Therefore, we fit our GP regression model directly on observed outputs from h, and refer to the distance between the output and the failure threshold to infer the probability of safety.
Let us denote X n as the D×n design matrix of [ x 1 · · · x n ], and y n and z n as the vector of n observations from f and h, respectively. With GP priors on f and h, the hyperparameters Θ f = {κ f , θ f } and Θ h = {κ h , θ h } can be estimated by maximizing the log marginal likelihoods, which are where K f,n and K h,n are covariance matrices comprised of every pair of x, x ∈ X n given Θ f and Θ h , respectively. Oncê f n andĥ n are obtained with maximizing their log marginal likelihoods, the predictive mean and variance off n at an unobserved design point x ∈ Ω can be derived as and so doesĥ n 's.

B. Safe Variance Reduction
Let us consider L 2 loss called the mean squared error (MSE), while it is not required in practice for our approach. The IMSE criterion, we also adopt partially in our method, can be derived from the risk with L 2 loss [33]. Suppose f * is an unbiased predictor of f with the minimum MSE (also called the best MSE predictor [2]) with respect to Ω S in the family of GP we consider. Then, the L 2 risk off n can be decomposed as which is the sum of the L 2 risk of f * and the variance off n . Since the L 2 risk of f * is negligible due to its unbiasedness, the generalization error with respect to the safe region can be reduced by focusing on the variance reduction inf n . The IMSE criterion seeks the design point that is expected to reduce the variance of our predictor the most based on the realizedf n . Let us denote the integrated variance of predictor in (2) as where s ∈ Ω S . Then, for an unobserved x ∈ Ω, the expected variance reduction over Ω S inf n given x is where X n+1 = [ X n x ], and K f,n+1 is the covariance matrix of X n+1 . Eq. (3) is always nonnegative (Theorem 1 and 2 in [34]) and V Ω S (f n ) is invariant to x, so it suffices to maximize the second term. For the GP, choosing the next design point that minimizes (3) is called active learning Cohn (ALC), and it is widely used along with the active learning Mckay (ALM) that chooses the next query point with maximum entropy [35]. In terms of computational cost, the IMSE criterion is dominated by the inversion of K f,n+1 , so we adopt the Cholesky update for IMSE criterion in [36] to alleviate the cost and improve the numerical stability. Unfortunately, (3) cannot be used directly, since the safe region is unknown a priori. Thus, we stick to our predictorĥ n to estimate the safe region as whereμ h n (x) andσ h n (x) are the mean and standard deviation ofĥ n (x), and Φ(β n ) = Pr(ĥ n (x) < ξ) = 1 − p n for which p n ∈ (0, 1). That is, β n is related to the probability of failure of x ∈ S n as shown in the following theorem.
Theorem 1 (Failure Probability). Suppose we have Nsampling budget andĥ n . For γ ∈ (0, 1), choosing design points x i 's from the safe region S i for i ∈ {n+1, . . . , n+N } has the failure probability as The proof is provided in Appendix A. By constraining our choice of next design point x ∈ S n , (3) can be written as where s ∈ S n . The first term of (6) indicates the irreducible variance induced by discrepancy between Ω S and S n . Meanwhile, the second term is the reducible variance in the estimated safe region. Eq. (6) shows that the choice of p n for S n affects the learnability of active learning and the safety. More explicitly, S n needs to be conservative to prevent failure by setting p n small enough, while too conservative setting of S n will increase the irreducible variance term in (6) and reduce the explorable region.
To extend the purview of variance reduction by x in (6) while keeping safety, we may consider another safe region called the progressive safe region that has a more generous safety level than S n as S + n = {x ∈ Ω|μ h n (x) + β + nσ h n (x) < ξ}, where Φ(β + n ) = Pr(ĥ n (x) < ξ) > 1 − p + of which p + n < p n . Straightforwardly, we have S n ⊂ S + n ⊆ Ω. By considering S + n as the reference set for the integrand (4), i.e., s ∈ S + n , we can reduce the irreducible variance in the first term. Consequently, we have the following acquisition function where x ∈ S n . While the increment of progressive safe region will reduce the unconsidered safe region Ω S \ S + n , it is noteworthy that it can also simultaneously increase the reducible variance in Ω F \ S + n , which is meaningless. Therefore, S + n also need to be defined not too generous.
Even though J f is designed carefully with appropriate safe regions, we can minimize the irreducible variance by minimizing the discrepancy betweenĥ n and Ω S . Generally speaking, J f only focuses on reducing the variance off n , and does not care about reducing the discrepancy. In order to efficiently expand the explorable region and improve the estimation accuracy of safe region, we need to incorporate the information from h as well as f in the information criterion. In the following section, we illustrate the safe region expansion that focuses on the estimation of safe region boundary.

C. Safe Region Expansion
The safe region expansion is required to reduce the error induced by the mismatch of S n and Ω S , and furnish higher confidence in exploration. It is metaphorically similar to that we can win when we know much more about the opposite. To expand the safe region without failure, we need to exploit the numerically informative output of h to approach to the boundary of safe region from inside thereof, and the expansion can be maximized when the design point is closest to the boundary [37]. Inspired by [38], we incorporate the uncertainty ofĥ n and closeness to the boundary with the following criterion: where η n (x) = ασ h n (x) of which α > 0. I(x) attains its maximum whenĥ n (x) = ξ, which is the case of x ∈ ∂Ω S . Otherwise, it gets additional scores whenĥ n (x) does not exceed the threshold within an acceptable interval. The role of α is magnifying the effect of uncertainty in I(x). Let the expected value of I(x) with respect toĥ n be J h (x), which is the acquisition function for the safe region expansion expressed as whereΦ h n andφ h n are the CDF and the PDF ofĥ n (x). Eq. (9) is composed of two terms: the first term is related to the uncertainty, and the second term is related to the closeness to the boundary. Consequently, maximizing (9) leads to sampling near the boundary with high uncertainty if such points exist, and the most uncertain point otherwise.
Obviously, interests of J f and J h are inherently different, since they are associated with different mechanisms, f and h. Also, they are formulated for different purposes. It implies that the safe approximation of target function and the safe region expansion have trade-off, thus we need to compromise both criteria to determine the most informative design point. We discuss the framework for addressing the balancing between two criteria in the next section.

D. Harmonizing Acquisition Functions
In this section, we integrate two acquisition functions to optimize (maximize) them judiciously to achieve the safe active learning. Two distinct acquisition functions are proposed to accomplish different objectives, and the optimization of two criteria is a multi-objective optimization (MOO) problem. Conceptually, we may think of a point x ∈ S n that achieves the maximum of each criterion simultaneously, called as the utopia point. However, MOO typically has no single optimal solution contrary to the usual single-objective optimization. Therefore, the Pareto optimality concept is mostly referred to define the optimality in this regime. Let Obviously, every Pareto optimal point is weakly Pareto optimal, while the reverse is not true.
It is common to scalarize the vector-valued objective functions in MOO, and the formulation of problem is critical for Pareto optimality of the solution. In this paper, we use the weighted sum, which is widely used, to scalarize our criteria with the integrated acquisition function: where p ∈ N. The weight parameter w in (10) exactly conveys the preference of decision maker between two objectives. For example, if one is more interested in the safe region expansion, the decision maker will weigh more on J h , and decrease w when the estimated safe region seems acceptable. Otherwise, some may begin with small w to see if the safe region expansion is necessary. It turns out that the integrated acquisition function guarantees the Pareto optimality of its solution given w as shown in the following proposition.
Proposition 1 (Sufficient Pareto Optimality, [39]). For any w ∈ (0, 1), a solution that maximizes the integrated criterion is Pareto optimal associated with w. When w = 0 or w = 1, a solution of the integrated criterion is weakly Pareto optimal associated with w.
However, J f and J h may be different in their scales, so the weight parameter w cannot be determined straightforwardly. It is common to transform component objective functions in the formulation of MOO, so we may normalize each criterion as where minimum and maximum of J i for i ∈ {f, h} stand for the minimum and maximum over S n . To normalize objective functions, we need their maxima and minima, so we provide two scaling options in this paper. The first scaling Evaluate S n , S + n over Ω 5: x * = arg max x∈Sn J(x) 6: Observe y * , z * at x * 7: Updatef ,ĥ with D 10: Update β, β + α, w (Optional) 11: end while method is for global searching in a dense-grid over the design space. By discretizing the design space into a dense-grid, we may evaluate all criteria over the grid. It can provide the heuristically global optimal solution, and make scaling more consistent. Another scaling method is to incorporate the lower and upper bounds of criteria. We already have that both criteria are nonnegative, so they are lower bounded by zero. For the upper bounds, let us consider the progressive safe region J f which is upper bounded by since the expected variance reduction in (7)  from its original formulation. In this way, we can scale both criteria by their tractable bounds that can be obtained prior to the evaluation of every candidate. The upper bounds of J f and J h can be referred to the asymptotic convergence of the integrated acquisition function in (10) as described in Theorem 2, which is an important property of the proposed method.
Theorem 2 (Asymptotic Convergence). The integrated criterion converges to zero as the number of samples n → ∞, which means the variance of target predictor is minimized, and the safe region is fully expanded.
The proof is given in Appendix B. Theorem 2 implies that the integrated criterion J leads to the best estimator of f and promising conservative estimation of h. Let us refer to our active learning as PhysCAL (Physics-Constrained Active Learning), and its pseudocode is provided in Algorithm 1. Note that PhysCAL can be terminated by not only the sampling budget, but also the prediction accuracy of the target model when the sampling budget is implicit or early stopping is reasonable. In order to do so, a separated testing dataset or cross-validation is required.

IV. SIMULATION STUDY
In this section, we apply our active learning to the approximation of a constrained 2-D simulation function. The response surface of the target function is shown in the left top of Fig. 3, and the constraint function's response surface and the failure region are in the left bottom thereof. We set the safe region as Ω S = {x ∈ [0, 1] 2 |h(x) < 0.7}, thereby the failure region ratio to the design space being approximately 0.28. Assuming we have no prior knowledge of the safe design settings, 10 initial samples are obtained over the design space using the maximin LHD, which yields 2-3 samples from the failure region in 10 replications. Observations from f and h are corrupted by Gaussian noise, and additional 20 samples are obtained with active learning. According to Theorem 1 and the sampling budget, we set p = 0.001/20 to retain the safety, and p + = 0.01 without update. As the benchmark method, the safe exploration for GP in [22] (referred as SEGP) is considered. For the other parameter settings, α = 2 and w = {0.0, 0.1, . . . , 0.9, 1.0} are used.
In the simulation study, setting the safety level in the benchmark method as high as ours is impossible due to the insufficient number of samples to make a nontrivial explorable region as they mentioned in their work. Moreover, the benchmark method fails to estimate a plausible failure region, thereby resulting frequent failures as shown in the right of Fig.  3. Meanwhile, our method estimated the failure region much better for every instance, and explored more safely as shown in the center of Fig. 3.
As a result, PhysCAL achieved the prediction error (MSE) 0.0023 with 0.5 additional number of failures on average, and there are five zero-failure cases of ten replications. Meanwhile, the benchmark method achieved a better prediction accuracy with an MSE of 0.001, but the averaged additional failures was 6.7 among 20 queries with no case of zero-failure. It is not surprising that the benchmark method did better in predictive accuracy, since exploration was unrestricted by underestimating the failure region. Hence, in the case that one failure is very crucial, our approach will be more suitable. Fig. 2 shows the performance of PhysCAL with different weight parameter settings. We can observe that a low weight parameter may improve the prediction accuracy, but does not necessarily make our approach more conservative due to lack of knowledge in failure region. Increasing w may compel our exploration being more aggressive with an increased number of failures, and induce low predictive accuracy due to indifference to the variance reduction in the target approximation. In this simulation, w = 0.4 was the best choice among considered values with promising predictive accuracy and the least number of failures.

V. CASE STUDY
In this section, the proposed method is applied to the predictive modeling of deformation of the composite fuselage assembly in aerospace manufacturing. Composite materials such as carbon fiber reinforced polymers are extensively applied to various domains (including aerospace, automotive, construction and energy) due to its versatility, high strengthto-weight ratio, and corrosion-resistance. However, composite materials are nonlinear and anisotropic due to their structural natures [40], so flexible models and adaptive design of experiments are required to predict its deformation precisely. Moreover, since they are also subject to structural failures, manufacturers should avoid unsafe load settings in the fabrication of composite structures.
In aircraft manufacturing, each fuselage section is subject to dimensional deviations at the joint rim due to the multi-batch manufacturing system. Consequently, to assemble fuselage sections, the shape control procedure is required to reshape them homogeneously. In the shape control procedure, the fuselage section is placed on the supporting fixture, and ten equispaced hydraulic actuators are introduced to reshape the fuselage as shown in the left figure of Fig. 4. A promising approach for the optimal shape control procedure is to employ a highly precise predictive model of composite fuselage's deformation instead of time-consuming physics-based model [41], [42]. However, the construction of such model validated with real experiment is very challenging due to the expensive cost of sampling, and the risk of structural failure due to unsafe loads.
There are different types of failures in composite materials such as fiber, matrix, and ply failures. Likewise, a number of composite material failure criteria (e.g., Tsai-wu, Tsai-hill, Hoffman, Hashin) are devised for different modes of failures [43]. In this paper, we considered the Tsai-wu criterion, which is one of the most widely used interactive failure criterion. Note that it is possible to consider more than one criterion simultaneously by taking the most parsimonious criterion, or considering the intersection of safe regions defined by multiple constraint GP models. Briefly, Tsai-wu criterion considers interactions between different stress components in addition to the principal stresses (in a homogeneous element). Using the principal material coordinate system on the cubic element of composite material, consider three directions: 1 is the fiber direction; and 2 and 3 are directions perpendicular to 1, respectively. Let σ T i and σ C i be the tensile failure stress and the compressive failure stress in i ∈ {1, 2, 3} direction, and τ F 12 be the shear failure stress in the 12 plain. The Tsai-wu criterion is defined as, where the left-hand side is the nonnegative criterion value, and the failure occurs when it exceeds one. The Tsai-wu criterion value induced by the shape adjustment solved by the FEM is shown in the right of Fig. 4. We can observe that failures are occurred at the bottom of the fuselage, since fixtures that sustain the fuselage are restricting its deformation. Not only limited to the Tsai-wu criterion, physics constraints have many assumptions such as homogeneity, absence of higher-order interactions, etc., although they are convincing apparatuses to consider the structural reliability. Hence, they are typically utilized with the safety-of-margin (the reciprocal of acceptable failure criterion) or UQ to prevent unexpected failures. Likewise, the safe shape control system should consider the failure criterion not only its value, but also the additional safety measures.

A. Experiment Settings
The proposed methodology is applied to the construction of predictive model of deformation. A well-calibrated FEM simulator of the procedure is referred to as our oracle considering the risk of real experiment. The target function's input is the vector of unidirectional forces (in lbf) of ten actuators, and the output is Y, Z-directional deformations (in microinch) of fuselage at one of 91 critical points around the rim. The maximum magnitude of actuator force is 1,000 lbs, which may cause failures in the structure (see Fig. 4). As our additional safety measure, the margin of safety with Tsai-wu criterion is set with 1.25 (i.e., the acceptable criterion is 0.8).
For the initial design, the maximin LHD is used to generate 20 observations for both deformation and failure criterion, and additional 20 samples are queried by different methods: random, ALM, ALC, SEGP, and PhysCAL. We adopted the pool-based scenario in this case study by providing the 400 size of candidate pool that uniformly spreads out the design space. Considering the variability in the initial design, we have generated ten initial dataset independently, and replicated the experiment.
For PhysCAL, we also considered different weight parameters as the simulation study, and set other parameters as α = 2, p = 0.001/20, and p + = 0.1. In SEGP, we reduced the safety level of which from the PhysCAL's until that SEGP induces a nonempty explorable space. For the model evaluation, we used 100 safe samples as the testing dataset that is independently generated with the candidate pool, and the mean absolute error (MAE) is used as the metric.

B. Result
The result is summarized in Table I. First, we can observe that PhysCAL outperforms other methods in the number of additional failures. It achieves nine zero-failures from ten cases. Also, we can observe that ALM, ALC, and SEGP incurred more failures than the random. The reason is that the design space is almost dominated by the safe region, but the failure region may be more interesting than elsewhere. Interestingly, SEGP is inapplicable in this case when the initial dataset does not contain failure samples, since the method uses the binary classifier. It implies that employing the regressor as the constraint model is more advantageous when the prior information has no failure.
In terms of prediction accuracy, PhysCAL performs much better than random sampling, and comparable to other active learning approaches considering the scale of metric. We can conjecture that other methods are able to observe from the failure region that may be informative, so their accuracy is the consequence of unsafe exploration. Furthermore, PhysCAL is more flexible than other methods, since we may update the weight parameter in PhysCAL during data acquisition to focus more on the variance reduction as well as other methods. Therefore, PhysCAL is more promising for this case considering the risk of failure in the system. 1) Weight Parameter: Different weight parameters are considered for PhysCAL in the case study, and the performance of different weight parameter settings is provided in Fig. 5. In terms of safety, increasing weight parameters can improve the safety of active learning in this case. That is, even though we set S n extremely safe as in this case study, the safe region expansion is critical to the safe exploration. However, we also observe that setting w = 1.0 also increases the number of failure slightly, so we should not make PhysCAL too aggressive by only focusing on the safe region expansion. For the predictive accuracy, we also observe that any extreme setting is not an optimal solution.
There are more related factors that may affect the significance of weight parameter. The failure region in this case is much smaller than the safe region, thus quite aggressive exploration may be acceptable (i.e., increasing w). In the simulation case, the failure region ratio is higher than this case, so w = 0.4, which is smaller than the best weight value of the case study, perform better. Also, the correlation between target and constraint functions can be a reason. More explicitly, we know that the safe region expansion term includes the predictive uncertainty term. Thus, when f and h are correlated, the safe region expansion term also can be helpful to improve the predictive accuracy of target model. In this case, increasing w is relevant.

C. Margin of Safety
In order to observe the effect of margin of safety, we increased the value from 1.25 to 1.5, which reduces the acceptable Tsai-wu criterion to 0.66. Although we may consider margin of safety higher than 1.5, such high value is irrelevant in practice. Increasing the margin of safety yields the increment of failure region, so it results with the increased numbers of failures. In average of 10 replications, PhysCAL got 0.5 additional failures (six zero-failures), which is the least among considered methods, and achieved the predictive accuracy of 3.74 microinches. Random sampling got 0.6 additional failures, and achieved the predictive error of 5.16 microinches. Meanwhile, ALM and ALC achieved 8.3 and 7.9 failures, respectively. Likewise, SEGP was applicable only for the last three cases, and induced 3.3 failures from those cases.

VI. SUMMARY
For physics-constrained systems that are expensive-toevaluate, failure-averse active learning is proposed in this paper. In order to achieve safe active learning under the regime of implicit physics constraints, GP priors are imposed on the target function and physics constraints, and two acquisition functions are developed for the safe variance reduction and the safe region expansion. For the safe variance reduction, two safe regions with different safety levels are employed in the IMSE criterion, thereby maximizing the safety and variance reduction over the underlying safe region. For the safe region expansion, the acquisition function is devised to sample near the safe region boundary considering uncertainty. Two acquisition functions are endowed with different objectives, so MOO framework with Pareto optimality is applied to integrate them into the flexible global criterion. The integrated acquisition function is sufficient for the Pareto optimality of the design point to be queried, and can be flexibly adjusted by decision maker's preference considering the trade-off between two acquisition functions. Also, it is shown that the integrated acquisition function asymptotically leads to the best estimation of system.
In the simulation study, the proposed approach showed promising performance with achievement of zero-failure, while the benchmark method fails to avoid failures in its learning process. Furthermore, with different parameters settings, we empirically observed that the safe variance reduction and the safe region expansion should be involved simultaneously for better predictive accuracy and higher safety. Our method has also shown remarkable performance in the case study, the predictive modeling of composite fuselage deformation considering its structural failure with Tsai-wu criterion. It achieved zero-failure in most cases, while other benchmark methods inducing more failures or inferior predictive accuracy. Our proposed method is adaptive, since it can incorporate domain knowledge and decision maker's preference with amenable parameters. Therefore, it is also applicable to other domains that are subject to implicit failures. APPENDIX A PROOF OF THEOREM 1 For any x i ∈ S i , the probability of failure is and it holds for ∀ i. Then, by Boole's inequality (union bound), we have Suppose we use p i = p for ∀ i consistently. Then, the failure probability can be simplified as i p i−1 = N p γ =⇒ p = γ N .
APPENDIX B PROOF OF THEOREM 2 First, we show that S n converges to a nonempty subset of Ω S an n → ∞ by the consistency of GP regression over Ω. Lemma 1. Let S 0 be the initial estimated safe region that contains a safe design point at least. With any dense sequence of design points in S n including PhysCAL, S n converges to a nonempty S * ⊆ Ω S .
Proof. Let us denote a union of estimated safe regions as which is bounded by Ω. Suppose any dense sequence of design points {x 1 , x 2 , . . . , x n } ⊂ T n . Since the sequence is monotone, and bounded by Ω, it converges almost surely to Ω as n → ∞. Note that the collection of dense sequences of design points includes PhysCAL, so the consistency of GP regression holds over Ω [44]. Thus,f n andĥ n converge almost surely to the best MSE predictorf * andĥ * associated with Ω. Sinceĥ * is unbiased and has zero-variance over Ω, for a large enough n, we have that S n = {x ∈ Ω|μ h * (x) < ξ} = {x ∈ Ω| h(x) < ξ} S * ⊆ Ω S , whereμ h * is the mean function ofĥ * .
Lemma 1 allows us to utilize the consistency of GP models to show that degeneracy of J(x). Recall that both J f and J h are nonnegative, and they are upper bounded by for any x ∈ Ω. Then, by the consistency of GP, we have Thus, by the squeeze theorem, J(x) converges to zero for every x ∈ Ω.