점탄성 데이터에 선형 스프링 모델 적용하기

Kelvin-Voigt 모델로 만든 점탄성 시편 데이터를 F = kx로 적합했을 때 무엇이 깨지는지 확인하고, 그 다음 단계인 댐퍼 항을 추가한 프로젝트입니다. 개인 학습으로 진행했습니다.

한눈에 보기

• 질문: 점성이 있는 재료의 힘-변위 데이터를 스프링만 있는 모델로 적합하면 어떤 한계가 나타나는가?
• 방법: F = kx + c·x_dot로 가상 데이터를 생성하고, F = kx로 회귀해 기준선을 잡은 뒤, 수치미분으로 속도 x_dot를 복원해 해석해와 비교한 뒤, [x, x_dot]설계행렬을 세워 최소제곱으로 k, c를 동시에 추정했다.
• 결과 1: 강성 k는 200.31 N/m로 정확히 추정됐지만 R²는 0.777로 떨어졌다. 모델이 틀려도 k가 맞는 이유는 왕복 운동에서 댐퍼 기여가 상쇄되기 때문이다.
• 결과 2: 중심차분 오차가 이론 절단오차의 3배였다. 원인은 알고리즘이 아니라 CSV 저장 정밀도(반올림 오차의 증폭)였고, 정밀도를 올리자 이론값과 일치했다.
• 결과 3: [x, x_dot] 설계행렬로 k, c를 동시에 추정하자 R²이 0.777에서 0.99995로 회복되었다. 과적합을 검증하기 위해 두 계수를 비교했고 두 계수 모두 참값에 근접했다(k 오차 0.008%, c 오차 0.268%).
• 진행 상태: [x, x_dot] 설계행렬로 k, c를 동시에 추정하는 다중선형회귀까지 4단계 완료.

배경

학기 전 연결학습에서 스프링 시편의 힘-변위 데이터로 다음을 했다.

• F = kx 회귀로 강성 추정
• 이상값 검출과 선형성 검사(R² >= 0.99를 합격 기준으로 설정)
• 여러 시편 CSV를 자동으로 분석해 판정

이 모델에는 시간이 없다. 실제 생체조직처럼 당기는 속도에 따라 힘이 달라지는 재료를 다루려면 속도 항이 필요하다. 그래서 스프링과 댐퍼를 병렬로 놓은 Kelvin-Voigt 모델 F = kx + c·x_dot로 확장하기로 했다. 정지 상태에서는 F = kx로 돌아가므로 기존 모델을 버리는 것이 아니라 포함한다.

진행 단계

단계	내용	파일	상태
1	Kelvin-Voigt 모델로 점탄성 가상 데이터 생성, 히스테리시스 확인	make_viscoelastic_data.py	완료
2	F = kx만으로 회귀해 기준선 확보	estimate_kx_only.py	완료
3	수치미분으로 x_dot 복원, 해석해와 오차 비교, 오차 원인 추적	check_numerical_derivative.py	완료
4	np.linalg.lstsq로 k, c 동시 추정, R² 회복 확인	estimate_full.py	완료

파일 구성

• make_viscoelastic_data.py: k=200 N/m, c=30 N·s/m, 진폭 0.03 m, 주기 2 s의 사인 변위로 힘을 계산하고 노이즈를 더해 저장한다.
• viscoelastic_test.csv: 생성된 테스트 데이터(101행, 시간·변위·힘).
• estimate_kx_only.py: F = kx 회귀와 R² 계산.
• check_numerical_derivative.py: np.gradient 중심차분으로 속도를 복원하고 해석해와 비교한다.
• estimate_full.py: [x, x_dot] 설계행렬로 k와 c를 최소제곱으로 동시에 추정하고, estimate_kx_only.py의 R²와 비교한다.
• 설계 노트/: 코드를 쓰기 전 손으로 작성한 목표 출력·변수 목록·의사코드(PDF). 수정 전후 두 판본을 모두 남겼다.

실행 방법

필요한 환경:

• Python 3.14
• NumPy

실행 순서:

python make_viscoelastic_data.py
python estimate_kx_only.py
python check_numerical_derivative.py
python estimate_full.py

결과

선형 스프링 모델의 한계

	k 추정	R²
F = kx	200.31 N/m (참값 200.00)	0.77688

모델이 데이터의 약 22%를 설명하지 못한다. 그 몫이 댐퍼 항이다. 주목할 점은 모델이 크게 틀렸는데도 강성 추정은 정확했다는 것이다. 왕복 운동에서 댐퍼의 기여가 갈 때와 올 때 상쇄되어 회귀 직선의 기울기를 밀지 못하기 때문이다. 즉 이 시편은 불량이 아니지만 이전 프로젝트에서 정한 검사 규격인 R² >= 0.99로는 불합격한다. 불합격한 것은 시편이 아니라 모델이다.

R²는 두 항의 변동 비율로 손계산할 수 있다. 3.410 / (3.410 + 0.999) = 0.773이고 실측 0.77688과 맞았다.

수치미분의 오차

x 저장 정밀도	Max abs error	Mean abs error
소수점 6자리	0.000024	0.000012
소수점 10자리	0.000008	0.000005

중심차분의 이론 절단오차는 (h²/6)|x'''| = 7.75e-6이다. 6자리 저장에서는 실측이 그 3배였는데, 원인은 반올림 오차 5e-7이 h = 0.02로 나뉘며 50배로 증폭된 것이었다. 정밀도를 높이자 이론값과 일치했다.

절단오차는 h가 작을수록 줄고 반올림오차는 커진다. 같은 반올림이 회귀 결과는 전혀 바꾸지 않았다. 미분은 이웃 점의 차이를 보고 회귀는 전체를 평균하기 때문이다. 연산의 성격이 오차 민감도를 정한다.

댐퍼 항을 추가한 결과

모델	k	c	R²
F = kx + b	200.31 N/m	—	0.77688
F = kx + c·x_dot	200.02 N/m	30.08 N·s/m	0.99995

설명하지 못하던 22%가 사라졌다. 남은 0.005%는 데이터에 넣은 측정 노이즈다.

R²만으로는 모델이 옳아진 것인지 표현력만 늘어난 것인지 구별할 수 없다. 여기서는 R²가 오른 동시에 두 계수가 참값을 맞혔다(k 200.02 대 200.00, c 30.08 대 30.00). 따라서 추가한 댐퍼 항이 데이터의 생성 과정과 잘 부합한다는 근거로 볼 수 있다.

두 계수의 정확도는 같지 않았다. 상대오차가 k는 0.008%, c는 0.268%였다. 두 항이 힘에 기여하는 크기가 다르기 때문으로 보인다. k·x의 진폭은 6.00 N, c·x_dot의 진폭은 1.41 N으로 4.2배 차이가 나는데, 노이즈는 양쪽에 같은 크기로 들어가므로 작은 항의 계수가 상대적으로 더 흔들린다.

이 프로젝트에서 배운 것

• 새 모델의 성능은 절대값이 아니라 기준선 대비로만 의미가 있다. 그래서 먼저 실패하는 모델을 만들었다.
• R²는 모델이 데이터를 얼마나 설명하는지를, k 오차는 추정값이 맞는지를 잰다. 둘은 독립이다.
• F = kx + c·x_dot는 미지수 k, c에 대해 여전히 1차식이다. 기저함수가 x와 x_dot, 가중치가 k와 c인 선형회귀이며, 머신러닝 수업의 y = Σ w_j f_j(x)와 같은 뼈대다.
• 다항 차수를 올리는 것과 물리에서 나온 항을 추가하는 것은 다르다. 전자는 표현력만 늘려 과적합 위험이 있고, 후자는 모델을 교정한다.
• 오차가 이론과 다르면 알고리즘을 의심하기 전에 데이터가 거쳐 온 경로(저장 정밀도)를 먼저 본다.
• R²가 오른 것만으로는 모델이 올바른지 확정할 수 없다. 계수들이 참값에 가까운지 함께 봐야 표현력만 늘린 것과 구별된다.

작업 방식

• 코드를 쓰기 전에 목표 출력을 먼저 적고, 빈칸에서 변수 목록을 뽑고, 한국어 의사코드를 쓴 뒤 한 줄씩 번역했다. 손 설계 노트는 설계 노트/에 PDF로 남겼다.
• README에는 직접 실행해 확인한 것만 단정해서 썼다.
• 매일 학습일지를 작성했다(로컬 보관, 저장소에는 포함하지 않음).

참고

코드와 설계는 직접 작성했고, 개념 설명과 결과 해석은 Claude Opus 5와의 대화를 통해 정리했다.