Yulong Qiao

Jinsheng Sun

Philippe Castagliola

Xuelong Hu

Optimal Design of One-sided Exponential EWMA Charts Based

Keywords: Exponential chart, Median run length, Expected median run length, Exponentially weighted moving average chart

come

Introduction

Quality of a product is important for manufacturing or service organizations to improve their competitiveness (see [START_REF] Oakland | Statistical Process Control[END_REF]). Statistical Process Control (SPC) plays an important role in the quality improvement of a product. As a useful tool in SPC, control charts are usually used to monitor the state of processes.

Control charts are usually designed under the assumption that the quality character-istic follows a normal distribution. This assumption is not appropriate in the monitoring of some particular processes, such as in monitoring the rate of occurrences of congenital malformations (see [START_REF] Lucas | Counted data CUSUM's[END_REF]), in monitoring the consecutive radiation pulses (see [START_REF] Luo | CUSUM analyses of time-interval data for online radiation monitoring[END_REF]), and in monitoring the medical errors (see [START_REF] Dogu | Monitoring time between medical errors to improve health-care quality[END_REF]), where the exponential distribution is usually assumed for the time between events (TBE) data. [START_REF] Lucas | Counted data CUSUM's[END_REF] first investigated the properties of the exponential chart. Following [START_REF] Lucas | Counted data CUSUM's[END_REF], many research works on exponential charts were performed. [START_REF] Gan | Design of optimal exponential CUSUM control charts[END_REF] developed an optimal design procedure of the exponential Cumulative Sum (CUSUM) control charts. Based on the average time to signal (AT S), [START_REF] Liu | A comparative study of exponential time between events charts[END_REF] compared the properties of several exponential charts.

The Phase I implementation of the exponential chart was further considered by [START_REF] Jones | Phase I control charts for times between events[END_REF]. More related research works on exponential charts can be found in [START_REF] Santiago | Control charts based on the exponential distribution: adapting runs rules for the t chart[END_REF], Yang et al. (2015), Zhang et al. (2005Zhang et al. ( , 2006)), to name a few. These researches are mainly focused on two-sided exponential charts. In fact, one-sided charts are more appropriate for processes for which the direction of the potential mean shift is anticipated (see [START_REF] Shu | A one-sided EWMA control chart for monitoring process means[END_REF]). Some existing researches focused on one-sided exponential charts can be found in [START_REF] Gan | Designs of one-and two-sided exponential EWMA charts[END_REF], [START_REF] Ozsan | Properties of the exponential EWMA chart with parameter estimation[END_REF], [START_REF] Pehlivan | Impact of model misspecification on the exponential EWMA charts: a robustness study when the time-between-events are not exponential[END_REF]. Moreover, if a two-sided chart is needed, a combined scheme of both the lower-and upper-sided charts could be designed instead, and there is practically no difference between the combined scheme and the two-sided case (see [START_REF] Reynolds | CUSUM charts with variable sampling intervals[END_REF]).

Exponentially Weighted Moving Average (EWMA) type charts gained extensive attention from many researchers for its advantage in detecting small to moderate shift sizes in a process (see [START_REF] Chandrasekaran | Modeling and analy-sis of EWMA control schemes with variance-adjusted control limits[END_REF], [START_REF] Gan | Computing percentage points of run length of an exponentially weighted moving average control chart[END_REF]), Waldmann (1986)). Much researches on EWMA type charts can be found in [START_REF] Haq | A new nonparametric synthetic EWMA control chart for monitoring process mean[END_REF], [START_REF] Knoth | Run length quantiles of EWMA control charts monitoring normal mean or/and variance[END_REF], [START_REF] Maravelakis | Run length properties of run rules EWMA chart using integral equations[END_REF] and so on. These researches are mainly based on the assumption of normally distributed observations. For non-normally distributed observations, researches have also been studied to investigate the performance of EWMA type charts. In high-quality processes, TBE data are sometimes assumed to be exponentially distributed. Based on this fact, [START_REF] Gan | Designs of one-and two-sided exponential EWMA charts[END_REF], [START_REF] Gan | Computing average run lengths of exponential EWMA charts[END_REF] developed the exponential EWMA control chart to monitor the rate of occurrences of rare events and presented a program to determine the average run length (ARL) profiles of the exponential EWMA chart. [START_REF] Liu | A study of EWMA chart with transformed exponential data[END_REF] proposed a new EWMA chart to monitor the exponentially distributed TBE data using the double square root (SQRT) transformation. [START_REF] Pehlivan | Impact of model misspecification on the exponential EWMA charts: a robustness study when the time-between-events are not exponential[END_REF] investigated the sensitivity of the exponential EWMA chart in detecting the process shift and they concluded that the exponential EWMA chart could be extremely robust to departures from the assumed exponential distribution.

The conditional expected values were employed to monitor small mean shifts in [START_REF] Raza | EWMA control chart for poisson-exponential lifetime distribution under type I censoring[END_REF]. By transforming the exponential data to approximately normal ones, [START_REF] Khan | A EWMA control chart for exponential distributed quality based on moving average statistics[END_REF] proposed an EWMA chart to monitor the exponentially distributed data. The effect of parameter estimation on the performance of the exponential EWMA chart has been analyzed in [START_REF] Ozsan | Properties of the exponential EWMA chart with parameter estimation[END_REF]. A survey of the related research works on exponential EWMA charts can be found in [START_REF] Ali | An overview of control charts for high-quality processes[END_REF].

Generally, the ARL is used as a criterion to evaluate the performance of control charts.

Nevertheless, the only use of ARL has been criticized by some researchers (see [START_REF] Khoo | Performance measures for the Shewhart control chart[END_REF], [START_REF] Khoo | Optimal designs of the multivariate synthetic chart for monitoring the process mean vector based on median run length[END_REF], [START_REF] Lee | Optimal designs of multivariate synthetic |s| control chart based on median run length[END_REF], Tang et al. (2019), Waldmann (1986) and so on).

The ARL is regarded as the representation of the "central tendency" of the run length (RL) distribution, however, it is clearly not suitable as its shape becomes asymmetric. Furthermore, it has been pointed out that the RL distribution of control charts may be skewed, especially when the process is in-control. Therefore, the interpretation of the ARL is sometimes confusing for practitioners to evaluate a control chart's performance with skewed RL distributions. For example (see Section 2, Table 1), a lower-sided exponential EWMA chart with (λ -, h -) = (0.4, 0.2045) has an in-control ARL of 500, while 50% of in-control RLs is smaller than 348. When the process is out-of-control (c = 0.6), the ARL is 57.1, and 50% of the out-of-control RLs is smaller than 41. Instead, the median run length (M RL), defined as the median number of samples that are plotted on a chart before an out-of-control signal is issued, has no such interpretation problem. In the above example, the in-control and out-of-control M RL values are 348 and 41, respectively, which indicate that half of the in-control and out-of-control RLs are smaller than 348 and 41, respectively. For practitioners, the interpretation based on the M RL is more intuitive than the one based on the ARL.

Hence, more and more researchers promote using the M RL as a better criterion to evaluate the performance of control charts (see Lee andKhoo (2006, 2017a), [START_REF] Nyau | Optimal statistical design of variable sample size multivariate exponentially weighted moving average control chart based on median run-length[END_REF]), Teoh et al. (2016)). Moreover, practitioners may not have enough historical knowledge of shift sizes in a process. To circumvent this problem, the expected M RL (EM RL) is recommended as a measure of a chart's performance when practitioners cannot specify the exact value of the shift size (see Tang et al. (2019), Teoh et al. (2016), You et al. (2016)).

In fact, RL distributions of one-sided exponential EWMA control charts are also skewed, especially when the process is in-control (see Section 2). The M RL and EM RL are good alternatives to evaluate the properties of one-sided exponential EWMA control charts. As far as we know, there is no research conducted on the M RL or EM RL properties of exponential charts yet. The purpose of this paper is to fill this gap and to develop the optimal design procedures for both the upper-and lower-sided exponential EWMA charts based on M RL and EM RL. The rest of this paper is organized as follows: Section 2 introduces the one-sided exponential EWMA charts and briefly discusses problems arising from the only use of ARL.

Then, an optimal design procedure is developed for the upper-and lower-sided exponential EWMA charts based on M RL and EM RL in Section 3. The performances of these onesided exponential EWMA charts are investigated in Section 4 and the effect of parameter estimation is briefly discussed in Section 5. Finally, some conclusions and recommendations are given in Section 6. The RL properties of the chart are presented using a Markov chain approach in the Appendix.

One-Sided Exponential EWMA Charts

As in most reported research works, let us assume that the TBE variable X follows an exponential distribution with a single parameter θ. The probability density function (p.d.f.) f (x) and cumulative distribution function (c.d.f.) F (x) of the exponential distribution are given as follows:

f (x) = 1 θ e -x θ , x ≥ 0 (1) F (x) = 1 -e -x θ . x ≥ 0 (2)
The process is said to be in-control when θ = θ 0 , where θ 0 is the known in-control mean of the process. Otherwise, the process is out-of-control when θ = cθ 0 (c = 1), where c = θ/θ 0 is a constant that reflects the shift in θ 0 .

Since this paper focuses on one-sided charts, we introduce the following upper-and lower-sided exponential EWMA charts:

• If the goal is to detect an increase in the process mean (c > 1), the upper-sided exponential EWMA chart is used. The charting statistic for the tth sample is computed as follows:

Q + t = max(B + , (1 -λ + )Q + t-1 + λ + X t ), t = 1, 2, . . . (3)
where 0 < λ + ≤ 1 and B + are the smoothing parameter and the reflecting boundary, respectively. The initial value Q + 0 is generally constrained by B + ≤ Q + 0 < h + , where U CL = h + is the upper control limit of the chart. An out-of-control signal is issued when

Q + t > h + .
• If the goal is to detect a decrease in the process mean (c < 1), the lower-sided exponential EWMA chart is used. The charting statistic for the tth sample is computed as follows:

Q - t = min(B -, (1 -λ -)Q - t-1 + λ -X t ), t = 1, 2, . . . ( 4
)
where 0 < λ -≤ 1 and B -are the smoothing constant and the reflecting boundary, respectively. The initial value

Q - 0 is generally constrained by h -< Q - 0 ≤ B -, where LCL = h -is the lower control limit of the chart. An out-of-control signal is issued when Q - t < h -.
As in [START_REF] Gan | Designs of one-and two-sided exponential EWMA charts[END_REF] and [START_REF] Pehlivan | Impact of model misspecification on the exponential EWMA charts: a robustness study when the time-between-events are not exponential[END_REF], without loss of generality, we assume that the in-control mean θ 0 = 1 and the initial values

Q - 0 = Q + 0 = θ 0 = 1 (other situations
with specific initial values can be analyzed similarly). Furthermore, the reflecting boundary

B + (B -
) is used to ensure that the charting statistic is at most a certain distance away from the control limit h + (h -), which can improve the sensitivity of the chart (see [START_REF] Pehlivan | Impact of model misspecification on the exponential EWMA charts: a robustness study when the time-between-events are not exponential[END_REF]).

In order to draw attention to the use of the M RL metric, some ARL and percentile values of the RL distribution are presented in Table 1 when λ -∈ {0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1} and the in-control ARL (denoted as ARL 0 ) roughly equals to 500 for the lower sided exponential EWMA chart. In this article, the RL properties of one-sided exponential EWMA charts can be obtained using the Markov chain approach given in the Appendix. Moreover, the values of the reflecting boundaries are set to be B + = 0.5 and B -= 2.0 (discussed in Section 3.1).

It is obvious that there is a large difference between ARL 0 and in-control M RL (denoted as M RL 0 ) when the process is in-control. However, the difference between ARL and M RL values decreases with the shift size increasing. For example, if (λ -, h -) = (0.05, 0.6861), when the process is in-control (c = 1), the ARL 0 and M RL 0 are 500.1 and 351, respectively.

When the process is out-of-control (c = 0.2), the ARL and the M RL values decrease down to 10.2 and 10, respectively. Similar results can also be drawn from other values of (λ -, h -)

shown in Table 1. It is noted that the shape of the RL distribution of the lower-sided exponential EWMA chart changes with the shift size, i.e., from highly right skewed when the process is in-control to nearly symmetric when the process is out-of-control.

On the other hand, for the upper-sided exponential EWMA chart, Table 2 also shows the ARL and some percentiles of the RL distribution with λ + ∈ {0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1}

and the ARL 0 roughly equals to 500. It can be seen that the difference between ARL and M RL also changes with the shift size. For example, if (λ + , h + ) = (0.2, 2.2378), ARL 0 and M RL 0 are 500 and 347, respectively, when the process is in-control (c = 1). When the mean shift c changes to c = 5, ARL and M RL values are 3.1 and 3, respectively.

(Please Insert Tables 1-2 to Here)

For a skewed RL distribution, the interpretation solely based on the ARL is different from the one with respect to a nearly symmetric distribution. A practitioner should be advised against viewing ARL as a symbol of "half of the time" because ARL only provided information about the expected RL (see [START_REF] Khoo | Optimal designs of the multivariate synthetic chart for monitoring the process mean vector based on median run length[END_REF], [START_REF] Nyau | Optimal statistical design of variable sample size multivariate exponentially weighted moving average control chart based on median run-length[END_REF]). It is not intuitive enough and can sometimes result in misunderstandings when the RL distribution changes with the shift size. Fortunately, the percentiles of the RL distribution provide more information on the performance of a control chart, no matter the RL distribution is skewed or symmetric. As a good alternative, the M RL has been suggested by some researchers to evaluate the properties of control charts (see [START_REF] Khoo | Optimal designs of the multivariate synthetic chart for monitoring the process mean vector based on median run length[END_REF], Lee and Khoo (2017a), [START_REF] Nyau | Optimal statistical design of variable sample size multivariate exponentially weighted moving average control chart based on median run-length[END_REF]), Tang et al. (2019), Teoh et al. (2017), to name a few). In this article, for one-sided exponential EWMA charts, the M RL is also suggested as the criterion to be investigated.

3 Optimal Design of One-sided Exponential EWMA

Control Charts

In this section, we mainly focus on the optimal design of one-sided (upper-and lower-sided) exponential EWMA charts based on the M RL and EM RL, respectively. The setting of the reflecting boundaries B + and B -are discussed first, and the procedure for the design of the smoothing parameter and control limits of the chart are developed lately. for different values of λ -∈ {0.1, 0.3, 0.5, 0.7} and mean shifts c ∈ {0.1, 0.2, . . . , 1.0} when M RL 0 ∈ {100, 500, 1000}, for the lower-sided exponential EWMA chart. The control chart limit h -is selected to produce the desired M RL 0 for different combinations of λ -and B -.

Owing to the space limitation, only the numerical simulation with M RL 0 = 100 is presented in Table 3. It can be noted that the control chart limit h -generally increases as the reflecting boundary B -increases. When B -is larger than 2.0, the value of h -converges to a constant. For example, when λ -= 0.3, the value of h -converges to 0.348 with the increase in B -. Moreover, it can also be seen that, for fixed values of c, the out-of-control M RL value (denoted as M RL 1 ) decreases as the reflecting boundary B -increases, which implies that the sensitivity of the control chart increases. For values of B -larger than 2.0, it can be noted that the detection ability of the lower-sided exponential EWMA chart generally remains the same. For example, when λ -= 0.3 and c = 0.8, if the value of B -increases from 1.0 up to 1.5, the M RL 1 values decrease from 44 down to 40. While if the value of B -increases up to 2.3 or 2.5, the M RL 1 values are both equal to 39. For other values of the desired M RL 0 , similar results can also be drawn. Consequently, as in [START_REF] Gan | Designs of one-and two-sided exponential EWMA charts[END_REF], a reasonable value of the reflecting boundary B -for the M RL-based lower-sided exponential EWMA chart is 2.0. In the rest of this article, the reflecting boundary B -is set to be 2.0.

For the upper-sided exponential EWMA chart, similar numerical simulations are also performed to investigate the choice of the reflecting boundary B + for different values of M RL 0 ∈ {100, 500, 1000}. Owing to the space limitation, these results are not presented here. It can be concluded from the results that a reasonable value of the reflecting boundary B + is 0.5, which leads to a good detection ability of the upper-sided exponential EWMA chart for different values of λ + and M RL 0 . In this article, the value of the reflecting boundary B + is set to be 0.5.

(Please Insert Table 3 to Here)

Design procedure of one-sided exponential EWMA charts based on the M RL

When the process is in-control (c = 1), it is expected that the control chart runs with a small false alarm rate, i.e., a large M RL 0 . While at the same time, the chart is also expected to detect a mean shift quickly, i.e., a small M RL 1 , when the process is out-of-control (c = 1).

The optimal design of the one-sided exponential EWMA charts based on the M RL consists in finding the smallest M RL 1 value for a specified mean shift c (denoted as c opt ) with the constraint on the desired M RL 0 , i.e., (λ * , h * ) = argmin

(λ,h) M RL 1 (λ, h, c opt ), (5)
subject to

M RL(λ * , h * , c = 1) = M RL 0 , (6)
where λ and h are the smoothing parameter and the control chart limit, respectively. It is noted that λ stands for λ + or λ -for the upper-or the lower-sided exponential EWMA charts, respectively. Similarly, h stands for h + or h -for the upper-or lower-sided exponential EWMA charts, respectively. The upper-sided exponential EWMA chart is used to detect an increase in the mean shift (c opt > 1). On the contrary, the lower-sided exponential EWMA chart is employed if the specified mean shift decreases (c opt < 1).

For a specified value of M RL 0 and the interested mean shift c opt , one can obtain the optimal parameter combinations of the one-sided exponential EWMA chart using Eqs. ( 5)

and ( 6). The procedure is summarized as follows:

(1) Specify the desired M RL 0 and the interested mean shift c opt .

(2) Initialize the smoothing constant λ = 0.05.

(3) Determine the control limit h for the corresponding λ, based on the desired M RL 0 specified in Step (1).

(4) Compute the M RL 1 for the corresponding combinations (λ, h) specified in Step (3) and c opt .

(5) Repeat Steps (3) to ( 4) for all values of λ varying from 0.05 to 1 with a step size 0.001.

(6) Record the combinations (λ, h) that produce the minimum M RL 1 for the specified c opt as the optimal parameter combinations.

Since the M RL is a discrete quantity, there may exist several values of control limit h for each smoothing constant λ, based on the constraint on the specified M RL 0 . Without loss of generality, the median of all the h values is selected for the corresponding λ. It is noted that there may exist several optimal parameter combinations (λ, h) that produce the same minimum value of M RL 1 for a specified mean shift c opt . In this situation, the median of all the optimal λ values and the corresponding h are selected in this article.

Design procedure of one-sided exponential EWMA charts based on the EM RL

The optimal M RL-based chart proposed in this paper is effective in detecting a shift of particular size and the computation of the M RL requires the exact value of the mean shift c. However, the practitioners are not always able to specify the exact shift size in advance, but the range of shift is roughly known in practice (see [START_REF] Celano | Performance of t control charts in short runs with unknown shift sizes[END_REF]). In this case, the shift size c should be considered as a random quantity that may vary according to a certain distribution and the EM RL can be used instead of the M RL, i.e.,

EM RL = Ω M RL(λ, h, c) × f c (c)dc, (7)
where f c (c) is the p.d.f. of the mean shift c over the range Ω. It is generally assumed explicitly (see [START_REF] Domangue | Some omnibus exponentially weighted moving average statistical process monitoring schemes[END_REF], Sparks ( 2000)) or implicitly (see [START_REF] Reynolds | Should observations be grouped for effective process monitoring[END_REF]) that all process shifts occur with an equal probability, for example, the shift in the The optimal parameter combinations (λ * , h * ) of the one-sided exponential EWMA chart based on the EM RL can be obtained as follows:

(λ * , h * ) = argmin (λ,h) EM RL 1 (λ, h, Ω), ( 8
) subject to M RL(λ * , h * , c = 1) = M RL 0 , (9)
where EM RL 1 is the out-of-control EM RL value computed using Eq. ( 7). The optimal design procedure of the EM RL-based one-sided exponential EWMA charts is similar to the steps in Section 3.2 except that the interested mean shift c opt and the M RL 1 are replaced by the interested mean shift range Ω and the EM RL 1 , respectively.

Performance Evaluation Based on M RL

In this section, the performances of both the upper-and lower-sided exponential EWMA charts are studied. The optimal charting parameter combinations used here are obtained following the procedures described in Sections 3.2 and 3.3.

Performance evaluation for the mean shift of known size

The M RL profiles for the lower-sided exponential EWMA chart and the lower-sided exponential Shewhart chart are displayed in Table 4, for different values of the specified mean shifts c opt ∈ {0.2, 0.4, 0.6, 0.8} and the desired M RL 0 ∈ {100, 200, 300, 500}. The optimal parameter combinations (λ -, h -) of the lower-sided exponential EWMA chart and the lower control limit (LCL) of the corresponding lower-sided exponential Shewhart chart are given in rows 4 and 5.

It is shown that the lower-sided exponential EWMA chart with a small optimal λ -value has better performance for small shift sizes. At the same time, the chart with a large optimal λ -value performs better in the cases of large shift sizes. For example, considering the case of Furthermore, as in [START_REF] Kumar | Design and implementation of q th quantile-unbiased tr-chart for monitoring times between events[END_REF], a comparative study is carried out to compare the charts designed for a fixed nominal ARL 0 and the corresponding M RL 0 . Table 6 presents the M RL and ARL values of the chart for different parameter combinations of (λ -, h -) and c when M RL 0 = 350 and ARL 0 = 500. It is noted that the parameter combinations (λ -, h -) are obtained by optimizing the M RL-based chart (the upper part in Table 6) and the ARL-based chart (the lower part in Table 6) for different values of (Please Insert Tables 4-6 to Here)

M RL 0 =

Performance evaluation for the mean shift of unknown size

When the mean shift is unknown, but the range of the shift size can be roughly specified, the optimal parameter combinations and properties of the lower-sided exponential EWMA chart and the lower-sided exponential Shewhart chart are provided in Table 7 for several uniform distributions of shift size Ω ∈ {[0.1, 0.3], [0.3, 0.5], [0.5, 0.7], [0.7, 0.9]} when M RL 0 ∈ {100, 200, 300, 500}. For example, if M RL 0 = 200, the optimal parameter combinations (λ -, h -) and the corresponding EM RL 1 value of the lower-sided exponential EWMA chart over the shift size range Ω = [0.3, 0.5] are (0.204, 0.406) and 10.3, respectively. Using the optimal parameter combinations (λ -, h -), the out-of-control M RL value of the exponential EWMA chart is 10 when c = 0.4.

From Tables 4 and7, it is noticed that the performance of the lower-sided exponential EWMA chart designed based on EM RL is close to the one designed based on M RL, if the specified mean shift c is in the mean shift range Ω. For example, when M RL 0 = 200 and c opt = 0.6 in Table 4, the M RL 1 values of the lower-sided exponential EWMA chart are 28, 19, 14 for shift sizes of 0.7, 0.6, 0.5, respectively. While in Table 7, the In order to investigate the effect of parameter estimation on the RL distribution of the chart, Tables 8 and9 show the effect of the amount of Phase I data (i.e., m) on the Phase II performance of the optimal lower-sided charts. Table 8 presents the in-control M RL values for the lower-sided exponential EWMA chart when m ∈ {5, 10, 20, 30, 50, 100, +∞}.

M
Note that the parameter combinations of (λ -, h -) are the optimal settings for the specific M RL 0 ∈ {100, 300} and c opt ∈ {0.2, 0, 4, 0.6, 0.8} when the process parameter is known (m = +∞). It is found that there is a difference between the in-control M RL values corresponding to the known parameter cases and the ones corresponding to the estimated parameter cases, especially for small values of λ -and m. For example, when M RL 0 = 100 and λ -= 0.487, the in-control M RL of the chart are 122 and 100 corresponding to m = 5 and m = +∞, respectively. If λ -= 0.06, the in-control M RL value of the chart are 174 and 100 corresponding to m = 5 and m = +∞, respectively.

Moreover, when M RL 0 = 300, for the same parameter combinations (λ -, h -) presented in Table 8, the out-of-control 5th percentile of RL distribution, M RL and 95th percentile of RL distribution calculated using Eqs. (A.11) and (A.12) in the Appendix are also presented in Table 9. It is shown that the difference of M RL 1 between the known and estimated parameters cases becomes large when m decreases. For example, when c opt = 0.8, the M RL 1 value equals to 78 when m = 5. If m increases from 5 up to 100, the M RL 1 value decreases from 78 down to 55, which is more close to the M RL 1 = 54 corresponding to the known parameter case. It is obvious that the M RL 1 value of the chart with estimated parameter tends to be close to their known parameter counterpart, as m increases. Furthermore, it can also be noted that as the c opt decreases, the difference of M RL 1 between the known and estimated parameters cases decreases. This fact implies that the effect of the parameter estimation becomes small when a large shift happens in the process. In addition, the 5th percentile of the RL distribution increases and the 95th percentile of the RL distribution decreases when m increases.

Due to the space limitation, the results for the upper-sided exponential EWMA chart are not presented here. Based on the results presented above, the RL distribution of the chart is obviously affected by the estimated parameter. To obtain a close performance to the chart with known parameter, the number m of Phase I samples should be large enough to minimize the effect of parameter estimation.

(Please Insert Tables 8-9 to Here)

Conclusions and Recommendations

Since the shape of the RL distributions of control charts changes with shift size, this paper studied the properties of the exponential EWMA chart based on M RL and EM RL when the process shift sizes are either known or unknown, respectively. For an increase or a decrease in the process mean, the optimal design procedures of the upper-and lower-sided exponential EWMA charts based on M RL have been presented. In practice, when the shift size is unknown in advance, based on EM RL, the optimal design procedures for both the upper-and lower-sided charts are also presented. Moreover, through the comparisons of the

Q =             p 1,1 p 1,2 • • • p 1,N p 2,1 p 2,2 • • • p 2,N . . . . . . . . . . . . p N,1 p N,2 • • • p N,N             , (A.1)
where p i,j is the probability of a transition from state i to state j in one step. Let q be the initial probabilities associated with the N transient states, i.e., q = (q 1 , q 2 , . . . , q N ) .

As stated in [START_REF] Latouche | Introduction to matrix analytic methods in stochastic modeling[END_REF] or [START_REF] Neuts | Matrix-geometric solutions in stochastic models: an algorithmic approach[END_REF], the RL of the one-sided exponential EWMA chart is a Discrete Phase-type (DPH) random variable of parameters (Q, q), the probability mass function (p.m.f.) f RL ( |Q, q) and the c.d.f. F RL ( |Q, q) of the RL distributions are given as follows:

f RL ( |Q, q) = q Q -1 r, (A.2) F RL ( |Q, q) = 1 -q Q 1. (A.3)
where r =(1 -Q1) with 1 = (1, 1, . . . , 1) . The ARL of the one-sided exponential EWMA chart is

ARL = q (I -Q) -1 1. (A.4)
Following [START_REF] Gan | An optimal design of EWMA control charts based on median run length[END_REF], the 100γth (0 < γ < 1) percentile of the RL distribution, denoted by γ is defined as

F RL ( γ -1|Q, q) γ and F RL ( γ |Q, q) > γ, (A.5)
Moreover, the M RL of the one-sided exponential EWMA chart can be computed by setting γ = 0.5,

F RL (M RL -1|Q, q) 0.5 and F RL (M RL|Q, q) > 0.5. (A.6)
For the upper-sided exponential EWMA chart, the interval between the reflecting boundary B + and the control limit h + is divided into N subintervals, each with width

∆ + = h + -B + N .
When the number N of subintervals is sufficiently large (N = 500 in this paper), this finite approach provides an effective method that allows the RL properties to be accurately evaluated. The statistic

Q + t is said to be in transient state j at time t if H + j -∆ + 2 ≤ Q + t < H + j + ∆ + 2
, where H + j = B + + (j -1 2 )∆ + (j = 1, 2, . . . , N ) represents the midpoint of the jth sub-interval. The statistic Q + t is in the absorbing state if Q + t falls outside the control limits.

Then, the transition probability p i,j is given as:

p i,j =        Pr(Q + t < H + 1 + ∆ + 2 |Q + t-1 = H + i ) for j = 1 Pr(H + j -∆ + 2 ≤ Q + t < H + j + ∆ + 2 |Q + t-1 = H + i ) for j = 2, . . . , N (A.7)
After some algebraic calculations, we can obtain

p i,j =                F B + +∆ + -(1-λ + )H + i cλ + for j = 1 F H + j + ∆ + 2 -(1-λ + )H + i cλ + -F H + j -∆ + 2 -(1-λ + )H + i cλ + for j = 2, . . . , N (A.8)
where F (•) is c.d.f. of the exponential distribution (see Eq. ( 2)).

Similarly, the transition probability p i,j of the lower-sided exponential EWMA chart can be obtained:

p i,j =                1 -F B --∆ --(1-λ -)H - i cλ - for j = 1 F H - j + ∆ - 2 -(1-λ -)H - i cλ - -F H - j -∆ - 2 -(1-λ -)H - i cλ - for j = 2, . . . , N (A.9)
where

H - j = B --(j -1 2 )∆ -(j = 1, 2, . . . , N ) represents the midpoint of the jth sub-interval and ∆ -= B --h - N
. By using Eq. (A.6), one can easily compute the M RL of the upper-and lower-sided exponential EWMA charts.

If the process parameter is estimated, let w = θ 0 / θ0 represents the ratio of the actual in-control mean to its estimator. It is known that w is distributed according to an inverse gamma distribution, denoted as IG(m, m), and the p.d.f. of w is,

f W (w) = m m (m -1)! w -m-1 e -m w . (A.10)
Since w is an independent random variable, the unconditional c.d.f. F RL ( | Q, q) and the M RL of the RL distribution can be obtained as follows:

F RL ( | Q, q) = +∞ 0 (1 -q Q 1)f W (w)dw. (A.11) F RL (M RL -1| Q, q) 0.5 and F RL (M RL| Q, q) > 0.5. (A.12)
where the conditional transition probability matrix Q = [p i,j ] N ×N is computed as follows:

• For the upper-sided exponential EWMA chart with estimated parameters pi,j =

               F B + +∆ + -(1-λ + )H + i wcλ + for j = 1 F H + j + ∆ + 2 -(1-λ + )H + i wcλ + -F H + j -∆ + 2 -(1-λ + )H + i wcλ +
for j = 2, . . . , N (A.13)

• For the lower-sided exponential EWMA chart with estimated parameters pi,j =

               1 -F B --∆ --(1-λ -)H - i wcλ - for j = 1 F H - j + ∆ - 2 -(1-λ -)H - i wcλ - -F H - j -∆ - 2 -(1-λ -)H - i wcλ - for j = 2, . . . , N (

3. 1

 1 Choices of the reflecting boundaries B + and B -for the M RLbased one-sided exponential EWMA charts In order to design the M RL-based one-sided exponential EWMA control charts, the reflecting boundary B + (B -) of the upper-sided (lower-sided) EWMA control chart should be selected first. For the ARL-based exponential EWMA charts, by doing many numerical analyses, values of the reflecting boundaries have been suggested to be B + = 0.5 and B -= 2.0 by[START_REF] Gan | Designs of one-and two-sided exponential EWMA charts[END_REF]. These values are also used in other researches conducted on the ARL-based exponential EWMA chart (see[START_REF] Ozsan | Properties of the exponential EWMA chart with parameter estimation[END_REF]). Now, we are wondering if these suggested values can also be used for the M RL-based one-sided exponential EWMA charts? Considering this issue, many numerical simulations have been performed

  mechanical parts produced by multi-purpose Computer Numerical Control (CNC) machines (see[START_REF] Celano | Performance of t control charts in short runs with unknown shift sizes[END_REF]), thus a uniform distribution for c is recommended. Whereas, other p.d.f. of the mean shift c, like for instance, truncated normal distribution (see[START_REF] Ryu | Optimal design of a CUSUM chart for a mean shift of unknown size[END_REF]), Rayleigh distribution (see Wu et al. (2004)), to name a few, can be considered according to the historic data of the out-of-control cases or the practitioner's knowledge about the process. Following[START_REF] Castagliola | Monitoring the coefficient of variation using EWMA charts[END_REF]), Tang et al. (2019), Teoh et al. (2017), a uniform distribution is assumed in this article.

  300, the M RL 1 values increases from 54 up to 77 when λ -increases from 0.052 up to 0.196 when the mean shift c = 0.8, which implies that the control chart with λ -= 0.052 detects the mean shift (c = 0.8) more quickly than the chart with λ -= 0.196. While, if the mean shift c decreases down to 0.2, the M RL 1 values are 10 and 7 for λ -= 0.052 and λ -= 0.196, respectively, which implies that the control chart with λ -= 0.196 detects the mean shift (c = 0.2) more quickly than the chart with λ -= 0.052.By comparing the performances of the lower-sided exponential EWMA chart with the lower-sided exponential Shewhart chart, it is also noted that the EWMA-type chart is superior to the Shewhart-type chart in most cases. The lower-sided exponential Shewhart chart only performs slightly better than the lower-sided exponential EWMA chart when the mean shift is large (c < 0.05). For instance, when M RL 0 = 200, the M RL 1 values of the lowersided exponential EWMA chart with optimal parameter combinations designed for c opt = 0.4 are all smaller than the ones of the lower-sided exponential Shewhart chart when c ≥ 0.03.When c = 0.01, the M RL 1 = 2 for the lower-sided exponential Shewhart chart is a little smaller than the M RL 1 = 4 for the lower-sided exponential EWMA chart.Moreover, the M RL profiles of the upper-sided exponential EWMA and Shewhart charts are also displayed in Table5for M RL 0 ∈ {100, 200, 300, 500}. In order to observe the effect of the shift size comprehensively, we chose several representative values from Gan (1998), i.e., c opt ∈ {1.04, 1.4, 2, 5}. The optimal parameter combinations (λ + , h + ) of the upper-sided exponential EWMA chart and the upper control limit (U CL) of the uppersided exponential Shewhart chart are displayed in rows 4 and 5. Similar to the lower-sided exponential EWMA chart, the upper-sided exponential EWMA chart with a small value of λ + gives better performance for smaller shift sizes compared with the upper-sided exponential Shewhart chart. From Table5, it is also clear that the upper-sided exponential EWMA chart consistently performs better than the upper-sided exponential Shewhart chart, for small to moderate shift sizes (1 < c < 3). The upper-sided exponential Shewhart chart only outperforms the upper-sided exponential EWMA chart when the mean shift is large (c > 3).

c

  opt . The M RL and ARL values corresponding to c = c opt are also bolded in the table. It can be concluded from the table that for the specific c = c opt in column 1, the M RL of the M RL-based chart performs a little better than the ARL-based chart and the ARL of the ARL-based chart performs a little better than the M RL-based chart. For example, when c opt = 0.6, the M RL = 23 for the M RL-based chart performs a little better than the M RL = 24 for the ARL-based chart when c = 0.6. At the same time, the ARL = 26.6 for the ARL-based chart performs a little better than the ARL = 27.1 for the M RL-based chart when c = 0.6 (see the bolded values in columns 4 and 5).

  RL 1 values are 28, 19, 14, respectively, when using the parameter combinations (λ -= 0.087, h -= 0.610) designed based on Ω = [0.5, 0.7]. This implies that the optimal parameter combinations obtained based on EM RL can be considered as a robust alternative to the parameter combinations obtained based on M RL if the specified mean shift c of the optimal M RL-based chart is in the range Ω of the optimal EM RL-based chart. Furthermore, Table7 alsoshows the superiority of the lower-sided exponential EWMA chart over the lowersided exponential Shewhart chart. For example, when M RL 0 = 200, the M RL 1 values of the lower-sided exponential EWMA chart are smaller than the ones of the lower-sided exponential Shewhart chart when c ∈ [0.05, 0.9]. The lower-sided exponential Shewhart chart is more sensitive than the EWMA type chart only when c < 0.05. Similar findings are also obtained for the upper-sided exponential EWMA chart. Due to the space limitation, these results are not presented here. (Please Insert Table 7 to Here) 5 The Effect of Parameter Estimation on the One-Sided Exponential EWMA Charts Parameter estimation is an important research issue in control charts because the selection of Phase I samples influences the performance of control charts via the degree of accuracy in the estimator. Suppose that {Z 1 , Z 2 , . . . , Z m } is a sample from an in-control Phase I sample, then a common estimator for θ 0 is the maximum likelihood estimator θ0 = 1

  M RL-based one-sided exponential EWMA chart and the Shewhart chart, it is noted that the EWMA type chart generally performs better than the Shewhart type chart, except for large shift sizes in the process. As this work is based on the univariate variable, future research works may be extended to the chart for monitoring the multivariate exponential distributed processes. interval into a finite set of sub-intervals and by considering the midpoint value within each sub-interval to approximate the value of the corresponding statistic. Assuming a discretetime Markov chain with N transient states, the transition probability matrix Q is given as follows:

Table 1 :

 1 A.14) You, H. W.,Khoo, M. B. C., Castagliola, P., and Qu, L. (2016). Optimal exponentially weighted moving average X charts with estimated parameters based on median run length The percentiles of the RL distribution and ARL for the lower-sided exponential EWMA chart when ARL 0 = 500 (λ -,h -)

	and expected median run length. International Journal of Production Research, 54(17):1-
	22.
	Zhang, C. W., Xie, M., and Goh, T. N. (2005). Economic design of exponential charts for time
	between events monitoring. International Journal of Production Research, 43(23):5019-
	5032.
	Zhang, C. W., Xie, M., and Goh, T. N. (2006). Design of exponential control charts using a
	sequential sampling scheme. IIE Transactions, 38(12):1105-1116.

Table 2 :

 2 The percentiles of the RL distribution and ARL for the upper-sided exponential EWMA chart when ARL 0 = 500 (λ + ,h + )

	c	ARL 5th 10th 20th 30th 40th 50th 60th 70th 80th 90th
	(0.05,1.4167) 1.00 500.0 31	58	116 182 258 348 459 601 801 1144
	1.04 306.7 23	39	74	114 160 215 282 368 489 696
	1.40 35.5	7	9	14	18	23	28	34	42	53	72
	5.00	3.4	1	1	2	2	3	3	4	4	5	6
	10.00 2.1	1	1	1	1	2	2	2	2	3	4
	(0.10,1.7141) 1.00 500.0 29	56	114 181 257 348 458 601 803 1147
	1.04 329.5 21	39	77	120 170 230 302 396 528 753
	1.40 39.6	5	8	13	18	23	30	37	47	61	84
	5.00	3.2	1	1	2	2	2	3	3	4	5	6
	10.00 1.9	1	1	1	1	1	2	2	2	3	3
	(0.20,2.2378) 1.00 500.0 27	54	113 180 256 347 458 602 804 1149
	1.04 353.1 20	39	81	127 181 245 324 425 567 810
	1.40 48.4	5	8	13	19	26	35	45	58	76	108
	5.00	3.1	1	1	1	2	2	3	3	4	4	6
	10.00 1.8	1	1	1	1	1	2	2	2	3	3
	(0.40,3.2213) 1.00 500.0 27	54	112 179 256 347 458 602 804 1150
	1.04 374.7 20	40	84	134 192 260 343 451 602 861
	1.40 63.0	4	8	15	23	33	44	58	76	101 143
	5.00	3.1	1	1	1	2	2	2	3	4	4	6
	10.00 1.8	1	1	1	1	1	1	2	2	2	3
	(0.60,4.2015) 1.00 500.0 26	53	112 179 256 347 458 602 804 1150
	1.04 384.8 20	41	86	138 197 267 353 463 619 885
	1.40 73.2	5	8	17	27	38	51	67	88	117 167
	5.00	3.2	1	1	1	2	2	2	3	4	5	6
	10.00 1.8	1	1	1	1	1	1	2	2	3	3
	(0.80,5.1986) 1.00 500.0 26	53	112 179 256 347 458 602 804 1151
	1.04 390.5 21	42	88	140 200 271 358 470 628 898
	1.40 80.1	5	9	18	29	41	56	73	96	129 184
	5.00	3.3	1	1	1	2	2	2	3	4	5	7
	10.00 1.8	1	1	1	1	1	1	2	2	3	3
	(1.00,6.2147) 1.00 500.0 26	53	112 179 256 347 458 602 804 1151
	1.04 393.7 21	42	88	141 201 273 361 474 633 906
	1.40 84.7	5	9	19	31	44	59	78	102 136 194
	5.00	3.5	1	1	1	2	2	3	3	4	5	7
	10.00 1.9	1	1	1	1	1	1	2	2	3	3

Table 3 :

 3 The M RL of the lower-sided exponential EWMA chart when M RL 0 = 100

					λ -= 0.1					λ -= 0.3
		B -	1.0	1.5	2.0	2.3	2.5	5.0	1.0	1.5	2.0	2.3	2.5	5.0
	c	h -0.591 0.629 0.631 0.631 0.631 0.630	0.323 0.344 0.348 0.348 0.348 0.348
	1.0		100	100	100	100	100	100	100	100	100	100	100	100
	0.9		59	53	52	52	52	52	67	63	62	62	62	62
	0.8		37	31	31	31	31	30	44	40	39	39	39	39
	0.7		24	20	20	20	20	20	29	26	25	25	25	25
	0.6		17	15	14	14	14	14	19	17	16	16	16	16
	0.5		13	11	11	11	11	11	13	12	11	11	11	11
	0.4		10	9	9	9	9	9	9	8	8	8	8	8
	0.3		8	7	7	7	7	7	7	6	6	6	6	6
	0.2		7	6	6	6	6	6	5	5	5	5	5	5
	0.1		6	6	5	5	5	5	4	4	4	4	4	4
					λ -= 0.5					λ -= 0.7
		B -	1.0	1.5	2.0	2.3	2.5	5.0	1.0	1.5	2.0	2.3	2.5	5.0
	c	h -0.190 0.201 0.204 0.205 0.205 0.205	0.106 0.110 0.112 0.113 0.113 0.113
	1.0		100	100	100	100	100	100	100	100	100	100	100	100
	0.9		72	70	69	69	69	68	77	75	75	74	74	74
	0.8		51	48	47	47	46	46	57	55	54	54	54	54
	0.7		35	32	31	31	31	31	42	40	39	38	38	38
	0.6		24	22	21	21	21	21	29	28	27	27	27	26
	0.5		16	14	14	14	14	14	20	19	18	18	18	18
	0.4		11	10	9	9	9	9	13	12	12	12	12	12
	0.3		7	7	6	6	6	6	9	8	8	8	8	8
	0.2		5	5	5	5	5	5	5	5	5	5	5	5
	0.1		4	3	3	3	3	3	3	3	3	3	3	3

Table 4 :

 4 The M RL of the lower-sided exponential EWMA chart and the lower-sided expo-

		nential Shewhart chart when M RL 0 ∈ {100, 200, 300, 500}
					M RL 0 = 100				M RL 0 = 200
				EWMA		Shewhart		EWMA		Shewhart
		c opt	0.8	0.6	0.4	0.2		0.8	0.6	0.4	0.2
		λ -0.060 0.101 0.258 0.487	LCL =	0.051 0.100 0.220 0.435	LCL =
	c	h -0.736 0.629 0.390 0.211	0.0069	0.716 0.579 0.386 0.207	0.0035
	1.00		100	100	100	100	100	200	200	200	200	200
	0.90		49	52	60	69	90	82	91	108	127	180
	0.80		29	31	37	46	80	43	48	60	80	160
	0.70		20	20	24	31	70	28	28	35	49	140
	0.60		14	14	16	21	60	20	19	21	30	120
	0.50		11	11	11	14	50	15	14	14	19	100
	0.40		9	9	8	9	40	12	11	10	12	80
	0.30		8	7	6	6	30	10	9	8	8	60
	0.20		7	6	5	5	20	9	7	6	5	40
	0.10		6	5	4	3	10	8	6	5	4	20
	0.05		6	5	4	3	5	7	6	5	4	10
	0.03		6	5	4	3	3	7	6	5	3	6
	0.01		5	5	4	3	1	7	6	4	3	2
					M RL 0 = 300				M RL 0 = 500
				EWMA		Shewhart		EWMA		Shewhart
		c opt	0.8	0.6	0.4	0.2		0.8	0.6	0.4	0.2
		λ -0.052 0.089 0.196 0.407	LCL =	0.051 0.079 0.180 0.375	LCL =
	c	h -0.688 0.580 0.392 0.206	0.0023	0.664 0.577 0.387 0.206	0.0014
	1.00		300	300	300	300	300	500	500	500	500	500
	0.90		111	124	150	183	270	161	179	229	288	450
	0.80		54	59	77	109	240	70	78	108	162	400
	0.70		33	34	42	64	210	40	41	55	90	350
	0.60		23	22	25	38	180	27	26	30	49	300
	0.50		17	16	16	22	150	20	18	19	28	250
	0.40		14	12	11	14	120	16	14	13	16	200
	0.30		11	10	9	9	90	13	11	10	10	150
	0.20		10	8	7	6	60	11	10	8	7	100
	0.10		8	7	6	4	30	9	8	6	5	50
	0.05		8	7	5	4	15	9	8	6	4	25
	0.03		8	7	5	4	9	9	7	6	4	15
	0.01		8	6	5	4	3	8	7	5	4	5

Table 6 :

 6 The M RL and ARL of the lower-sided exponential EWMA chart when M RL 0 = 350 and ARL 0 = 500

				M RL-based with M RL 0 = 350
	c opt		0.8		0.6		0.4		0.2
	λ -	0.0510	0.0830	0.1480	0.3950
	h -	0.6826	0.5857	0.4535	0.2074
	c	M RL ARL	M RL ARL	M RL ARL	M RL ARL
	1.0	350	499.2	350	500.4	350	501.6	350	502.6
	0.9	124	170.5	137	191.9	159	225.0	209	299.8
	0.8	58	75.0	64	85.4	77	106.6	123	175.0
	0.7	34	41.2	36	44.7	41	54.7	71	100.2
	0.6	24	26.6	23	27.1	25	30.9	41	56.6
	0.5	18	19.2	17	18.5	16	19.3	24	32.0
	0.4	14	14.8	13	13.7	12	13.3	14	18.3
	0.3	12	12.0	10	10.8	9	9.8	9	10.9
	0.2	10	10.1	9	8.9	7	7.7	6	6.9
	0.1	9	8.8	8	7.6	6	6.4	5	4.8
				ARL-based with ARL 0 = 500.0
	c opt		0.8		0.6		0.4		0.2
	λ -	0.0501	0.0530	0.1446	0.3519
	h -	0.6859	0.6755	0.4592	0.2363
	c	M RL ARL	M RL ARL	M RL ARL	M RL ARL
	1.0	350	499.6	350	499.5	349	499.9	348	500.0
	0.9	123	169.9	125	172.0	157	222.8	202	288.6
	0.8	58	74.8	59	75.7	76	105.2	115	163.5
	0.7	34	41.1	34	41.4	41	53.9	65	91.6
	0.6	24	26.6	24	26.6	24	30.5	37	51.1
	0.5	18	19.2	18	19.1	16	19.2	22	28.9
	0.4	14	14.9	14	14.7	12	13.2	13	16.9
	0.3	12	12.1	12	11.9	9	9.8	9	10.4
	0.2	10	10.2	10	10.0	7	7.7	6	6.9
	0.1	9	8.8	9	8.7	6	6.4	5	4.9

Table 7 :

 7 The M RL and EM RL of the lower-sided exponential EWMA chart and the lowersided exponential Shewhart chart when M RL 0 ∈ {100, 200, 300, 500}

					M RL 0 = 100				M RL 0 = 200
				EWMA		Shewhart		EWMA		Shewhart
		Ω	[0.7,0.9] [0.5,0.7] [0.3,0.5] [0.1,0.3]		[0.7,0.9] [0.5,0.7] [0.3,0.5] [0.1,0.3]
		λ -	0.050	0.093	0.244	0.510	LCL =	0.050	0.087	0.204	0.414	LCL =
	c	h -	0.768	0.647	0.405	0.198	0.0069	0.720	0.610	0.406	0.220	0.0035
		EM RL 1	30.3	14.7	8.3	4.7	-	46.7	19.6	10.3	5.6	-
	0.90		48	52	60	69	90	82	89	106	126	180
	0.80		29	30	37	47	80	43	46	58	78	160
	0.70		19	20	23	32	70	27	28	34	48	140
	0.60		15	14	16	21	60	19	19	21	29	120
	0.50		12	11	11	14	50	15	14	14	18	100
	0.40		10	9	8	9	40	12	11	10	12	80
	0.30		8	7	6	6	30	10	9	8	8	60
	0.20		7	6	5	5	20	9	8	6	5	40
	0.10		6	6	4	3	10	8	7	5	4	20
	0.05		6	5	4	3	5	7	6	5	4	10
	0.03		6	5	4	3	3	7	6	5	4	6
	0.01		6	5	4	3	1	7	6	5	3	2
					M RL 0 = 300				M RL 0 = 500
				EWMA		Shewhart		EWMA		Shewhart
		Ω	[0.7,0.9] [0.5,0.7] [0.3,0.5] [0.1,0.3]		[0.7,0.9] [0.5,0.7] [0.3,0.5] [0.1,0.3]
		λ -	0.050	0.085	0.179	0.370	LCL =	0.050	0.070	0.157	0.330	LCL =
	c	h -	0.695	0.590	0.415	0.230	0.0023	0.667	0.602	0.421	0.237	0.0014
		EM RL 1	59.5	22.8	11.6	6.2	-	79.5	27.1	13.3	6.9	-
	0.90		111	122	146	178	270	160	173	219	276	450
	0.80		54	59	74	104	240	70	75	101	150	400
	0.70		33	33	41	60	210	40	40	51	81	350
	0.60		23	22	24	35	180	27	26	29	44	300
	0.50		17	16	16	21	150	20	19	18	25	250
	0.40		14	12	11	13	120	16	14	13	15	200
	0.30		11	10	9	9	90	13	12	10	10	150
	0.20		10	8	7	6	60	11	10	8	7	100
	0.10		9	7	6	5	30	9	8	6	5	50
	0.05		8	7	5	4	15	9	8	6	4	25
	0.03		8	7	5	4	9	9	8	6	4	15
	0.01		8	6	5	4	3	8	8	6	4	5

Table 8 :

 8 The in-control M RL values for the lower-sided exponential EWMA chart with estimated parameters when m ∈ {5, 10, 20, 30, 50, 100, +∞} and M RL 0 ∈ {100, 300}

				M RL 0 = 100			M RL 0 = 300
		c opt	0.8	0.6	0.4	0.2	0.8	0.6	0.4	0.2
		λ -0.060 0.101 0.258 0.487	0.052 0.089 0.196 0.407
	m	h -0.736 0.629 0.390 0.211	0.688 0.580 0.392 0.206
	5		174	159	135	122	650	554	453	387
	10		127	123	114	109	421	392	356	332
	20		111	109	105	104	346	334	319	311
	30		106	105	103	102	326	318	310	305
	50		103	103	101	101	312	307	303	302
	100		101	101	101	100	303	301	300	300
	+∞		100	100	100	100	300	300	300	300

Table 9 :

 9 out-of-control 5th, 50th and 95th percentiles of RL distribution for the lower-sided exponential EWMA chart with estimated parameters when m ∈ {5, 10, 20, 30, 50, 100, +∞} and M RL 0 = 300

	c opt		0.8			0.6			0.4			0.2
	λ -		0.052			0.089			0.196			0.407
	h -		0.688			0.580			0.392			0.206
	m	5th M RL 95th	5th M RL 95th	5th M RL 95th	5th M RL 95th
	5	13	78	> 10 4	9	26	87787	6	12	392	4	6	32
	10	14	64	> 10 4	10	24	994	6	12	74	4	6	18
	20	15	58	5566	10	23	181	6	12	40	4	6	14
	30	16	57	1524	10	23	115	6	11	34	4	6	13
	50	16	56	595	10	22	83	6	11	30	4	6	13
	100	17	55	312	10	22	67	6	11	27	4	6	12
	+∞	17	54	174	11	22	55	6	11	25	4	6	12

Acknowledgments

The authors are grateful to the anonymous reviewers for their useful and detailed comments which contributed to the improvement of this article. This work are supported by National Natural Science Foundation of China No. 71802110, Natural Science Foundation of JiangSu Province No. BK20170894, Humanity and Social Science Foundation of Ministry of Education of China No. 19YJA630061, Key Research Base of Philosophy and Social Sciences in Jiangsu-Information Industry Integration Innovation and Emergency Management Research Center.

Appendix

The RL properties of one-sided exponential EWMA charts can be obtained using the Markov chain approach. The states of the Markov chain are defined by partitioning the control Sparks, R. S. (2000). CUSUM charts for signalling varying location shifts. Journal of Quality Technology, 32(2):157-171. Tang, A. A., Castagliola, P., Sun, J. S., and Hu, X. L. (2019). Optimal design of the adaptive EWMA chart for the mean based on median run length and expected median run length.

Quality


H Borouchaki
email: fhouman.borouchaki@univ-troyes.fr

P Lafon

P Laug

P L George

MINIMAL VARIATIONAL SURFACES AND QUALITY MESHES

Keywords: minimal variational surfaces, quality meshes, equilateral meshes, regular meshes, adaptative meshing, nite elements, convex quadratic optimization

Many p h ysical phenomena in science and engineering can be modelled by partial di erential equations (PDEs) and solved using the nite element method (FEM). Such a method uses as computational spatial support a mesh of the domain where the equations are formulated. The mesh quality i s a k ey-point for the accuracy of the numerical solution. This paper describes a methodology to construct a quality mesh of the domain from a given discretization of its boundary. W e show that the size map related to such a mesh constitutes a minimal variational surface supported by a given contour. This surface can be constructed, from its boundary, using the nite element method or by the resolution of a simple discrete optimization problem. The quality mesh of the domain is then a mesh conforming to the size map given by this surface. A numerical example is given to demonstrate the method.

INTRODUCTION

The resolution of a physical problem formulated in terms of PDEs using the nite element method (FEM) is based on a spatial discretization, or mesh, of the computational domain. The convergence of such a method, as well as the quality of the resulting solution, is strongly related to the element shape quality (the ideal shape is that of an equilateral or regular element). Many methods (cf. 2]) can be used to construct a mesh of the domain, in general from the data of its boundary. Historically, the rst methods have been based on the \local" regularity of the mesh. One can realize that such methods can fail in the case of a signi cant size variation in the boundary discretization of the domain. The present methods establish a continuous size eld in the domain from the size given by the speci ed boundary discretization, and propose various algorithms to generate a mesh that conforms to this size eld.

This paper discusses what can be an optimal size eld which is related to a quality mesh of the domain (i.e. a mesh which is as regular as possible). Such a s i z e eld presents a \minimal size variation" among all the size elds verifying the speci ed size on the boundary of the domain. The size variation can be quanti ed in two di erent w ays which c haracterize the gradation of the related meshes. Let h be the size function in the domain, the rst approach consists in considering the \usual" gradient rh of h. As for the second approach, the variation is measured by the \logarithmic" gradient r logh of h. Let us denote by (h) the variation of h irrespective of the type of variation. A surface with minimal variation de ned in a domain supported by a given smooth contour i s a s u rface in which the sum of the square of the modulus of the variation on is minimal. Mathematically, i f h corresponds to such a surface, it is characterized by: h(;) = and J (h) =min u2H 1 ( ) J (u)

with J (u) = Z jj (u)jj 2 dX
where ; is the boundary of and

H 1 ( ) = fu 2 L 2 ( ) 8i i (u) 2 L 2 ( )g i (u) denot- ing the i-th component o f (u).
T o solve this problem, we propose a collocation method which allows us to construct approximatively the minimal variational surfaces. In Section 1, the concept of a regular mesh is recalled. In Section 2, the minimal variational surfaces are introduced and their relation with regular meshes is indicated. In Section 3, an algorithm which creates approximatively minimal variational surfaces is given. In Section 4, the method is applied to build regular meshes of a domain whose boundary discretization in prescribed.

REGULAR MESH

Let be a closed bounded domain in R 2 de ned by its boundary. A quality triangular mesh, equilateral or regular, o f is a mesh whose elements are equilateral or regular. The existence of such a mesh is not guaranteed in general. Indeed, it depends, to some degree, on the domain boundary discretization. Therefore, we will call a regular triangular mesh the \best" triangular mesh which can be completed. As the issue of constructing a regular mesh for an arbitrary domain is still an open problem, there are various methods which allow the construction of \almost"open regular meshes.

As an illustration, let us consider the domain de ned by a circular ring whose outer boundary is a circle with radius R = 1 centered at the origin, and the inner boundary is a circle with radius r = ( 1 ; p 3 80 ) [START_REF] Borouchaki | Generation de maillages de qualite, la marche unite, traite MIM, methodes numeriques, maillage[END_REF] also centered at the origin. First, the outer (resp. inner) boundary is discretized with a constant step which equals 2 R 80 (resp. 2 r 80 ), so that the dis- cretization of each boundary has 80 segments. For this discretization, let us notice that the value r of the inner radius is calculated in order to identify the equilateral mesh of the domain (in the sense of the best possible mesh). This mesh, illustrated by F i g u r e 1 , is actually made up of 12 layers of nearly equilateral triangles. Then, a slightly di erent discretization of the inner boundary is considered, taking a step which equals 2 r 80 + 0 :007. Figure 2 shows an almost equilateral mesh of the domain corresponding to this new discretization. As shown in this gure, non regular triangles appear around the fth and sixth layers.

In a classical context, two main types of discretization of the domain boundary can be considered. The rst one involves uniform discretizations with a constant stepsize. The advantage of this type of discretization is its possibility, in general, to construct equilateral meshes. However, it cannot guarantee a good representation of the domain boundary for a given stepsize. Figures 3 and4 show respectively two uniform discretizations of the contour of a wheel cross-section (French railway company), each circle representing the stepsize along any direction. Given a uniform discretization of the domain boundary, a regular mesh of the domain can only be a mesh whose element size equals to the boundary stepsize. So, the desired size of the elements is a priori known at each v ertex of the regular mesh of the domain. Figures 5 and 6 illustrate the two uniform (with constant stepsize), almost equilateral, meshes of the wheel corresponding to the two g i v en discretizations of its contour. The desired size (here constant) along all the directions at each v ertex of these meshes are drawn on Figures 7 and8. We can notice that, at each v ertex of these meshes, the circle representing the desired size of the elements almost passes through the adjacent v ertices. This is a characteristic, among others, of almost equilateral meshes. Let us also notice that the domain geometry is better represented by the second mesh, conceding that the number of elements is greater.

Let us consider the case where the domain boundary is made up of several connected components, discretized with di erent stepsizes (this is a frequent case for domains involved in uid dynamics problems). A regular mesh of such a domain is a mesh whose element s i z e in the neighborhood of each boundary is close to the stepsize of the boundary discretization. As for the element size at other points of the domain, it must remain locally almost constant. Figure 9 shows a domain whose boundary contains ve connected components, one representing the in nite, and the four others different parts of an aircraft wing. The in nite component is discretized with a relatively large stepsize, and the others with much smaller and di erent stepsizes. Figures 9 and11 show almost equilateral mesh of the domain. Figures 10 and12 show respective enlargements of these meshes around the wing body, t h us illustrating the boundary discretizations and the corresponding meshes.

The second type of discretization is called \geometric", for it is adapted to the geometry of the domain boundary. In that case, it can be shown that the step-Figure 9. Almost equilateral mesh of a domain in computational uid dynamics (CFD). size must be proportional to the radius of curvature.

The drawback o f t h i s t ype of discretization is that it can give rise to large variations of the stepsize. To reduce this phenomenon, di erent methods of stepsize smoothing can be used.

Figures 13 and14 show, respectively, t wo geometric discretizations of the wheel contour. The rst one is de ned so that any edge of the discretization is suciently close to the contour (the stepsizes are locally proportional to the radii of curvature). This condition guarantees, among others, that the angle between each edge and the tangents at its extremities remains smaller than a given threshold value (here 4 ). Regarding the second discretization, the stepsizes of the rst one are smoothed so that the ratio between the lengthes of two consecutive edges is bounded. As can be seen on the gures, the circles representing the desired sizes violate in the rst case, and respect (more or less) in the second case, the discretization obtained. Figure 16. Controlled g e ometric mesh. Figures 15 and 16 show t wo almost equilateral meshes corresponding to the two previous geometric discretizations (in the rst case, the stepsizes are locally proportional to the radii of curvature, while in the second case the stepsizes are locally averaged, starting from the rst discretization). As can be noticed on these gures, the second mesh is more regular. This observation can be con rmed by examining Figures 17 and18, which s h o w the sizes at the vertices of the meshes.

As an indication, Figures 19 and20 show the surfaces de ned by the Cartesian functions of the sizes associated with these two meshes.

Remark. Apart from the discretization of the domain boundary, the geometric characteristics can have an e ect on the generation of the equilateral mesh of the domain. These may b e o f d i e r e n t natures, amongst which are angular corners and very close boundary segments (or thin domains).

In the case of a geometric discretization (generally non uniform) of the domain boundary, it is not obvious to a priori establish the desired size of the elements at the vertices of the regular mesh of the domain. Let us recall that a regular mesh is a mesh whose element size variation is locally almost null. Accordingly, t h e idea consists in nding, among all the continuous size functions which conform to the prescribed sizes on the boundary, the least varying one. Therefore, this function represents a kind of \minimal variational surface".

The size variation can be quanti ed in two di erent ways which c haracterize the gradation of the related meshes. Let h be the size function in the domain, the rst approach consists in considering the \usual" gradient rh of h. As for the second approach, the variation is measured by the \logarithmic" gradient r logh of h. In the next sections, we de ne the minimal variational surfaces based on these two t ypes of variations and we g i v e a simple algorithm to construct these surfaces.

MINIMAL VARIATIONAL SURFACES

Let us denote by (h) the variation of h irrespective of the type of variation considered. A surface with minimal usual variation de ned in a closed bounded domain of R 2 , supported by a g i v en contour (sufciently smooth), is a surface in which the sum of the square of the modulus of the variation on is minimal. Mathematically, i f h corresponds to such a surface, it is characterized by: h(;) = and J (h) =min u2H 1 ( ) J (u)

with J (u) = Z jj (u)jj 2 dX where ; is the boundary of and H 1 ( ) = fu 2 L 2 ( ) 8i i (u) 2 L 2 ( )g L 2 denoting the space of functions having a summable square on , and i (u) the i-th component o f (u). In the case where (h) = rh and assuming that ; is regular enough, this surface is nothing other than the solution of the non homogeneous Dirichlet problem: h(;) = and h = 0in which c haracterizes harmonic functions. Using a classical result of numerical analysis, it can be proved that h is a solution of the problem:

h(;) = and 8v 2 H 1 0 ( ) Z < rh rv > d X = 0
where H 1 0 ( ) is the set of fonctions in the Sobolev space of order 1, H 1 ( ), h a ving a null value on the boundary ; of domain . Likewise, if (h) = r log h, it can be shown that h is a solution of the problem:

h(;) = and 8v 2 H 1 0 ( ) Z ; 2 h 3 jjrhjj 2 v + 1 h 2 < rh rv > dX = 0 :
These two problems can be solved using the nite element method. Considering piecewise P 1 solutions, in the rst case, the problem reduces to solving a linear system and, in the second case, a non linear system. Although these solutions give the desired exact surfaces, we nd it simpler to use a collocation method, as will be seen in the following, to approximate these minimal variational surfaces.

APPROXIMATE CONSTRUCTION ALGORITHM

Like in the nite element method, we s e a r c h for minimal variational surfaces by using a given mesh of .

Then, the problem consists in determining the elevations at the vertices of this mesh. Let T be a uniform mesh of domain , including N vertices, the M rst vertices being on the boundary. Knowing the contour (which represents the sizes associated with the discretization of the boundary ; of ), we h a ve t h e elevations at the M rst vertices of T . W e m ust now determine the elevations at the other vertices.

Let D (h e) be a discrete measure of the variation of h along the edge e of the mesh, we consider the following optimization problem:

Minimize J (h) with J (h) = 1 2 X e2T 2 D (h e)
where e sweeps all the edges of the mesh. This problem consists in minimizing the size variation along the directions given by the mesh considered. We nd that this simpler formulation is more e cient in practice to approximate the minimal variational surfaces. Let us denote by h i the elevation (or the desired size) at vertex i of the mesh, and by d ij the (\theoretically" constant) size of edge ij of the mesh. Considering the rst ty p e o f v ariation, we h a ve:

D (h ij) = h i ; h j d ij and M+1 i j N h i ; h j d ij 2
and, in this case, we obtain a simple problem of convex quadratic optimization which reduces to solving a linear system whose matrix is de nite positive. It can be shown that this system of equations represents, indeed, the discretization at order 2 of a non homogeneous and non isotropic transport-di usion equation, of the type:

(X) @h @x + (X) @h @y ;a(X) @ 2 h @ 2 x ; b(X ) @ 2 h @x @y ; c(X ) @ 2 h @ 2 y = 0 : with h(;) = (X = ( x y) representing the generic point of the domain). Similarly, analyzing the second ty p e o f v ariation, we can consider:

D (h ij) = 1 d ij h i ; h j h i + h j and J (h) = 1 2 X M+1 i j N 1 d ij h i ; h j h i + h j 2
which corresponds to a more complex optimization problem which can be solved by an approach based on a quasi-Newton method (algorithm \BFGS", cf. 3]).

Analyzing this new discrete criterion would exhibit a non linear partial di erential equation, having the same anisotropic feature as the rst type of variation. Considering a uniform mesh of the domain would partly avoid this shortcoming.

EXAMPLES

In this section, we present the academic example of a circular ring which is centered at the origin and whose outer boundary has a radius equal to 1. We s u p p o s e that the boundary of this domain is discretized conforming to a size function with a sinusoidal variation. This discretization is illustrated in gure 21. Besides, we consider a uniform mesh of the domain (Figure 22), on which w e will construct the minimal surfaces supported by the curve de ned by the size function on the domain boundary. Figures 23 and24 show, respectively, the two t ypes of minimal variational surfaces, usual and logarithmic, obtained by solving the optimization problem introduced in the previous section. These two surfaces can be di erentiated by examining the isovalues corresponding to constant e l e v ations. To solve both problems, the iterative algorithm BFGS is used. In the rst case, the searched function h is initialized on the domain boundary by the speci ed size function and, for each v ertex inside the domain, by h min + ( h max ; h min ) i N : In the latter expression, h min and h max are respectively the smallest and the largest size prescribed on the boundary, i is the vertex number, and N is the number of vertices of the uniform mesh considered (this initialization voluntarily brings a perturbation of rJ(h) in the domain).

The initial values obtained are J (h) = 4.43e04 and jjrJ(h)jj = 4.9e04. The convergence of the method is obtained after 119 iterations in 6 seconds (on a HP9000-J7000-440Mhz) and corresponds to J (h) = 9.12e02 and jjrJ(h)jj = 1.2e-02. In the second case, the searched function h is initialized by the size function obtained in the rst case, and the initial values are J (h) = 1.21e05 and jjrJ(h)jj = 2.9e05. The convergence of the method is obtained after 129 iterations in 8 seconds and corresponds to J (h) = 1.07e05 and jjrJ(h)jj = 1.7.

The meshes on gures 25 and 26 represent respectively the meshes conforming to the size elds given by t h e minimal variational surfaces obtained, and are made up of respectively 1638 and 2460 elements. The average shape quality of these meshes (using a chosen monotonous measure varying from 0 ( at triangle) to 1 (regular triangle)) is in the order of 0.96 and 0.97. These meshes are generated by a method based on a combined frontal-Delaunay approach 1], 4]. Figure 26. Logarithmic quality mesh.

Figure 1 .

 1 Figure 1. Equilateral mesh.

Figure 2 .

 2 Figure 2. Almost equilateral mesh.

Figure 3 .

 3 Figure 3. Uniform discretization.

Figure 4 .

 4 Figure 4. Another uniform discretization.

Figure 5 .

 5 Figure 5. Uniform mesh.Figure6. Another uniform mesh.

Figure 7 .

 7 Figure 7. Uniform sizes.Figure8. Other uniform sizes.

Figure 8 .

 8 Figure 7. Uniform sizes.Figure8. Other uniform sizes.

Figure 10 .

 10 Figure 9. Almost equilateral mesh of a domain in computational uid dynamics (CFD).Figure10. Partial enlargement around the wing body.

Figure 11 .

 11 Figure 11. Another almost equilateral mesh of the domain.Figure12. Partial enlargement around the wing body.

Figure 13 .

 13 Figure 13. Geometric discretization.

Figure 14 .

 14 Figure 14. Controlled g e ometric discretization.

Figure 15 .

 15 Figure 15. Geometric mesh.Figure16. Controlled g e ometric mesh.

Figure 17 .Figure 19 .Figure 20 .

 171920 Figure 17. Geometric sizes.Figure18. Controlled g e ometric sizes.

Figure 21 .

 21 Figure 21. Speci ed b oundary discretization.

Figure 22 .Figure 24 .

 2224 Figure 22. Uniform mesh.

Figure 25 .

 25 Figure 25. Usual quality mesh.Figure26. Logarithmic quality mesh.


Fanjuan Meng
email: fanjuan.meng@utt.fr

Carl Labergere

Pascal Lafon

Metal Forming Parameters Optimization of Metal Forming Process

Keywords: Metal Forming Process, Numerical Simulation, Optimization, The Initial Billet Dimension, Response Surface Method (RSM)

In metal forming process, the initial shape of the billet is usually simple. Even for some complicated parts, the cylinder is always be used. But defining the dimension of the initial cylinder billet is the most important step for product quality control. Reasonable initial billet dimension can not only reduce raw material cost but also improving mechanical properties of the part.

In this paper, a metal forming parameters optimization system was developed. It arrives at the optimum values by using finite element analysis numerical simulation and response surface method (RSM). To develop this optimization system, several numerical tools are required: geometric modelling (CATIA V5™), FE analysis (ABAQUS®), and optimization computation (ModeFRONTIER®). A two-step axisymmetric forging project was studied as an example by using proposed method.

Introduction

Nowadays, in metal forming process the numerical simulation is used to avoid long and expensive experiments and many research works have been made in this aspect [START_REF] Im | A computer-aided-design system for forming processes[END_REF][START_REF] Kobayashi | Metal Forming and the Finite Elements Methods[END_REF]. In this paper, the main objective is to simulate metal forming process coupled with an optimization procedure in order to obtain the optimal parameters. In the Mechanical System and Concurrent Engineering Laboratory (ICD-LASMIS), an extensive work has been developed since ten years. A finite element package has been developed to solve thermo-visco-elasto-plastic problems with ductile damage in large deformation. This package allows realizing automatically numerical simulation of metal forming processes. In order to avoid large mesh distortions and achieve properly the convergence of calculation, meshing and re-meshing procedure are also included [START_REF] Labergere | Improving numerical simulation of metal forming processes using adaptive remeshing technique[END_REF][START_REF] Labergere | Numerical design of extrusion process using finite thermo elastoviscoplasticity with damage. Prediction of chevron shaped cracks[END_REF][START_REF] Saanouni | Application to Metal Forming, Computational Damage Mechanics, Numerical and Computational methods[END_REF]. The whole model of metal forming process optimization system is shown in Figure 1. In order to connect all these operations together automatically in ModeFRONTIER®, the different macro commands are used to realize it. First, the CATIA TM macro command is used to create the geometry of the tools and define the parameters that we want to optimizer. Then, the ABAQUS® macro command is used to input the shape and processing parameters into a software package based on ABAQUS® to simulate the metal forming process. Finally, the advanced optimization algorithm is performed to reach the optimum variable values on the response surface established.

Numerical Simulation of Metal Forming Process

The numerical simulations of metal forming processes by using the finite element package developed in our laboratory are presented in this paragraph.

THERMO-VISCO-ELASTO-PLASTIC BEHAVIOUR MODEL

A thermo-elato-viscoplastic behaviour model is used for the simulation of forging operation imposed different multi-physic coupling (temperature, plastic strain, and plastic strain rate). This advanced constitutive equations proposed a non linear isotropic and kinematic hardenings strongly coupled with ductile isotropic damage. The coupling between the ductile damage and the elasto-viscoplastic constitutive equations is formulated in the framework of the thermodynamics of irreversible processes together with the Continuum Damage Mechanics (CDM) theory [START_REF] Labergere | Numerical design of extrusion process using finite thermo elastoviscoplasticity with damage. Prediction of chevron shaped cracks[END_REF][START_REF] Saanouni | Application to Metal Forming, Computational Damage Mechanics, Numerical and Computational methods[END_REF][START_REF] Lestriez | Numerical prediction of damage in metal forming process including thermal effects[END_REF].

NUMERICAL ASPECT

The fully coupled thermomechanical constitutive equations presented above have been implemented into Abaqus/Explicit FE code using the Vumat user subroutine. The dynamic explicit resolution procedure has been used in order to solve the thermomechanical problem based on both weak forms relative to the equilibrium and thermal equations. In the case of an updated Lagrangian formulation, the high non-linearity and the large plastic deformation can distort the element of the mesh. To avoid this numerical problem and modulate the fracture propagation, a new numerical methodology based on different tools (DIAMESH2D: 2D triangular and quadrangular mesher, error estimates, transfer fields, automatic sequence load or re-meshing) has been developed in ICD/LASMIS in collaboration with UTC/Roberval [START_REF] Labergere | Improving numerical simulation of metal forming processes using adaptive remeshing technique[END_REF][START_REF] Labergere | Numerical design of extrusion process using finite thermo elastoviscoplasticity with damage. Prediction of chevron shaped cracks[END_REF].

Optimization Problems

VARIABLES

There are so many factors associated with the metal forming process. In this paper, we deal with the example of an aluminum forged part which is obtained by two different operations (platting and forging). The platting operation (Figure 2) is realized with a mechanical press, and the upper tool is translated with a constant velocity V tool =50 s mm / . The forging operation (Figure 3) is realized with a hammer, two blows are needed to completely deform the part. The kinematic energy imposed to the upper tool is the same for blow 1 and blow 2, which is equal to 30KJ. We defined that the initial dimension of the cylinder billet R initial and H initial and the final height of the billet after platting process H final are the variables to be optimized which are shown in Figure 2 and the final shape of the part after forging operation is shown in Figure 3.

OPTIMIZATION OBJECTIVE FUNCTIONS

The forging energy is strongly related to the over heating. From an economic point of view, it is important to minimize forging energy [START_REF] Vieilledent | Optimisation des Outils en Forgeage à Chaud par Simulation Éléments Finis et Méthode Inverse Applications à des Problèmes Industriels[END_REF]. Another criterion that should be taken into account during the optimization is the metallurgic quality of forged parts. It's related directly to the uniformity of distribution of plastic deformation [START_REF] Poursina | Optimum pre-form dies in two-stage forging[END_REF]. So in our work, the objective functions are defined as the minimization of forging energy, improving metallurgic quality of forged parts and reducing the maximum forging force during the operations. And it can be summarized as the following functions (1-3):

    ) . : . : min( 2 0 2 1 0 1 2 1 dt dV dt dV F t V t V I obj forging platting             (1)   p avg p i Nb i T i II obj ele V V F      1 min with    ele Nb i p i T i p avg V V 1   (2) max min F F III obj  (3)
Where  is the Cauchy stress matrix;  is the (elastic and plastic) strain rate;

( ) (       i i n n i   (5) max ) ( S i gap S    (6)
In function (4), ) (i D is the damage value of the element i inside the machining part of the piece which is shown in Figure 4 and the maximum value of the damage D max is set to 0.1. If the damage value outside this part is larger than D max , we also think that the forged part is available for the following operations. In function [START_REF] Saanouni | Application to Metal Forming, Computational Damage Mechanics, Numerical and Computational methods[END_REF],  is the angle value of the nearby two contour elements as shown in Figure 5 and the minimum value of the angle min

 is set to 3  . If ) (i  is larger
than min  , we can get the conclusion that there is no wrinkle in the forged part. In function [START_REF] Lestriez | Numerical prediction of damage in metal forming process including thermal effects[END_REF], S is the total value of the gap distance between the contour element nodes of the achieved part and the required part and the maximum value of the total distance S max is set to 6mm. This value is to verify whether the part fills well or not in the tools. As shown in Figure 6, in order to make the value more precise, we make a box at the end of the flash part, because when calculating the gap, this part is useless to our results, so we just accumulate the gap value outside this box.

Optimization System Model

There are several steps for this metal forming optimization system, including input part, output part and the simulation running part. The workflow in the Mode-FRONTIER software is shown in Figure 7.

Response Surface Method

A simulation of the forging process requires a large amount of calculation. In this situation, a meta-model based on the responses surfaces can be built with few simulation processes. A meta-model can be defined as a collection of statistical and mathematical techniques useful for developing, improving, and optimizing processes [START_REF] Tang | Preform tool shape optimization and redesign based on neural network response surface methodology[END_REF]. So to rapidly get the optimum design values, an appropriate optimizer algorithm is used and all evaluations on the objective and constraint functions are evaluated with a specific meta-model.

Applications

To illustrate the methodology described above, an initial billet dimension optimization system based on response surface method is proposed. In our application, it is assumed that the starting shape is cylinder, and we set the value of variable H initial are 110, 115, 121, 128 and 135mm. The volume of the initial billet (V initial ) is fixed which is equal to 552640 3 mm , so R initial can be calculated.

Another variable H final are defined as 70, 75, 80 and 85mm. There are 20 simulations have been down and the Simulation parameters are as follows: material: 2017A, forging temperature: 360°. F is equal to 8.5722E5N.

I obj F 70 75 80 85 110 1,584E+07 1,602E+07 1,869E+07 1,917E+07 115 1,985E+07 1,578E+07 1,441E+07 2,162E+07 121 1,485E+07 1,664E+07 1,655E+07 1,454E+07 128 1,539E+07 1,509E+07 2,212E+07 1,430E+07 135 1,902E+07 1,416E+07 1,649E+07 1,563E+07

Conclusions

In this paper, the advanced simulation technique is applied to a two-operation metal forming process. The optimization module uses response surface method to search the optimal parameters by changing the variables within the specified limits.

In the future, it is desired that we should incorporate more global optimization algorithms with the surface repose method. Because the response surface approach is a method of the Meta-model, it does not have the function of improving the response surface. It will be better that meta-modelling combined to an evolutionary strategy, or to a genetic algorithm, as in the hybrid approach.

Figure 1 .

 1 Figure 1. The metal forming process optimization system.

Figure 2 .Figure 3 .

 23 Figure 2. Platting operation. (a: H=Hinitial=110mm b: H=90mm c: H=Hfinal=70mm)

ele

  Nb is the total number of the elements; i V is the volume of element i; T V is the total volume; p i  is the equivalent plastic strain computed at the unique Gauss point of the element i; p avg  is the average equivalent plastic strain; max F is maximum force value during the platting operation. 2.3 OPTIMIZATION CONSTRAINTS Several constraints should be considered, including damage, wrinkle and achieving the right configuration at the end. And it can be summarized as the following functions (4

Figure 4 .

 4 Figure 4. The constraint of Damage condition.

Figure 5 .

 5 Figure 5. The constraint of Wrinkle condition.

Figure 6 .

 6 Figure 6. The constraint of Filling condition.

Figure 7 .

 7 Figure 7. The optimization system workflow in ModeFRONTIER.

F

  (Newton) After finishing all the simulations, we can get the interesting values that we need, including the objective functions results and the constraints results. The Tables(1-3) show the three objective functions results. Based on all these results, the surface response model are generated which are shown in Figures (8-10). For each figure, the left 3D picture is the response surface model. It represents the relationship between the two variables and the objective functions values. We get it by using classical two degree polynomial model which is easy to implement in solving our problem. The right one shows the contours plots of the variables values and the objective function results.

Figure 8 .

 8 Figure 8. The response surface representing the relationships between the objective function I obj F and the variables.

Figure 9 .

 9 Figure 9. The response surface representing the relationships between the objective function II obj F and the variables.

Figure 10 .

 10 Figure 10. The response surface representing the relationships between the objective function III obj F and the variables.

Figure 8

 8 Figure 8 shows that the distribution of the forging energy is nonlinear. The objective function I obj F reaches the maximum value when H final is equal to 85mm. So during the forging process, the larger of the value of H final , the more energy is needed. From Figure 9 it can be learned that II objF is more influenced by the value of H final and decreasing H final leads to increase the deformation uniformity. When the value of H final is small, the equivalent plastic strain propagates uniformly in all part which may benefit the forging part. But for a large value of H final , the

Table 1

 1 The simulation results of

	I obj F (Mille Joules)

Table 2

 2 The simulation results of

	II obj F

Table 3

 3 The simulation results of

	III
	obj

Acknowledgement

This study was supported by the UTT/LASMIS, Conseil Régional de Champagne Ardenne, FRANCE and CSC, CHINA, which are gratefully acknowledged.