A Segmented Model with Upside-Down Bathtub Shaped Failure Intensity

Upside-Down 욕조 곡선 형태의 고장 강도를 가지는 세분화 모형

Abstract

In this study, a segmented model with Upside-Down bathtub shaped failure intensity for a repairable system are proposed under the assumption that the occurrences of the failures of a repairable system follow the Non-Homogeneous Poisson Process. The proposed segmented model is the compound model of S-PLP and LIP (Segmented Power Law Process and Logistic Intensity Process), that fits the separate failure intensity functions on each segment of time interval. The maximum likelihood estimation is used for estimating the parameters of the S-PLP and LIP model. The case study of system A shows that the S-PLP and LIP model fits better than the other models when compared by AICc (Akaike Information Criterion corrected) and MSE (Mean Squared Error). And it also implies that the S-PLP and LIP model can be useful for explaining the failure intensities of similar systems.

1. 서론

무기체계를 포함하여 수리 가능한 체계 (Repairable System)의 고장 강도 함수 (Failure Intensity Function)는 일반적으로 단조 (Monotone) 함수의 형태를 보인다. 그러나 개발시험 단계의 특정 기동화력 무기체계나 항공기의 에어컨 [1], 증기 배기장치 [2]와 같은 체계의 경우에는 체계의 고장 강도 함수가 비단조 (Non-Monotone)형태인 Upside-Down 욕조 곡선 형태로 나타나는 경우가 있다.

현재까지 Upside-Down 욕조 곡선 형태의 고장 강도에 대한 연구는 주로, NHPP (Non-Homogeneous Poisson Process)를 가정한 단일고장 강도 (Single Failure Intensity) 모형이나 복합 고장 강도 (Compound Failure Intensity) 모형 중 세분화 (Segmented) 모형을 대상으로 진행되었다.

Upside-Down 욕조 곡선 형태의 고장 강도 모형은 소프트웨어 분야에서 많이 활용되었으며 S-Shaped를 가지는 단일 고장 강도 모형을 고려한 연구들이 주로 수행되었다. Yamada et al. [3]과 Goel [4]이 각각 제시한 Yamada 모형과 GGO (Generalized Goel Okumoto) 모형이 대표적이다. 수리 가능한 체계에 관한 연구로는 최근 Jiang [1]이 Guo et al. [5]의 모형을 변형한 NN (New NHPP) 모형을 단일 고장 강도 모형으로 제시하였다. 그리고 Du et al. [2]은 고장 강도의 기울기가 변하는 지점을 기준으로 두 구간을 세분화하여 구간별로 PLP (Power Law Process) 모형 [6]을 추정하는 Two Segment 모형을 통해 단조 증가, 감소, 욕조 곡선, Upside-Down 욕조 곡선을 모두 적합할 수 있는 복합 고장 강도 모형을 제안하였다.

이 연구에서는 일부 국내 무기체계의 시험평가 단계에서 나타난 Upside-Down 욕조 곡선 형태의 고장 강도 함수를 설명하기에 적합한 S-PLP and LIP (Segmented Power Law Process and Logistic Intensity Process) 복합 모형을 제안하고, 최우추정법 (Maximum Likelihood Estimation)을 이용하여 S-PLP and LIP 모형의 모수를 추정하는 방법을 제시하였다. 그리고 국내 기동화력 무기체계 중 하나인 A체계의 사례를 통해 이 연구에서 제안한 S-PLP and LIP 모형의 적용 가능성을 확인하였다. A체계의 사례연구에서 S-PLP and LIP 모형의 분석 결과와 더불어 기존에 많이 활용되는 PLP 모형, 소프트웨어 신뢰성 분야의 GGO 모형과 Yamada 모형, 최근 연구에서 제시된 NN 모형과 Two Segment 모형의 분석 결과를 서로 비교 분석하였다.

2. 기존 고장 강도 모형

수리 가능한 체계가 주어진 [0, T] 동안 고장이 발생하는 경우를 가정했을 때, 고장 강도 함수 (Failure Intensity Function)를 λ(t), 누적 강도 함수 (Cumulative Intensity Function)를 Λ(T) = ∫^T₀λ(x)dx 라고 하면 기존 연구에서 다루어진 PLP [6], GGO [4], Yamada [3], NN [1], Two Segment [2] 모형의 고장 강도 함수와 누적 강도 함수 식은 Table 1과 같다.

Table 1. Failure Intensity and Cumulative Intensity Functions of Failure Intensity Models

PLP 모형 [6]에서 λ는 척도 모수, β는 형상 모수이며, 고장 강도 함수의 형태는 β>1인 경우 단조 증가, β<1인 경우 단조 감소, β=1인 경우 상수로 나타나게 된다. Yamada 모형 [3]은 시험의 시작 후 고장 발생 비율이 특정한 시점까지 증가하다가 지수 형태로 감소하는 형태를 가진다. GGO 모형 [4]은 Goel이 Goel-Okumoto 모형 [7]을 수정하여 소프트웨어 신뢰성 시험 초기에 고장 발생 비율이 증가하였다가 감소하는 형태를 모형화한 것이며 척도 모수 a와 b, 형상 모수 c를 갖는다. NN 모형 [1]은 Jiang이 Guo et al. [5]의 모형을 변형한 모형이며, β>1일 때 Guo et al. 의 모형보다 고장 강도 함수의 감소 구간에서 기울기가 완만하게 감소하는 것이 특징이다. Du et al. [2]의 Two Segment 모형은 t₀를 기준으로 구간을 구분하여 구간별로 각각 고려한 PLP 모형의 모수를 추정하고 두 구간의 연속성을 만족시키기 위해 조정 계수를 추가한 모형으로 단조 증가, 감소, 욕조 곡선, Upside-Down 욕조 곡선 형태의 고장 강도를 모두 모형화할 수 있다. 일반적으로 Two Segment 모형은 t₀를 기준으로 구간을 둘로 구분(Segment) 하고 구간별로 동일한 또는 서로 다른 고장 강도 모형을 적합시킨 모형으로 단조 증가, 단조 감소, 욕조 곡선 형태 및 Upside-Down 욕조 곡선 형태를 모두 나타낼 수 있는 것이 특징이다.

3. S-PLP and LIP (Segmented Power Law Process and Logistic Intensity Process) 모형

S-PLP and LIP 모형의 고장 강도 함수는 Upside-Down 욕조 곡선의 변화점 T_c를 기준으로 구분된 각 구간에 대해 각각 PLP 모형과 Logistic 함수 형태의 고장 강도 모형이 적합 된 형태를 갖는다. 이 때 S-PLP and LIP 모형의 고장 강도 함수 λ(t)는 식 (1), 누적 강도 함수 Λ(T)는 식 (2)와 같다. S-PLP and LIP 모형에서 첫 번째 구간 (0<t≤T_c)의 고장 강도는 증가하는 형태로 나타나므로 PLP 모형의 β는 1보다 커야 한다. Du et al. [2]이 제안한 Two Segment 모형의 경우 두 구간의 고장 강도 함수가 PLP 모형을 따르므로 두 번째 구간에서 PLP 모형의 특성에 따라 아래로 볼록한 (Convex) 형태로 고장 강도가 감소하게 된다. 그러나 일반적으로 수리 가능한 체계에 대한 Upside-Down 욕조 곡선 형태의 고장 강도는 두 번째 구간에서 완만하게 감소하는 경우가 많으므로, 이런 경우 두 번째 구간의 고장 강도 함수 형태는 PLP 모형보다 Logistic 함수 형태 (a/1+e^bt))가 더 적합할 수 있다.

\(\begin{align}\lambda(t)=\left\{\begin{array}{l}\lambda \beta t^{\beta-1}, 0<t \leq T_{c} \\ \frac{a}{1+e^{b t}}, T_{c}<t<\infty\end{array}\right.\end{align}\),       (1)

( λ > 0, β > 1, a, b > 0, t > 0 )

\(\begin{align}\begin{aligned} \Lambda(t) & =\int_{0}^{t} \lambda(x) d x \\ & =\left\{\begin{array}{cl}\lambda t^{\beta} & , 0<t \leq T_{c} \\ \frac{-a\left[b t-\ln \left(e^{b t}+1\right)\right]}{b}, & T_{c}<t<\infty\end{array}\right.\end{aligned}\end{align}\)       (2)

S-PLP and LIP 모형에서 모수 값의 변화에 따른 고장 강도는 Fig. 1과 같다. PLP 모형에서 λ는 척도 모수이고 β는 형상 모수이므로 첫 번째 구간의 고장 강도는 λ에 비례하며, β가 1에 가까울수록 첫 번째 구간의 곡선 기울기가 증가하고 β가 커질수록 아래로 볼록한 (Convex) 형태를 가진다. 두 번째 구간에서 a는 척도 모수이고 b는 형상 모수이므로, 고장 강도는 a에 비례하며, b는 작을수록 두 번째 구간의 감소 곡선 기울기가 완만해지고 b가 커질수록 기울기가 커지게 된다.

Fig. 1 Behavior of S-PLP and LIP Model

관측 중단 시간 T까지 발생한 n개의 고장 시간이 0 < t₁ < t₂ <⋯< t_τ < T_c <⋯< t_n < T 일 때 S-PLP and LIP 모형의 우도 함수는 식 (3), 대수우도 함수는 식 (4)와 같다.

\(\begin{align}\begin{aligned} L= & \prod_{i=1}^{\tau} \lambda \beta t_{i}^{\beta-1} \times e^{-\lambda T_{e}^{\beta}} \\ & \times \prod_{i=\tau+1}^{n}\left(\frac{a}{1+e^{b t_{i}}}\right) \times e^{\frac{a\left[b T-\ln \left(e^{b T}+1\right)\right]}{b}}\end{aligned}\end{align}\)       (3)

\(\begin{align}\begin{aligned} \ln L= & \sum_{i=1}^{\tau} \ln \lambda \beta t_{i}^{\beta-1}-\lambda T_{c}^{\beta} \\ & +\sum_{i=\tau+1}^{n}\left(\frac{a}{1+e^{b_{i}}}\right)+\frac{a\left[b T-\ln \left(e^{b T}+1\right)\right]}{b}\end{aligned}\end{align}\)       (4)

세분화 모형의 경우 구간별로 고장 강도 함수의 형태가 다르므로 변화점 T_c에서 서로 다른 고장 강도 값을 가질 수 있다. 그러나 T_c에서 고장 강도 함수의 연속성을 만족해야 하는 경우 t = T_c에서 두 고장 강도 함수의 값이 같아야 하므로 λβT^β-1_c=a/(1+e^bT_c)가 된다. 따라서 λ=a/(1+e^bT_c)(βT^β-1_c)]가 되며, 이를 식 (4)에 대입하면 대수우도 함수는 식 (5)와 같게 된다.

\(\begin{align}\begin{aligned} \ln L & =\sum_{i=1}^{\tau} \ln \frac{a}{1+e^{b T_{c}}}\left(\frac{t_{i}}{T_{c}}\right)^{\beta-1}-\frac{a T_{c}}{\left(1+e^{b T_{c}}\right) \beta} \\ & +\sum_{i=\tau+1}^{n} \ln \left(\frac{a}{1+e^{b t_{i}}}\right)+\frac{a\left[b T-\ln \left(e^{b T}+1\right)\right]}{b}\end{aligned}\end{align}\)       (5)

S-PLP and LIP 모형의 모수 추정은 최우추정법을 활용하며, β, a, b의 최우 추정량은 식 (5)의 대수우도 함수를 최대로 하는 \(\begin{align}\hat \beta \end{align}\), \(\begin{align}\hat {a}\end{align}\), \(\begin{align}\hat {b}\end{align}\)가 된다. 그러나 \(\begin{align}\hat \beta\end{align}\), \(\begin{align}\hat {a}\end{align}\), \(\begin{align}\hat {b}\end{align}\)에 대한 Closed-form 형태의 해를 구할 수 없으므로 모수 β, a, b에 대한 최우 추정은 Newton-Raphson 방법과 같은 수치해석 방법을 활용하여 구한다. λ에 대한 최우추정량 \(\begin{align}\hat{\lambda}\end{align}\)은, 모수 β, a, b에 대한 최우추정량 \(\begin{align}\hat {\beta}\end{align}\), \(\begin{align}\hat {a}\end{align}\), \(\begin{align}\hat {b}\end{align}\)와 T_c를 이용하여 식 (6)과 같이 구한다.

\(\begin{align}\hat{\lambda}=\hat{a} /\left[\left(1+e^{\hat{b} T_{c}}\right)\left(\hat{\beta} T_{c}^{\hat{\beta}-1}\right)\right]\end{align}\)       (6)

T_c 값을 모르는 경우, 다양한 방법을 이용하여 T_c 값을 구할 수 있으나 이 연구에서는 Guo et al.[8]이 제안한 휴리스틱 방법을 이용하여 T_c 값을 구하였다. Guo et al. [8]의 휴리스틱 방법은 5단계로 이루어져 있으며, T_c를 포함하여 β, a, b의 최우 추정량을 구하는 방법은 다음과 같다.

단계 1 : 시간에 따른 고장 수 그래프 작성

단계 2 : 단계 1 그래프에서 첫 번째 고장을 T_c,min, 마지막 고장을 T_c,max로 하여 T_c의 범위를 산출하고, 이를 바탕으로 적합한 ∆T_c 값을 선정

단계 3 : i번째 T_c 를 T_c,i = T_c,min + i × ΔT_c로 놓고 T_c,i에서 \(\begin{align}\hat{\lambda}\end{align}\), \(\begin{align}\hat{\beta}\end{align}\), \(\begin{align}\hat{a}\end{align}\), \(\begin{align}\hat{b}\end{align}\)을 추정

단계 4 : T_c,i에서 대수우도 값 계산

단계 5 : i=i+1, T_c,i가 T_c,max 될 때까지 단계 3과 단계 4를 반복

단계 1∼단계 5에 따라 T_c,i별 \(\begin{align}\hat{\lambda}\end{align}\), \(\begin{align}\hat{\beta}\end{align}\), \(\begin{align}\hat{a}\end{align}\), \(\begin{align}\hat{b}\end{align}\) 추정 값과 대수우도 값을 계산하였을 때, 대수 우도값을 최대로 하는 T_c,i와 모수 추정 값들을 S-PLP and LIP 모형의 모수에 대한 최우추정의 결과로 한다.

4. 적용 사례

적용 사례에서 다룬 대상 체계는 기동화력 무기체계 중 하나인 A체계이며, 시험평가 단계 가운데 개발시험 (DT, Developmental Test) 단계에서 수집된 고장 데이터를 사용하였다. 그리고 주어진 개발시험 기간에 관측된 A체계의 고장은 NHPP에 따라 발생한 것으로 가정하였다.

수집된 고장 데이터를 S-PLP and LIP 모형에 적합하여 구한 분석 결과와 기존의 모형들 (PLP, GGO, Yamada, NN, Two Segment 모형)을 적용했을 때의 분석 결과를 비교하여 S-PLP and LIP 모형의 적용 가능성과 유효성을 확인하였다. 모형 간의 비교 기준은 AICc [9] (Akaike Information Criterion corrected)와 평균제곱오차인 MSE (Mean Squared Error)이며, 모형의 적합도는 CVM (Cramér-von Mises) 검정을 이용하여 확인하였다.

4.1 대상 체계

적용 사례에서 다루고 있는 A체계는 후방 지역작전 간 기동 타격 및 중요 시설 방호를 목적으로 개발된 체계이다. 시험평가 단계 중 개발시험 단계에서 수집된 고장 데이터는 거리로 기록되며 단위는 km이다. A체계의 시험평가 기간에 발생한 총 고장 수는 22건이며, 측정된 누적 고장 거리 데이터는 Table 2와 같다. Table 2에서 두 번째 및 세 번째 누적 고장 거리가 같은 것은 서로 다른 고장 모드가 거의 동시에 연속적으로 발생하였기 때문이다. 그리고 시험평가 종료 시점 즉, 정시 관측 중단 (Type I Censoring) 시점에서 A체계의 총 누적 주행거리는 35,081km이다.

Table 2. Cumulative Failure Time of System A (km)

(+: Type I Censoring Time)

4.2 모형별 분석 결과

A체계 개발시험 단계 고장 데이터를 대상으로 모형 적합을 고려한 고장 강도 모형들의 대수우도값 산출과 모수 추정은 R 소프트웨어의 maxLik 패키지[10]를 활용하였다. A체계의 Upside-Down 욕조 곡선 형태의 고장 강도에 가장 적합한 모형을 비교 분석하는 비교 기준으로 사용한 AICc와 MSE는 식 (7), 식 (8)과 같다. 식 (7)에서 L은 우도 값, k는 모수 개수, n은 고장 수이고 식 (8)에서 \(\begin{align}\hat{\lambda}(t_i)\end{align}\)는 t_i에서 모형별로 추정된 고장 강도 추정 값이며 Y_i는 실제 데이터 t_i에서 관측된 고장 강도 값이다. AICc의 경우 값이 낮을수록 모형의 적합도가 좋은 모형을 의미하며, MSE의 경우에도 값이 낮을수록 주어진 데이터를 해당 모형이 더 잘 설명함을 의미한다.

AICc = -2ln(L) + 2k + (2k² + 2k)/(n - k - 1)       (7)

\(\begin{align}MSE=\sum_{i=1}^{n}\left(Y_{i}-\hat{\lambda}\left(t_{i}\right)\right)^{2} / n\end{align}\)       (8)

이 연구에서 제안한 S-PLP and LIP 모형의 모수 추정 결과와 PLP, GGO, Yamada, NN, Two Segment 모형들의 모수 추정 결과 및 모수 추정 결과를 바탕으로 작성한 모형별 고장 강도 그래프는 Fig. 2∼Fig. 7과 같다.

Fig. 2 PLP Model

Fig. 3 GGO Model

Fig. 4 Yamada Model

Fig. 5 NN Model

Fig. 6 Two Segment Model

Fig. 7 S-PLP and LIP Model

그리고 AICc와 MSE를 기준으로 모형들을 비교한 결과는 Table 3과 같다. AICc와 MSE는 값이 작을수록 해당 모형이 주어진 데이터를 잘 설명한다고 말할 수 있는데, S-PLP and LIP 모형의 경우 AICc, MSE 값이 각각 284.93, 4.83×10^-8로 다른 고장 강도 모형의 AICc, MSE 값들보다 작음을 알 수 있다. 즉, S-PLP and LIP 모형이 A체계 개발시험 평가 단계 데이터의 고장 강도를 가장 잘 설명한다고 말할 수 있다. 실제로 S-PLP and LIP 모형의 경우 변화점을 기준으로 구간을 세분화하여 구간별로 모형을 적합하고 연속성 조건을 만족하는 형태로 이루어졌기 때문에 기존의 Upside-Down 욕조 곡선 형태의 고장 강도 모형보다 더 잘 적합할 수 있다고 판단된다.

Table 3. MLE, AICc, MSE of Failure Intensity Models

최우추정법을 이용한 S-PLP and LIP 모형의 모수 추정값은 각각 \(\begin{align}\hat{\lambda} = 1.88 \times 10^{-6}\end{align}\), \(\begin{align}\hat{\beta} = 1.6903\end{align}\), \(\begin{align}\hat{a}=0.0040\end{align}\), \(\begin{align}\hat{b}=6.13\times 10^{-5}\end{align}\)이고, 변화점 T_c값은 7,797이다. 이를 이용하여 S-PLP and LIP 모형의 고장 강도 함수식을 구하면 식 (9)와 같다.

\(\begin{align}\hat{\lambda}(t)=\left\{\begin{array}{ll}\left(1.88 \times 10^{-6}\right)(1.6903) t^{(1.6903-1)}, & 0<t \leq 7,797 \\ \frac{0.0040}{1+e^{\left(6.13 \times 10^{-5}\right)},} & 7,797<t<\infty\end{array}\right.\end{align}\)       (9)

A체계의 개발시험 단계 고장 데이터에 대한 S-PLP and LIP 모형의 적합도 검정은 CVM (Cramér-von Mises) 검정 방법을 이용하였다. T시간까지 관측된 N개의 고장 데이터가 정시 중단 데이터인 경우 CVM 검정통계량 C²은 식 (10)과 같다. 여기서 t_i는 i번째 고장 시간, Λ(t_i)는 t_i에서의 누적 고장 강도, Λ(T)는 T에서의 누적 고장 강도이다.

\(\begin{align}C^{2}=\frac{1}{12 N}+\sum_{i=1}^{N}\left[\frac{\Lambda\left(t_{i}\right)}{\Lambda(T)}-\frac{2 i-1}{2 N}\right]^{2}\end{align}\)       (10)

유의 수준 α에 따른 모형의 CVM 통계량에 대한 기각 값 (Critical Value)의 산출은 Park and Kim [11]이 제안한 몬테카를로 시뮬레이션을 이용한 방법을 활용하였으며 그 결과는 Table 4와 같다.

Table 4. Critical Value of Cramér-von Mises Statistic

A체계 개발시험 단계의 고장 데이터의 경우, S-PLP and LIP 모형에 대한 C²는 0.2684이며 유의수준이 0.05일 때 기각 값은 Table 4로부터 0.4946이다. 따라서 S-PLP and LIP 모형의 C²값이 기각 값보다 작으므로 유의수준 0.05에서 주어진 데이터가 해당 모형을 따른다고 말할 수 있다.

5. 결론

이 연구에서는 수리 가능한 체계의 고장 발생이 NHPP를 따른다는 가정하에서, Upside-Down 욕조 곡선 형태의 고장 강도를 갖는 S-PLP and LIP 모형을 제시하고 최우추정법을 이용하여 모수에 대한 추정량을 구하는 방법을 함께 제시하였다. S-PLP and LIP 모형은 변화점을 기준으로 구간을 구분하고, 고장 강도가 증가하는 첫 번째 구간의 고장 강도 적합은 PLP 모형을, 고장 강도가 감소하는 두 번째 구간의 적합은 LIP 모형을 고려하였다.

적용 사례에서는 기동화력 무기체계인 A체계의 시험평가 단계 고장데이터를 이용하여 S-PLP and LIP 모형을 적합시키고 모수를 추정하였다. CVM 검정을 이용한 적합도 검정 결과 S-PLP and LIP 모형의 C² 값은 0.2684이며, 유의수준 0.05일 때 기각 값 0.4946 보다 작으므로 A체계의 고장 데이터를 S-PLP and LIP 모형으로 적합시켜도 문제가 없음을 확인하였다. 또한 기존 모형들의 추정 결과와 비교할 때, S-PLP and LIP 모형의 AICc 는 284.93이고 MSE는 4.83×10^-8으로 다른 모형들에 비해 충분히 작은 값을 나타냄으로 S-PLP and LIP 모형이 주어진 A체계의 고장 데이터를 상대적으로 더 잘 적합할 수 있음을 확인하였다.

이 연구에서 제안한 S-PLP and LIP 모형은 고장 강도가 Upside-Down 욕조 곡선 형태를 따를 때 고장 강도 모형의 대안으로 사용할 수 있으며, 수리 가능한 체계에 대한 보다 정확한 신뢰도 추정과 예측에 활용할 수 있을 것으로 기대한다.