Sheng Di

Member, IEEE Yves Robert

Ieee ; Fellow

Vivien Frédéric

Member, IEEE Franck Cappello

Frédéric Vivien

Toward an Optimal Online Checkpoint Solution under a Two-Level HPC Checkpoint Model

Keywords: High-Performance Computing, Fault Tolerance, Optimization, Multilevel Checkpoint

) it shows that periodic patterns are optimal and determines the best pattern. We evaluate the proposed solution and compare it with the most up-to-date related approaches on an extreme-scale simulation testbed constructed based on a real HPC application execution. Simulation results show that our proposed solution outperforms other optimized solutions and can improve the performance significantly in some cases. Specifically, with the new solution the wall-clock time can be reduced by up to 25.3% over that of other state-of-the-art approaches. Finally, a brute-force comparison with all possible patterns shows that our solution is always within 1% of the best pattern in the experiments.

INTRODUCTION

The execution scale of the largest high-performance computing (HPC) applications steadily increases over time, driven by the demand of solving the problems that were so far out of reach. In the near future, extreme-scale HPC environments with several millions of cores will be the norm for scientific simulations. An 83,000-processor supercomputer, for example, can match only 1% of the human brain [START_REF] Feinberg | An 83,000-processor supercomputer can only match 1% of your brain[END_REF].

In an extreme-scale HPC environment, however, users often experience a number of hardware crashes or transient errors during execution. Indeed, the most powerful existing computers suffer on the order of one failure per day [START_REF] Casanova | On the impact of process replication on executions of large-scale parallel applications with coordinated checkpointing[END_REF]. Without fault-tolerance mechanisms protecting their execution, HPC applications (such as MPI programs) would have to be restarted from scratch after any kind of failure or error, leading to huge resource waste and resulting in extremely low performance.

The main approaches used to protect HPC executions against unexpected interruptions or failures are based on replication, either space-based or time-based. With space-based replication, a number of replicated processes (i.e., replicas) are launched in order to process the same work as the original processes. This model, called the replica execution model [START_REF] Fiala | Detection and correction of silent data corruption for large-scale high-performance computing[END_REF], [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF], [START_REF] Walters | A Scalable Asynchronous Replicationbased Strategy for Fault Tolerant MPI Applications[END_REF], [START_REF] Walters | Replication-Based Fault Tolerance for MPI Applications[END_REF], leads to a huge waste of resources. The reason is that the total number of cores or processes used to process an application will double or triple, depending on the fault tolerance capability demanded by users. With time-based replication, the replication happens only in the case of failures and takes the form of a re-execution. This is also known as the checkpoint/recovery model. It protects the execution by periodically setting checkpoints (i.e., saving the runtime memory or the key variable states into stable storage devices such as disks). In the case of failures, the saved states are read back, and the execution restarts from the snapshot of the application that they provide. Such a model avoids some of the resource waste created by the replica execution model (e.g., a doubling of required computational resources). Nevertheless, it suffers from several overhead sources: time required to write checkpoints to stable storage, time lost in recovery (read-back of checkpoints), and time elapsed between the last checkpoint and the time the failure struck (the corresponding computations must be re-executed). Furthermore, the time required to store the runtime memory state into the parallel file system increases with the application's execution scale [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], [START_REF] Di | Optimization of a multilevel checkpoint model with uncertain execution scales[END_REF], and it is expected to be prohibitive in extreme-scale HPC environments.

Multilevel checkpointing [START_REF] Gomez | Low-overhead diskless checkpoint for hybrid computing systems[END_REF], [START_REF] Gomez | Distributed diskless checkpoint for large scale systems[END_REF], [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF], [START_REF] Moody | Design, modeling, and evaluation of a scalable multi-level checkpointing system[END_REF], [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF], [START_REF] Kulkarni | The design and implementation of a multi-level content-addressable checkpoint file system[END_REF] has recently been proposed to resolve the huge checkpoint/recovery overhead issue. Multilevel checkpointing allows the use of different levels/types 1 of checkpoints during the execution. Obviously, a multilevel checkpoint model is much more flexible than the traditional single-level checkpoint model that uses only one type of storage device or mechanism (typically stable but also slow) to perform checkpoints. However, these different checkpoint levels correspond to different resilience and recovery capabilities and thus suffer from different checkpoint/recovery overheads.

One of the critical issues with the multilevel checkpoint model, then, is how to optimize it: What should be the checkpointing period on each different checkpoint level for a given HPC application? This problem is challenging because of the following factors:

(1) the different levels protect the execution against different types of errors or failures, and each type of failure follows its own distribution; and (2) the different levels cannot be optimized independently because their relative performance impacts each other, so that the whole multilevel checkpointing system must be globally optimized.

In this paper, we formulate the research as a two-level (or two-type) checkpoint model, in comparison with the traditional single-level checkpoint model. The two-type checkpoint model is able to cover most of the checkpoint demands and is easy to use in practice. Many of the existing multilevel checkpoint toolkits provide two levels or types of checkpoints, such as the double in-memory checkpoint scheme [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF], two-level incremental checkpoint recovery scheme [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF], and Fault Tolerance Interface (FTI) [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF]. On the other hand, based on the in-depth analysis in our previous work [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], the two-level model (e.g., one level is storing checkpoints onto the local memory/disk and the other one is using partner-copy technology) is already able to deal with a vast majority of failure cases, since the probability with many components being crashed simultaneously is relatively low. Under the twolevel checkpoint model, we derive an optimal solution (i.e., the optimal checkpoint intervals for every checkpoint level) based on the checkpoint/recovery overheads and the failure probabilities for each level. The key contributions are listed below.

• We derive an online optimal solution (where the fault-free job length 2 is unknown a priori) to determine the optimal checkpoint intervals for the two checkpoint levels/types. Compared with our previous work [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], we rebuild the checkpoint model and remove some assumptions or conditions that had to be provided in the previous optimization work. For example, our previous work is an offline solution because it requires knowing the fault-free job length before the optimization. It also assumed that the failures would not occur during the checkpoint and that the total wall-clock length must be close to the fault-free job length. Such assumptions simplified the problem but also introduced skewness in the optimization such that the derived solution deviates from the best experimental solution significantly in some cases, as confirmed by our experiments. On the contrary, we do not rely on any firstorder approximation in this paper. • We theoretically prove that the optimal solution is to set the checkpoints periodically with the same checkpoint intervals at each particular level, and we analytically determine the optimal intervals. To confirm the optimality of our online solution, we also derive an offline optimal solution by assuming that the job length is known already. We prove that the level-1 (i.e., type-1) checkpoint intervals are the same for the online solution and offline solution.

• We evaluate our solution by using an elaborate HPC simulation testbed, which is constructed according to the real execution progress of an MPI program. We compare our new solution with the solutions of the most up-to-date related research such as [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF]. The accuracy of the simulator is also confirmed by running a MPI scientific benchmark (called Heat Distribution) on a real cluster environment with up to 1024 cores. We highlight three significant findings: (1) our online optimal solution leads to exactly the same results with the offline optimal solution, which confirms the optimality of our online solution; (2) our online optimal solution always outperforms other optimized solutions and improves the performance significantly in some cases. Specifically, with the new solution, the wall-clock time can be reduced by up to 25.3% compared with that of other state-of-the-art approaches; and (3) the optimal solution derived in this paper is close to the best experimental solution obtained via bruteforce search, with the difference (in wall-clock time) less than 1% in most cases. The rest of the paper is organized as follows. In Section 2, we formulate the two-level checkpoint problem and propose a flexible pattern-based checkpoint model, in which any HPC execution can be split into multiple consecutive patterns, each containing a number of low-level checkpoint intervals (also called chunks in the paper). In Section 3, we compute the expected execution time for a pattern and show that it is minimized when all the chunks in the pattern have the same size. This result is key to deriving the optimal solution, which is presented in Section 4 for online scheduling and in Section 5 for offline scheduling. Although the whole derivation is technically involved, the solution is still simple to implement. In Section 6, we provide an in-depth analysis of how to use our derived solution in practice. In Section 7, we evaluate the proposed solution on an extreme-scale simulation environment and compare it with other state-of-the-art approaches. In Section 8, we discuss related work and highlight the differences with our new solution. In Section 9, we provide concluding remarks and directions for future work.

PROBLEM FORMULATION

In this section, we formally present the problem under study. We consider a failure-prone platform where two types of faults can strike. Type-1 faults are independent and identically distributed (IID) and follow an exponential distribution of failure rate λ 1 . Type-2 faults are IID and follow an exponential distribution of failure rate λ 2 . Faults of different types are independent. Two types of checkpoints are possible.

• A type-2 checkpoint takes time C 2 . Recovery from a type-2 checkpoint takes time R 2 . A type-2 checkpoint enables recovery from both type-1 and type-2 faults. • A type-1 checkpoint takes time C 1 . Recovery from a type-1 checkpoint takes time R 1 . A type-1 checkpoint enables recovery only from type-1 faults. In other words, type-1 checkpoints are destroyed by type-2 faults. After a fault or failure 3 occurs, there is a downtime of length D followed by a recovery of the same type (e.g., a recovery from a type-2 checkpoint after a type-2 fault). In our analysis, for simplicity we assume that faults can strike during checkpoints but not during recoveries. Consider the situation in which a type-2 fault strikes during a type-1 recovery following a type-1 fault. Then the actual recovery should be from a type-2 checkpoint, rather than from a type-1 checkpoint. On the one hand, taking into account all the possible scenarios of failures striking during 3. Fault often refers to the transient error, whereas failure means a fail-stop error. We will use the two terms interchangeably in the following text, in that our generic model is suitable for both of them. recoveries would lead to a complex analysis. On the other hand, the probability of a failure striking the system during a recovery is low because recoveries are expected to be faster than checkpoints 4 . Therefore, in the analysis we assume that faults or failures do not strike during recoveries, since we expect that this approximation will have only a limited impact. In fact, this expectation will be corroborated by the performance evaluation (Section 7) because in our simulations the failures can strike at any time, including during recoveries.

Intuitively, type-2 faults are more dramatic than type-1 faults because they delete all type-1 data. As a motivating example, consider a two-level memory system, with main memory and hard disks. Type-1 checkpoints are kept in main memory and enable recovery form soft errors [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF], while type-2 checkpoints are kept in stable storage and enable recovery from situations where the content of the main memory is lost. With existing multilevel checkpoint tools (such as FTI [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF]), a type of checkpoint can be set based on two modes: according to either the workload processed (i.e., the productive time) or the number of other types of checkpoints taken during the execution. For example, the type-2 checkpoints can be taken every 100 seconds (mode-A) or every five type-1 checkpoints (mode-B) in the execution. Mode-A checkpointing is also called interval-based checkpointing, and mode-B is also known as pattern-based checkpointing. Our checkpoint solution is suitable for both modes. In what follows, we will formulate the problem based on mode-B and extend it to mode-A later. We define a checkpoint pattern as follows (see Figure 1). Definition 1. A pattern PATTERN(K, W, w 1 , ..., w K ) (or simply PATTERN(K, W ) when there is no ambiguity) is a sequence of K computational chunks. The ith chunk is of size w i and K i=1 w i = W . Each chunk is followed by a type-1 checkpoint. The last type-1 checkpoint is followed by a type-2 checkpoint.

Throughout the paper, we assume unit-speed execution, so we can speak of time or work size interchangeably. A key technical contribution of this paper (provided in Section 3) is to compute the expected time E(PATTERN(K, W, w 1 , ..., w K )) to execute a given pattern PATTERN(K, W, w 1 , ..., w K ). Then we envision two scenarios, depending on whether the total fault-free job length is known a priori or not. Both scenarios use the value of E(PATTERN(K, W, w 1 , ..., w K )) computed in Section 3.

Online scheduling

The first scenario, where the total job length is not known, applies to online scheduling problems. We can also think of jobs running a very long (even infinite) time on the platform. In this scenario, we seek the best pattern, namely, the one whose overhead is minimal. The overhead of a pattern is defined as follows. [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF]. Recovery/restart time is expected to be smaller than checkpoint time. Our experiments with the real-world simulation benchmark protected by FTI toolkits confirmed this point (as shown in Figure 2). Moreover, the difference between recovery time and checkpoint time would be even larger with data lossy compression techniques, because recent evaluation [START_REF] Di | Fast error-bounded lossy hpc data compression with sz[END_REF] indicates that the decompression time is generally only 1/5 of the compression time for an error-bounded lossy compressor.

Definition

2. The overhead of the pattern PATTERN(K, W, w 1 , ..., w K ) is the ratio

OVERHEAD(PATTERN(K, W, w 1 , ..., w K )) = E(PATTERN(K, W, w 1 , ..., w K )) W -1. (1)
One main contribution of this paper is to provide an algorithm to compute the optimal pattern PATTERN(K online , W online ). First we compute the best pattern with a fixed number of chunks K, and then we show how to determine the optimal value K online of K. While the derivation is technical, the approach is fully constructive and determines the total length of the pattern W online and the size of each of the K online chunks inside the pattern.

Determining the optimal pattern, and then repeatedly using it until job completion, is the optimal approach with exponential failure distributions and infinite jobs. Indeed, once a pattern is successfully executed, the optimal strategy is to re-execute the same pattern. This is because of the memoryless property of exponential distributions: the history of failures has no impact on the solution, so if a pattern is optimal at some point in time, it stays optimal later in the execution, because we have no further information about the amount of work still to be executed.

Offline scheduling

The second scenario assumes that the total work length W total is known before execution, either exactly or through an accurate estimation based on benchmarking and sampling (this is the case for some scientific applications such as numerical linear algebra libraries [START_REF] Choi | The design and implementation of the ScaLAPACK LU, QR, and Cholesky factorization routines[END_REF]). The optimal strategy for this scenario is more complicated to derive. Indeed, after executing a first pattern PATTERN(K 1 , W 1 ), the remaining work to execute is W total -W 1 , and the second pattern might have a different length W 2 . We derive an optimal solution when we relax the problem and allow patterns to have a noninteger number of chunks. In this case, we show that all patterns have the same length, and we fully characterize the optimal solution.

For very long jobs (i.e., for large values of W total ), the optimal online strategy (repeating the pattern PATTERN(K online , W online )) also can be used, at the price of a shorter (suboptimal) pattern at the end of the computation when W total is not a multiple of W online . In Section 7, we experimentally compare both approaches.

PERFORMANCE OF PATTERNS

In this section, we analyze the expected execution time of a pattern and prove that the best patterns have same-size chunks.

Expected execution time of a pattern

We first determine the expectation E(PATTERN(K, W, w 1 , ..., w K )) of the execution time for an arbitrary pattern composed of K chunks, where the ith chunk has size w i , with K i=1 w i = W . The derivation is lengthy and technical, and the reader may want to skip the proof. We need a few definitions before stating the main result. For convenience, main definitions are summarized in Table 1. Optimal chunk size for a pattern with K chunks Definition 3. Let λ = λ 1 + λ 2 , L = λ2 λ and R = 1+λ1R1+λ2R2 λ + D. For 1 ≤ j ≤ K, let e j = e λ(wj +C1) -1. Let F 0 = 0 and, for 1 ≤ j ≤ K, let

F j = j k=1 L k-1 σ k (e 1 , ..., e j )
where σ k is the symmetric function of order k.

In Defnition 3, λ is the cumulated failure rate, which includes both failure types, while L is the fraction of type-2 failures among all failures. Therefore, 0 ≤ L ≤ 1, with L = 0 if there are no type-2 failures, and L = 1 if there are no type-1 failures. Finally, R is some averaged overhead, which simplifies to R = 1 λ1 + D + R 1 if there are no type-2 failures. Also, here is an example for the F j 's with K = 3: we have σ 1 (e 1 , e 2 , e 3 ) = e 1 + e 2 +e 3 , σ 2 (e 1 , e 2 , e 3 ) = e 1 e 2 +e 1 e 3 +e 2 e 3 and σ 3 (e 1 , e 2 , e 3 ) = e 1 e 2 e 3 . We obtain E(PATTERN(K, W, w 1 , ..., w K ))

F 1 = e 1 , F 2 = e 1 + e 2 + L(e 1 e 2 ) = F 1 + (1 + LF 1
= RF K + R (1 + LF K ) e λC2 -1 (2)
Proof. We introduce some simplified notation to ease the reading of the proof. Let E(PATTERN) denote E(PATTERN(K, W, w 1 , ..., w K )). Let E(w 1 ) be the expected time to successfully execute the first chunk w 1 and the following type-1 checkpoint C 1 . For i ∈ [2, K], let E(w i ) be the expected time to successfully execute chunk w i and the following type-1 checkpoint, after the completion of the (i -1)th type-1 checkpoint. Let E(C 2 ) be the expected time to successfully execute the final type-2 checkpoint in the pattern, after the Kth (and last) type-1 checkpoint. By the linearity of expectations,

E(PATTERN) = K i=1 E(w i ) + E(C 2 ). (3)
Computing E(w i ) Consider the beginning of the execution of the ith chunk w i , right after the successful completion of the (i -1)th type-1 checkpoint. Let X 1 be the time at which the first type-1 fault strikes the system. X 1 is thus a random variable following an exponential distribution of failure rate parameter λ 1 . Let X 2 be the time at which the first type-2 fault strikes the system. X 2 is thus a random variable following an exponential distribution of parameter λ 2 . Let X be the time at which the first fault (of any type) strikes the system. We consider three cases:

1) No fault strikes the system before time

w i + C 1 (expected completion time E 0 )
2) A fault strikes the system before time w i + C 1 , and the first fault to strike the system is a type-1 fault (expected completion time E 1 ). In this case, we can recover from the last type-1 checkpoint and re-execute only the current chunk. This holds true even if i = 1: for the first chunk of the pattern the most recent checkpoint is a type-2 checkpoint, but it was preceded by a type-1 checkpoint, which is cheaper to recover from. 3) A fault strikes the system before time w i + C 1 , and the first fault to strike the system is a type-2 fault (expected completion time E 2 ). In this case, we have to recover from the last type-2 checkpoint, and we have to re-execute all the chunks of the pattern, together with their type-1 checkpoints, starting again from w 1 . By the law of total probability, E(

w i ) = E 0 + E 1 + E 2 .
We compute E 0 , E 1 , and E 2 as follows:

Case 1. Without fault, we have X = min{X 1 , X 2 } ≥ w i + C 1 .
Because the minimum of two exponential distributions is an exponential distribution whose parameter is the sum of the parameters, this happens with probability e -(λ1+λ2)(wi+C1) . Hence, E 0 = e -(λ1+λ2)(wi+C1) (w i + C 1 ). Case 2. When a type-1 fault strikes first, we have X 1 < w i + C 1 and X 1 ≤ X 2 . Then the system performs a type-1 recovery and retries processing the same amount of work. The re-execution has expectation E(w i ) by the memoryless property. Hence, we have the following (see the extended version [START_REF] Di | Toward an Optimal Online Checkpoint Solution under a Two-Level HPC Checkpoint Model[END_REF] for intermediate steps in the computation of E 1 ):

E 1 = X1=wi+C1 X1=0 X2=+∞ X2=X1 (X 1 +D+R 1 +E(w i )) × λ 1 e -λ1X1 λ 2 e -λ2X2 dX 2 dX 1 = X1=wi+C1 X1=0 (X 1 +D+R 1 +E(w i ))λ 1 e -(λ1+λ2)X1 dX 1 = (D + R 1 + E(w i )) X1=wi+C1 X1=0 λ 1 e -(λ1+λ2)X1 dX 1 + X1=wi+C1 X1=0 X 1 λ 1 e -(λ1+λ2)X1 dX 1
We derive that

E 1 = (D + R 1 + E(w i )) λ1 λ1+λ2 1 -e -(λ1+λ2)(wi+C1) -λ1 λ1+λ2 (w i + C 1 )e -(λ1+λ2)(wi+C1) + λ1 (λ1+λ2) 2 1 -e -(λ1+λ2)(wi+C1) .
Case 3. When a type-2 fault strikes first, we have X 2 < w i + C 1 and X 2 ≤ X 1 . Then we have to recover from the last type-2 checkpoint and re-execute all chunks from w 1 up to w i . Hence, we have the following.

E 2 = X2=wi+C1 X2=0 X1=+∞ X1=X2 (X 2 +D+R 2 +E(w 1 )+...+E(w i )) = (D+R 2 +E(w 1 )+...+E(w i )) X2=wi+C1 X2=0 λ 2 e -(λ1+λ2)X2 dX 2 + X2=wi+C1 X2=0 X 2 λ 2 e -(λ1+λ2)X2 dX 2
Transposing the integral computations of Case 2, we obtain the following.

E 2 = (D + R 2 + E(w 1 ) + ... + E(w i )) λ2 λ1+λ2 × 1 -e -(λ1+λ2)(wi+C1) -λ2 λ1+λ2 (w i + C 1 )e -(λ1+λ2)(wi+C1) + λ2 (λ1+λ2) 2 1 -e -(λ1+λ2)(wi+C1) Since E(w i ) = E 0 + E 1 + E 2 , we derive: E(w i ) = e -(λ1+λ2)(wi+C1) (w i + C 1 ) +(D + R 1 + E(w i )) λ1 λ1+λ2 × 1 -e -(λ1+λ2)(wi+C1) -λ1 λ1+λ2 (w i + C 1 )e -(λ1+λ2)(wi+C1) + λ1 (λ1+λ2) 2 1 -e -(λ1+λ2)(wi+C1) +(D + R 2 + E(w 1 ) + ... + E(w i )) λ2 λ1+λ2 × 1 -e -(λ1+λ2)(wi+C1) -λ2 λ1+λ2 (w i + C 1 )e -(λ1+λ2)(wi+C1) + λ2 (λ1+λ2) 2 1 -e -(λ1+λ2)(wi+C1) = 1+λ1R1+λ2(R2+E(w1)+...+E(wi-1)) λ1+λ2 + D + E(w i ) × 1 -e -(λ1+λ2)(wi+C1)
This leads to

E(w i ) = 1 + λ 1 R 1 + λ 2 (R 2 + i-1 j=1 E(w j )) λ 1 + λ 2 + D × e (λ1+λ2)(wi+C1) -1 .
With the notations of Definition 3, we obtain

E(w i ) =   R + L i-1 j=1 E(w j )   e i . (4)
Simplifying E(w i ) Lemma 1.

E(w i ) = R(1 + LF i-1 )e i for 1 ≤ i ≤ K (5)
Proof. Recall the notation given in Definition 3. We prove the result by induction. It holds for i = 1 because Equation ( 5) leads to E(w 1 ) = Re 1 , since F 0 = 0. Assume that the result holds up to rank i, and consider E(w i+1 ). From Equation (4), we have

E(w i+1 ) =   R + L i j=1 E(w j )   e i+1 = R   1 + L i j=1 (1 + LF j-1 )e j   e i+1 , (6)
the last equality coming from the induction hypothesis. We further detail the term inside the summation:

(1 + LF j-1 )e j = e j + Le j j-1 k=1 L k-1 σ k (e 1 , .
.., e j-1 ) = e j + j-1 k=1 L k σ k (e 1 , ..., e j-1 )e j = e j + j-1 k=1 L k (σ k+1 (e 1 , ..., e j ) -σ k+1 (e 1 , ..., e j-1

)) = e j + j-1 k=1 L k σ k+1 (e 1 , ..., e j ) -j-2 k=1 L k σ k+1 (e 1 , ..., e j-1 )
Plugging back into Equation ( 6), we obtain the following.

E(w i+1 ) = R 1+L i j=1 (1+LF j-1 )e j e i+1 = R 1+L i j=1 e j + j-1 k=1 L k σ k+1 (e 1 , ..., e j ) - j-2 k=1 L k σ k+1 (e 1 , ..., e j-1 ) e i+1 = R 1+Lσ 1 (e 1 , ..., e i )+ i j=1 j-1 k=1 L k+1 σ k+1 (e 1 , ..., e j ) -i j=1 j-2 k=1 L k+1 σ k+1 (e 1 , ..., e j-1 ) e i+1 = R 1+Lσ 1 (e 1 , ..., e i )+ i j=2 j-1 k=1 L k+1 σ k+1 (e 1 , ..., e j ) -i j=3 j-2 k=1 L k+1 σ k+1 (e 1 , ..., e j-1 ) e i+1 = R 1+Lσ 1 (e 1 , ..., e i )+ i j=2 j-1 k=1 L k+1 σ k+1 (e 1 , ..., e j ) -i-1 j=2 j-1 k=1 L k+1 σ k+1 (e 1 , ..., e j ) e i+1 = R 1+Lσ 1 (e 1 , ..., e i )+ i-1 k=1 L k+1 σ k+1 (e 1 , ..., e i ) e i+1 = R 1 + i-1 k=0 L k+1 σ k+1 (e 1 , ..., e i ) e i+1 = R 1+L i k=1 L k-1 σ k (e 1 , ..., e i ) e i+1 = R(1+LF i ) e i+1
Using Equation ( 5), we simplify the expression j i=1 E(w i ). Indeed, comparing with Equation ( 4), we obtain

i j=1 E(w j ) = RF i for 0 ≤ i ≤ K. (7)
Computation of E(C 2 )

We compute E(C 2 ), the expected time to complete the final type-2 checkpoint of the pattern, similar to the way we computed E(w i ) before. We consider three cases: 1) No fault strikes the system before time C 2 . Then the expectation is E 0 = e -(λC2) C 2 . 2) A fault strikes the system before time C 2 , and the first fault to strike the system is a type-1 fault. Then we perform a type-1 recovery and re-execute C 2 . The expectation is

E 1 = X1=C2 X1=0 X2=+∞ X2=X1 (X 1 + D + R 1 + E(C 2 ))
× λ 1 e -λ1X1 λ 2 e -λ2X2 dX 2 dX 1 .

3) A fault strikes the system before time C 2 , and the first fault to strike the system is a type-2 fault: X 2 < C 2 and X 2 ≤ X 1 . Then we recover from the last type-2 checkpoint and re-execute the whole pattern. The expectation is

E 2 = X2=C2 X2=0 X1=+∞ X1=X2 (X 2 +D+R 2 + K i=1 E(w i )+E(C 2 )) × λ 2 e -λ2X2 λ 1 e -λ1X1 dX 1 dX 2 .
We compute the integrals as before and use

K i=1 E(w i ) = RF K to derive the final result. E(C 2 ) = R (1 + LF K ) e λC2 -1 (8)
We have from Equations ( 7) and ( 8) that

K i=1 E(w i ) + E(C 2 ) = RF K + R (1 + LF K ) e λC2 -1 ,
which concludes the proof of Proposition 1.

Best patterns have equal-size chunks

Proposition 1 enables us to prove a key result: when both the total size W and the number of chunks K of a pattern are known, then the expectation of the execution time E(PATTERN(K, W, w 1 , ..., w K )) is minimized when all chunks have same size W/K. Note that the overhead of the pattern is also minimized in that case, because W is given.

Proposition 2. The expectation E(PATTERN(K, W, w 1 , ..., w K )) of a pattern with K chunks and total size W = K i=1 w i is minimized when w i = W K (same-size chunks). Letting E(PATTERN(K, W )) denote this minimal value, we have

E(PATTERN(K, W )) = α + β L (1 + L(e λ( W K +C1) -1)) K , ( 9
)
where α = R(e λC2 -1) -β L and β = R(1 + L(e λC2 -1)). Equation ( 9) simplifies when λ 2 = C 2 = 0, which corresponds to type-1 failures only. In that case,

E(PATTERN(K, W )) = K( 1 λ1 + D + R 1 )(e λ W K +C1 -1)
, in accordance with [START_REF] Bougeret | Checkpointing strategies for parallel jobs[END_REF]. Similarly when λ 1 = C 1 = 0 (only type-2 failures), the pattern reduces to a single chunk of length W and E(PATTERN(K, W )) = ( 1λ2 + D + R 2 )(e λ(W +C2) -1). Proof. Consider a pattern whose number of chunks K and total size W = K i=1 w i are given. From Equation ( 2), the expected execution time is minimal when F K is minimal. Recall the definition of F j :

F j = j k=1 L k-1 σ k (e 1 , ..., e j )
, where e i = e λ(wj +C1) -1. Note that F j is a symmetric function of e 1 , ..., e j . Consider any two indices i and j between 1 and K, and fix the sum w i + w j = ω as well as all the w k 's except w i and w j . Then we can rewrite F K as a sole function of w i . Then there exist two constants γ 1 and γ 2 , with γ 1 > 0, such that

F K (w i ) = γ 1 e λ(wi+C1) + γ 1 e λ(ω-wi+C1) + γ 2 .
By convexity, F K (w i ) is minimal when w i = w j = ω 2 . By induction, F K is minimized when all the w j 's are equal.

We now compute the value of F K when all chunks have same size w = W K . We show by induction on j that

F j = 1 L 1 + L e λ(w+C1) -1 j -1 , (10)
when w 1 = ... = w K = w. For j = 1, Equation ( 10) is

F 1 = e λ(w+C1) -1, which is the definition of F 1 when w 1 = w.
Assume that Equation ( 10) holds up to some index i included.

From Equations ( 5) and ( 7) we have F i = i j=1 e j (1 + LF j-1 ). Then

F i+1 = i+1 j=1 e j (1 + LF j-1 ) = i+1 j=1 e λ(w+C1) -1 1+L 1 L 1+L e λ(w+C1) -1 j-1 -1 = i+1 j=1 e λ(w+C1) -1 1+L e λ(w+C1) -1 j-1 = i+1 j=1 (-1+1+L(e λ(w+C 1 ) -1)) L 1 + L e λ(w+C1) -1 j-1 = -1 L i+1 j=1 1 + L e λ(w+C1) -1 j-1 + 1 L i+1 j=1 1 + L e λ(w+C1) -1 j = 1 L 1 + L e λ(w+C1) -1 i+1 -1 ,
which proves the result for F K . Plugging this value into Equation (2) leads to the result.

We give an example to describe how to compute the expected wall-clock time by Proposition 2 (or Equation ( 9)). Suppose the checkpoint/recovery overheads on the two levels are C 1 =R 1 =20 seconds and C 2 =R 2 =50 seconds respectively. The failure rates on the two levels are 24/86,400=2.78×10 -4 and 4/86,400=4.63×10 -5 respectively, which means that there are 24 (4) failures occurring at level 1 (level 2) per day. The downtime is set to 0 for simplicity. Consider a pattern with K=4 same-size chunks and the workload size W is 1,472 seconds, then the chunk size is equal to 368 seconds in length, and the expected wall-clock length can be computed as 1773.2 seconds based on Equation (9).

ONLINE SCHEDULING

In this section, we characterize the optimal pattern for jobs whose length is not known in advance. We first determine the optimal pattern when the number K of chunks is given, and we then derive the best value of K.

Optimal pattern with given number K of chunks

Proposition 3. When the number of chunks, K, is given, the optimal pattern has same-size chunks w i = w opt (K), where w opt (K) is the unique solution of the equation

βλKwe λ(w+C1) (1 + L(e λ(w+C1) -1)) K-1 = α + β L (1 + L(e λ(w+C1) -1)) K (11)
The pattern overhead is then

OVERHEAD(K) = α + β L (1+L(e λ(wopt (K)+C1) -1)) K Kw opt (K) -1.
Proof. Consider a pattern whose number of chunks K is given.

If the total amount of work W is given, Proposition 2 shows that the overhead of the pattern is minimal when chunks have samesize W/K. But what is the optimal pattern size W opt (K) for a pattern with K chunks? Or equivalently, what is the optimal size w opt (K) of the chunks for a pattern with K chunks (and then the total size of the pattern will be W opt (K) = Kw opt (K))? Consider a pattern with K chunks of same size w. We express its overhead as a function f (w) of the chunk size w. Using Proposition 2 (with the constants α and β defined there), we have

f (w) = α + β L 1 + L e λ(w+C1) -1 K Kw -1. (12)
To simplify the derivation, let E(w) = e λ(w+C1) and N (w) = 1 + L(E(w) -1). Note that E(w) ≥ 1 and N (w) ≥ 1 for all w ≥ 0. We have f (w) =

α+ β L (N (w)) K Kw -1.
Differentiating, we obtain that f ′ (w) = 0 if and only if w is a solution of the equation

βλKwE(w)(N (w)) K-1 = α + β L (N (w)) K , (13)
which is identical to Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF].

To conclude the proof, we show that Equation ( 13) has a unique positive solution w opt (K) for any given value of K. Let g(w) = βλKwE(w)(N (w)) K-1 -α-β L (N (w)) K . Differentiating again, we get g ′ (w) = βλ 2 KwE(w)(N (w)) K-2 (N (w)+ (K -1)LE(w)), which is always positive. Hence g(w) is increasing. We have

g(0) = -α -β L (N (0)) K and N (0) ≥ 1; hence g(0) ≤ -(α + β L ) = -R(e λC2 -1) < 0 . Furthermore,
we easily see that lim w→+∞ g(w) = +∞. Hence g(w) has a unique positive zero, which shows that Equation ( 11) has a unique minimum.

Remark 1. Equation ( 11) simplifies when K = 1 and λ 2 = C 2 = 0, which corresponds to the optimal single-chunk pattern for type-1 failures of rate λ = λ 1 . In that case, L = 1, α = -β.

After some easy transformations, we obtain

(λw opt (1) -1)e λwopt (1)-1 = L(-λC 1 -1). ( 14
)
Here L is the Lambert function, defined as

L(z) = x if xe x = z.
We retrieve a result of [START_REF] Bougeret | Checkpointing strategies for parallel jobs[END_REF] that, when approximated to the first order, is similar to a result of [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF].

Optimal pattern

In this section we show how to determine the number K online of chunks in the optimal pattern. From Proposition 3 we know that each chunk will have length w opt (K online ) and that the total length of the optimal pattern will be W online = K online × w opt (K online ).

Theorem 1. The number of chunks K online in the optimal pattern is either max{1, ⌊K * ⌋} or ⌈K * ⌉ (whichever leads to the pattern with smallest overhead), where K * is the unique solution of the equation

N (w * ) ln(N (w * )) = λLw * E(w * ) (15)
Here w * = w opt (K * ) is the unique solution of Equation ( 11)

when K = K * .
Proof. From Proposition 3, the overhead of the optimal pattern with K chunks is

OVERHEAD(K) = α+ β L (1+L(e λw opt (K)+C 1 -1)) K Kwopt (K) -1,
where w opt (K) is the unique solution of Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF]. We note that we can still define w opt (K) as the unique solution of Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF] for arbitrary values of K, not just integer ones. This allows us to define OVERHEAD(K) for (positive) real values of K. We study this function OVERHEAD(K) and will show that it has a unique minimum K * , which is the unique solution of Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF]. To simplify notations, let f (K) = OVERHEAD(K) for K > 0. To differentiate f (K), we need to differentiate w opt (K) as a function of K. We know w opt (K) only implicitly, through Equation ( 11), but we can differentiate both sides of this equation to determine w ′ opt (K). We will compute w ′ opt (K) later when proving existence and unicity of the solution of Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF]. We first differentiate f (K). We have

         (E(w opt (K))) ′ = λE(w opt (K))w ′ opt (K) (N (w opt (K))) ′ = L(E(w opt (K))) ′ (N K (w opt (K))) ′ = (e K ln(N (wopt (K))) ) ′ = (ln(N (w opt (K))) + K (N (wopt (K))) ′ N (wopt (K)) )N K (w opt (K))
from which we derive

f ′ (K) = 1 K 2 (wopt (K)) 2 Kw opt (K) β L (N K (w opt (K))) ′ -(α + β L N K (w opt (K)))(Kw ′ opt (K) + w opt (K)) .
The numerator of f ′ (K) is of the form Aw ′ opt (K) + Bw opt (K), where

A = Kw opt (K) β L N K (w opt (K)) KLλE(wopt (K)) N (wopt (K)) -K(α + β L N K (w opt (K))) B = K β L N K (w opt (K)) ln(N (w opt (K))) -(α + β L N K (w opt (K))
). From Equation (13), we obtain A = 0, which means that the term in w ′ opt (K) cancels out, and

B = K β L N K-1 (w opt (K)) [N (w opt (K)) ln(N (w opt (K))) -Lλw opt (K)E(w opt (K))]
which shows that the sign of f ′ (K) is that of

g(K) = N (w opt (K)) ln(N (w opt (K))) -λLw opt (K)E(w opt (K)).
In particular, f ′ (K) = 0 if and only if Equation ( 15) is satisfied.

Here, g(K) uses only the value w opt (K), but we can study the function

h(w) = N (w) ln(N (w)) -λLwE(w)
for arbitrary positive values of w. We have g(K) = h(w opt (K)).

Thus, if we show that h(w) > 0 for all w, we will derive that g(K) > 0 for all K. We differentiate h(w) as follows.

h ′ (w) = λLE(w)[ln(N (w)) -λw]
We need to study yet another auxiliary function to determine the sign of h ′ (w). Let h 1 (w) = ln(N (w)) -λw for w > 0. Recall that N (w) = 1 + L(E(w) -1) and E(w) = e λ(w+C1) . We get h ′ 1 (w) = λ(L-1) N (w) < 0 because L < 1; hence, h 1 (w) is decreasing. We have h 1 (0) > 0 and lim w→+∞ h 1 (w) = ln(L)+ λC 1 . We study two cases: (1) when the latter limit is positive, then h 1 (w) is always positive; and (2) when it is negative, then h 1 (w) is positive for w > w 0 and negative for w < w 0 , for some w 0 whose exact value does not matter for this proof.

For case (1): h ′ (w) > 0 for all w > 0; hence, h(w) is increasing. Recall that N (0) > 1; hence h(0) = N (0) ln(N (0)) > 0, and h(w) is positive for all w, which implies that g(K) is positive for all K. We conclude that f ′ (K) is positive for all K; hence f (K) is increasing, and the optimal pattern is obtained for K = 1, with a single chunk.

For case (2): h ′ (w) > 0 for w < w 0 , and h ′ (w) < 0 for w > w 0 ; hence h(w) is increasing up to w = w 0 and then decreasing after that, down to some limit ℓ. We have h(0) > 0 as before; so if ℓ ≥ 0, then h(w) is positive for all w, and we get the same conclusion as for case [START_REF] Feinberg | An 83,000-processor supercomputer can only match 1% of your brain[END_REF]. On the contrary, if ℓ < 0, then h(w) has a unique zero w 1 > w 0 , and h(w) is positive for w < w 1 and negative for w > w 1 . This implies that g(K) is positive if w opt (K) < w 1 and negative if w opt (K) > w 1 . In the following we prove that w opt (K) is a decreasing function of K and that lim K→∞ w opt (K) = 0. This result was expected: the more chunks in a pattern, the shorter these chunks to trade off failure-free overhead and re-execution time. Assuming these results hold true, we again have two cases. If w opt (1) < w 1 , we have w opt (K) < w 1 for all K; hence g(K) and f ′ (K) are always positive. Then the best pattern is obtained with K = 1, as before. But if w opt (1) > w 1 , then g(K) remains negative while w opt (K) > w 1 and becomes positive when w opt (K) < w 1 (which will eventually happen, because w opt (K) is decreasing down to its limit 0). Then f (K) is decreasing to a unique minimum K * obtained when g(K * ) = 0, which is the solution of Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF]. The best pattern requires an integer number of chunks; hence it is obtained with either max{1, ⌊K * ⌋} or ⌈K * ⌉ chunks.

To conclude the proof, we need to show that w opt (K) is a decreasing function of K and that lim K→∞ w opt (K) = 0. We first show that w ′ opt (K) < 0 for all K ≥ 1. To do so, we differentiate Equation ( 13) and derive that P w ′ opt (K) = Q, where (after some simplification)

P = λ 2 βKE(w opt (K))N K-2 (w opt (K)N (w opt (K)) +L(K -1)) Q = -α ln(N (w opt (K))) -1 K (α + β L N K (w opt (K))
). We have P > 0 (obvious) and will show that Q < 0, so that w ′ opt (K) < 0 for all K. To obtain the latter inequality Q < 0, we rewrite Q as

Q = -(α + β L )(ln(N (w opt (K))) + 1 K ) + β L (ln(N (w opt (K))) + 1-N K (wopt (K)) K
).

We have α + β L = R(e λC2 -1) > 0. We will show that ln(N (w opt (K)))+

1-N K (wopt (K))) K ) < 0, which leads us to con- clude that Q < 0. Consider the function h 2 (x) = ln(x) + 1-x K K for x ≥ 1 (remember that N (w opt (K)) > 1). We have h ′ 2 (x) = 1-x K x < 0 and h 2 (1) = 0, hence h 2 (x) < 0 for x > 1.
To determine lim K→∞ w opt (K), we consider Equation ( 11) and divide both sides by N (w opt (K)) K . We obtain

βλKw opt (K) E(w opt (K)) N (w opt (K)) = α N (w opt (K)) K + β L .
Since we know that E(w) and N (w) are increasing functions of w, for all w > 0 we have

1 L + 1-L E(0) ≤ E(w) N (w) = 1 L + 1-L E(w) ≤ 1 L β L ≤ α N (w) K + β L ≤ α N (0) K + β L .
Hence Kw opt (K) is bounded above and below for all K by two constants. This shows that lim K→∞ w opt (K) = 0, and concludes the proof. We here use the example presented in the end of Section 3.2 to describe how to obtain the optimal number of chunks K * and the optimal chuck size w * by combining Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF] and Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF]. First, w * can be computed by conducting Newton's method on Equation (15) because w is the unique variable in this equation. For the example presented in Section 3.2, w * is computed as 368 seconds. Then, we can obtain K * by running Newton's method for Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF] with the computed w * . In that example, the optimal number of chunks K * would be computed as 3.51 by Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF]. Note that the value of K * computed here is not an integer, so we need to tune the solution for fitting the practice, either by selecting a rounded integer number to stick with the pattern-based checkpoint model or by extending the solution to the intervalbased checkpoint model with rational number of K * , which will be discussed in more details later.

OFFLINE SCHEDULING

In this section we discuss the best strategy for offline scheduling, that is, when the total amount of work W total (i.e., fault-free job length) is known in advance. In fact, one can always use the online approach and repeat the optimal pattern (of length W online ) until there remains a final piece of work of size smaller than W online to execute. When W total is not an exact multiple of W online , however, the last piece of work might introduce some extra overhead. For small values of W total , say requiring only a few patterns W online , this extra overhead may be significant. To this end, if the job length is already known, one may prefer to further improve the performance by adopting an offline solution instead of always using the online approach.

The question then is how to optimize the offline solution by determining the optimal number and size of the patterns. We will show that the optimal solution is to use the same-size patterns, and we will compute their optimal number (or equivalently, their optimal size). First, however, we address the following question: Given a pattern size W , what is the optimal number of chunks to execute it? This is the dual question of the one addressed in Section 4.1, where the number of chunks was given and we were searching for the optimal work size. Due to space limitations, we do not provide proofs in this section, but instead refer to the extended version [START_REF] Di | Toward an Optimal Online Checkpoint Solution under a Two-Level HPC Checkpoint Model[END_REF].

Optimal pattern with given size W

In this section we derive the best pattern when the work size W is given.

Proposition 4. When W is given, the optimal pattern of size W has K opt = max{1, ⌊K * ⌋} or K opt = ⌈K * ⌉} same-size chunks, where K * = W wopt and w opt is the unique solution (if it exists) of the following equation:

ln 1 + L e λ(wopt +C1) -1 = λw opt 1 + L -1 1 + L e λ(wopt +C1) -1 (16)
We observe that Equation ( 16) can be rewritten as Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF]. More precisely, if we let w * = w opt in Equation ( 15), we retrieve Equation [START_REF] Di | Fast error-bounded lossy hpc data compression with sz[END_REF]. Hence the optimal size w opt of the chunks in a pattern is always obtained in the same way, by solving Equation ( 16), and this result is true for any value of the pattern length W , not only when W = W online (as was shown by Equation ( 15)).

Solution with rational number of chunks

We are ready to determine the optimal offline solution when allowing rational numbers of chunks. Given a job whose total size is W total , we divide it into p patterns of size W 1 , W 2 , . . . , W p , where W total = p i=1 W i . From Proposition 4, once we know the pattern sizes, we know how to subdivide them into chunks. We will accept a fractional number of patterns and a fractional number of chunks for each pattern in our solution. In Section 6, we explain how to round up the solution. The following result gives the optimal number of patterns and their sizes: Proposition 5. The optimal rational offline solution for a job of total size W total is to divide it into p * same-size patterns, where p * = W total ln(N (w opt ))

L αL βe + 1 w opt , ( 17
)
where w opt is given by Equation ( 16), and letting N (w) = 1 + L(e λ(w+C1) -1) as before. Here L is the Lambert function (see Equation ( 14)).

DISCUSSION

In this section, we discuss how to implement the above-mentioned optimal solutions. In practice, we recommend using the intervalbased optimal online checkpoint solution, because it exhibits the best performance in our simulations (as shown in the next section).

As mentioned in Section 2, our derivation is based on a patternbased checkpoint mode (the mode-B of Section 2), which can be extended easily to an interval-based checkpoint mode (the mode-A). For this purpose, one simply needs to further compute the optimal level-2 checkpoint interval (denoted by w * opt2 ) based on the optimal chunk size w * opt and on the optimal number (denoted by K * ) of chunks in each pattern, as follows:

w * opt2 = K * • w * opt ( 18
)
Under a pattern-based mode, type-2 checkpoints are taken every K * type-1 checkpoints, and K * must be an integer. Under an interval-based mode, a new type-2 checkpoint is taken when the amount of workload successfully processed since the last type-2 checkpoint was taken is equal to w * opt2 . In the latter case, K * can take non-integer values in Equation [START_REF] Di | Toward an Optimal Online Checkpoint Solution under a Two-Level HPC Checkpoint Model[END_REF].

We have a total of four solutions based on either interval-or pattern-based checkpoint modes, as listed below.

• Online-Pattern-OPT: The first one is the optimal pattern solution without knowing the productive job length before execution. With this solution, the system will set the level-1 checkpoints periodically, based on the optimal level-1 checkpoint interval (denoted by w * opt ), and set level-2 checkpoints periodically too, based on the number of level-1 checkpoints taken (denoted by K * ) in the pattern. In our implementation, the values of w * opt and K * are obtained by solving Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF] and Equation ( 11) using Newton's method. Note that the optimized checkpoint frequency at level 2 is supposed to be an integer, so K * will be rounded to the nearest nonzero integer number in this solution. This solution is represented by {w * opt ,K * }, where K * is a positive integer. • Online-Interval-OPT: This solution is the optimal online checkpoint solution, in which the checkpoint positions (for both level 1 and level 2) during the execution are always set based on checkpoint intervals. Specifically, the level-1 checkpoint interval is the same as that of Online-Pattern-OPT, while the level-2 checkpoints will not be set based on the number of level-1 checkpoints taken but will be set periodically based on a level-2 checkpoint interval (denoted by w * opt2 ). In our implementation, we first compute w * opt and K * by solving Equation [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF] and Equation [START_REF] Li | Two-level incremental checkpoint recovery scheme for reducing system total overheads[END_REF], respectively, and then obtain w * opt2 based on Equation [START_REF] Di | Toward an Optimal Online Checkpoint Solution under a Two-Level HPC Checkpoint Model[END_REF]. This optimal solution is represented by {w * opt ,w * opt2 }. • Offline-Pattern-OPT: This solution is an offline solution, which requires to know the productive job length in advance. Similarly to the Online-Pattern-OPT solution, Offline-Pattern-OPT sets the level-1 checkpoints based on checkpoint intervals (i.e., the productive time used for processing the workload) and sets the level-2 checkpoints every K * level-1 checkpoints, where K * is a positive integer number. Implementing this method involves three steps: (1) given a known job length W total , use Proposition 5 to compute p * (i.e., the optimal number of patterns in the execution);

(2) compute the optimal level-1 checkpoint interval w * opt (i.e., chunk size) by solving Proposition 4 using Newton's method; and (3) compute the optimal number K * of chunks in each pattern by W total p * •w * opt . Similar to the Online-Pattern-OPT the value of K * computed will likely not be an integer number, so we will use its nearest nonzero integer number instead in practice.

• Offline-Interval-OPT: Similarly to Offline-Pattern-OPT, Offline-Interval-OPT requires knowing the job length beforehand. The difference is that level-1 checkpoints and level-2 checkpoints both will be set based on checkpoint intervals under this solution. In particular, w * opt is the same as the one in Offline-Pattern-OPT, and w * opt2 is equal to w * opt •K * .

PERFORMANCE EVALUATION

In this section, we present the evaluation results, by comparing our proposed optimal solution with other optimized multilevel checkpoint methods. We first describe the experimental setting and then present the evaluation results.

Experimental Setup

For our experimental evaluation, we have access to the Argonne FUSION cluster [START_REF]Fusion cluster[END_REF]. However, there are at most 128 physical nodes (for a grand total of 1,024 cores) available in that cluster with a limited resource usage quota. This is far too limited a scale to validate our results, since our fault-tolerance research is designed for exascale applications. Therefore, we have to use simulations to evaluate our multilevel checkpoint solutions. In addition, we validate the accuracy of our simulator, by performing practical experiments deployed with FTI [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF] and scientific simulation benchmark based on a real-world diffusion problem (2D heat distribution) over the Fusion cluster. The application used in our benchmark is a well-known MPI program, called Heat Distribution, whose MPI communication algorithms (such as the ghost array design between adjacent blocks) are commonly adopted in real scientific projects such as parallel ocean simulation [START_REF] Smith | The Parallel Ocean Program (POP) reference manual: Ocean component of the Community Climate System Model (CCSM)[END_REF], [START_REF]Community Earth System Model (CESM)[END_REF]. The key reason we adopted this application in our evaluation is that it involves many MPI function calls, including MPI_Bcast, MPI_Barrier, MPI_Recv, MPI_Send, MPI_Irecv, MPI_Isend, MPI_Waitall, and MPI_Allreduce. The checkpoint and recovery overheads are both dependent on two factors: the program memory size determined by the problem size, and the execution scale (i.e., the number of processes).

Our simulations are performed as follows. The simulator emulates the running of some MPI program by processing a fixed amount of parallel workload. The parallel productive time is fixed for a particular test case. Each simulation is driven by ticks (one tick is set to one second in the simulation), simulating the whole procedure of running the MPI program. Our simulator takes into account any possible unexpected events during the execution, such as the occurrence of failures during the recovery.

The simulator emulates the checkpoint/recovery procedure based on the real-world multi-level checkpoint toolkit FTI [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF], so there are a total of four checkpoint levels implemented in our simulator, including local disk storage, partner-copy technique, RS-encoding technique, and PFS file system. Local disk storage level allows the checkpoints to be stored on local disks of the execution nodes. Partner-copy technology makes each checkpoint file have two copies stored in the local storage device and another partner-node respectively. Thus, upon a failure event with multiple simultaneous hardware crashes, the whole execution can still be recovered via partner-copy as long as there are no adjacent/partner nodes crashed. Reed-solomon encoding is a more advanced technology, which allows the application to be recovered with M missing checkpoints. PFS is the top level, which means the checkpoints are to be stored in the parallel file system. Obviously, local storage corresponds to no hardware failure, partner-copy corresponds to nonadjacent node failure, Reed-solomon encoding corresponds to adjacent node failures with at most M failed nodes, and PFS corresponds to any other failure cases. In the evaluation of our propsoed solution, we take two out of the four optional checkpoint levels with different overheads.

We verify the accuracy of our simulation by comparing its performance results to that of the real runs on the Argonne FUSION cluster [START_REF]Fusion cluster[END_REF] with Heat Distribution benchmark using up to 1024 cores. All of the failures were injected by using the command "kill -9" to terminate some of the processes running by the MPI program. As suggested by the developer of FTI, we just need to remove some checkpoint files to emulate the missing nodes. For example, in order to emulate the single node failure, we select one node and remove its local checkpoint file and then restart the application using FTI. FTI automatically recovers the missing checkpoint file before moving on. In order to simulate the multiple node crashes simultaneously striking the application, we just need to delete multiple checkpoint files. The FTI would also automatically recover the multiple missing checkpoint files, which may suffer from longer recovery overhead.

The verification results are shown in Figure 2, which confirms very similar results between simulation and the real experiments. The failure rates injected follow a Poisson distribution, with the mean failure intervals being set to 861 seconds, 1722 seconds, 4305 seconds and 8610 seconds respectively. In absolute terms, with the same setting on the productive length and checkpoint/restart overheads, the difference between the simulation results and experimental results is less than 4%. Note that the model proposed in this work involves two levels instead of four, so the users can select any two of the checkpoint levels based on the characterization or analysis of the failure rates in practice, as we did in our evaluation. Our evaluation involves nine test cases, regarding various checkpoint/restart overheads, different failure rates, and different application execution lengths. All of the test cases respect the observations based on our experiments with the benchmark running on the Argonne FUSION cluster under the protection of FTI. Specifically, we observe that the higher the checkpoint level is, the higher checkpoint/recovery overhead is, while the checkpoint overhead does not change a lot with different scales on the same level. Moreover, the higher-level failure rate should be lower than lower-level failure rate, as analyzed in our previous work [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF].

The settings of the nine test cases are presented in Table 2, where C 1 and C 2 refer to checkpoint overheads (in seconds), R 1 and R 2 refer to restart overheads (in seconds), and λ 1 and λ 2 refer to the failure rates (number of failures per day) at level 1 and level 2, respectively. We conservatively overestimate recovery times and set them equal to checkpoint times, because this represents the worst-case scenario for the simulations. The failures are injected randomly in terms of the Poisson process, which is a standard approach [START_REF] Jin | Optimizing HPC Fault-Tolerant Environment: An Analytical Approach[END_REF], [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF]. That is, the failure intervals follow exponential distribution with the failure rates λ 1 and λ 2 at the two levels. We compare five methods. The first two methods (namely, Online-Pattern-OPT and Online-Interval-OPT) are online solutions without knowing fault-free job length. The third and fourth solutions are the offline solutions that require to know the job length beforehand. These two offline solutions are derived in Section 5. The fifth method (namely, Approximate-OPT) is the approach proposed in our previous work [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], and it is also an offline solution because it requires knowing the job length beforehand. In our previous work, Approximate-OPT was deemed the best solution compared with other state-of-the-art approaches such as Young's formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF]. In that work [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], an approximate optimal algorithm was proposed, by leveraging the fixed-point iteration method. Such a solution is subject to a few assumptions; for example, the failures are assumed not to occur during the checkpoint/restart period, and the final wall-clock length is not extended significantly compared with the fault-free productive length. As such, the output of the algorithm may deviate from the practical optimal solution, especially when the checkpoint/restart overhead is relatively large compared with the length (as is confirmed in our evaluation, shown later).

Evaluation Results

First of all, we compare the outputs (including interval-based model and pattern-based model) of the optimal online solution versus the optimal offline solution. We also compare these solutions to our previous work [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], as shown in Table 3.

From the table we observe two significant findings. We can see that our online optimal solution leads to the same results as the offline optimal solution does. In fact, the difference between the online optimal solution and offline optimal solution is tiny. For case 1, for instance, the offline optimal solution is derived as {w * opt =368.64474109273493, K * =3.5134717932162443}, while the online optimal solution is {w * opt =368.64474109270884, K * =3.513471793216452}. Hence, the optimality of our online solution (optimization without knowing the job length) is equal to that of the offline optimal solution (the optimization with knowledge of job length). However, we can also see that the three approaches (Online-Interval-OPT, Online-Pattern-OPT, and Approximate-OPT) output different checkpoint solutions. As for the Online-Pattern-OPT (i.e., {W * opt ,K * }), the level-2 checkpoint frequency has to be an integer number in practice, so its final value will be rounded to its nearest integer for checkpoint level 2. In the table, the K * column shows the rounded integer value (outside parentheses) and the original real number (inside parentheses) computed by our solution. Such a rounding operation will degrade the optimality of checkpoint solution, as will be shown later.

W opt2 =1215 W opt2 =1235 W opt2 =1255 W opt2 =1275 W opt2 =1295 W opt2 =1315 W opt2 =1335 W opt2 =1355 W opt2 =1375 (
W opt2 =693 W opt2 =713 W opt2 =733 W opt2 =753 W opt2 =773 W opt2 =793 W opt2 =813 W opt2 =833 W opt2 =853 (
W opt2 =631 W opt2 =651 W opt2 =671 W opt2 =691 W opt2 =711 W opt2 =731 W opt2 =751 W opt2 =771 W opt2 =791 (
W opt2 =406 W opt2 =426 W opt2 =446 W opt2 =466 W opt2 =486 W opt2 =506 W opt2 =526 W opt2 =546 W opt2 =566 (
W opt2 =239 W opt2 =259 W opt2 =279 W opt2 =299 W opt2 =319 W opt2 =339 W opt2 =359 W opt2 =379 W opt2 =399 (
W opt2 =420 W opt2 =440 W opt2 =460 W opt2 =480 W opt2 =500 W opt2 =520 W opt2 =540 W opt2 =560 W opt2 =580 (
W opt2 =333 W opt2 =353 W opt2 =373 W opt2 =393 W opt2 =413 W opt2 =433 W opt2 =453 W opt2 =473 W opt2 =493 (
W opt2 =370 W opt2 =390 W opt2 =410 W opt2 =430 W opt2 =450 W opt2 =470 W opt2 =490 W opt2 =510 W opt2 =530 (h) Case 8
W opt2 =370 W opt2 =390 W opt2 =410 W opt2 =430 W opt2 =450 W opt2 =470 W opt2 =490 W opt2 =510 W opt2 =530 (i) Case 9
We validate the optimality of our proposed solution (Online-Interval-OPT) to see how close it is to the best experimental solution 5 . In this evaluation, different solutions are generated by combining different level-1 checkpoint intervals and level-2 checkpoint intervals, to evaluate whether our solution indeed leads to the best performance (with minimum wall-clock length). For each test case and each checkpoint solution, we run the simulation 1,000 times each with different random arrival time points of failure events but with the same settings on job length (i.e., productive time), failure rate, checkpoint overhead and recovery overhead. We then compute the mean values for every portion of times consumed in the execution, including productive time, checkpoint overhead, restart overhead, and rollback time.

5. The best experimental solution is obtained by brute-force searching/traversing different combinations of checkpoint intervals (regarding the two levels/types of failures), which are increased gradually by a tiny increment (5 seconds in experiments) from a small checkpoint interval (20 seconds).

Figure 3 presents the wall-clock times consumed with different checkpoint intervals set in various failure cases, where the wall-clock time is the sum of the above four time portions. To make Figure 3 easy to read, we highlight our solutions by solid curves. We present the results only for the solutions with a particular checkpoint interval range, because the solutions with other checkpoint intervals lead to much higher wall-clock lengths (worse results). In Figure 3 we see that the wall-clock time with our optimal checkpoint solution is really (or fairly close to) the minimum wall-clock length among all other solutions. In Table 4, we present the wall-clock time of our optimal checkpoint solution versus that of the best experimental solution selected by traversing all the checkpoint intervals. It shows that the wall-clock time difference between our solution and the best experimental solution is always within 0.7% for the first seven test cases, and it can be limited within 7.7% for the other two test cases in which the checkpoint/recovery overheads are large (up to 5 minutes for each level-2 checkpoint) and failure rates are also high (about one failure every three minutes on average). The relatively high difference between our solution and the best experimental solution in the two test cases is due to the assumption that no faults occur during the recovery period, whereas the failure-overrecovery would occur often in the last two test cases. For these two test cases, we later show that our optimal solution achieves significant performance gains compared with the results of the previous optimized solution proposed in [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF]. In Figure 4, we compare the total wall-clock time and different portions of time consumed by our proposed online checkpoint solution (Interval-OPT) versus that of the other two optimized approaches (Pattern-OPT and Approximate-OPT). For the first six test cases, we observe that the three solutions lead to very close wall-clock times: the difference is smaller than 2% in wall-clock length. For the last three test cases, our optimal online solution (Interval-OPT) outperforms the other two solutions significantly. In test case 8, the wall-clock time of our solution is reduced by about 11% and 25.3%, compared with Pattern-OPT solution and Approximate-OPT [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], respectively. In test case 9, the wall-clock time can be reduced by 12.5% and 23.6%, respectively. The key reason for the poor performance of Pattern-OPT is that the level-2 checkpoint frequencies adopted are always conducted strictly based on the number (i.e., rounded integer) of level-1 checkpoints, definitely deviating from the best experimental solution. The key reason for the degraded performance of the Approximate-OPT algorithm proposed in [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF] is that it is subject to a few significant assumptions such as the relatively low failure rates and notmuch-extended wall-clock length in comparison to the fault-free productive time. In particular, when the failure rate is not high or the checkpoint/restart overhead is relatively low (i.e., the first 6 test cases), the performance of our previous Approximate-OPT algo-rithm is close to that of the optimal solution proposed in this paper. However, with high failure rates and heavy checkpoint/restart overhead, the performance degradation of Approximate-OPT is clear, as confirmed by the last three test cases in Figure 4.

RELATED WORK

The optimization issue of finding the optimal checkpoint intervals for HPC applications has been studied for decades. In 1973, Young [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF] proposed a formula to optimize the checkpoint interval with the first-order approximation. In 2006, Daly [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] extended his solution to a higher-order approximation and took recovery overhead into account. These two works consist of the fundamental and classic optimization strategies. Thereafter, more checkpoint solutions were proposed to resolve more complicated cases. Chen and Ren [START_REF] Chen | Adaptive optimal checkpoint interval and its impact on system's overall quality in soft real-time applications[END_REF], for example, proposed an optimal checkpoint solution from the perspective of the overall system performance. Liu et al. [START_REF] Liu | An optimal checkpoint/restart model for a large scale high performance computing system[END_REF] proposed a reliability-aware optimal checkpoint/restart strategy for a large-scale HPC system. Compared with other checkpoint/restart optimization work, their model can deal with a varying checkpoint interval with different failure distributions. Jin et al. [START_REF] Jin | Optimizing HPC Fault-Tolerant Environment: An Analytical Approach[END_REF] optimized the checkpoint intervals and execution scales simultaneously. The optimal checkpoint interval is derived by assuming exponential failure distribution and using first-order approximation. They derived the expected wall-clock length regarding various overheads and rollback loss, but they did not prove that it is a convex function with respect to the variables before using Newton's method to approximate the optimal solution. In that situation, the converged result may not be globally optimized; and even worse, the algorithm may not converge given inappropriate initial values. Most important, all these works are restricted to the single-level checkpoint model, which cannot work efficiently in an extreme-scale HPC environment because of inevitable huge checkpoint/recovery overheads.

Multi-level checkpoint model was introduced in order to address the limitation of single-level checkpoint strategies. In 1995, Vaidya et al. proposed a two-level checkpoint model [START_REF] Vaidya | A case for two-level distributed recovery schemes[END_REF], which is similar to the pattern-based checkpoint model in our paper. By numerical analysis, they confirmed that multi-level checkpoint model is able to improve the performance compared to singlelevel checkpoint model. To this end, some toolkits (such as [START_REF] Moody | Design, modeling, and evaluation of a scalable multi-level checkpointing system[END_REF], [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF]) already support multilevel checkpoint protocols. That is, the users are allowed to set various types of checkpoints with different overheads for protecting have been developed to support applications against different types of failures. However, all of these works did not investigate what the checkpoint intervals are supposed to be for a particular HPC application.

How to optimize the checkpoint intervals for multilevel checkpointing has been studied in the recent years. Hakkarinen and Chen [START_REF] Hakkarinen | Multilevel diskless checkpointing[END_REF] proposed a multilevel diskless checkpoint model. This model is similar to our proposed multi-level checkpoint model, although it focuses on a diskless checkpoint environment. Because of the complexity of the problem, the authors mainly derived the representation formula for the expected wall-clock time regarding different checkpoint overheads and failure rates and narrowed the value range for the optimal checkpoint intervals but did not provide exact optimal checkpoint intervals for the multilevel checkpoint model. Moreover, they evaluated their solutions by a hypothetical simulation based on the estimated wall-clock length with numerically summed overheads and rollback loss. In contrast, we carefully implement a fundamental simulation testbed, and we confirm the accuracy of our simulation using a real-world HPC application (called Heat Distribution). In our previous work, we proposed an approximate solution to optimize the checkpoint intervals for different checkpoint levels [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], [START_REF] Di | Optimization of a multilevel checkpoint model with uncertain execution scales[END_REF]. In [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], we formulated the multilevel checkpoint model based on the existing multilevel checkpoint toolkits [START_REF] Moody | Design, modeling, and evaluation of a scalable multi-level checkpointing system[END_REF], [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF]. In that work, to make the problem tractable, we made two important assumptions: (1) failures do not occur during the checkpoint period, and (2) the total wall-clock time is not expanded significantly compared with the fault-free execution length. With such two assumptions, an iterative algorithm was proposed to get an approximate optimal solution based on convex-optimization theory. An extended work was proposed in [START_REF] Di | Optimization of a multilevel checkpoint model with uncertain execution scales[END_REF] by taking into account the environment with uncertain execution scales. In comparison to such previous works, the new solution proposed in this paper completely removes the two assumptions mentioned above and provides exact solutions instead of first-order approximations or fixed-point iteration approximations. Specifically, we reconstruct the multilevel checkpoint model and derive a new optimal checkpoint solution that is much closer to the best experimental solution. Our elaborate simulation shows that the new optimal solution outperforms the previous solutions by up to 25.3% in some cases.

CONCLUSION AND FUTURE WORK

In this paper, we optimize the checkpoint solution based on a two-level checkpoint model, in which the system is allowed to set two different types of checkpoints for protecting the HPC applications against two types of failures. Optimizing the two-level checkpoint solution is difficult in that each type of failure follows its own distribution and the overall performance is synthetically determined by the two types of checkpoints and failures. To address the issue, we first formulate the research to be a patternbased checkpoint model and then extend it to the interval-based checkpoint model. We derive the optimal two-level solutions based on both online scheduling and offline scheduling and prove that the optimal solution must adopt the equal-size checkpoint intervals on particular checkpoint levels. The key theoretical conclusion about how to conduct the optimal checkpoint strategies in practice is also discussed in Section 6. To evaluate the proposed solutions, we build a simulation testbed and confirm its accuracy by running a real-world MPI program on a real cluster with 1,024 cores. We highlight three significant findings from our evaluation.

• Our online optimal solution leads to the same results as the offline optimal solution, thus confirming the optimality of our online solution. • Our online optimal solution always outperforms other optimized solutions, and improves the performance significantly in some cases. Specifically, the wall-clock time under the new solution can be reduced by up to 25.3% compared with the results of other state-of-the-art methods. • The optimal solution derived in this paper is close to the best experimental solution (obtained via brute-force searching), with the difference (in the wall-clock time) less than 1% in most cases.

A possible future project could involve investigating how to optimize the two-level checkpoint solution by removing the assumption that the failures may not occur during the recovery period. In addition, we plan to evaluate the optimal solutions proposed using additional real-world applications.

Fig. 1 :

 1 Fig. 1: Pattern including K chunks of respective sizes w 1 , ..., w K .

)e 2

 2 and F 3 = e 1 + e 2 + e 3 + L(e 1 e 2 + e 1 e 3 + e 2 e 3 ) + L 2 (e 1 e 2 e 3 ) = F 2 + (1 + LF 2 )e 3 . Proposition 1.

Fig. 2 :

 2 Fig. 2: Confirming the Effectiveness of the Simulator

Fig. 3 :

 3 Fig. 3: Verification of the optimality of our online interval-OPT solution.

Fig. 4 :

 4 Fig. 4: Performance comparison of our solution with others.

TABLE 1 :

 1 Notations.

	PATTERN(K, W )	Pattern of size W , with K chunks
	wi	Size of i-th chunk in pattern
	λ1, C1, R1	Fault-rate, checkpoint time, recovery time for level-1 faults
	λ2, C2, R2	Fault-rate, checkpoint time, recovery time for level-2 faults
	λ	Cumulated rate of faults: λ = λ1 + λ2
	L	Fraction of level-2 faults: L = λ2/λ
	R	Averaged overhead: R = D + (1 + λ1R1 + λ2R2)/λ
	wopt (K)

TABLE 2 :

 2 Simulation setting of 9 test cases.

		C1	C2	R1	R2	λ1	λ2	Length
	Case 1	20	50	20	50	24	4	86400
	Case 2	20	50	20	50	50	10	86400
	Case 3	20	100	20	100	100	20	86400
	Case 4	10	40	10	40	100	20	86400
	Case 5	10	40	10	40	200	40	86400
	Case 6	10	100	10	100	200	40	43200
	Case 7	40	200	40	200	300	60	21600
	Case 8	50	300	50	300	400	60	21600
	Case 9	50	300	50	300	400	60	10800

TABLE 3 :

 3 The optimized solutions under different approaches.

			Online Solution				Offline Solution
		Pattern-OPT	Interval-OPT	Pattern-OPT	Interval-OPT	Approximate-OPT
		w * opt	K *	w * opt	w * opt2	w * opt	K *	w * opt	w * opt2	w * opt	w * opt2
	Case 1	368.6	4 (3.51)	368.6	1295.2	368.6	4 (3.51)	368.6	1295.2	385.7	1433
	Case 2	252.7	3 (3.06)	252.7	773	252.7	3 (3.06)	252.7	773	269.6	896.9
	Case 3	175.9	4 (4.04)	175.9	711.3	175.9	4 (4.04)	175.9	711.3	195.2	885
	Case 4	126.4	4 (3.85)	126.4	486.1	126.4	4 (3.85)	126.4	486.1	135.7	567.3
	Case 5	88.0	4 (3.63)	88.0	319	88.0	4 (3.63)	88.0	319	97.1	395.8
	Case 6	88.0	6 (5.68)	88.0	499.9	88.0	6 (5.68)	88.0	499.9	99.5	626.5
	Case 7	134.4	3 (3.07)	134.4	412.7	134.4	3 (3.07)	134.4	412.7	168.8	682.4
	Case 8	124.1	4 (3.62)	124.1	449.5	124.1	4 (3.62)	124.1	449.5	166.5	815.1
	Case 9	124.1	4 (3.62)	124.1	449.5	124.1	4 (3.62)	124.1	449.5	166.5	815.1

TABLE 4 :

 4 Our derived optimal solution vs. best experimental solution (in seconds).

		{w * opt ,w * opt2 }	Best Experimental Solution	Difference
	Case 1	104024	103788	0.23%
	Case 2	115220	114890	0.28%
	Case 3	144883	144461	0.29%
	Case 4	119451	119144	0.26%
	Case 5	140029	139801	0.16%
	Case 6	84884	84517	0.43%
	Case 7	126407	125523	0.7%
	Case 8	389354	362493	6.9%
	Case 9	190764	175701	7.7%

Job length here refers to the fault-free productive length. For simplicity of presentation, we will interchangeably use different terms (such as work size and productive length) in the following text. They all mean job length.

ACKNOWLEDGMENTS

We would like to thank the reviewers for their comments and suggestions, which greatly helped improve the final version of the paper. This research was funded in part by the European project SCoRPiO, and by the PIA ELCI project. Yves Robert is with Institut Universitaire de France. This research was also based upon work supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research Program, under Contract DE-AC02-06CH11357; and by the INRIA-Illinois-ANL-BSC-JSC-Riken Joint Laboratory on Extreme Scale Computing, and Center for Exascale Simulation of Advanced Reactors (CESAR) at Argonne.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory ("Argonne"). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide to the public, and perform publicly and display publicly, by or on behalf of the Government.

Sheng Di Sheng Di received his Ph.D. degree from the University of Hong Kong in 2011. He is currently a postdoctoral researcher at Argonne National Laboratory. Dr. Di's research interest involves resilience on high-performance computing (such as silent data corruption, optimization checkpoint model, and in-situ data compression) and broad research topics on cloud computing (including optimization of resource allocation, cloud network topology, and prediction of cloud workload/hostload). He is working on multiple HPC projects, such as detection of silent data corruption, characterization of failures and faults for HPC systems, and optimization of mutlilevel checkpoint models. He is the author of 13 papers published in international journals and 33 papers published at international conferences. He served as programming committee member 14 times for different conferences, and served as an external conference/journal reviewer over 50 times. Contact him at sdi1@anl.gov.

Yves Robert


Iñigo Bediaga

Christian Ramirez
email: cramirez@ideko.es

Kepa Huerta
email: khuerta@ekin.es

IN-PROCESS TOOL WEAR MONITORING IN INTERNAL BROACHING

Keywords: Broaching, monitoring, tool wear

In this work, an -in process‖ tool wear prediction system for broaching process is presented based on a cutting force analysis. The aim is to link the changes in the measured force with the tool wear. Once the broaching cycles have been performed, different strategies have been analyzed including the evolution of the maximum peak force and the evolution in the used broaching energy per section. Finally, the results of the analysis and the main conclusions are presented.

INTRODUCTION

Machining is an essential process in most industrial manufacturing procedures, because it can create a variety of part shapes with tight tolerances over a wide range of materials.

However, due to the aggressive nature of metal cutting, the presence of chatter and/or excessive wear on the tool can dramatically impact on the product quality and productivity [START_REF] Altintas | Chatter stability of metal cutting and grinding‖[END_REF][START_REF] Zatarain | Estudio comparativo de los modelos matemáticos de chatter en fresado: monofrecuencia, multifrecuencia y simulación en el tiempo‖[END_REF][START_REF] Mancisidor | Receptance coupling for tool point dynamic prediction by fixed boundaries approach‖[END_REF].

There has been a significant amount of research in monitoring of machining processes such as turning, milling and drilling [START_REF] Dimla | Sensor signals for tool-wear monitoring in metal cutting operations-a review of methods‖[END_REF]. However, less effort has been dedicated to the monitoring of broaching and similar special operations. Broaching is a high productivity operation that permits mass production of parts with complex geometries. Considering the significant technical and economic implications of the broaching processes in strategic industrial sectors like aerospace, the aim of this work is to develop a system to monitor the broaching process to maintain the quality of the parts.

Broaching is a shaping process that is used to manufacture massively parts with complex internal and/or external profiles, which are not possible to machine using more flexible metal cutting methods. One of the main challenges is to monitor tool wear, taking into account that the condition of the tool has crucial importance for the surface quality of the machined components. In broaching, the cost of the tool named broach can be really high compared with other processes. Therefore, it is important to maximize the productivity reducing the tool cost by programing properly the change of the broach, but avoiding tool malfunctions that can damage severely the final component.

Nowadays, the tool is often inspected visually by an expert that determines the situation. The final goal of this research is to implement a continuous and automatic monitoring system to define the condition of the tool including the wear.

STATE OF THE ART

Process monitoring is based on the utilization of observable and measurable process and product parameters to evaluate the state of a system including different mechanical elements [START_REF] Fernandez | Feasibility study on diagnostic methods for detection of bearing faults at an early stage‖[END_REF]Bediaga et al, 2013]. These parameters include product features such as; surface roughness or geometric dimensions and process variables, such as, machine speeds and forces. There is a need to adapt and refine current statistical process monitoring and quality control methods to take full advantage of the available data. In the last years, the monitoring has increased its impact thanks to the introduction of new embedded sensors and the tremendous enhancement in the data storage possibilities. This revolution has open new possibilities for big data treatment by means of smart analytics.

However, broaching processes are not a common target for this new techniques and the state of the art in the monitoring processes for that operations is limited. First of all, it is important to consider the articles related to broaching process modeling. The geometry of the tool is especially important in broaching, therefore, the broaching process models have been The objective of this work is to perform the first step to develop an industrial broaching monitoring system based on smart analytics. In this first step, the life of a broaching tool has been monitored during the life of the tool in industrial conditions by means of different sensors including acceleration and force transducers.

DESIGN AND VERIFICATION OF THE EXPERIMENTAL SET-UP

An experimental set-up was prepared to measure the force and the acceleration during the internal broaching process of 2000 parts to capture the evolution of the wear of the broaching tool.

Experimental set-up

Taking into account the big amount of pieces to be machined, a test layout was designed in order to facilitate the tests. The set-up consists on a Kistler dynamometric plate, a signal amplifier for the plate, an acquisition card and an embedded PC.

The layout is shown in Figure 1

Broaching parameters

The selected tool is a 30-spline broach with 58 teeth with a pitch of 8mm. In Figure 2, the basic geometry of the cutting flutes of the broach is described. The selected tool has a tooth thickness of 2.5mm, a depth of tooth of 2.7mm, a radius of 1.6mm, a rake angle of 15 degrees, and a clearance angle that varies from 0 degree to 2 degrees. The depth of cut per tooth is 0.043mm, and it is maintained constant for each tooth. In this case, the broach has not any coating to accelerate the wear process on it. The cutting speed is fixed on 50m/min without any coolant.

The target parts are made of C45 steel having a height of 40mm and a diameter of 50mm (See Figure 3). Initially, the raw signal (See Figure 4) is acquired with a sampling frequency of 4096 samples per second. If the time domain signal of the axial broaching force is observed, the force measurement looks noisy but it is possible to appreciate a clear trend during the machining process. However, one of the objectives of the monitoring system is to identify a failure in the broach and to determine the damaged flute. In order to avoid the high frequency components a low pass filter can be applied considering a frequency higher than the tooth passing frequency (see Figure 5). In this case the force signal becomes clearer, and it is easy to distinguishing the force produced by each tooth checking the time domain filtered signal (see Figure 6).

With the objective of analyzing separately the different parts of the broaching process, the force signal was divided in two portions (See Figure 6). The initial portion starts with the force component of the first teeth, while the final portion includes the last teeth. Then, a spectral analysis of the unfiltered or raw force signal was carried out by means of a Fast Fourier Transformation (FTT). The spectrum of the first portion shows (see Figure 7) the appearance of the maximum peak at 792 Hz, which is not a harmonic of the tooth passing frequency. However, the spectral analysis of the last time portion shows that the previous peak disappears. This peak could exist due to different causes, and with the aim of finding the real cause, a dynamic analysis was performed. This analysis consisted in taking frequency response functions (FRF) at different points of the system.

After analyzing all the system, it was found out that the dynamometer plate mounting block (See Figure 8) had a natural frequency close to the peak shown in the spectrum, as is shown in the FRF (See Figure 8). This FRF clearly shows a natural frequency in the same range than the peak mentioned before. Several tests at different cutting speeds were carried out appearing always a peak close to 800 Hz. Therefore, the highest peak was concluded to be chatter, and the source of the chatter problem comes from the flexibility of the torsional mode of the dynamometer plate block.

Finally, the dynamometric plate overhang was reduced and stiffen, getting almost disappeared the torsional mode as result and removing chatter vibrations from the broaching cycle.

BROACHING FORCE EVOLUTION CONSIDERING WEAR

A procedure has been developed to monitor tool wear by means of different indicators in every broaching cycle. The filtered time domain signal is divided in three sections and different thresholds are defined. The main parameters are described in Figure 9. The analysis of the complete life of the broach requires the monitoring of a large number of parts. Therefore, it is important to automatize the measuring process of the broaching process. The designed program incorporates several thresholds and time lapses to ensure that force signal is adequate, and is completely acquired. There is a threshold fixed at the beginning of the signal warning when the broaching force starts. When this threshold is achieved a pre trigger time lapse is also activated that permits to capture a portion of the signal before the trigger.

On the other hand, there is another threshold at the last part of the signal that is in charge of defining the end of the measurement. This threshold is used to corroborate that the machining is finished, and to fix post broaching time lapse that adds final data to consider the entire signal.

More than 2000 broaching tests were carried out, and different indicators were obtained for each part including the maximum force peak, the energy consumption and the average force value per tooth. Finally the real tool wear was checked visually every 500 test.

Evolution of the maximum force peak:

The first indicator is the maximum value of the filtered signal of the measured broaching force. The maximum force peaks fluctuate from one measurement to another, but it shows a clear trend. Logically, the maximum force rises with the wear of the flutes. In order to obtain the mentioned trend a linear regression is applied.

The trend of the maximum force peak values is increasing from an initial value of 27230 N to a final value of 32520 N (See Figure 2). The mathematical principle used to calculate the energy (N*m or J) consists in multiplying the force (N) by the cutting speed (m/min) and by the sampling period (s).

         j N K K K s i F F F f e 1 1

120

(Equation 1)

Where, e : energy in the i th test number [J],

f : cutting speed [m/min], s F : sampling frequency [Hz],

i: test number, K : sample, j N : number of samples in the j th portion.

The broaching tool is composed by a high number of flutes, and it is important to identify where the wear is growing more. Therefore, the whole process is divided in three sections to In this particular case, the trend of the energy value in the first portion is increasing from an initial value of 3023J to a final value of 3334J. The increase can be established in a 9.3% (See Figure 11(a)). In the second portion the energy increase is lower during these 2000 broaching operations reaching a 3.4% of increase from 2537 Nm to 2626 Nm (See Figure 11(b)). The trend for the third period maintains almost constant, even is slightly decreasing (See Figure 11(c)). According to this indicator, wear grows basically in the first section of the broaching tool.

Evolution of mean value of each tooth:

This analysis gives the chance to know the situation of each tooth. The previous indicator showed that the first teeth suffer more wear, and therefore, this analysis was focused on the first flutes (See Figure 12). The analysis of the average value of the force determined that the bigger difference in absolute values was appreciated in the fifth tooth (See Table 1).

CONCLUSIONS

This paper presents a force analysis of the broaching process during whole life of the broaching tool. This analysis has been made using different indicators including the evolution of the maximum peak force, the evolution of the consumed energy and the average value of the force per flute. A general conclusion of the research is that the axial force peak increases with the wear, and that the wear of the broach does not evolve constantly through all the surface of the broach because of a non-linear distribution of the broaching forces that are higher on the first 10 teeth.

The first section of the broach shows clearly an increasing trend for the consumed energy.

This fact is validated by the conventional analysis performed taking pictures of every 500 tests.

  applied for tool design and optimization [Konkurk and Budak, 2004; Ozelkan et al, 2010]. Ozlu et al. (2010) simulated broaching operations, proposing analytical approaches for both mathematical and geometrical modeling that result in very fast solution times and accurate predictions. Mandrile et al. (2012) implemented a methodology to simulate cutting forces in fir tree broaching. Recently, Chopaldi and Kuttan (2014) implemented a mechanistic model for the cutting force in rectangular slotting operation, validating this model with experimental test. An energy based analysis of broaching operation was performed by Kishawy et al. (2012), where the cutting force simulation was applied to quantify the contribution of different power components of the cutting system. These efforts reached the market when MAXIMA launched a broaching process simulator named BOSS ® . Similarly, Klocke et al. (2015) implemented a mathematical model employing successfully the Kienzle-force model to get the cutting force of a non-tested broach [Klocke et al 2015; Munoa et al 2006]. On the other hand, there are some articles related to the experimental monitoring of the broaching process: For instance, Axinte and Gindi (2003) showed that acoustic emission, vibration and cutting force signals were sensitive to tool condition, and a correlation could be made between the broaching tool situation and sensor signals. This paper reveals the high affectivity of cutting forces detecting the tool conditions like chipped edges, broken flute and worn tool. The same research group reports an attempt to use a combination of acoustic emission (AE) and force signals to detect surface anomalies broaching parts made in TI6-4 alloy [Axinte et al, 2005]. They emphasized that acceleration signals, due to their low sensitivity to material deformation, are not capable to detect surface anomalies, relying this purpose on AE and cutting force signals. The same author continued working in this area, studying the use of probabilistic neural networks for automated classification of tool malfunctions. The aim of this work was to define a strategy to process the output signals acquired during broaching in order to detect and clarify the condition of the teeth of broaching tools [Axinte, 2006]. In this work, Axinte presents a method to automatically detect and clarify tool malfunctions based on specific pattern of the force. He also carried out the monitoring of experimental broaching tests, analyzing simultaneously the acceleration and cutting force signals reporting for the first time aspects related to dynamics of broaching when complex geometries were generated [Axinte, 2007]. Bound and Gilti (2006) monitored broaching process and analyzed the measurement of acoustic emission, cutting force and vibration signal for the identification of surface anomalies and tool condition. Finally, Rathinam (2013) implemented also a methodology for in process monitoring.

  where a Kistler dynamometric plate is connected with the signal amplifier and the acquisition card. The acquisition card records the measured force and sends it to the PC. The start and stop of the measurement is established by the PLC of the machine, which is connected to the PC and signal amplifier (triggering signal).

Figure 1 .

 1 Figure 1. Test layout for multiple signal recording.

Figure 2 .

 2 Figure 2. Geometry of the broach.

Figure 3 .

 3 Figure 3. Test workpiece

Figure 4 .

 4 Figure 4. Original force signal Figure 5. Filtered force signal (180Hz).

Figure 6 .Figure 7 .

 67 Figure 6. Kistler plate measured force divided in two stretches

Figure 8 .

 8 Figure 8. Dynamic analysis of dynamometric plate (excitation point green arrow, and measurement point indicated by a blue circle)

Figure 9 .

 9 Figure 9. Presentation of the parameters and indicators used for describing the broaching force.

Figure 2 . 5 . 2

 252 Figure 2. Evolution of the force peak

Figure 11 .

 11 Figure 11. Evolution of the energy expended in the first (a), second (b) and third (c) sections.

Figure 12 . 5 . 3

 1253 Figure 12. Average value per tooth

Figure 13 .

 13 Figure 13. Pictures of the 5 th , 20 th , 35 th and 55 th for every 500 tests

Table 1 .

 1 Difference per tooth ,51 33,30 25,73 30,20 9,41 14,47 12,83 6,16 6,28 5,83

	Nº Tooth	1º	2º	3º	4º	5º	6º	7º	8º	9º	10º	11º
	Difference (N)	2540 2480 5230 8580 8870 2820 4210 3520 1530 1590 1370
	Difference (%)	49,29 29

ACKNOWLEDGEMENT

The authors would like to acknowledge the support of the European Commission under research activities carried out inside FoF.NMP.2013-8 EASE-R3 project (grant number 608771).