컨베이어벨트 수명분포 추정 — 고장 시각을 모르는 데이터로 수명을 추정한다

무엇을 물었나 — 같은 설비에 '언제 고장 나는가'

고장의 원인이 아니라 고장까지 걸리는 시간을 확률분포로 답하는 문제였다.

대용량 스토리지 시스템(High Rack Storage System)의 컨베이어 벨트 고장 데이터를 가지고, 최적화된 방식으로 작동하는 벨트와 최적화되지 않은 방식으로 작동하는 벨트의 수명분포를 각각 추정하는 것이 이 프로젝트의 목적이었다.

같은 설비의 로그를 두고 ‘어느 구간의 컨베이어가 고장의 원인인가’를 분류로 물었던 적이 있다. 이번에는 질문이 달랐다. 원인이 아니라 시간이었고, 답도 라벨이 아니라 분포였다.

수명분포를 물으면 곧바로 걸리는 문제가 하나 있다. 수명을 알려면 고장 시각을 알아야 하는데, 이 데이터에는 고장 시각이 없다. 있는 것은 ‘언제 확인해 보니 이미 고장이었다’뿐이다.

두 가지 운전 방식 — 무엇이 '최적화'인가

물체가 아래로 내려가는 동안 옆으로 미는 컨베이어를 함께 돌리느냐 마느냐가 두 방식을 가른다.

최적화된 방식은 물체가 아래로 내려갈 때 옆으로 이동시키는 컨베이어 벨트도 같이 동작하는 운전이다. 최적화되지 않은 방식은 아래로 이동하는 중에는 옆으로 이동하는 컨베이어가 작동하지 않는 운전이다. 두 방식은 같은 설비에서 같은 일을 하지만, 한 사이클 안에서 동시에 움직이는 부분의 수가 다르다.

설비는 가로로 움직이는 네 개의 컨베이어(BLO, BHL, BHR, BRU)와 위아래로 움직이는 두 부분(HR, HL)으로 이루어져 있다. 데이터는 이 부분들 각각의 power·전압·전력·거리를 일정 간격의 타임스탬프마다 기록한다.

따라서 이 프로젝트의 비교는 '두 제품의 수명 비교'가 아니라 같은 설비를 다르게 운전했을 때의 수명 비교에 가깝다. 운전 방식이 수명분포의 모양 자체를 바꾸는지가 궁금한 지점이었다.

데이터 — 두 개의 로그, 20개 열

고장 여부는 0/1 라벨로만 들어 있고, 라벨이 바뀌는 지점이 사실상 유일한 정보였다.

두 운전 방식에 대해 각각 하나의 로그가 있었다. 열은 20개이고, 첫 열이 타임스탬프, 둘째 열이 고장 라벨(정상 0, 고장 1)이다. 나머지는 각 컨베이어 부분의 전기·거리 신호다.

측정 횟수는 타임스탬프가 0인 지점의 개수로 셌다. 타임스탬프가 0으로 되돌아가는 지점이 새 측정의 시작이므로, 이 개수가 곧 관측 회차의 수가 된다. 행 수를 회차 수로 나누면 회차 하나가 평균 221행(비최적화)·177행(최적화)이다. 원자료 화면에서 타임스탬프 간격이 0.05 안팎인 것과 곱해 보면 한 회차의 길이는 10 안팎이 되고, 뒤에서 볼 고장 확인 시각의 최대값 13.460999과도 대체로 맞는다.

수명분포에 쓸 수 있는 정보는 결국 둘째 열이 0에서 1로 바뀌는 지점 하나뿐이었다. 나머지 18개 신호 열은 이 프로젝트에서는 쓰지 않았다.

왜 관측이 구간으로만 주어지는가

고장 시각을 본 것이 아니라, 고장이 이미 일어났다는 사실을 나중에 확인한 것이다.

이 데이터는 일정 시간 간격으로 고장 여부를 확인한 구간관측중단(interval censored) 데이터다. 어떤 벨트가 t = 3.956001에 확인했을 때 정상이었고 t = 4.011002에 확인했을 때 고장이었다면, 우리가 아는 것은 다음 두 가지뿐이다.

• 수명이 3.956001보다는 길다 → 이 시점 기준으로는 우측관측중단
• 수명이 4.011002보다는 짧다 → 이 시점 기준으로는 좌측관측중단

두 정보를 합치면 수명 T가 (3.956001, 4.011002] 안 어딘가에 있다는 것까지만 말할 수 있다. 원본 보고서가 각 구간의 하한을 '우측관측중단 시각', 상한을 '좌측관측중단 시각'이라고 부른 것도 이 분해를 그대로 옮긴 이름이다. 전처리는 이 두 시각을 따로 저장하는 일이었다.

이 구조를 우도에 반영하지 않고 고장이 확인된 시각을 정확한 고장 시각처럼 쓰면 수명을 한쪽으로 밀어서 추정하게 된다. 다만 정직하게 적어 두면, 이 데이터에서 그 편향의 크기는 구간 폭을 넘지 못한다. 보고서에 예시로 남은 네 개의 구간에서 폭은 0.045~0.075였고 뒤에서 추정될 평균 수명은 6~8 부근이므로, 구간관측중단 처리가 결과를 크게 바꾸는 요소는 아니었다. 그래도 우도를 구간으로 쓰는 것이 이 과제의 요구였고, 편향의 방향과 크기를 스스로 확인할 수 있다는 점에서 의미가 있었다.

관측 구간 추출 — 라벨이 바뀌는 순간을 잡는 알고리즘

고장 상태를 기억해 두고, 새 사이클이 시작되면 그 상태를 되돌리는 규칙으로 구간을 뽑았다.

전처리의 핵심은 23,645행과 19,634행을 훑으면서 고장 라벨이 0에서 1로 바뀌는 지점을 찾아, 그 직전 행의 타임스탬프와 그 행의 타임스탬프를 구간의 하한과 상한으로 각각 쌓는 것이다. 규칙은 셋이었다.

1. 라벨이 1인데 아직 고장으로 표시해 두지 않은 상태라면 → 고장 상태로 바꾸고, 상한에 현재 시각을, 하한에 직전 행(마지막 정상 관측)의 시각을 기록한다.
2. 라벨이 0인데 고장 상태로 표시되어 있다면 → 새 사이클이 시작된 것이므로 상태를 정상으로 되돌린다.
3. 그 외에는 넘어간다.

상태를 기억해야 하는 이유는 고장 라벨이 한 행에만 찍히는 것이 아니라 여러 행에 걸쳐 이어지기 때문이다. 상태를 기억하지 않으면 같은 고장 하나가 수십 개의 구간으로 중복 기록된다. 반대로 새 사이클에서 상태를 되돌리지 않으면 이후 사이클의 고장을 전혀 잡지 못한다.

지금 다시 보면 이 절차는 결과를 냈을 뿐 견고하지는 않았다. 훑을 행 수를 데이터에서 읽어 오는 대신 23,645·19,634라는 수를 직접 적어 넣었고, 두 벨트를 이어서 처리할 때 앞 벨트에서 남은 고장 상태가 그대로 넘어가면 뒤쪽 벨트의 첫 고장을 놓칠 수 있는 구조였다. 결과가 맞았던 것과 절차가 견고한 것은 다른 문제였다.

이상값 제거 — IQR 규칙을 수명 데이터에 적용하면 생기는 일

규칙을 그대로 적용했지만, 한쪽 벨트에서는 실제로 제거된 관측이 하나도 없었다.

추출한 구간의 하한 시각에 IQR 규칙을 적용해 이상값을 걸렀다. 비최적화 벨트의 사분위수는 다음과 같았다.

IQR은 10.007500 − 5.903748 = 4.103752, 1.5배는 6.155628이므로 경계는 −0.25188과 16.163128이 된다. 최적화 벨트는 1사분위수 3.956001, 3사분위수 8.796997, IQR 4.840996이어서 경계가 −3.305493과 16.058491이었다.

여기서 두 가지가 눈에 걸린다. 첫째, 비최적화 벨트의 최소값 3.462997과 최대값 13.460999은 모두 경계 안에 있다. 즉 이 규칙으로 실제로 제거된 관측은 하나도 없었다. 둘째, 두 경우 모두 하한이 음수다. 수명은 정의상 양수이므로 하한 조건은 애초에 걸릴 수 없고, 규칙은 사실상 상한 하나로만 작동한다.

대칭을 전제하는 IQR 규칙을 오른쪽으로 치우친 수명 데이터에 그대로 옮긴 결과였다. 당시에는 '전처리 절차를 밟았다'는 데 만족했지만, 규칙을 적용한 뒤 몇 개가 제거되었는지를 세어 보지 않은 것이 이 절의 실질적인 공백이다.

우도함수 — 구간을 그대로 우도에 넣는다

정확 관측이 하나도 없으므로 밀도함수 f(t)는 등장하지 않는다. 전부 CDF의 차이다.

수명 T의 누적분포함수를 F라 하면, 하한 a와 상한 b 사이에서 고장이 확인된 관측 하나가 우도에 기여하는 값은 F(b) − F(a)다. 관측이 독립이라고 보면 전체 우도는 이 값들의 곱이고, 로그를 취하면 합이 된다. 최대화 대신 음의 로그우도를 최소화하는 목적함수로 두고 수치 최적화로 모수를 찾았다.

세 분포를 후보로 두었고, 누적분포함수만 바꿔 끼우면 되도록 같은 형태의 목적함수를 세 번 썼다.

와이블의 척도모수 초기값을 상한 시각의 63.2% 분위수로 둔 것은 우연이 아니다. 와이블에서 t = η이면 F(η) = 1 − exp(−1) ≈ 0.632이므로, 척도모수는 대략 관측의 63.2% 분위수 근처에 있다. 초기값을 이렇게 잡으면 수치 최적화가 엉뚱한 곳에서 출발하지 않는다.

목적함수는 결국 한 줄이었다 — 관측마다 상한의 누적확률에서 하한의 누적확률을 빼고, 그 로그를 모두 더한 뒤 부호를 뒤집은 것. 이 한 줄이 ‘구간관측중단 데이터의 우도’라는 말 전체에 해당한다는 것을 손으로 써 보고 나서야 납득했다.

직접 구현한 결과를 기성 도구와 맞춰 보았다

직접 구현한 결과를 믿으려면 비교 대상이 필요했다.

같은 데이터를 관측중단 데이터를 최대우도추정으로 분포에 적합시키는 기성 도구에도 넣었다. 세 분포 모두 같은 방식으로 적합하고 요약과 적합도 그림을 확인했다.

결과는 두 벨트 모두에서 직접 최적화해 얻은 모수와 기성 도구의 모수가 같았다. 이것이 이 프로젝트에서 가장 얻고 싶었던 확인이다. 기성 도구를 부르면 결과는 한 줄로 나오지만, 그 한 줄이 각 관측에 대해 누적분포함수의 차이를 어떻게 곱해 우도를 만드는지는 보이지 않는다. 직접 구현해 두면 같은 값이 나오는지 검증할 수 있고, 결과를 남의 말이 아니라 제 말로 설명할 수 있게 된다.

여기서 지금은 판정할 수 없는 대목이 하나 남는다. 내가 구간의 양 끝에 붙인 이름은 관측중단의 방향을 기준으로 한 것이어서, 기성 도구가 기대하는 하한·상한의 뜻과 정반대였다. 그대로 넘겼다면 상한과 하한이 뒤바뀌어 누적확률의 차이가 음수가 되고 로그를 취하는 순간 계산이 깨졌을 것이므로, 문서에 남지 않은 정리 단계가 있었다고 보는 편이 자연스럽다. 다만 지금 자료로는 확인되지 않는다.

그럼에도 이 대목에서 배운 것은 남는다. 내가 도메인 용어로 붙인 이름과 도구가 기대하는 이름의 뜻이 다를 수 있다는 것, 그리고 그 어긋남은 통계의 문제가 아니라 자료를 넘기는 규약의 문제라는 것이다.

추정 결과 — 두 벨트 모두 와이블

AIC로 고른 분포는 같았지만, 모수는 뚜렷하게 달랐다.

지수·와이블·대수정규 세 분포를 각각 적합하고 AIC를 비교했다. 두 벨트 모두에서 AIC가 가장 작은 것은 와이블분포였고, 적합도 그림으로도 같은 결론이 나왔다.

형상모수가 두 경우 모두 1보다 뚜렷하게 크다는 점은 그 자체로 읽을 것이 있다. 와이블의 형상모수가 1보다 크면 고장률이 시간에 따라 증가한다는 뜻이고, 이는 초기 불량이나 우발 고장이 아니라 마모에 의한 고장에 가깝다는 해석으로 이어진다. 지수분포(형상모수를 1로 고정한 경우에 해당)가 선택되지 않은 것도 같은 이야기다.

다만 AIC와 BIC의 실제 수치는 보고서 본문에 서술로만 남아 있고 표의 그림 안에 있어 지금은 확인할 수 없다. 방법론에는 AIC와 BIC를 모두 쓰겠다고 적어 두었지만, 결과 서술에서 실제로 언급한 것은 AIC뿐이다.

두 운전 방식의 비교 — '최적화'가 수명을 늘리지는 않았다

척도모수가 9.36에서 7.02로, 형상모수가 3.67에서 2.76으로 내려갔다.

추정된 모수만 놓고 보면 최적화된 방식으로 운전한 벨트 쪽의 수명이 더 짧고, 고장 시점의 흩어짐은 상대적으로 더 크다. 아래 값은 보고서에 없는 값으로, 위의 추정 모수로부터 내가 계산한 것이다. 와이블 평균은 η·Γ(1 + 1/β), 중앙값은 η·(ln 2)^(1/β), B10 수명은 η·(−ln 0.9)^(1/β)로 구했다.

이렇게 보면 아래로 내려가는 동안 옆 컨베이어까지 함께 돌리는 운전은 처리량 관점에서는 최적화일 수 있어도, 설비 수명 관점에서는 그렇지 않다고 읽힌다. 동시에 움직이는 부분이 많아질수록 단위 시간당 부하가 커지기 때문이라는 설명이 가능하다.

정직하게 덧붙이면, 원본 보고서는 이 비교를 여기까지 밀고 가지 않았다. 두 벨트의 결과를 나란히 적고 ‘둘 다 와이블이 가장 잘 적합한다’는 데서 결론이 끝난다. 두 데이터의 시간 단위가 같은지, 두 벨트가 같은 기간·같은 부하에서 관측되었는지도 보고서로는 확인되지 않는다. 위의 비교는 그 조건들이 같다는 가정 위에서만 성립한다.

확인되지 않은 채 남은 것들

결과를 다시 열어 보려니 정작 필요한 숫자가 문서에 남아 있지 않았다.

이 프로젝트를 다시 정리하면서 확인하지 못한 것들을 그대로 적어 둔다.

• AIC·BIC의 실제 값 — 보고서에서는 세 분포의 결과를 표에 스크린샷으로 넣고 본문에는 ‘와이블이 가장 작다’는 서술만 남겼다. 텍스트로 남은 수치가 없어 세 분포의 격차가 얼마나 컸는지 알 수 없다.
• 적합도 그림 — 추정된 분포의 누적분포함수를 실측과 겹쳐 그린 비교 그림이 문서에 있었으나 지금 남은 자료에서는 복원되지 않는다.
• 시간 단위 — 타임스탬프의 단위(초·분·시)가 기록되어 있지 않다. 그래서 ‘평균 수명 8.44’가 물리적으로 얼마인지 말할 수 없다.
• 이상값 제거 이후의 표본 수 — 제거 전후의 개수를 세어 두지 않았다. 처리할 관측 수를 측정 회차 수(107·111)로 고정해 두었기 때문에, 추출된 구간의 개수가 그와 다르면 뒤쪽 구간이 조용히 빠질 수 있는 구조였다.
• 교체 주기 — 목적은 수명분포 추정까지였고, 추정된 분포를 보전 정책(예방 교체 주기, 비용 최소화)으로 옮기는 단계는 하지 않았다.

수치를 그림으로만 남긴 것이 가장 아쉽다. 결과를 표로 옮겨 적는 습관이 없으면, 몇 년 뒤의 자신은 자기 결과를 인용할 수 없다.

여기서 남은 것

도구 한 줄로 끝날 일을 굳이 손으로 쓴 이유가 이후의 기준이 되었다.

이 프로젝트에서 실제로 어려웠던 부분은 통계가 아니라 관측 방식을 우도로 번역하는 일이었다. 데이터가 ‘언제 고장 났는가’가 아니라 ‘언제 확인해 보니 고장이었는가’로 기록되어 있다는 사실을 인정하면, 우도에는 밀도함수 대신 CDF의 차이가 들어간다. 모델을 고르기 전에 관측이 어떻게 만들어졌는지를 먼저 봐야 한다는 것을 이때 처음 몸으로 배웠다.

기성 도구를 부르면 한 줄로 끝나는 일이었다. 그럼에도 로그우도를 직접 짜고 수치 최적화로 같은 값에 도달하는지 확인한 뒤에야, 결과를 내 말로 설명할 수 있게 되었다. ‘패키지의 결과라도 원리를 한 번은 직접 구현해 확인한다’는 원칙이 여기서 생겼다.

동시에 남은 반성도 분명하다. 전처리 절차는 결과가 맞았을 뿐 견고하지 않았고, IQR 규칙은 절차를 밟았다는 만족만 남긴 채 실제로는 아무것도 걸러 내지 않았을 가능성이 높다. 결과 수치를 그림으로만 보관한 탓에 지금은 자기 결과의 절반을 인용하지 못한다. 방법을 제대로 이해하는 것과, 그 방법을 재현 가능하게 남기는 것은 별개의 일이었다.