Regularized Operating Envelope with Interpretability and Implementability Constraints

Let $Y$ denote the response variable, which represents the KPI of interest in our operating envelope problem. The state variables that affect the value of $Y$, i.e., the input in machine learning language, are denoted by a $p-$dimensional vector $\mathbf{X}=[X_{1},...,X_{j},...,X_{p}]^{T}$. Denote the underlying density distribution of $\mathbf{X}$ as $g(\mathbf{X})$, which is unknown. To simplify our notations, we mainly focus on scenarios where the $p$ state variables are numerical variables in this paper. The proposed problem formulation and the corresponding solutions can be extended to handle categorical state variables with mild modifications. The space spanned by $\mathbf{X}$ is denoted as $\Omega(\mathbf{X})$. Suppose that there is an underlying mapping from $\mathbf{X}$ to the response variable $Y$, i.e.,

where $\epsilon$ represents the combination of any other factors that is not in $\mathbf{X}$ and the random disturbance. We assume that the mapping has been normalized such that the expectation of $\epsilon$ is 0 and the variance is $\sigma_{\epsilon}^{2}$, and $\epsilon$ is independent of $\mathbf{X}$. Neither the structure of $f(\cdot)$ nor $\sigma_{\epsilon}^{2}$ is known beforehand. Instead, we have access to $n$ independent pairs of samples, $(y_{i},\mathbf{x}_{i})$, with $\mathbf{x}_{i}=[x_{1}^{(i)},...,x_{j}^{(i)},...,x_{p}^{(i)}]^{T}$ for $i=1,...,n$.

The ultimate goal of our operating envelope problem is to identify an optimal region in the state variable space $\Omega(\mathbf{X})$ such that the response variable is maximized, and the detected region is interpretable and implementable by making use of the observed data points $(y_{i},\mathbf{x}_{i})$, $i=1,...,n$. Before presenting the proposed formulation for our new operating envelope problem, we first describe the concepts of interpretability and implementability of any region $R$ in the state variable space.

To enhance the interpretability of the achieved envelope, we propose to define eligible candidates as regions of certain simple geometrical shapes. Specifically, we assume that the qualified regions are $p$-sided hyper-rectangles with sides being parallel to the $p$ state variable axes or unions of $L$ ($L<\infty$) disjoint such hyper-rectangles. The mathematical representations of this type of hyper-rectangles is provided in Definition 1. Hyper-rectangles that are parallel to the orthogonal state variable axes are essentially the Cartesian product of intervals on each axis. The identified optimal hyper-rectangle explicitly indicates within which interval each state variable should be, increasing the explainability of data-driven results to operators.

For any region $R$ defined above, we propose to use the probability that state variable $\mathbf{X}$ falling into it to quantify its implementability, i.e. $P_{R}=P(\mathbf{X}\in R)$. A larger value of $P_{R}$ indicates that the corresponding region has a better implementability. To ensure the implementability in the output region, we introduce a lower bound $\beta\in[0,1)$ and define implementable regions as all regions that satisfy the coverage constraint

In Eq. (2), higher values of $\beta$ place higher requirements on the implementability. If any region size is implementable, $\beta$ can be set as 0. $\beta$ is a hyperparameter that needs to be specified in advance according to the practical needs.

Based on the described interpretability and implementability concepts, the proposed operating envelope problem is formally defined in Definition 1 below.

In Definition 1, Eq. (3), the main objective function, indicates that the goal of our proposed operating envelope is to maximizing the region-wise response variable over $2pL$ region-characterizing factors, i.e, $r_{l1}^{l},r_{l1}^{u},...,r_{lp}^{l},r_{lp}^{u}$, for $l=1,...,L$. To the best of our knowledge, this is the first formal definition for the operating envelope identification problem that directly considers the expected value of response variable $Y$ within regions as the objective function. In previous literature on data-driven operating envelope problems, machine learning classifiers treat the accuracy measures for predicting the labels ‘$Y\geq$ cutoff’ and ‘$Y<$ cutoff’[18, 19] as the objective function. The proposed definition also originally incorporates customizable interpretability and implementability constraints by Eq. (4) and Eq. (5) to account for important practical requirements. The desired operating envelope that is a union of $L=3$ disjoint rectangles within a space spanned by $p=2$ covariates is given in Figure 1.

Note that our proposed definition can be generalized to handle candidate regions with other shapes that can be characterized by a finite set of parameters, such as hyperspheres. In Definition 1, $L$ is a hyperparameter that indicates the number of desired disjoint hyper-rectangles in the output. Specifying multiple disjoint sub-regions in the output is sometimes important. For instance, the operations may be conducted at multiple different operating modes. For each mode, we need to specify an optimal operating region accordingly. This hyperparameter is optional for our constrained operating envelope problem. When it is not specified, we can output the operating envelope and the corresponding optimal mean response variable under different $L$, and the user can choose among the results based on their own needs.

In the next section, we describe step by step how we use the sample data $(y_{i},\mathbf{x}_{i})$, $i=1,...,n$, to obtain the optimal envelope defined in Definition 1.

A natural way of evaluating the objective function as well as the probabilistic coverage for any region $R$ could be that we first estimate $f(\mathbf{X})$ by $\hat{f}(\mathbf{X})$, for example, through regression models, and estimate $g(\mathbf{X})$, for example, via the kernel density estimation $\hat{g}(\mathbf{X})$. Then we estimate the objective function and the probabilistic coverage by plugging in $\hat{f}(\mathbf{X})$ and $\hat{g}(\mathbf{X})$ to the objective function and the probabilistic coverage, respectively. Specifically, the estimates are

and

where $\hat{\int_{R}}$ represents the numerical integration using techniques reviewed by [20]. However, this function approximation first and then numerical integration method is computationally heavy, especially in the high dimensional data cases.

To overcome the calculation challenge, we proposed an alternative model-free approach. Based on the observed samples $(y_{i},\mathbf{x}_{i})$, $i=1,...,n$, for any given region $R$, we propose to use the sample average approximation (SAA) to estimate $E[Y|\mathbf{X}\in R]$. Let $I(U)$ be an indicator function of event $U$, i.e.,

then the proposed estimator for $E[Y|\mathbf{X}\in R]$ is

For any given region $R$, we propose to estimate the probabilistic coverage $P_{R}=P(\mathbf{X}\in R)$ by the sample proportion estimator. Mathematically,

The model-free sample average estimator and the sample proportion estimator in Eq. (9) and Eq. (10) are efficient to calculate. Our simulation studies also indicate that they yield accurate estimates for $E[Y|\mathbf{X}\in R]$ and $P_{R}=P(\mathbf{X}\in R)$, respectively, as long as the sample data is representative over the entire space $\Omega(\mathbf{X})$.

When the raw data $(y_{i},\mathbf{x}_{i})$, $i=1,...,n$ is not representative over the entire space. We propose to improve data’s representability by first generating more data from the learned models $\hat{f}(\mathbf{X})$ and $\hat{g}(\mathbf{X})$ before using Eq. (9) and Eq. (10) to estimate the required components. We focus on the model-free estimation approach given in Eq. (9) and Eq. (10) in this paper.

Based on the estimation module in the previous section, empirically, for any given $L$ and $\beta$, we need to calculate the optima for the following constrained optimization problem,

subject to

where $R=R_{1}\bigcup...R_{l}...\bigcup R_{L}$, with $R_{l}=[r_{l1}^{l},r_{l1}^{u}]\times...\times[r_{lj}^{l},r_{lj}^{u}]\times...\times[r_{lp}^{l},r_{lp}^{u}]$, $r_{lj}^{l},r_{lj}^{u}\in\mathbb{R}$ and $0<r_{lj}^{l}<r_{lj}^{u}<\infty$. $\forall l,l^{\prime}\in\{1,...,L\}$ and $l\neq l^{\prime}$, there exist $j\in\{1,...,p\}$ such that

Note that $\hat{E}[Y|\mathbf{X}\in R]$ and $\hat{P}_{R}$ are given by Eq. (9), (10).

Note that the above constrained optimization problem involves $2Lp$ decision variables (i.e., $r_{lj}^{l},r_{lj}^{u},l\in\{1,...,L\}$ and $j\in\{1,...,p\}$) and the disjoint constraint Eq. (13) actually involves $L\choose 2$$*p$ individual conditions. To further boost the computational efficiency, we define a new objective function with the disjoint constraint Eq. (13) considered as $O(R)$. Let the set consist of all regions that satisfy the requirements regarding region $R$ in Eq. (13) be $\Theta_{\mathbf{X}}$. The new objective function $O(R)$ is defined as

with $\eta$ being a numerical number that is smaller than the minimum value of $y_{i}$, $i=1,...,n$. The problem in Eq. (11), Eq. (12) and Eq. (13) is simplified to

subject to

Lagrangian relaxation is a popular approach to solve the constrained optimization problem like our problem expressed in Eq. (15) and Eq. (16). However, the Lagrangian relaxation approach is known to suffer the problem of the duality gap under the general setting. It is not guaranteed that we always obtain the global or even local optima. In our numerical studies, we experimented the Lagrangian relaxation approach but failed to find the optimal solutions in some situations. In this paper, we propose to use an augmented Lagrangian method, called the penalty method. To be more specific, the original constrained optimization problem in Eq. (15) and Eq. (16) is equivalently transformed into

where $(v)^{+}=\max(v,0)$, $c^{*}$ is the minimal point of function $K(c)$, $c>0$, defined as

The next step is to build a module to solve the single objective optimization problem in Eq. (18) for any given $c>0$. Due to the inaccessibility to the closed form of the objective function and the high dimensional nature (total number of parameters for any given $R$ is $2pL$), we propose to use the genetic algorithm (GA), one of the evolutionary algorithms, to conduct the optimization task [21]. GA does not require information or approximation on derivatives, which is not available in our problem setting. Moreover, the GA is more robust than coordinate cycling directed search methods as it maintains a population of potential solutions in each iteration of searching and bears less risk of being trapped in a local minimum.

GA tries to mimic the natural process and consists of three components: selection, based on the fitness score, determining which individual candidates are chosen for the mating process and how many offspring each candidate produces; crossover, an exchange is performed between some variables of the parents, producing two new individuals; and mutation, for a few selected offspring, one variable is altered by a small perturbation[21]. When using GA to optimize over continuous variables, the real coding is used and the fitness of individual candidates is evaluated by the magnitude of the objective function[16].

To solve the optimization problem in Eq. (18) using GA, for any given value of $c,L,\beta$, we proceed as follows:

• •

  Step 1: Initiate a candidate population of size $W$: for $w=1,...,W$, denote the $w$-th initialization as a $2Lp$-dimensional vector $[r_{l,j,w}^{l},r_{l,j,w}^{u}]^{T}$, for $l=1,...,L;j=1,...,p$. The region spanned by these $2Lp$ values is denoted as $R_{w}$.

• •

  Step 2: Evaluate the fitness (i.e., how good it is in terms of maximizing the objective function in Eq. (18) ) of each initialization by evaluating the objective function $O(R_{w})-c(\beta-\hat{P}_{R_{w}})^{+}$.

• •

  Step 3: Check whether certain stopping criteria is met. For instance, the criteria can be the increment in the fitness, compared to the previous population, is smaller than some threshold. If yes, output the best individual as the identified region $R^{*}$. Otherwise, move to the next step to produce a new population.

• •

  Step 4: There are several steps in the new population producing process:

  • –

    Select which individuals among the current population will produce offspring based on their fitness, i.e., the fitted values of the objective function $O(R)-c(\beta-\hat{P}_{R})^{+}$.Individuals with higher fitness values have higher chances to be selected.

  • –

    Some of the $2Lp$ variables are exchanged between the selected parents to produce offspring.

  • –

    The produced offspring are mutated by certain perturbation with a certain probability.

  • –

    The produced new population is fed into step 2.

In this paper, we used the R package ‘GA’ [22] to numerically implement the above procedures. The built-in parallel running option is adopted to increase the efficiency.

In supervised machine learning problems, generalization refers to the ability of an algorithm being effective for unseen data [23]. The generalization of machine learning algorithms is important as it ensures that a similar accuracy level can be achieved when the model is deployed to make predictions for new data sets.

Analogously, in our operating envelope problem, it is also important to ensure that the achieved mean value of the response variable (i.e., the KPI) when applying the identified region $R^{*}$ to unseen data does not deviate much from the mean value of the response variable achieved in the training phase due to sampling error. To ensure the generalization of the identified operating envelope, we propose the following regularized objective function

where $\text{SD}(\hat{E}[Y|\mathbf{X}\in R])$ is the sample standard deviation of the achieved estimate $\hat{E}[Y|\mathbf{X}\in R]$, $\gamma>0$ represents the size of regularization on the sample standard deviation part. The intuition is that if the sample standard error (i.e., variation across different data sets) of the achieved optimal average KPI within the detected hyper-rectangle is small, we have better chance to maintain the same optimal KPI level for any new data set.

There is usually no closed form for $\text{SD}(\hat{E}[Y|\mathbf{X}\in R])$ in Eq. (19). We propose to use the bootstrapping re-sampling technique to non-parametrically estimate this quantity [24]. The detailed estimation procedure is summarized as follows. Let the total number of bootstrapping be $M$,

• •

  Step 1: For iteration $m=1$ to $M$:

  • –

    Sample with replacement from the raw index {1,…,n}, and get another index set of length $n$, $\{d_{m,1},...,d_{m,i},...,d_{m,n}\}$.

  • –

    The new data set of size $n$ is $(y_{d_{m,i}},\mathbf{x}_{d_{m,i}})$ for $i=1,...,n$.

  • –

    Compute

    $$\hat{E}[Y|\mathbf{X}\in R]_{m}=\frac{\sum_{i=1}^{n}y_{d_{m,i}}I(\mathbf{x}_{d_{m,i}}\in R)}{\sum_{i=1}^{n}I(\mathbf{x}_{d_{m,i}}\in R)}\stackrel{{\scriptstyle\text{def}}}{{=}}\hat{E}_{m}.$$

    (20)

• •

  Step 2: Compute and output

  $$\hat{\text{SD}}_{\text{bs}}(\hat{E}[Y|\mathbf{X}\in R])=\sqrt{\frac{1}{M-1}\sum_{m=1}^{M}(\hat{E}_{m}-\bar{\hat{E}}_{m})^{2}}.$$

  (21)

Let the estimated sample standard deviation be $\hat{\text{SD}}_{\text{bs}}(\hat{E}[Y|\mathbf{X}\in R])$. The empirical regularized objective function for any given region $R$ is then

The counterparts of the bias and variance concepts in supervised machine learning problems exist in the proposed operating envelope problem in Definition 1. Specifically, we would like to achieve as high KPI value as possible in our setting. Hence, we define the bias as the gap between the expectation of the achieved KPI values for an unseen data set and the true optimal KPI value. The true optimal KPI value refers to the optimal objective value in Definition 1 assuming that we could have had access to the whole data population. In informal notation, $\text{Bias}=\text{True Optimal KPI}-E[\hat{E}^{*}_{\text{test}}[Y|\mathbf{X}\in R]]$. Since we would not have the true optimal KPI value, we propose to use another non-increasing function of the expected KPI value of unseen data to represent bias. In particular,

The variance concept here is similar to its counterpart in machine learning. In particular, we define variance as the variability of achieved KPI values of an unseen data set when applying the identified regions by solving the problem in Definition 1 using different training data sets. It is formally defined as follows:

For any given regularization parameter $\gamma$, the bias and variance defined above can be estimated by the cross-validation technique [25]. Similar to supervised machine learning, there exists a trade-off between the bias and variance, which requires judgments of the real needs to decide how to balance between bias and variance in the operating envelope problem. This fact is demonstrated by the simulation studies and the real-world data analysis in Section IV.

All the components discussed in Section III-A to III-D works together and forms the proposed regularized ‘GA + penalty’ algorithm summarized in Algorithm 1.

In this subsection, we consider a set of simulations with $p=1$ state variable to demonstrate the effectiveness of our proposed regularized ‘GA + penalty’ algorithm in Section III for the operating envelope with interpretability and implementability constraints problem defined in Section II. The first two experiments focus on investigating the impacts of different components in Definition 1 on the region detection results. In these two experiments, the variability is simulated to be the same across different sub-regions of $\Omega(\mathbf{X})$ and the ‘GA + penalty’ algorithm without regularization is utilized. In the latter two experiments, we investigate the benefit of adding regularization. For all the simulation studies, we assume that data $(y_{i},x_{i})$, for $i=1,...,n$ are generated from model:

where $x_{i}\sim\text{i.i.d }g(X)$ and $\epsilon_{i}\sim\text{i.i.d }N(0,\sigma_{\epsilon}^{2})$.

In this section, we consider a simulation study with $p=2$ state variables. We illustrate that it is infeasible to conduct a multidimensional operating envelope identification by separately learning the lower and upper limits on each state variable dimension when the state variables are correlated. We also show the results for different numbers of the disjoint regions.

The data are simulated as follows. The joint probability density function of the state variable vector is $g(X_{1},X_{2})=\frac{1}{4}\sum_{s=1}^{4}N({\mbox{\boldmath${\mu}$}}_{s},{\mbox{\boldmath${\Sigma}$}}))$, where ${\mbox{\boldmath${\Sigma}$}}=\left(\begin{smallmatrix}1.2&1\\ 1&1.2\end{smallmatrix}\right)$, the means are ${\mbox{\boldmath${\mu}$}}_{1}=(\pi,\pi)$, ${\mbox{\boldmath${\mu}$}}_{2}=(\pi,3\pi)$, ${\mbox{\boldmath${\mu}$}}_{3}=(3\pi,\pi)$ and ${\mbox{\boldmath${\mu}$}}_{4}=(3\pi,3\pi)$. The mapping from $(X_{1},X_{2})$ to $Y$ is specified as $f(X_{1},X_{2})=\sqrt{f_{0}(X_{1})f_{0}(X_{2})}$, with $f_{0}(Z)$ being $-2.25(\cos(Z)-4.5)$ when $0\leq Z\leq 2\pi$ and $-2(\cos(Z)-4.5)$ when $2\pi<Z\leq 4\pi$. The standard deviations of random errors are respectively 0.11, 0.15, 0.05, 0.5 in the four sub-regions ($S_{1}$, $S_{2}$, $S_{3}$, $S_{4}$) in Fig. 9(a). The contour plot of the generated data can be found in Fig. 9.

Under this simulation setting, it can be seen that we can achieve a better expected region-wise response variable by bounding both state variables. Let’s first consider to calculate a single rectangle shaped operating envelope (i.e., $L=1$) under the probabilistic coverage constraint $P(X_{1}\in[r_{1}^{l},r_{1}^{u}],X_{2}\in[r_{2}^{l},r_{2}^{u}])>\beta$. To calculate the 2-dimensional operating envelope, one natural thinking is to bound each state variable separately and then construct a rectangle shaped envelope from the Cartesian product of the detected intervals. However, it is infeasible to appropriately specify the required probabilistic coverage on each state variable dimension, due to the fact that relation in Eq. (30) no longer holds when the state variables are not independent.

In this study, we use simulated data to justify this argument. Specifically, suppose that the required coverage for the targeting rectangle shaped operating envelope is $\beta=0.2$. We calculate the interval on $X_{1}$ and $X_{2}$ separately with probabilistic coverage $\sqrt{\beta}=\sqrt{0.2}$. The identified intervals are shown in Fig. 9(b) and Fig. 9(c). However, the rectangle determined by the Cartesian product of these two intervals only covers $14.67\%$ of data, which is smaller than the required coverage threshold $\beta=0.2$. A more reasonable approach is to utilize our proposed ‘GA + penalty’ algorithm in Section III to jointly search over the two-dimensional state variable space. The detected region with a valid empirical coverage ($20.02\%>\beta=0.2$) is visualized in Fig. 9(d).

Next, we evaluate the performance of the proposed regularized ‘GA + penalty’ algorithm when target results consist of $L=2$ disjoint regions. According to our simulation setting, sub-region $S_{3}$ has the largest mean and the smallest variance (best sub-region). Sub-regions $S_{1}$ and $S_{4}$ have the second largest mean values, and the variability within sub-region $S_{1}$ is smaller than that within sub-region $S_{4}$. These facts indicate that the detected $L=2$ disjoint envelopes should be within $S_{3}$ and $S_{1}$, such that the achieved high mean response variable is more probably to be reached in unseen data. We first set the required probabilistic coverage as $\beta=0.25$, the bias and variance terms as a function of the regularization parameter $\gamma$ are given in Fig. 10(a). Based on these plots, we choose $\gamma=20$ and the resulted envelope is shown at the bottom of Fig. 10(a). The bias and variance curves, and the achieved operating envelope for $\beta=0.35$ are given in Fig. 10(b). The regularization term introduced in Section III-D helps to push the second envelope into the sub-region $S_{1}$ with less variability, instead of $S_{4}$, even though the mean response is the same within these two sub-regions. All the simulation results in Section IV-A and Section IV-B demonstrate the effectiveness of our proposed regularized operating envelope identification algorithm.

In this section, we apply the proposed regularized ‘GA + penalty’ operating envelope identification algorithm to real industrial data on Kaggle.com. The data set is uploaded to Kaggle by Eduardo Magalhães Oliveira and can be downloaded through https://www.kaggle.com/edumagalhaes/quality-prediction-in-a-mining-process. The data set includes the quality of produced ore in a flotation plant, variables indicating the quality of the raw ore, as well as state variables that quantify the operations along the mining processes. Among all the state variables, the data disposer mentioned that five of them, including starch flow, amina flow, ore pulp flow, ore pulp pH, and ore pulp density, are highly related to the quality of the produced ore. Our objective is to jointly identify an implementable set of upper and lower limits for the five important state variable such that the resulted ore quality is higher or highest given the quality of raw ore fed into the mining processes. Totally, we used 36872 records, and $50\%$ of them are used for training and the remaining $50\%$ data are used for testing.

First, to move the effect of raw ore’s quality on the quality of the produced ore, we propose to normalize the quality of the produced ore by subtracting its fitted value from the random forest regression given the raw ore’s quality. The random forest regression is learned from the training data. Next, we apply the proposed ‘GA + penalty’ approach both with and without regularization to identify the region spanned by the top five important state variables such that the corresponding mean normalized quality score of the produced ore is high. The results for coverage threshold $\beta=0.2$, the number of disjoint regions $L=1,2$, and the regularization parameter $\gamma=5,10$ are summarized in Table I. As shown by Table I, the results match with our intuition. First, when learning without regularization, we achieve a higher optimal mean response when $L=2$ compared to $L=1$. This is because we search over a larger set of candidate solutions when $L$ is higher. Second, given the number of disjoint regions $L$, the bias increases with the increase of regularization parameter, while the variance is smaller for larger regularization parameter values. This phenomenon justifies the bias and variance trade-off discussed in Section III-D. In practice, users can choose the combinations $L$, $\beta$, $\gamma$ that best meets their needs.

We notice that the computational cost of our GA based algorithm increases with the increment in dimensions (i.e., the number of state variables considered), due to the nature of GA. For the real data analysis, it took us about 5 hours to get the results. However, as discussed in the previous sections, GA is the most reasonable solution for the considered operating envelope with interpretability and implementability constraints problem. From the application perspective, the operating envelope is not needed to be learned in a streaming mode, hence the computational cost is not a big issue for the proposed ‘GA + penalty’ approach.