Process optimization via confidence region: a case study from micro-injection molding

In industrial research, experiments are designed to determine the optimal factor levels of the process parameters. Typically, experimental data are used to fit empirical models (for example, regression models) to derive one set of optimal conditions that maximize (or minimize) the response. Unfortunately, the optimization rarely provides a Confidence Interval for the location of the optimal solution, even though the optimal solution itself is subjected to variability. From a practitioner's point of view, identifying a region of possible optimal values provides high operational flexibility to adjust process parameters online during production. This paper provides a procedure for computing a confidence region for the optimal point based on experimental data, bootstrapping, and data depth. The procedure is validated using a case study from micro-injection molding, where the part weight is maximized under a constraint of the probability of flash formation. The proposed method considers that the objective function (part weight) and the constraint (probability of flash formation) are estimated from experimental data and subjected to sampling variability.


Introduction
The objective of experimental design is to build empirical models of the process to identify the relevant process factors and, eventually, the optimal levels of the processing conditions. However, the processes often considered in technological applications are "noisy," meaning that the same input parameters will not result in identical response values. This is called sampling variability. Usually, sampling variability is not considered in the optimization phase, and only one set of process parameters is selected as optimal. Still, the process cannot repeat itself means that the optimal value is also subjected to sampling variability. Therefore, it is possible to derive a Confidence Region, CR, of the location of the optimal point. The interest in building Confidence regions for the location of the optimal point From a practitioner's point of view, knowing only one set of parameters might be too conservative. Boundary conditions might change, and machine performances vary during their lifetime when running in a steady state. Identifying a region of possible optimal values could be more interesting and provide higher flexibility in a manufacturing environment. Often, machine operators change the process parameters online during production to slightly adjust the settings based on their experience. Leaving the choice entirely to the operators could be dangerous. To this extent, it would be favorable to provide a region of optimal processing conditions that could be used to adjust the process.
Building a Confidence Region for the location of the optimal solution is a well-known problem in statistics. Early work by Box and Draper (1954) focused on the Confidence region of a stationary point for a response surface function. In the case of a polynomial model and normality assumption, procedures were proposed in the literature (Cahya et al., 2004;Peterson et al., 2002;Wan et al., 2016), and these methods were applied in drug testing (Carter et al., 1982) and biology (Brooks et al., 2005). Del Castillo et al. (2020) generalized the previous results providing a framework to identify CR for the location of the optimal point for distribution-free models.
However, multiple objectives must be optimized simultaneously in the manufacturing process, such as costs, performances, and manufacturability. In this case, we refer to Multi-Objective Decision Making (MODM) problems (Nuno Ricardo Costa & Lourenço, 2016). In the literature, several approaches are available to solve MODM: desirability functions, grey analysis, and goal programming to cite the most important ones. In desirability analysis (Nuno R. Costa et al., 2011;Derringer & Suich, 1980), each response is converted into a desirability function that varies over the range [0,1]; if the response is close to the target, the desirability function approaches 1, and if the response is outside an acceptable region, the desirability function is 0. The desirability approach was applied to different manufacturing processes (Camposeco-Negrete, 2020;Gupta et al., 2019). Grey relational analysis (Ju-Long, 1982) associates each experimental data with a grey relational coefficient and grey relational grade based on how each data behaves compared to the expected value for each response. A greater grey relational grade corresponds to better performance. Grey analysis was often applied in combination with the Taguchi design of the experiments (Palanikumar, 2011) (Kuram & Ozcelik, 2013). Desirability and grey analysis transform a multidimensional problem into a one-dimensional optimization problem using the desirability function and the grey relation grade. In Goal programming, all the objectives are assigned target values and a relative priority on achieving these values. Goal programming treats these values as "targets to aspire for" (Reklaitis et al., 1983). The algorithm attempts to find a solution as close as possible to all target values following the priorities given by the experimenter. All these methods, and others, are viable solutions to solve MODM problems. However, to the author's knowledge, none of them provides information on the Confidence Region of the optimal solution found but only one set of optimal process parameters.
Recently Del Castillo (2020) addressed the problem of constructing the CR of the optimal point considering a Multi-Objective problem under uncertainty. However, only the objective function was subjected to error, while the constraint was deterministic. In an industrial environment, the constraint could be estimated from experimental data. For example, it could represent a probability of defect formation as a function of the process parameters and, therefore, subject to sampling variability.
The literature review showed no methodologies involving the evaluation of CR of the optimal point in MODM problems where all the objectives are empirical functions subject to sampling variability. For this reason, this work focuses on providing an optimization procedure that, starting from an experimental approach, allows to determine a Confidence Region for the location of the optimal set of parameters when multiple responses need to be optimized at the same time, and all of them are estimated from experimental data.
The proposed procedure is validated on a case study considered regarding the micro-injection molding process (μIM). μIM represents a manufacturing process that allows the production of micro-components with low tolerances at a relatively low cost. It finds applications in the biomedical, automotive, and microfluidic sectors. The μIM is similar to the standard injection molding process. Nevertheless, major criticalities arise when miniaturizing the process. The microinjection molding process is challenging to optimize because it operates in extreme conditions (Eladl et al., 2018). Furthermore, the quality of a μIM component is usually associated with the ability to fill the cavity during the process, resulting in the final component's weight. Researchers tried to optimize the μIM process by maximizing the weight of the part, as shown in several works (Attia & Alcock, 2011;Eladl et al., 2018). The high weight of the part implies the lack of voids and shrinkage phenomena. The process parameters that affect the part weight are known, and they have been investigated. The most studied parameters were mold temperature, melt temperature, injection speed, holding pressure, and cooling time. However, if not well-calibrated, the conditions that ensure a complete cavity filling produce parts with a well-known flash defect. Flash is the formation of unwanted additional material between the two molds due to the pressure of the injection material that exceeds the clamping force of the mold (Eladl et al., 2018). The process parameters that describe the formation of flash in μIM are the same that define the part weight, so one experimental campaign can be designed to study both part weight and flash formation simultaneously. The μIM is a suitable case study for investigating the advantage of using Confidence regions to optimize industrial-technological problems. As already described, μIM works in extreme conditions, and therefore, the optimization process is challenging. The knowledge of a set of optimal parameters, rather than just one condition, is helpful in an industrial environment to obtain higher operational flexibility. The optimization is constrained to the formation of the flash defect. The optimal solution might fall on or near the boundary when solving a constrained problem, resulting in dangerous conditions. Using a confidence region provides a broader knowledge of the feasible conditions that satisfy the constraint.
This work provides a procedure for evaluating the confidence region of the optimal solution for MODM problems. Contrary to the currently available methods for process optimization in MODM problems, the proposed methodology provides a region of the parameters containing the optimal solutions of the problem with a determined probability. Among all the parameters belonging to the CR, the practitioner can select the most suitable one based on the specific needs, which makes the proposed approach useful primarily in an industrial environment. We propose two versions of the same procedure used in an industrial application. One procedure is called the Improved Standard procedure (ISp), and the second is the Generalized Optimal Region procedure (OGRp). The two procedures use the same input data to obtain a CR for the location of the optimal solution; however, they provide a different approach to selecting the optimal set of processing conditions. The procedures exploit experimental data, bootstrap techniques (Adjei & Karim, 2016), and data depth (Pokotylo et al., 2016) to build the Confidence Region. The paper is organized as follows. In section "Case study: micro injection molding" the μIM is described, and the experimental design is used to estimate the regression equations for weight and flash probability as a function of the process parameters. In section "Description of the optimization procedure", the optimization methodology is described. The results of the two procedures are presented in section "Results" and a discussion of the methodologies is proposed in section "Discussion".

Case study: micro injection molding
The machine used for experimentation is the DESMA Tec Formica Plast 1 K ( Fig. 1 a), and the material was POM (Polyoxymethylene). The part identified for the experimentation is a benchmark for micro injection moulding, a double thin plate component with a thickness of 500 μm and an aspect ratio of about 8, Fig. 1b).
The process parameters considered were melting temperature T melt, holding pressure P hold , mold temperature T mold , injection speed v inj, and holding time t hold . For simplicity, in the analysis, we consider x as the vector of the coded process parameters (x 1 , x 2 , x 3 , x 4 , x 5 ), while natural variables are indicated with the name of the factor (T melt , P hold , T mold , v inj , t hold ). A summary of the variables and the nomenclature is reported in Table 1. The levels of the process parameters are reported in Table 2. The response variables are part weight and the presence of flash on the part (1 flash, 0 no flash); part weight is a continuous variable while the flash is a binary variable.
A central composite design CCD (Montgomery, 2019) was carried out to determine the influence of the process parameters on the part weight and the flash formation. The CCD design was selected because curvature was expected in the responses. A CCD design is composed of a 2 5 factorial points, with the addition of 10 axial points. Additionally, 20 replicates of the center point were added to obtain a precise estimate of the process variability. The CCD was replicated two times, resulting in 104 experiments, (2 5 + 10) × 2 + 20 104. The significant factors affecting the probability of flash formation are reported in Table 3. A logistic regression model was fitted to describe the probability of flash formation as a function of the process parameters. The model in coded units is reported in Eq. (1).
(1) The conditions that did not result in flash formation were used to estimate the weight of the parts as a function of the process parameters. Out of 104 experimental runs, 33 parts showed the flash defect. So, 71 conditions were used to fit the part weight as a function of the process parameters. The significant factors for weight are reported in Table 4.
The regression equation for the weight in coded units is reported in Eq. (2).
The objective function depends only on melting temperature and holding pressure, while the constraint is a function of four process parameters (mold temperature, melt temperature, injection speed, and hold pressure). Hold time does not affect the weight or the flash formation. The selection of the optimal solution is achieved through the definition of a utility function, defined as If the probability of flash formationP F L AS H is large (Eq. (1) is close to 1), then the value of the utility function U x, β, γ becomes close to 0 whatever the value of the weight of the part, that is precisely what we are aiming at: avoid regions of the parameters where the probability of flash formation is large independently from the value of the objective function. On the contrary, when the probability of flash formation is close to 0, the utility function is a proxy of the weight of the part.
The condition that maximizes Eq. (3) is reported in Table  5.
The solution of the standard optimization states that x 3 and x 4 should be selected at their lowest level, as they influence only the constraint. Factor x 5 (holding time) does not influence both the process and the flash, so it is selected at the lowest value to reduce the processing time. Factor x 1 (melting temperature) is selected at the highest level, while x 2 (holding pressure) is set at a coded value of 0.7, which corresponds to 1350 bar in natural variables.

Description of the optimization procedure
We aim to build a confidence region for a constrained optimization problem where both the objective function and the constraint are estimated from experimental data. Two key points are proposed in this work: 1. The constraint of the problem is non-deterministic. Its coefficients are subject to bootstrap, so the proposed procedure can be used in cases where the objective function and the constraint are estimated from experimental data. In principle, the number of constraints that can be considered is not limited. 2. The trimming is carried out on the objective function and constraint coefficients rather than the optimal solutions. As both the objective function and the constraint are generated by bootstrapping, some combinations of a realization could result in optimal solutions with extreme combinations of process parameters that are unfeasible in an industrial environment.

Improved standard procedure (ISp)
We start by describing the Improved Standard procedure.
(ISp). The five steps of the ISp are described here, and they are schematized in Fig. 2a). In Step 1, the objective and constraint functions are estimated based on experimental data. All 104 experimental conditions were used to estimate the probability of flash formationP F L AS H x,γ as in Eq. (1). A subset of the experimental conditions was used to estimate the part weight, Eq.
(2). The subset consisted of all experimental conditions that did not result in flash formation.
In step 2, bootstrapping is used. Bootstrapping is a statistical resampling method used to create simulated samples from a single data set (Efron & Tibshirani, 1994).
We implemented the bootstrapping procedure proposed in (Adjei & Karim, 2016) for logistic regression using the original dataset, obtaining B sets ofγ * i . For each bootstrap iteration (i 1,…,B), the conditions that resulted in a lack of flash defect were considered for bootstrapping the objective function; the result is B sets ofβ * i . Thus, the subset of conditions that resulted in a lack of flash changed based on the simulation result for each bootstrap iteration. This approach was selected because it replicates the experiment and how Eqs. (1) and (2) were estimated in step 1. The output of this step is the bootstrap parameter region B defined as In Step 3, the trimming procedure is carried out using the projection data depth (Pokotylo et al., 2016) method. Data depth technique were initially introduced in multivariate data analysis to measure the centrality of a point x ∈ R d in a data set. The data depth algorithm assigns to any point in the data set a real number between [0,1] which defines its degree of centrality. The result is a non-parametric ordering of the points, where points near the center are assigned high data depths values. Data depth found application in detecting outliers in multivariate data analysis (Hubert et al., 2015). In this paper, the concept of data depth was used to eliminate the combinations of bootstrapped coefficients γ * i ,β * i that resulted in large values of data depth. The B coefficients γ * i ,β * i were ordered by trimming the α% outermost. This yields the trimming parameter region I, and the set of indexes of the coefficients that belong to the region is called In step 4, the utility function is maximized. A utility function is defined as the product between the weight of the part and the complement to 1 of the probability of flash formation. The utility function is the following: The process parameters that maximize the utility in the range − 1 < x < 1 are obtained ∀i ∈ I B .
The optimization problem is the following: The results are (1-α) B optimal solutions x * opt, i , so x *

opt
(1 − α)B In Step 5, the parameter confidence region is built. The convex hull that contains all the optimal solutions x * opt, i determines the Confidence Region.

Generalized optimal region procedure (GORp)
It is possible to generalize the ISp procedure to provide additional information on the realization of the utility function in the space of the parameters x; a schematization is illustrated in Fig. 2b. The Confidence Region C R 1−α resulting from ISp did not provide practical information regarding selecting the parameters inside the region itself, as shown in the next Section. For this reason, a modified procedure was selected Fig. 2 Workflow of the Improved Standard procedure IS procedure (a) and generalized optimal region GOR procedure (b) to help the practitioner identify an optimal region of process parameters. The Generalized Optimal Region procedure (GORp) aims to provide information regarding the combinations of process parameters x that ensure a specific value of utility, thus supplying technological-based criteria for selecting the region of the process parameters.
For generalizing the procedure, Steps 4 and 5 were modified as follows. In Step 4, the utility function was evaluated over a grid of points, therefore N 41. As described in Sect. "Bootstrap results", the optimization problem can be simplified considering only two out of 4 process variables (x 5 is not significant); for this reason, the grid is defined in two dimensions. However, it is possible to generalize the procedure by considering a multi-dimensional grid. We identify x jk as a single point on the grid, where j 1,…,N and k 1,…,N. The utility function is computed as follows: U is therefore a matrix with size N × N × (1 − α)B. In Step 5, the optimal region is identified. The utility function must be maximized; however, B(1-α) utility functions result in (1-α)B optimization problems. This means that for each grid value x jk , the utility has B(1-α) different values because the bootstrapped coefficients γ * i ,β * i are different at each iteration i. For the selection of the optimality region, the median value of the median utilityũ jk at each point x jk is used: Therefore, the optimal region R is the set of the grid points x jk resulting in a median utilityũ jk higher than a specific value q.

R
x jk ũ jk > q , j 1, ..., N , k 1, ..., N The value q represents the median value of the expected weight in a specific operating condition. So, the optimal region is the set of process parameters that ensure that 50% of the utility is higher than the selected q value. The practitioner can select the value q based on the specific application, providing a more flexible tool than ISp. The resulting combination of process parameters is not technically a confidence region but rather an optimal region that guarantees a minimum value of median utility.
The selection of the median as an index will be discussed in Sect. "Results of the standard procedure".
The following parameters were set B 100,000 and α 5% for the numerical simulation. Bootstrap procedures and utility functions evaluations were carried out using Matlab®; data depth analysis was done using the "ddalpha" R-package (Pokotylo et al., 2016).

Results
In this section, the results of the two proposed procedures are reported. Firstly, a general overview of the results of the Bootstrapped coefficients is carried out. Later, the results of ISp and GORp are illustrated.

Bootstrap results
By inspecting Eq. (1), we find that ifγ * 4 (coefficient of x 3 ) is positive and ifγ * 6 (coefficient of the interaction x 24 ) is negative, then the solution of the optimization problem requires that x 3 − 1 and x 4 − 1. By looking at the B bootstrapped coefficients, these two conditions are always met: the bootstrap did not result in negativeγ 4 and positiveγ 6 coefficients. So, to minimize the probability of flash formation x 3 (mold temperature) and x 4 (injection speed) should be set at their lowest level. So, from now on all graphs and considerations are based on setting (x 3 , x 4 ) (−1, − 1). This solution is coherent with the one found in the standard optimization problem described in Sect. "Case study: micro injection molding".
The optimization problem requires only the optimization of x 1 (melt temperature) and x 2 (holding pressure) because they both influence the objective function and the constraint. In detail, as x 1 and x 2 increase, the weight increases, but at the same time also, the defect probability increases.

Results of the standard procedure
In the standard procedure, an optimization problem is solved for each realization i of the trimmed bootstrap coefficients, i.e., ∀i ∈ I B . The optimization problem in Eq. (5) involves Some optimal solutions are more frequent than others, as shown in the histogram in Fig. 4b. The condition (x 1 , x 2 ) (− 1,1) is the most frequent optimal solution, followed by (x 1 , x 2 ) (1,1). In Fig. 4b the histogram shows the probability that a condition is selected as optimal by the ISp. For graphical purposes, the optimal solutions were grouped. For example, in the yellow bar are summed all the solutions that fall in a small range of parameters (− 1 ≤ x 1 < -0.95, 0.95 ≤ x 2 < 1). The conditions with low x 1 and high x 2 are the most likely to be selected as optimal. The solutions in the upper-right corner, i.e., x 1 x 2 1, indicate that in some realizations of the bootstrap procedure, the constraint did not belong to the region of the parameters (− 1 ≤ x ≤ 1). The combination of parameters (x 1 , x 2 ) (1,1) was selected as optimal 7505 times out of 95,000. In this condition, the utility maximization is equal to maximizing the weight as the probability of flash formation is zero. As a matter of fact, by looking at Eq. (2), the weight increases as both x 1 and x 2 increase. In general, the conditions where x 2 is set near its maximum value show a larger probability of being optimal. The histogram in Fig. 5 shows the utility value of the optimal solutions u * opt . The utility varies in a small range of values (90 to 101 mg) despite obtaining a wide set of optimal values x * opt .
In conclusion, the confidence region C R x * opt 1−α for the location of the optimal solution is the convex hull of the points in Fig. 4 a). The convex hull contains the optimal solution to the standard problem, (x 1 , x 2 ) (1, 0.7). This solution is reported in Table 5 and marked as a green triangle in Fig. 4a. This combination is optimal in 435 cases out of 95,000 according to ISp. The CR seems too large to be useful in an industrial environment. Based on this procedure, the practitioner cannot derive information on the utility for each point that belongs to the CR; this is because the utility of each point depends on the specific realization i of the bootstrap (for each i ∈ I B the values of the coefficientsγ * i andβ * i change). Especially for high-frequency optimal conditions (for example, x 1 x 2 − 1), it is impossible to derive information on the utility's distribution at a specific value of x. Moreover, it is not advisable to compare the utility estimation for low-frequency solutions and high-frequency solutions. An additional criticality of the ISp solution is that the condition with the highest probability of flash formation (x 1 , x 2 ) (1,1) belongs to the confidence region C R x * opt 1−α . This is because, by solving the optimization problem individually, 10% of the time (x 1 , x 2 ) (1,1) is the optimal solution, as shown in Fig. 4b. This result allows the practitioner to select a parameter combination that leads to a high flash formation probability.
In conclusion, the ISp presents criticalities in the resulting CR for the proposed technological problem. For these reasons, a second procedure was derived, GORp. The GORp provides an iso-utility region, i.e., a sub-region of the parameters that result in the same median utility.

Results of the GORp
Compared to the previous case, for each point of the grid x jk , (1 − B) α utility values are available and used to build the optimal region. In this case, no optimization problem is directly solved, but each point's median utility is considered for building the optimal region R.
For clarity purposes, the distribution of the weight and probability of flash formation for two selected points are shown in Fig. 6 (the points are x 1 x 2 0.8 and x 1 x 2 1). For condition x 1 x 2 0.8, the probability of flash formation varies between 0 and 0.4. However, around 8500 simulations resulted in a lack of flash formation, less than 10% of the generated data. The weight distribution in Fig. 6b can be approximated by a gaussian distribution centered on 97 mg. The resulting utility function for this specific combination of parameters is shown in Fig. 6c. The utility has a left-skewed distribution due to the shape of the flash probability distribution. For the upper right condition at x 1 x 2 1, the weight distribution is similar to the one shown in Fig. 6a. However, the mean weight is 100.2 mg, higher than 97 mg (the part weight increase is expected from Eq. (2)). The flash probability increases greatly compared to x 1 x 2 0.8; the probability varies between 0 and 0.7 because the higher values of melting temperature and holding pressure result in a higher probability of flash formation. The result is a large range of utility values for this condition (x 1 x 2 1), as visible in Fig. 6f). This results in a heavily skewed utility distribution that indicates a high-risk condition from an industrial point of view. A low utility indicates a defective part because the utility is reduced in the presence of flash.
The median is selected as an index to be maximized for deriving the optimal region of the parameters R. As the utility distribution is highly skewed, the median is the best representation of the utility for each combination of parameters. For each point of the grid, the median utilityũ jk was evaluated, and Fig. 7 represents its distribution in the region of the parameters.
At low levels of both parameters, the median utility is low because of the low weight of the part; however, in this region, the probability of flash formation is close to 0. As x 2 increases, the utility increases, reaching its maximum when x 2 is larger than 0.8. However, in the upper right corner (x 1 x 2 1), there is a sharp decrease in the median utility, mostly due to the high probability of defect formation, as seen in Fig. 6 e). According to the GORp, when both x 1 and x 2 move toward their highest values, the probability of flash formation increases (see Eq. (1). This behavior is more penalized in the optimization phase in the GORp because it considers all the realization of the bootstrap (after trimming) simultaneously. In only 10% of the cases, the flash constraint does not belong to the optimal region, and therefore the optimal solution is (x 1 , x 2 ) (1,1); in all the other cases, the utility of this combination of parameters is reduced by the high probability of flash formation. When the generalized procedure evaluates the median distribution of the utility in (x 1 , x 2 ) (1,1), this value is affected by the presence of a heavier left tail due to the high probability of flash formation, as illustrated in Fig. 6f. The maximum value of the median utilityũ jk is equal to 97.12 mg. So, we select 96 mg as the reference value q for constructing the optimal region. All the conditions x that show a median utility higher than 96 mg belong to the optimal region, R and these points are colored in yellow in Fig. 7b. This region is characterized by a very high value of x 2 [0.9,1] and a quite large interval for x 1 [− 1, 0.4]. The optimal region of the GORp is contained in the confidence region obtained for the ISp. GORp determines a much smaller region of the parameters, and it also provides information on the utility values expected in this specific region. In the ISp case, the final output of the procedure was a set of points (x 1 , x 2 ) resulting from an optimization problem. The utility values of the optimal points in the ISp are not important in the final assessment because they depend on Fig. 6 Examples of distributions of weight (a), flash probability (b), and utility function (c) for a specific point in the region of the parameters. In this case T melt 0.8 and P hold 0.8 the specific iteration of the bootstrap, i.e., on the specific values of the parameters γ * i ,β * i . On the contrary, the GORp determines the optimal region by looking at the utility values distribution in all the bootstrap realizations. It is important to note that utility is directly related to the quality of the part, as it represents a compromise between the part weight and the probability of flash formation. Therefore, the GORp is a tool that can be used to determine a confidence region for industrial applications because the value q has important technological implications. The optimal solution to the standard problems does not belong to the optimal region R obtained by setting q 96 mg. The standard optimal solution, shown as a green triangle in Fig. 7b, is characterized by a median utility equal to 95.1 mg. This results from the fact that the standard problem results from only one set of coefficients γ ,β , while the GORp considers all the bootstrap results after trimming.

Discussion
In this work, we have shown the possibility of exploiting the concept of CR to solve optimization problems in the industrial environment. The procedure relies on experimental data and the estimation of regression models to identify the objective function and process constraints. Two methods were developed: ISp and GORp. The ISp procedure generates a real confidence region, that is, a region that has a (1 − α)% probability of containing the true optimum of the problem. The region obtained from this procedure, however, cannot be used in a production environment substantially for two reasons: 1. it is too large. The bounding box of all (1 − α)B solutions results in a range of melting temperature [− 1,1] and a range of holding pressure between [0.6,1]. 2. it contains a process condition (x 1 , x 2 ) (1,1) that has a very large probability of producing a defective part For this reason, a generalization of ISp was derived using the same experimental data and the same regression relationships. The new procedure has been defined as GORp, and it cannot be strictly defined as a Confidence Region because it does not provide information as "there is a (1 − α)% probability that this region contains the optimal solution". GORp exploits technological considerations to obtain a region of parameters that allows obtaining a median value of utility at least equal to a value specified by the experimenter. GORp avoids recommending non-optimal conditions from a technological point of view and has the advantage of indicating the median value of the response in the optimal region.
Since both procedures are based on experimental data and the estimated regression models, the quality of the region inevitably depends on the quality of the experiment. It is recommended to carry out an experimental design that gives the possibility of estimating quadratic effects of the parameters and not only linear ones.
Compared to the literature, the methodology proposed in this work allows obtaining a region of process parameters instead of a single combination. The standard methods used in MODM, such as desirability, grey analysis, and goal programming, provide only one set of optimal process parameters. No information is provided on the effect of a slight variation of these parameters on the quality of the final parts.
In this regard, the benefits of the CR approach are remarkable as the result provides the practitioner an indication of how much process parameters can be varied before a variation in the process becomes important. For example, selecting (x 1 , x 2 ) (0, 0.9) we expect a median utility of 96 mg (Fig. 7), while (x 1 , x 2 ) (0.5, 0.8) results in a median utility of 95 mg. Depending on the quality requirements of the specific component, GORp provides all possible combinations of process parameters that meet the requirement.

Conclusions
In this work, a new framework to build confidence region for the optimal solution of a technological problem is provided.
The procedure starts with experimental data used to estimate two functions: an objective function and a constraint. These two equations are used in combinations with bootstrap and data depth to build the confidence region. The procedure was applied to an industrial problem in micro-injection moulding, where the weight of the part must be maximized under the constraint that a defect, called flash, occurs. For the first time, such a procedure was applied to a technological problem, and, additionally, the procedure was generalized also to consider a stochastic constraint. Two alternative procedures were proposed (ISp and GORp). The ISp resulted in a large confidence region which is not tolerable in an industrial environment. Additionally, high-risk combination of parameters (x 1 , x 2 ) (1, 1) belonged to the optimal region. For these reasons, the GORp was implemented. The GORp determines the combinations of parameters that ensure the maximum median utility.
The confidence region of the GORp indicates that holding pressure should be set at high values (> 1250 bar) while melting temperature should be preferably selected between 190 and 220°C. The two alternative procedures proposed select almost similar confidence regions for the optimal value; however, the generalized one penalizes the large probability of flash formation when both holding pressure and melting temperature are set at their highest values. The GORp also provides precise information on the utility value for each combination of the parameters, which could not be obtained in the ISp because the frequency of the optimal solutions varied inside the confidence region. The GORp procedure has a great advantage for the practitioner because it provides a region of the parameters where the expected utility is known based on the choice of the value q.

Symbol
Definition Valuê P F L ASH x,γ Logistic regression equation describing the probability of flash formation based on experimental data Equation (1) w x,β Least square estimate of the weight of the part based on experimental data Equation (2) x Vector of the process parameters: x (x 1 , x 2 , x 3 , x 4 , x 5 )  (4) U(x jk ,γ * i ,β * i ) Utility function as defined in the GORp Equation (6) C R

Declarations
Conflict of interest All co-authors have seen and agree with the contents of the manuscript and there are no conflicts of interest and financial disclosures to report. The authors also claim that none of the material in the paper has been published or is under consideration for publication elsewhere.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.