Compensating Data Shortages in Manufacturing with Monotonicity Knowledge

Abstract

Systematic decision making in engineering requires appropriate models. In this article, we introduce a regression method for enhancing the predictive power of a model by exploiting expert knowledge in the form of shape constraints, or more specifically, monotonicity constraints. Incorporating such information is particularly useful when the available datasets are small or do not cover the entire input space, as is often the case in manufacturing applications. We set up the regression subject to the considered monotonicity constraints as a semi-infinite optimization problem, and propose an adaptive solution algorithm. The method is applicable in multiple dimensions and can be extended to more general shape constraints. It was tested and validated on two real-world manufacturing processes, namely, laser glass bending and press hardening of sheet metal. It was found that the resulting models both complied well with the expert’s monotonicity knowledge and predicted the training data accurately. The suggested approach led to lower root-mean-squared errors than comparative methods from the literature for the sparse datasets considered in this work.

1. Introduction

Systematic decision making in manufacturing—such as finding optimal parameter settings for a manufacturing process—requires appropriate models for that process. In particular, such models have to be sufficiently accurate, and at the same time, sufficiently quick at evaluating. In principle, for many industrial processes, precise simulation models based on detailed physical modeling can be built. Yet, these so-called white-box models are typically too slow to be of any practical use in exploring the process parameter space and in eventually finding optimal process parameters (online and offline). In this respect, machine learning models can be very useful surrogates with short runtimes.

Conventional machine learning models are purely data-based (so-called black-box models). Accordingly, the predictive power of such models is generally bad if the underlying training data $\mathcal{D}=\{({\mathit{x}}_l,t_l):l\in \{1,\dots ,N\}\}$ are insufficient. Unfortunately, such data insufficiencies occur quite often, and they can come in one of the following forms: On the one hand, the available datasets can be too small and have too little variance in the input data points ${\mathit{x}}_1,\dots ,{\mathit{x}}_N$. This problem frequently occurs in manufacturing [1] because varying the process parameters beyond well-tested operating windows is usually costly. On the other hand, the output data $t_1,\dots ,t_N$ can be too noisy.

Aside from potentially insufficient data, however, one often also has additional knowledge about the relation between the input variables and the responses to be learned. Such extra knowledge about the considered process is referred to as expert knowledge in the following. In [2], the interactions of users with software was tracked to capture their expert knowledge in a general form as training data for a classification problem. In [3], expert knowledge was used in the form of a specific algebraic relation between input and output to solve a parameter estimation problem with artificial neural networks. Such informed machine learning [4] techniques beneficially combine expert knowledge and data to build hybrid or gray-box models [5,6,7,8,9,10,11,12], which predict the responses more accurately than purely data-based models. In other words, by using informed machine learning techniques, one can compensate data insufficiencies with expert knowledge.

An important and common type of expert knowledge is prior information about the monotonicity behavior of the unknown functional relationship $\mathit{x}\mapsto y\left(\mathit{x}\right)$ to be learned. A large variety of concrete application examples with monotonicity knowledge can be found in ([13], Section 4.1) and ([14], Section 1), for instance. The present article exclusively deals with regression under such monotonicity requirements. For classification under monotonicity constraints, see, e.g., [14,15]. Along with convexity constraints, monotonicity constraints are probably the most intensively studied shape constraints [16] in the literature, and correspondingly, there exist plenty of different approaches to incorporate monotonicity knowledge in a machine learning model. See [17] for an extensive overview. Very roughly, these approaches can be categorized according to when the monotonicity knowledge is taken into account: in or only after the training phase. In the terminology of [4], this corresponds to the distinction between knowledge integration in the learning algorithm or in the final hypothesis.

A lot of methods—especially from the mathematical statistics literature, such as [18,19,20,21,22,23,24]—incorporate monotonicity knowledge only after training. These articles start with a purely data-based initial model, which in general does not satisfy the monotonicity requirements, and then monotonize this initial model according to a suitable monotonization procedure, such as projection [18,19,20,24], rearrangement [22,23,25] or tilting [21]. Among other things, it is shown in the mentioned articles that, in spite of noise in the output data, the arising monotonized models are close to the true relationship for sufficiently large training datasets. Summarizing, these articles show that for large datasets, noise in the output data can be compensated by monotonization to a certain extent.

In contrast to that, in some works, such as [13,17,26,27,28,29], monotonicity knowledge was incorporated already in training. In these articles, the monotonicity requirements were added as constraints—either hard [17,26,28,29] or soft [13,26]—to the data-based optimization of the model parameters. In [13,28], probabilistic monotonicity notions are used. In [26,27,28,29], support vector regressors in the linear-programming or the more standard quadratic-programming form, Gaussian process regressors or neural network models were considered, and the monotonicity of these models was enforced by constraints on the model derivatives at predefined sampling points [26,28,29] or on the model increments between predefined pairs of sampling points [27].

A disadvantage of the projection- and rearrangement-based approaches [22,23,24] from the point of view of manufacturing applications is that these methods are tailored to large datasets. Another disadvantage of these approaches is that the resulting models typically exhibit distinctive kinks, which are almost always unphysical. Additionally, the models resulting from the multidimensional rearrangement method by [23] are not guaranteed to be monotonic when trained on small datasets. A drawback of the tilting approach from [21] is that it is formulated and validated only for one-dimensional input spaces (intervals in $\mathbb{R}$). Accordingly, naively extending the non-adaptive discretization scheme from [21] to higher dimensions would result in long computation times. A downside of the in-training methods from [26,28,29] is that the sampling points at which the monotonicity constraints are imposed have to be chosen in advance (even though they need not coincide with the training data points).

With the method proposed in the present article, we address the aforementioned issues and shortcomings. In Section 2, our methodology for monotonic regression using semi-infinite optimization is introduced. It incorporates the monotonicity knowledge during training. Specifically, polynomial regression models are assumed for the input–output relationships to be learned. Since there is no after-training monotonization step in the method, our models are smooth, and in particular, do not exhibit kinks. Also, due to the employed adaptive discretization scheme, the method is computationally efficient also in higher dimensions. To our knowledge, such an adaptive scheme has not been applied to solve monotonic regression problems before, especially not in situations with sparse data. In Section 4, the method is validated by means of two applications to real-world processes which are both introduced in Section 3, namely, laser glass bending and press hardening of sheet metal. It turns out that the adaptive semi-infinite optimization approach to monotonic regression is better suited for the considered applications with their small datasets and the resulting models are more accurate than those obtained with the comparative approaches from the literature.

2. Semi-Infinite Optimization Approach to Monotonic Regression

In this section, our semi-infinite optimization approach to monotonic regression is introduced. It will be referred to as the SIAMOR method later on for brevity.

2.1. Semi-Infinite Optimization Formulation of Monotonic Regression

In our approach to monotonic regression, polynomial models

$$\begin{array}{c}\hfill \mathit{x}\mapsto {\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)=\sum _{\left|\alpha \right|\le m}{w}_{\alpha }{\mathit{x}}^{\alpha }\in \mathbb{R}\end{array}$$

(1)

are used for all input–output relationships $\mathit{x}\mapsto y\left(\mathit{x}\right)$ to be learned. In the above relation (1), the sum extends over all d-dimensional multi-indices ([30], Section 1) $\alpha =({\alpha }_1,\dots ,{\alpha }_d)\in {\mathbb{N}}_0^d$ with degree $\left|\alpha \right|:={\alpha }_1+\cdots +{\alpha }_d$ less than or equal to some total degree $m\in \mathbb{N}$. The terms ${\mathit{x}}^{\alpha }:=x_1^{{\alpha }_1}\cdots x_d^{{\alpha }_d}$ are the monomials in d variables of degrees less than or equal to m, and $w_{\alpha }$ are the corresponding model parameters to be tuned by regression. Since there are exactly

$$\begin{array}{c}\hfill N_m=\sum _{k=0}^m\left(\genfrac{}{}{0pt}{}{k+d-1}{d-1}\right)=\left(\genfrac{}{}{0pt}{}{m+d}{m}\right)\end{array}$$

(2)

d-dimensional monomials of degrees less than or equal to m, the polynomial regression model (1) can be equivalently written as

$$\begin{array}{c}\hfill {\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)=\sum _{i=1}^{N_m}{w}_i{\varphi }_i\left(\mathit{x}\right)=w^{\top }\mathbf{\varphi }\left(\mathit{x}\right),\end{array}$$

(3)

where the basis functions ${\varphi }_1,\dots ,{\varphi }_{N_m}$ constitute any enumeration of the d-dimensional monomials of degrees less than or equal to m, while $\mathit{w}:={(w_1,\dots ,w_{N_m})}^{\top }$ and $\mathbf{\varphi }\left(\mathit{x}\right):={({\varphi }_1\left(\mathit{x}\right),\cdots ,{\varphi }_{N_m}\left(\mathit{x}\right))}^{\top }$.

Standard polynomial regression without regularization ([31], Section 2.1) is about solving the unconstrained optimization problem:

$$\begin{array}{c}\hfill \underset{\mathit{w}\in W}{min}\frac{1}{2}\sum _{l=1}^N{\left({\widehat{y}}_{\mathit{w}}\left({\mathit{x}}_l\right)-t_l\right)}^2,\end{array}$$

(4)

or in other words, about optimally adapting the model parameters $w_i\in [-r,r]$ of the polynomial model (3) to the available dataset $\mathcal{D}=\{({\mathit{x}}_l,t_l):l\in \{1,\dots ,N\}\}$ containing N points. In the above relation, the monomial coefficients are allowed to vary in the compact hyperbox

$$\begin{array}{c}\hfill W=\{\mathit{w}\in {\mathbb{R}}^{N_m}:-r\le w_i\le r\mathrm{for}\mathrm{all}i\in \{1,\dots ,N_m\}\}\end{array}$$

(5)

with some large but finite $r>0$. Since W is compact and non-empty, and since the mean-squared error objective function of (4) is continuous, the standard polynomial regression problem (4) for any given dataset $\mathcal{D}$ has a solution $w$ (which is unique if, for instance, an ${\ell }^2$-regularization term is added).

In general, however, the resulting model $\mathit{x}\mapsto {\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)$ does not necessarily exhibit the monotonicity behavior an expert expects for the underlying true physical relationship $\mathit{x}\mapsto y\left(\mathit{x}\right)$. In order to enforce the expected monotonicity behavior, the following constraints on the signs of the partial derivatives ${\partial }_{x_j}{\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)$ are added to the unconstrained standard regression problem (4):

$$\begin{array}{c}\hfill {\sigma }_j·{\partial }_{x_j}{\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)\ge 0\mathrm{for}\mathrm{all}j\in J\mathrm{and}\mathit{x}\in X.\end{array}$$

(6)

The numbers ${\sigma }_j\in \{-1,0,1\}$ indicate the expected monotonicity behavior for each coordinate direction $j\in \{1,\dots ,d\}$:

• ${\sigma }_j=1$ and ${\sigma }_j=-1$ indicate that $\mathit{x}\mapsto y\left(\mathit{x}\right)$ is expected to be, respectively, monotonically increasing or decreasing in the jth coordinate direction;
• ${\sigma }_j=0$ indicates that one has no monotonicity knowledge in the jth coordinate direction.

Also, $J:=\{j\in \{1,\dots ,d\}:{\sigma }_j\ne 0\}$ is the set of all directions for which a monotonicity constraint is imposed, and the vector

$$\begin{array}{c}\hfill \mathbf{\sigma }:=({\sigma }_1,\dots ,{\sigma }_d)\end{array}$$

is referred to as the monotonicity signature of the relationship $\mathit{x}\mapsto y\left(\mathit{x}\right)$. Finally, $X\subset {\mathbb{R}}^d$ is the (continuous) subset of the input space on which the polynomial model (1) is supposed to be a reasonable prediction for $\mathit{x}\mapsto y\left(\mathit{x}\right)$. In this work, X was chosen to be identical with the range covered by the input training data points ${\mathit{x}}_1,\dots ,{\mathit{x}}_N$. I.e., X is the compact hyperbox

$$\begin{array}{c}\hfill X=[a_1,b_1]\times \cdots \times [a_d,b_d]\end{array}$$

(7)

with $a_j:={min}_{l=1,\cdots ,N}{x}_{l,j}$ and $b_j:={max}_{l=1,\cdots ,N}{x}_{l,j}$ and with $x_{l,j}$ denoting the jth component of the lth input data point $x_l$. Writing

$$\begin{array}{c}\hfill f\left(\mathit{w}\right):=\frac{1}{2}\sum _{l=1}^N{({\widehat{y}}_{\mathit{w}}\left({\mathit{x}}_l\right)-t_l)}^2\mathrm{and}{g}_j(\mathit{w},\mathit{x}):={\sigma }_j·{\partial }_{x_j}{\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)\end{array}$$

(8)

for brevity, our monotonic regression problem (4)–(6) takes the neat and simple form

$$\begin{array}{c}\hfill \underset{\mathit{w}\in W}{min}f\left(\mathit{w}\right)\mathrm{s}.\mathrm{t}.g_j(\mathit{w},\mathit{x})\ge 0\mathrm{for}\mathrm{all}j\in J\mathrm{and}\mathit{x}\in X.\end{array}$$

(9)

Since the input set X is continuous and hence contains infinitely many points $\mathit{x}$, the monotonic regression problem (9) features infinitely many inequality constraints. Consequently, (9) is a semi-infinite optimization problem [32,33,34,35,36] (or more precisely, a standard semi-infinite optimization problem, as opposed to a generalized one). It is well-known that any semi-infinite problem—and, in particular, the monotonic regression problem (9)—can be equivalently rewritten as a bi-level optimization problem [35,37,38], namely,

$$\begin{array}{c}\hfill \underset{\mathit{w}\in W}{min}f\left(\mathit{w}\right)\mathrm{s}.\mathrm{t}.\underset{\mathit{x}\in X}{min}{g}_j(\mathit{w},\mathit{x})\ge 0\mathrm{for}\mathrm{all}j\in J.\end{array}$$

(10)

Commonly, the minimization subproblems in the constraints of (10) are called lower-level problems of (9).

2.2. Adaptive Solution Strategy

Since the feasible set of (9) is compact (by the finiteness of the parameter r in (5)) and non-empty (it contains $\mathit{w}:=0\in {\mathbb{R}}^{N_m}$), the monotonic regression problem (9) has a solution by virtue of the Weierstraß extreme-value theorem. In order to compute such a solution of (9), a variant of the adaptive, iterative discretization algorithm by [39] is used. In a nutshell, the idea is the following: the infinite index set X of the original regression problem (9) is iteratively replaced by discretizations, that is, finite subsets $X^k\subset X$, and these discretizations are adaptively refined from iteration to iteration. In that manner, in every iteration k one obtains the ordinary (finite) optimization problem

$$\underset{\mathit{w}\in W}{min}f\left(\mathit{w}\right)\mathrm{s}.\mathrm{t}.g_j(\mathit{w},\mathit{x})\ge 0\mathrm{for}\mathrm{all}j\in J\mathrm{and}\mathit{x}\in X^k$$

(11)

featuring only finitely many inequality constraints. As usual, we refer to (11) as the kth discretized problem. In each iteration k, two steps are performed, namely, an optimization step and an adaptive refinement step. In the optimization step, a solution ${\mathit{w}}^k$ of the kth discretized problem (11) is computed. In the refinement step, for each direction $j\in J$, a point ${\mathit{x}}^{k,j}\in X$ is computed at which the jth monotonicity constraint at $\mathit{w}={\mathit{w}}^k$ is violated most. In more precise terms, for every $j\in J$, an approximate solution ${\mathit{x}}^{k,j}$ of the global optimization problem

$$\begin{array}{c}\hfill \underset{\mathit{x}\in X}{min}{g}_j({\mathit{w}}^k,\mathit{x})\end{array}$$

(12)

is computed. All the points ${\mathit{x}}^{k,j}$ for which a monotonicity violation occurs are then added to the current discretization $X^k$ in order to obtain the new discretization $X^{k+1}$. If no more monotonicity violations occur, the iteration is stopped. As usual, (12) is referred to as the $(k,j)$th lower-level problem in the following.

With regard to the practical implementation of the above solution strategy, it is important to observe that the discretized problems (11) are standard convex quadratic programs [40]. Indeed, by inserting (3) into (8) and using the design matrix $\mathsf{\Phi }$ with entries ${\mathsf{\Phi }}_{li}:={\varphi }_i\left({\mathit{x}}_l\right)$, one obtains

$$\begin{array}{c}\hfill f\left(\mathit{w}\right)=\frac{1}{2}{\parallel \mathsf{\Phi }\mathit{w}-t\parallel }_2^2=\frac{1}{2}{\mathit{w}}^{\top }{\mathsf{\Phi }}^{\top }\mathsf{\Phi }\mathit{w}-t^{\top }\mathsf{\Phi }\mathit{w}+\frac{1}{2}{t}^{\top }t.\end{array}$$

Consequently, the objective function of (11) is indeed quadratic and convex with respect to $\mathit{w}$. It is not strictly convex, though, in the sparse-data case $N<N_m$ considered in this paper. (Indeed, the kernel of the matrix ${\mathsf{\Phi }}^{\top }\mathsf{\Phi }\in {\mathbb{R}}^{N_m\times N_m}$ is equal to the kernel of $\mathsf{\Phi }\in {\mathbb{R}}^{N\times N_m}$, and an $N\times N_m$ matrix with $N<N_m$ has a non-trivial kernel, of course.)

Also, in view of

$$\begin{array}{c}\hfill g_j(\mathit{w},\mathit{x})={\sigma }_j·{\partial }_{x_j}{\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)={\sigma }_j·{\mathit{w}}^{\top }\left({\partial }_{x_j}\mathbf{\varphi }\left(\mathit{x}\right)\right),\end{array}$$

(13)

the constraints of (11) are indeed linear with respect to $w$.

With regard to the practical implementation, it is also important to observe that the objective functions $\mathit{x}\mapsto g_j({\mathit{w}}^k,\mathit{x})$ of the lower-level problems (12) are non-convex polynomials and therefore in general have several local minima. Consequently, (12) needs to be solved numerically with a global optimization solver.

2.3. Algorithm and Implementation Details

In the following, our adaptive discretization algorithm is described in detail. As has already been pointed out above, it is a variant of the general algorithm developed by ([39], Section 2), and it is explained after Algorithm 1 how our variant differs from its prototype [39].

Algorithm 1. Adaptive discretization algorithm for monotonic regression.

Choose a coarse (but non-empty) rectangular grid $X^0$ in X. Set $k=0$ and iterate over k. Solve the kth discretized problem (11) to obtain optimal model parameters ${\mathit{w}}^k\in W$. Solve the $(k,j)$th lower-level problem (12) approximately for every $j\in J$ to find approximate global minimizers ${\mathit{x}}^{k,j}\in X$. Add those of the points ${\mathit{x}}^{k,j}$ for which substantial monotonicity violations occur, i.e., for which $$g_j({\mathit{w}}^k,{\mathit{x}}^{k,j})<-{\epsilon }_j,$$ to the current discretization $X^k$ and go to Step 1 with $k=k+1$. If for none of the points ${\mathit{x}}^{k,j}$ substantial monotonicity violations occur, go to Step 3. Check for monotonicity violations on a fixed, fine rectangular reference grid $X_{\mathrm{ref}}\subset X$. If there are no such violations, that is, if $$g_j({\mathit{w}}^k,\mathit{x})\ge -{\epsilon }_j$$ for all $j\in J$ and $\mathit{x}\in X_{\mathrm{ref}}$, then terminate. Else, for every direction j with violations, add the reference grid point ${\mathit{x}}_{\mathrm{ref}}^{k,j}$ with the largest violation to $X^k$ and go to Step 1 with $k=k+1$.

In contrast to [39], the algorithm above does not require exact solutions of the (non-convex) lower-level problems. Indeed, Step 2 of Algorithm 1 only requires finding an approximate solution. Also, slight constraint violations are tolerated (Step 2 and 3), and a feasibility check on a reference grid (Step 3) is performed before termination. Without the feasibility check on the reference grid, it could happen that the algorithm—because of the merely approximate solutions of the lower-level problems—terminates at models ${\widehat{y}}_{{\mathit{w}}^k}$ which do not satisfy the imposed monotonicity constraints sufficiently well. Another difference to the algorithm from [39] is that there there are several lower-level problems in each iteration in this work and not just one, because monotonicity is enforced in multiple coordinate directions in general.

In our specific applications, the parameters inherent in Algorithm 1 were chosen as follows. The degrees m of the polynomial models in this work were chosen as the largest possible values that did not result in an overfit, because increasing m enhances the model’s accuracy in general. In this respect, the number of model parameters was allowed to exceed the number of data points ($N_m\ge N$), since the constraints represent additional information supplementing the data. As for the parameter r in (5), one only has to make sure that it is so large that the resulting box constraints in the discretized problems (11) actually do not restrain the solutions ${\mathit{w}}^k$ (Step 1 of Algorithm 1). In other words, r should be so large that relaxing or even dropping the pertaining box constraints does not improve the minimizer computed for (11) anymore. In the specific applications considered here, $r={10}^5$ turned out to meet this requirement. As for the tolerances ${\epsilon }_j$ (Steps 2 and 3 of Algorithm 1), a monotonicity violation of 1% of the ranges covered by the in- and output training data was allowed for:

$$\begin{array}{c}\hfill {\epsilon }_j=0.01\frac{{max}_{l=1,\dots ,N}{t}_l-{min}_{l=1,\cdots ,N}{t}_l}{{max}_{l=1,\dots ,N}{x}_{l,j}-{min}_{l=1,\cdots ,N}{x}_{l,j}}.\end{array}$$

(14)

And finally, concerning the reference grid $X_{\mathrm{ref}}$ in the finalization step (Step 3 of Algorithm 1), twenty values per input dimension, equidistantly distributed from the lower to the upper bound along each direction, were used.

Algorithm 1 was implemented in Python and the package sklearn was used for the numerical representation of the models. Since the discretized problems (11) are standard convex quadratic programs, a solver tailored to that specific problem class was used, namely, quadprog [41]. It can solve quadratic programs with hundreds of variables and thousands of constraints in just a few seconds because it efficiently exploits the simple structure of the problem. Since, on the other hand, the lower-level problems (12) are global optimization problems with possibly several local minima, a suitable global optimization solver was required. We chose the solver scipy.optimize.shgo [42], which employs a simplicial homology strategy, and which, in our applications, turned out to be a good compromise between speed and reliability. For the problems considered in this article, shgo’s internal local optimization was configured to occur in every iteration, to multi-start from a Sobol set of 100$·d$ points and to be executed using the algorithm L-BFGS-B with analytical gradients.

3. Applications in Manufacturing

In this section, two real-world manufacturing applications are described to which our monotonic regression algorithm was applied.

3.1. Laser Glass Bending

A first application example is laser glass bending. In the industrial standard process of glass bending [43], a flat glass specimen is placed in a furnace with furnace temperature $T_{\mathrm{f}}$, and then the heated specimen is bent at a designated edge driven by gravity. As an additional feature, a laser can be added to the industrial standard process in order to specifically heat up the critical region of the flat glass around the bending edge, and thus to speed up the process and achieve smaller bending radii [44,45]. The laser can generally scribe in arbitrary directions. In the process considered here, however, the laser path is restricted to three straight lines parallel to the bending edge. While the middle line defines the bending edge, the two outer lines are at a fixed distance ${\Delta }_{\mathrm{l}}/2=5.75$ mm in each direction to it. The laser spot moves along this path in multiple cycles with the number of cycles denoted by $n_{\mathrm{c}}$. The scribing speed and the power of the laser are kept constant. A mechanical stop below the bending edge guarantees that the bending angle does not exceed ${90}^{°}$. An illustration of the laser glass bending process is shown in Figure 1.

The goal of the glass bending process considered here is to obtain bent glass parts with a prescribed bending angle. In order to achieve such a pre-defined bending angle, the process operator has to find suitable combinations of the two process parameters $T_{\mathrm{f}}$ and $n_{\mathrm{c}}$, which is usually done based on experience. A more systematic approach, however, is to set up an appropriate model of the bending angle $y:=\beta$ as a function of the process variables

$$\begin{array}{c}\hfill \mathit{x}:=(x_1,x_2):=(T_{\mathrm{f}},n_{\mathrm{c}})\in X,\end{array}$$

(15)

where $X\subset {\mathbb{R}}^2$ is the rectangular set with the bounds specified in

Table 1. Ranges for the process variables of laser glass bending.

Variable	Min	Max	Phys. Unit
$T_{\mathrm{f}}$	480	560	${}^{°}$C
$n_{\mathrm{c}}$	40	50	—

. In particular, such a model should allow sufficiently precise real-time predictions in order to support the process operator in quickly searching the parameter space X for optimal process parameter settings. Since SIAMOR models are polynomial by construction, they perfectly satisfy this real-time requirement. In contrast, a repeated evaluation of finite-element simulation models of the glass-bending process is too time-consuming to be of any practical use in quickly exploring the parameter space.

As generating experimental training data from the real process is cumbersome, a two-dimensional finite-element model was set up to generate data numerically. The simulation of the process was based on a coupled thermo-mechanical problem with finite deformation. Since the CO${}_2$ laser used in the process operates in the opaque wavelength spectrum of glass, the heat supply was modeled as a surface flux into the deforming sheet. In this two-dimensional setting, the heat was assumed to be deposited instantaneously along the thickness direction and also instantaneously on all three laser lines. Radiation effects were ignored, and heat conduction inside the glass was described by the classical Fourier law with the heat conductivity obtained experimentally via laser flash analysis. In view of the relevant relaxation and process time scales for the applied temperature range, the mechanical behavior of the glass was described by a simple Maxwell-type visco-elastic law. The deformation due to gravity is heavily affected by the pronounced temperature dependence of the viscosity above the glass transition, which was described in our models using the Williams–Landel–Ferry approximation [46]. The generation of the simulated data was conducted using the commercial finite-element code Abaqus^©. It was used to create a training dataset comprising 25 data points sampled on a 2D rectangular grid. The values used for the two degrees of freedom (five for $T_{\mathrm{f}}$ and five for $n_{\mathrm{c}}$) were placed equidistantly from the lower to the upper bounds given in Table 1.

Within these ranges and for the laser configuration described above, process experts expect the following monotonicity behavior: the bending angle $y=\beta$ should increase monotonically with increasing glass temperature in the critical region, and thus with increasing $T_{\mathrm{f}}$ and $n_{\mathrm{c}}$. In other words, the monotonicity signature $\mathbf{\sigma }$ of the bending angle y as a function of the inputs $x$ from (15) is expected to be

$$\begin{array}{c}\hfill \mathbf{\sigma }=({\sigma }_1,{\sigma }_2)=(1,1).\end{array}$$

(16)

3.2. Forming and Press Hardening of Sheet Metal

Another application example is press hardening [47]. In the experimental setup considered here, a blank is placed in a chamber furnace with a furnace temperature $T_{\mathrm{f}}$ above $900{}^{°}$C. After heating the blank, an industrial robot transports it with handling time $t_{\mathrm{h}}$ into the cooled forming tool. In the following, the extra handling time $\Delta t_{\mathrm{h}}=t_{\mathrm{h}}-10$s is used instead, with $10$s being the minimum time the robot takes to move the blank from the furnace to the press. The final combined forming and quenching step allows for the variation of the press force $F_{\mathrm{p}}$ and the quenching time $t_{\mathrm{q}}$. Afterwards, the formed part is transferred by the industrial robot to a deposition table for further cooling. An illustration of the process chain is shown in Figure 2.

The goal of the press hardening process considered in this work is to obtain a formed metal part with a prescribed hardness level, where the hardness is measured in units of the Vickers hardness number (unit symbol HV). In order to achieve such a pre-defined hardness, the process operator has to find suitable combinations of the four process parameters $T_{\mathrm{f}},\Delta t_{\mathrm{h}},F_{\mathrm{p}},t_{\mathrm{q}}$. And for that purpose, in turn, an appropriate model is needed for the hardness y of the formed part (at distinguished measurement points on the surface of the part) as a function of the process variables

$$\begin{array}{c}\hfill \mathit{x}:=(x_1,\dots ,x_4):=(T_{\mathrm{f}},\Delta t_{\mathrm{h}},F_{\mathrm{p}},t_{\mathrm{q}})\in X,\end{array}$$

(17)

where $X\subset {\mathbb{R}}^4$ is the hypercuboid set with the bounds specified in

Table 2. Ranges for the process variables of press hardening.

Variable	Min	Max	Phys. Unit
$T_{\mathrm{f}}$	871	933	${}^{°}$C
$\Delta t_{\mathrm{h}}$	0	4	s
$F_{\mathrm{p}}$	1750	2250	kN
$t_{\mathrm{q}}$	2	6	s

. In particular, such a model has to allow sufficiently accurate real-time predictions in order to help the process operator in quickly searching the parameter space X for optimal process parameter settings. Since SIAMOR models are polynomial by construction, they perfectly satisfy this real-time requirement. In contrast, due to the four-dimensional parameter space, already a single evaluation of a representative finite-element simulation model of the press-hardening process is prohibitively time-consuming to be of any practical use in quickly exploring the parameter space.

As in the case of glass bending, experiments for the press hardening process are expensive because they usually require manual adjustments, which tend to be time-consuming. Additionally, the local hardness measurements at the chosen measurement points on the surface of the quenched part are time-consuming as well. This is why the training data base we used is rather small. It contains 60 points resulting from a design of experiments with the four process variables $T_{\mathrm{f}}$, $\Delta t_{\mathrm{h}}$, $F_{\mathrm{p}}$ and $t_{\mathrm{q}}$ ranging between the bounds in Table 2, along with the corresponding hardness values at six local measurement points (referred to as MP1, …, MP6 in the following).

In order to compensate this data shortage, expert knowledge is brought into play. An expert for press hardening expects the hardness to decrease monotonically with $\Delta t_{\mathrm{h}}$ and to increase monotonically with $T_{\mathrm{f}}$ as well as with $t_{\mathrm{q}}$. In other words, the monotonicity signature $\mathbf{\sigma }$ of the hardness y (at any given measurement point) as a function of the inputs $x$ from (17) is expected to be

$$\begin{array}{c}\hfill \mathbf{\sigma }=({\sigma }_1,\dots ,{\sigma }_4)=(1,-1,0,1).\end{array}$$

(18)

In fact, a press hardening expert expects even a bit more, namely that the hardness should grow in a sigmoid-like manner with $T_{\mathrm{f}}$ and that it grows concavely towards saturation with increasing $t_{\mathrm{q}}$. All these requirements result from qualitative physical considerations and are supported by empirical experience.

4. Results and Discussion

In this section, we describe the results of our semi-infinite adaptive optimization approach to monotonic regression (SIAMOR) in the industrial processes described in Section 3, and compare them to the results of other approaches to incorporating monotonicity knowledge, which are well-known from the literature.

4.1. Informed Machine Learning Models for Laser Glass Bending

To begin with, the SIAMOR method was validated on a 1D subset of the data for laser glass bending, namely, the subset of all data points for which $n_{\mathrm{c}}=50$. This means that, out of the 25 data points, five points remained for training. First of all, ordinary unconstrained regression techniques were tried (see Figure 3a). A polynomial model of degree $m=3$ (solid line) and a Gaussian process regressor [31] (GPR, dashed line) did not comply with the monotonicity knowledge at high $T_{\mathrm{f}}$. A radial basis function (RBF) kernel was used for the GPR. This non-parametric model is always a reasonable choice for simulated data because it accurately reproduces the data themselves if the noise-level parameter is kept small. For all GPR models in this work, that parameter was set to 10${}^{-5}$. Next, the polynomial model was regularized in a ridge regression (dotted line), where the squared ${\ell }^2$-norm ${\lambda \parallel \mathit{w}\parallel }_2^2$ with a regularization weight $\lambda$ was added to the objective function in (4). $\lambda =0.003$ was chosen, which was roughly the minimum necessary value to achieve monotonicity. However, the resulting model does not predict the data very well. Thus, all three models from Figure 3a were unsatisfactory.

As a next step, the monotonicity requirement with respect to $T_{\mathrm{f}}$ was brought to bear, and monotonic regression with the SIAMOR method ($m=5$) was used (see Figure 3b) and compared to the rearrangement [22] and to the monotonic projection [24] of the Gaussian process regressor from Figure 3a. As mentioned before, both comparative methods are based on a non-monotonic pre-trained reference predictor. This makes them fundamentally different to the SIAMOR method, which imposes the monotonicity already in the training phase. The projection was calculated as described in the Appendix A with $\left|G\right|=80$ grid points. For the rearrangement method, the R package monreg was invoked from Python using the package rpy2. The degree m of the polynomial ansatz (1) used in the SIAMOR method was chosen as described in Section 2.3. For the specific case considered here, the curve started to vary unreasonably (albeit still monotonically) between the data points for $m\ge 6$, and therefore $m=5$ was chosen. The SIAMOR algorithm was initialized with five equidistant constraint locations in $X^0$, and it converged in iteration 5 with a total number of nine constraints. The locations of the constraints are marked in Figure 3b by the gray, vertical lines. The adaptive algorithm automatically places the non-initial constraints in the non-monotonic region at high $T_{\mathrm{f}}$. In terms of the root-mean-squared error

$$\begin{array}{c}\hfill \mathrm{RMSE}=\sqrt{\frac{1}{N}\sum _{l=1}^N{\left(\widehat{y}\left({\mathit{x}}_l\right)-t_l\right)}^2}\end{array}$$

(19)

on the training data, the SIAMOR model fits the data best; see

Table 3. Root-mean-squared deviations (RMSE) of the monotonic regression models from the training data for laser glass bending (1D).

Monotonic Regression Type	RMSE [${}^{°}$]
projection [24]	1.3822
rearrangement [22]	1.8432
SIAMOR	1.1598

. Another advantage of the SIAMOR model is that it is continuously differentiable, whereas the rearrangement and projection models exhibit (slight) unphysical kinks, which are typical for these methods [22]. And finally, the rearrangement and the projection models, for temperatures higher than $540{}^{°}$C, both predict bending angles larger than ${90}^{°}$, which is unphysical due to the mechanical stop used in the glass bending process. In contrast, the predictions of the SIAMOR model do not (significantly) exceed ${90}^{°}$.

After these calculations on a 1D subset, the full 2D dataset of the considered laser glass bending process with its 25 data points was used. The results are shown in Figure 4. Again, part (a) of the figure displays an unconstrained Gaussian process regressor for comparison. The RBF kernel contained one length scale parameter per input dimension, and sklearn correctly adjusted these hyperparameters using L-BFGS-B. I.e., the employed length scales maximize the log-likelihood function of the model. Nevertheless, the model is unsatisfactory because it exhibits a bump in the rear right corner of the plot, contradicting the monotonicity knowledge.

Figure 4b shows the 2D monotonic projection of the GPR with the monotonicity requirements (16) with respect to $T_{\mathrm{f}}$ and $n_{\mathrm{c}}$. It was calculated according to the Appendix A on a rectangular grid G consisting of 40${}^2$ points (40 values per input dimension). The resulting model looks generally reasonable, and in particular, satisfies the monotonicity specifications, but it exhibits kinks and plateaus. The most conspicuous kink starts at about $T_{\mathrm{f}}=546{}^{°}$C, $n_{\mathrm{c}}=50$ and proceeds towards the front right. The rearrangement method by [23] was not used for comparison here because for small datasets in $d>1$, it does not guarantee monotonicity.

Figure 4c displays the corresponding response surface of a polynomial model of the form (1) with degree $m=7$ trained with SIAMOR. For $m=7$ there are $N_m=36$ model parameters. The discretization $X^0$ was initialized with a rectangular grid using five equidistant values per dimension. The algorithm converged in iteration 11 with 69 final constraints. The resulting model is smoother than the one in Figure 4b, and it predicts the training data more accurately. Indeed, the corresponding RMSE values are 1.2518${}^{°}$ for projection and 0.6607${}^{°}$ for SIAMOR.

4.2. Informed Machine Learning Models for Forming and Press Hardening

As in the glass bending case, the SIAMOR method was first validated on a 1D subset of the data for the press hardening process. Namely, only those data points with $F_{\mathrm{p}}=2250$kN, $\Delta t_{\mathrm{h}}=4$s and $t_{\mathrm{q}}=2$s were considered. These specifications are met by six data points, and these were used to train the models shown in Figure 5. The data are not monotonic due to experimental noise. However, they reflect the expected sigmoid-like behavior mentioned in Section 3.2, and this extends to the monotonized models. An unconstrained polynomial with $m=3$ was chosen as the reference model to be monotonized for the comparative methods from the literature. Degrees lower than that resulted in larger deviations from the data and degrees higher than that resulted in overfitting. Thus, out of all models of the form (1), the hyperparameter choice $m=3$ yielded the lowest RMSE values for projection and rearrangement. For the monotonic regression with SIAMOR, $m=6$ and five equidistant initial constraint locations in $X^0$ were chosen. It converged in iteration 8 with a total of 12 monotonicity constraints. In terms of the root-mean-squared error, the SIAMOR model predicts the training data more accurately, as can be seen in

Table 4. Root-mean-squared deviations (RMSE) of the monotonic regression models from the training data for forming and press hardening of sheet metal (1D).

Monotonic Regression Type	RMSE [HV]
projection [24]	5.0893
rearrangement [22]	4.8346
SIAMOR	3.3583

. The reason is that the rearrangement- and projection-based models are dragged away from the data by the underlying reference model, especially at high $T_{\mathrm{f}}$.

After these 1D considerations, the SIAMOR method was validated on the full 4D dataset of the press hardening process. Polynomial models with degree $m=3$ were used for unconstrained regression and monotonic projection, and polynomials with $m=6$ were used for the SIAMOR method. The resulting models are visualized in the surface plots in Figure 6. The unconstrained model from Figure 6a clearly shows non-monotonic predictions with respect to $\Delta t_{\mathrm{h}}$. Furthermore, the hardness slightly decreases with the furnace temperature at $T_{\mathrm{f}}$ close to 930${}^{°}$C, which is not the behavior expected by the process expert either. Figure 6b shows the monotonic projection of the unconstrained model. It was computed according to the Appendix A on a grid G consisting of 40${}^4$ points. The monotonic projection exhibits the kinks that are characteristic of that method, and it yields an RMSE of 28.84 HV on the entire dataset.

With an overall RMSE of 10.14 HV, the model resulting from SIAMOR is more accurate for this application. A corresponding response surface is displayed in Figure 6c. In keeping with (18), monotonicity was required with respect to $T_{\mathrm{f}}$ (increasing), $\Delta t_{\mathrm{h}}$ (decreasing) and $t_{\mathrm{q}}$ (increasing). As $m=6$, there are $N_m=210$ model parameters and the dicretization $X^0$ was initialized with a grid using four equidistant values per dimension. The algorithm converged in iteration 246 with 1372 final constraints. Our first try was with only two monotonicity requirements (namely, with respect to $T_{\mathrm{f}}$ and $t_{\mathrm{q}}$). We observed, however, that the final number of iterations decreased when the third monotonicity requirement was added. Thus, the monotonicity requirements in each direction promoted each other numerically within the algorithm and for the used data. This reduction in the number of iterations was not accompanied by a decrease in total calculation time because more lower-level problems have to be solved when there are more monotonicity directions.

With SIAMOR, monotonicity was achieved in all three input dimensions where it was required. See, e.g., Figure 6c, which is the monotonic counterpart of Figure 6a. A comparison of Figure 6a–c clearly shows how incorporating monotonicity expert knowledge helps compensate data shortages. Indeed, taking no monotonicity constraints into account at all (Figure 6a), we obtained an unexpected hardness minimum with respect to $\Delta t_{\mathrm{h}}$ at $\Delta t_{\mathrm{h}}\approx 2.5$s and small $T_{\mathrm{f}}$. This also resulted in unnecessarily low predictions of the monotonic projection for small $T_{\mathrm{f}}$ and $\Delta t_{\mathrm{h}}⪆1.5$s in Figure 6b. The SIAMOR model (Figure 6c), by contrast, predicted more reasonable hardness values in this range without needing additional data, because it integrated the available monotonicity knowledge in the training phase.

For the SIAMOR plots in Figure 7, $\Delta t_{\mathrm{h}}$ was reduced to 0 s. This figure shows that monotonicity is also achieved with respect to $t_{\mathrm{q}}$. Without having explicitly demanded it, the hardness y shows the expected concave growth towards saturation with respect to $t_{\mathrm{q}}$ in Figure 7a. An additional increase in $F_{\mathrm{p}}$ leads to Figure 7b, where the sign of the second derivative of y with respect to the quenching time $t_{\mathrm{q}}$ changes along the $T_{\mathrm{f}}$-axis. I.e., the model changes its convexity properties in this direction and increases convexly instead of concavely with $t_{\mathrm{q}}$ at high $T_{\mathrm{f}}$, $\Delta t_{\mathrm{h}}=0$s and $F_{\mathrm{p}}=2250$kN. This contradicts the process expert’s expectations. A possible way out is to measure additional data (e.g., in the rear left corner of Figure 7b), which is elaborate and costly, however. Another possible way out is to add the concavity requirement ${\partial }_{x_4}^2{\widehat{y}}_{\mathit{w}}\left(\mathit{x}\right)\le 0$ for all $\mathit{x}\in X$ with respect to the $x_4=t_{\mathrm{q}}$ direction to the monotonicity constraints (6) used exclusively so far. In order to solve the resulting constrained regression problem, one can use the same adaptive semi-infinite solution strategy, which was already used for the monotonicity constraints alone.

5. Conclusions and Outlook

In this article, a proof of concept was conducted for the method of semi-infinite optimization with an adaptive discretization scheme to solve monotonic regression problems (SIAMOR). The method generates continuously differentiable models, and its use in multiple dimensions is straightforward. Polynomial models were used, but the method is not restricted to this type of model, even though it is numerically favorable because polynomial models lead to convex quadratic discretized problems. The monotonic regression technique was validated by means of two real-world applications from manufacturing. It resulted in predictions that complied very well with expert knowledge and that compensated for the lack of data to a certain extent. At least for the small datasets considered here, the resulting models predicted the training data more accurately than models based on the well-known projection or rearrangement methods from the literature.

While the present article is confined to regression under monotonicity constraints, semi-infinite optimization can also be exploited to treat other types of shape constraints such as concavity constraints, for instance. In fact, the shape constraints can be quite arbitrary, in principle. Additionally, this is only one of several aspects in the field of potential research on the method opened up by this work. Others are the testing of SIAMOR in combination with different model types, datasets or industrial processes. When using Gaussian process regressors instead of the polynomial models employed here, one can try out and compare various kernel types. Additionally, the SIAMOR method can be extended to locally varying monotonicity requirements (i.e., ${\sigma }_j={\sigma }_j\left(\mathit{x}\right)$).

Another possible direction of future research is to systematically investigate how to speed up the solution of the global lower-level problems. When more complex models or shape constraints are used, this will be particularly important. The solution of multiple lower-level problems and the final feasibility test on the reference grid can be parallelized to reduce the calculation time, for example. A rigorous investigation of the convergence properties and the asymptotic properties of the SIAMOR method and its possible generalizations is left to future research as well.