제3부 동역학

3.12 접촉·충돌 동역학: 상보성(LCP)·마찰 원뿔·충격 — 헤르츠·헌트-크로슬리 연성 접촉과 하이드로엘라스틱 모델

3.12.1 이 장에서 배우는 것 / 선수 지식 / 특히 중요한 독자

이 장은 로봇이 환경과 닿고 부딪힐 때 접촉력을 어떻게 정하는지를 동역학 모형으로 정리한다. 모형은 두 계열이다. 강체 계열은 접촉력을 "거리가 0이거나 힘이 0"이라는 상보성 조건과 쿨롱 마찰 원뿔로 규정하고 다음 시간 간격의 속도를 선형 상보성 문제(linear complementarity problem, LCP)로 푼다. 연성 계열은 침투 깊이의 함수로 접촉력을 직접 정하며 헤르츠 스프링, 헌트-크로슬리 감쇠, 하이드로엘라스틱 압력장이 여기에 속한다. 이 장은 두 계열의 식을 세우고, 레믹 알고리즘으로 LCP를 풀고, 다면체 마찰 원뿔의 이방성·강체 접촉력의 비유일성·충격의 유효 질량·연성 충돌의 반발계수·하이드로엘라스틱 접촉력을 계산으로 확인한 뒤 SciPy·Pinocchio·MuJoCo·Drake와 같은 입력에서 비교한다. 시간 적분 방식과 엔진별 구현·IPC는 16.1, 접촉이 있는 파지 역학은 13.2, 관절 마찰은 3.8, 구속 동역학과 라그랑주 승수는 3.11, 부유기저와 다접촉 보행은 3.13·12.7, 충돌 안전의 생체역학 한계는 26.9과 이어진다.

배우는 것:

• 접촉의 규정 — 법선 상보성, 쿨롱 마찰 원뿔(고착·미끄러짐·이탈), 접촉 야코비안과 Delassus 연산자(접촉 공간 역관성).
• LCP 정식화와 풀이 — 다면체 원뿔과 시간 간격 이산화, 레믹 피벗과 투영 가우스-자이델, 마찰 없는 경우의 QP 동치.
• 다면체 근사의 대가 — 슬립 임계의 방향 이방성 $\cos(\pi/m_d)$와 미끄러짐 방향 오차 $\pi/m_d$.
• 충격 — 뉴턴 반발 법칙, 유효 질량 $\Lambda_n$과 충격량 $(1+e)\Lambda_n s$, 에너지 손실 $\tfrac12(1-e^2)\Lambda_n s^2$.
• 연성 접촉 — 헤르츠 접촉력과 접촉 시간, 켈빈-포크트의 결함, 헌트-크로슬리 감쇠, 하이드로엘라스틱 압력장과 구-평면 닫힌 식.

선수 지식: 3.1(강체 동역학), 3.4(운동방정식과 질량행렬), 3.11(구속 동역학), 1.13(QP), 1.18(ODE 적분기).

특히 중요한 독자: 시뮬레이터와 접촉 모형을 선택·설정하는 연구자, 다접촉 제어와 궤적 최적화를 구현하는 개발자, 충돌 안전 해석에서 유효 질량과 충격량을 다루는 안전 엔지니어.

출처 방침. 기준 시점은 2026-11-05이다. 직접 읽은 자료: Modern Robotics 사전 인쇄본(2017-05-03) 12.2.1절(쿨롱 마찰과 내접 다면체 원뿔)·12.3절(접촉 모드 열거와 모호성·불일치 서술)·장 끝 문헌 안내(상보성 문헌 지목), MuJoCo 공식 문서 "Computation" 장(볼록 소프트 접촉, 피라미드·타원 원뿔, 정규화 $R$, 솔버), Drake 문서 "Hydroelastic Contact User Guide"와 "Modeling Compliant Contact"(점 접촉 식, 헌트-크로슬리 소산, 압력장·접촉면), Pinocchio 4.1.0의 impulseDynamics 설명. 서지만 확인한 것(본문은 읽지 못했다): Stewart·Trinkle(1996), Anitescu·Potra(1997), Pang·Trinkle(1996), Stewart(2000), Hertz(1882), Hunt·Crossley(1975), Marhefka·Orin(1999), Stronge(1990), Todorov(2014), Elandt 외(2019), Masterjohn 외(2022), Castro 외(2023), Johnson(1985), Featherstone(2008). 이들 문헌의 식은 이 장이 직접 세워 계산으로 검증했으며 원문 식과의 일치는 확인하지 못했다 (확인 필요). 헌트-크로슬리 힘의 형식 $F=kx^{3/2}(1+\tfrac32\alpha\dot x)$와 헤르츠 계수 $\tfrac43E^*\sqrt{R^*}$는 이 장이 정의로 둔 것이다. 구의 압력장 $p=E_h(1-r/R)$은 Drake 1.58.0의 계산과 일치하는 것으로 확인했다(문서에서 식을 읽은 것은 아니다). 강성·감쇠·마찰계수·질량은 모두 이 장이 정한 예시이며 어떤 제품의 제원도 아니다. 수치는 순수 파이썬 계산이다.

3.12.2 개념 — 접촉은 부등식과 불연속이다

결론부터 적으면 다음과 같다. 접촉은 연속 방정식이 아니라 부등식 조건이다. 접촉력은 압축만 가능하고(법선 힘 $\ge0$), 두 물체는 파고들지 않으며(거리 $\ge0$), 둘 중 하나는 반드시 0이다(상보성). 마찰력은 법선 힘의 $\mu$배를 넘지 못하고, 미끄러질 때는 미끄럼 방향과 반대이며 에너지를 소산한다. 충돌은 속도의 불연속이다. 강체로 보면 이 조건들이 모든 시간 간격에서 접촉마다 성립해야 하므로 시간 간격 하나가 LCP 하나이다. 연성으로 보면 침투를 허용하고 침투와 침투 속도로 접촉력을 정하므로 상미분방정식이 되지만, 강한 접촉일수록 강성이 커져 적분이 어려워진다(수치 비교는 풀이 예제 ⑤).

두 접근은 연속적이다. 강성 $k\to\infty$인 연성 접촉이 강체 접촉이고, 반대로 MuJoCo는 엄격한 상보성을 버리고 정규화 $R$을 둔 볼록 문제로 풀어 강체 모형과 연성 모형의 중간에 해당한다. Drake는 점 접촉에 헌트-크로슬리 형 소산을 두고, 접촉면 전체에 압력을 적분하는 하이드로엘라스틱 접촉을 더한다. 이 장의 계산 결과를 미리 적으면 다음과 같다. ① 레믹 알고리즘이 마찰 없는 LCP 500개와 마찰 접촉 LCP 300개를 광선 종료 없이 풀었다. ② 내접 다면체 원뿔의 슬립 임계가 모서리 중점 방향에서 $\cos(\pi/m_d)$배(4변 $0.7071$)로 줄고 미끄러짐 방향이 $\pi/m_d$까지 어긋난다. ③ 상자 모서리 4점의 법선 힘은 풀이기마다 다르게 나오며 MuJoCo의 피라미드 모형은 최소 노름 해와 $2\times10^{-4}$ N 이내로 같다. ④ 2R 팔의 말단 충격은 유효 질량 $\Lambda_n$이 $0.597\sim20.06$ kg으로 변하고 Pinocchio impulseDynamics와 $3.2\times10^{-14}$ N·s 안에서 같다. ⑤ 헌트-크로슬리 반발계수가 $\alpha v_0\le0.04$에서 $1-\alpha v_0$와 $1.6\times10^{-3}$ 이내이고 켈빈-포크트는 접촉 종료 때 인장력 $-1108$ N을 낸다. ⑥ 하이드로엘라스틱 구-평면 접촉력이 $\pi E_h\delta^2(1-2\delta/3R)$이고 Drake의 계산이 메시를 세분할수록 이 식에 수렴한다($\delta=4$ mm에서 $0.8420\to0.9974$).

3.12.3 원리

3.12.3.1 접촉 기하와 Delassus 연산자

접촉 $k$의 부호 있는 거리를 $\phi_k(q)$라 하면 법선 방향 상대 속도는 $\dot\phi_k=N_kv$이고 $N_k$는 $1\times n_v$ 접촉 야코비안 행이다. 접선 평면의 단위 방향 $d_{k,i}$에 대한 속도 행을 $D_{k,i}$로 적는다. 접촉력은 일반화 힘 $N_k^{T}f_{n,k}+\sum_iD_{k,i}^{T}\beta_{k,i}$로 운동방정식 $M(q)\dot v+h(q,v)=\tau+\cdots$(3.4)에 들어간다. 접촉 행을 쌓은 $G$에서 $A=GM^{-1}G^{T}$를 Delassus 연산자라 부른다(MuJoCo 문서와 Pinocchio의 DelassusOperatorDense가 사용하는 이름이며 이 장은 이 이름을 그대로 사용한다). 이것은 접촉 공간에서의 역관성이다. 단일 법선 접촉의 $A_{nn}=N_kM^{-1}N_k^{T}$의 역수 $\Lambda_n=1/A_{nn}$가 그 접촉에서의 유효 질량(작업공간 관성의 법선 성분)이며 충격량과 충돌 안전 해석의 핵심 양이다(3.4, 26.9).

3.12.3.2 법선 상보성과 쿨롱 마찰 원뿔

법선 방향의 규정은 상보성이다.

$$0\le\phi\ \perp\ f_n\ge0,\qquad \phi=0\ \Rightarrow\ 0\le\dot\phi\ \perp\ f_n\ge0.$$

마찰은 쿨롱 법칙이다(MR 12.2.1절). 접선 힘 $f_t$는 $\lVert f_t\rVert\le\mu f_n$을 만족하고 고착($v_t=0$)이면 원뿔 내부 어디든 가능하며 미끄러지면 $f_t=-\mu f_n\,v_t/\lVert v_t\rVert$로 원뿔 표면에 있다. 두 경우를 하나의 식으로 나타내면 최대 소산 원리 $f_t=\arg\min_{\lVert f\rVert\le\mu f_n}f^{T}v_t$이다. 법선 힘이 0이면 마찰력의 상한도 0이다. MR은 접촉 모드(고착·미끄러짐·이탈)를 열거해 힘과 가속도가 일관된 모드를 찾는 방법을 제시하고, 이 모형이 해가 없거나(불일치) 둘 이상인(모호) 문제를 만들 수 있다고 적는다. 시간 간격 이산화는 이 어려움을 피해 항상 해가 있는 LCP로 바꾼다.

3.12.3.3 다면체 원뿔과 시간 간격 LCP

원뿔을 $m_d$개 방향의 내접 다면체(꼭짓점이 $\mu d_i$)로 바꾸면 마찰이 선형 부등식이 된다. $m_d$를 키울수록 원뿔에 가까워지고, 내접이므로 마찰 용량을 과소평가한다(MR 12.2.1절). 시간 간격 $h$ 동안의 충격량 $c_n=hf_n$, 마찰 충격량 $\beta_i\ge0$(방향 $d_i$)으로 속도는 $v^+=v^*+M^{-1}(N^{T}c_n+D^{T}\beta)$, 자유 속도는 $v^*=v+hM^{-1}(\tau-h(q,v))$이다. 상보성 변수 $z=(c_n,\beta,\lambda)$와 $w=Az+b$는 다음이다.

$$0\le z\ \perp\ w\ge0,\quad A=\begin{pmatrix}NM^{-1}N^{T}&NM^{-1}D^{T}&0\\ DM^{-1}N^{T}&DM^{-1}D^{T}&E\\ \operatorname{diag}(\mu)&-E^{T}&0\end{pmatrix},\quad b=\begin{pmatrix}Nv^*+\phi/h\\ Dv^*\\ 0\end{pmatrix}.$$

첫 행은 법선 상보성(틈 $\phi/h$를 포함)이고, 둘째 행은 슬립 속도의 하한 $Dv^++E\lambda\ge0$(여기서 슬랙 $\lambda_k$는 접촉 $k$의 다면체 의미 슬립 속도), 셋째 행은 원뿔 조건 $\mu c_n-E^{T}\beta\ge0$이다. $E$는 접촉마다 방향 수만큼 1을 모은 행렬이다. 이 구조는 Stewart·Trinkle과 Anitescu·Potra가 발표한 시간 간격 LCP의 구조를 따라 이 장이 다시 세운 것이며 원문 식과의 일치는 확인하지 못했다 (확인 필요). 대신 이 장은 이 LCP의 결과가 쿨롱 법칙의 닫힌 식(슬립 임계)과 MuJoCo의 거동에 맞는지를 아래에서 확인한다.

3.12.3.4 LCP의 성질과 풀이

마찰이 없으면 $A=NM^{-1}N^{T}$는 반정부호이고 LCP는 볼록 QP $\min_{c\ge0}\tfrac12c^{T}Ac+b^{T}c$의 최적 조건과 같다(1.13). 이때 사후 속도 $v^+$는 유일하다(운동 에너지 노름으로 $v^*$에서 가장 가까운 허용 속도). 접촉 수가 자유도 수보다 많으면 $A$가 특이하고 접촉 충격량 $c$는 유일하지 않다. 마찰이 있으면 $A$는 대칭이 아니어서 QP와 같지 않고, 이 장은 레믹 피벗(사전식 비율 검정)으로 푼다. 해의 존재는 Anitescu·Potra가 풀 수 있는 LCP로 정식화했다고 알려져 있으나 증명은 읽지 못했다 (확인 필요). 이 장은 임의 LCP 300개에서 광선 종료가 없었음을 확인했을 뿐 존재 증명이 아니다. 투영 가우스-자이델(PGS)은 대칭 $A$에서 한 변수씩 $x_i\leftarrow\max(0,x_i-w_i/A_{ii})$로 갱신하며 구현이 간단하나 수렴이 느리고 풀이가 시작값과 순서에 의존한다. MuJoCo는 엄격한 상보성을 비용으로 바꾼 볼록 문제 $\min\tfrac12\lambda^{T}(A+R)\lambda+\lambda^{T}(a_u-a_{ref})$를 $\lambda\in\Omega$에서 푼다(문서 식 (10)). $R$은 대각 양정부호이며 $R_{ii}=\tfrac{1-d_i}{d_i}\hat A_{ii}$로 정해지므로 문제가 엄격히 볼록해져 해가 유일하다. 문서는 이 모형에서 법선 속도와 힘이 동시에 양일 수 있고 마찰이 최대 소산이 아닐 수 있음을 인정하며 소프트 접촉에서는 오히려 물리적이라고 설명한다.

3.12.3.5 충격과 반발 법칙

충돌 순간의 속도 불연속을 충격량으로 다룬다. 법선 접촉 하나에서 충격량 $c\ge0$이면 $v^+=v^-+M^{-1}N^{T}c$이다. 뉴턴 반발 법칙 $Nv^+=-eNv^-$에서 닫힘 속도를 $s=-Nv^->0$로 두면 다음이 나온다.

$$c=(1+e)\,\Lambda_n\,s,\qquad T^--T^+=\tfrac12(1-e^2)\,\Lambda_n\,s^2.$$

충격량은 유효 질량에 비례하고 에너지 손실은 $e=1$에서 0, $e=0$에서 $\tfrac12\Lambda_ns^2$이다. 두 물체의 질량이 아닌 접촉 방향의 유효 질량이 결정한다는 점이 로봇 팔의 충돌 안전에서 중요하다. 푸아송 법칙은 압축 충격량 $c_c=\Lambda_ns$(법선 속도가 0이 될 때까지)에 $e$를 곱한 복원 충격량 $c_r=ec_c$를 더하며, 마찰 없는 단일 접촉에서는 뉴턴 법칙과 총충격량이 같다. 마찰이 있거나 접촉이 여럿이면 뉴턴·푸아송·에너지 반발계수(Stronge, 서지만 확인)가 서로 다른 결과를 낼 수 있으며 이 장은 비교하지 않았다. 여러 접촉의 뉴턴 충격은 LCP $0\le c\ \perp\ NM^{-1}N^{T}c+(1+e)Nv^-\ge0$으로 풀며 이 장의 impact_lcp가 이를 구현한다.

3.12.3.6 연성 접촉 — 헤르츠, 켈빈-포크트, 헌트-크로슬리

구가 강체 평면을 누를 때 헤르츠 이론의 접촉력은 $F=\tfrac43E^*\sqrt{R^*}\,x^{3/2}$이며 $1/E^*=(1-\nu_1^2)/E_1+(1-\nu_2^2)/E_2$, $1/R^*=1/R_1+1/R_2$이다(식은 서지만 확인, 이 장은 정의로 둔다). 질량 $m$이 속도 $v_0$로 접근하는 에너지 보존 충돌에서 최대 침투와 접촉 시간은 다음이다.

$$\delta_{\max}=\Bigl(\frac{5mv_0^2}{4k}\Bigr)^{2/5},\qquad T_c=c_T\,\frac{\delta_{\max}}{v_0},\quad c_T=2\int_0^1\frac{dx}{\sqrt{1-x^{5/2}}}=2.9433.$$

$k=\tfrac43E^*\sqrt{R^*}$이다. 감쇠를 더한 가장 단순한 모형이 켈빈-포크트(선형 스프링-댐퍼) $F=kx+c\dot x$이다. 감쇠비 $\zeta=c/(2\sqrt{km})$일 때 $e=\exp(-\pi\zeta/\sqrt{1-\zeta^2})$이고 접촉 시작에서 $F=cv_0$로 힘이 불연속하며 접촉이 끝날 때 $x\to0$이면서 $\dot x<0$이므로 $F=c\dot x<0$인 인장력이 나온다. 헌트-크로슬리 형식은 감쇠를 침투의 거듭제곱에 곱해 이 결함을 없앤다.

$$F=k\,x^{3/2}\Bigl(1+\tfrac32\,\alpha\,\dot x\Bigr),\qquad x>0.$$

$\alpha$의 단위는 s/m이고 접촉 시작과 끝($x=0$)에서 $F=0$이다. 작은 $\alpha v_0$에서 반발계수는 $e\approx1-\alpha v_0$이며 이 관계는 풀이 예제 ⑤의 계산으로 확인한 것이다. Drake는 점 접촉에 $f=f_e(x)(1-d\,v_n)_+$ 형식(탄성력에 감쇠 인자를 곱함, $v_n$은 접근할 때 음수, $d$는 s/m)을 사용하며 반발 속도가 $1/d$로 제한된다고 문서에 적는다. 연성 접촉의 적분은 강성 때문에 어렵다. 침투 최대점의 접선 강성 $k_{\mathrm{eff}}=\tfrac32k\sqrt{\delta_{\max}}$에서 각진동수 $\omega=\sqrt{k_{\mathrm{eff}}/m}$이고, 명시적 적분기의 안정 한계는 대체로 $\Delta t\lesssim2/\omega$이다.

3.12.3.7 하이드로엘라스틱 접촉

점 접촉은 "최대 침투" 하나로 힘을 정하므로 접촉 면적과 압력 분포를 알 수 없다(Drake 문서). 하이드로엘라스틱 모형은 변형을 풀지 않고 물체 내부에 압력장 $p(x)$를 미리 정한다. 표면에서 0이고 안쪽으로 갈수록 커진다. 두 물체가 겹칠 때 두 압력장이 같은 면(평형 면) $\Gamma$를 접촉면으로 잡고 압력을 적분해 합력과 모멘트를 얻는다.

$$F=\int_\Gamma p\,\hat n\,dA,\qquad M=\int_\Gamma r\times p\,\hat n\,dA.$$

한 물체가 강체이면 접촉면은 강체의 표면에 놓인다. 접촉 패치의 가장자리에서 압력이 0이므로 합력은 상대 자세에 연속이다(Drake 문서). 소산은 압력에 $(1-d\,v_n)_+$를 곱해 응력 수준에서 넣는다. 압력장이 $p=E_h(1-r/R)$인 반지름 $R$의 구가 강체 평면에 깊이 $\delta$로 눌릴 때 평면은 중심에서 $R-\delta$ 떨어지고 패치 반지름은 $a=\sqrt{2R\delta-\delta^2}$이다. 평면 위 점의 중심 거리가 $r=\sqrt{(R-\delta)^2+\rho^2}$이므로 $\rho\,d\rho=r\,dr$로 치환하면 닫힌 식이 나온다.

$$F=2\pi E_h\int_{R-\delta}^{R}\Bigl(1-\frac rR\Bigr)r\,dr=\pi E_h\,\delta^2\Bigl(1-\frac{2\delta}{3R}\Bigr).$$

힘이 $\delta^2$에 비례하고(헤르츠는 $\delta^{3/2}$) 패치 반지름이 $\sqrt{2R\delta}$로 헤르츠의 $\sqrt{R\delta}$보다 $\sqrt2$배 크다. 하이드로엘라스틱 계수 $E_h$는 재료 성질이 아니라 재료와 기하의 결합이며 영률에서 시작하되 단단한 재료에서는 $10^5\sim10^7$ Pa의 훨씬 작은 값을 사용한다(Drake 문서, 강의 영률은 약 $2\times10^{11}$ Pa). 이 모형의 근거 문헌은 Elandt 외(2019)와 Masterjohn 외(2022)이며 이 장은 원문을 읽지 못했다.

3.12.4 알고리즘

시간 간격 접촉 풀이(q, v, τ, h):
  충돌 검사 → 접촉 {φ_k, N_k, D_k(방향 m_d개), μ_k}          (11.2)
  v* ← v + h M⁻¹ (τ − h(q, v))
  A, b ← 03-03 식으로 구성                                     (Delassus 연산자 W = G M⁻¹ Gᵀ)
  z ← 레믹(A, b)            (광선 종료이면 PGS·타원 원뿔 솔버로 대체)
  v⁺ ← v* + M⁻¹ (Nᵀ c_n + Dᵀ β);   q ← q ⊕ h v⁺
충격(v⁻, e):   c ← 레믹(N M⁻¹ Nᵀ, (1+e) N v⁻);   v⁺ ← v⁻ + M⁻¹ Nᵀ c
연성 충돌(m, F, v0):   RK4로 x'' = −F(x, x')/m;  x가 다시 0이 되면 e ← −x'/v0
하이드로엘라스틱(E_h, R, δ):   F ← ∫ p dA   (구-평면은 π E_h δ² (1 − 2δ/3R))

3.12.5 풀이 예제 — LCP 풀이기, 피라미드 이방성, 힘의 비유일성, 팔의 충격, 연성 충돌, 하이드로엘라스틱

결론부터 적으면 다음과 같다. ① 레믹 풀이기가 500+300개 LCP를 풀고 NNLS와 $9.5\times10^{-13}$ 안에서 같다. ② 피라미드의 슬립 임계가 닫힌 식 $\cos(\pi/m_d)$와 $7\times10^{-8}$ 안에서 같고 MuJoCo 피라미드가 같은 이방성을(타원은 이방성이 없음을) 보인다. ③ 상자 모서리 힘이 풀이기마다 다르고 MuJoCo가 최소 노름 해와 같다. ④ 팔 충격이 Pinocchio와 같고 유효 질량이 최대 $20.06$ kg이다. ⑤ 연성 충돌의 $e$·접촉 시간이 SciPy와 같고 켈빈-포크트의 인장력을 확인했다. ⑥ 하이드로엘라스틱 힘이 SciPy quad와 $2\times10^{-15}$, Drake와 $0.1\sim2\,\%$ 안에서 같다.

구현. lemke·pgs·contact_lcp·impact_lcp·f_hc·f_kv·collide·hydro_sphere·pyramid_cap. 아래 코드는 1.1의 matmul·transpose를 사용한다. 임의 LCP 생성, 질점·상자·2R 팔의 입력, 그림은 스크립트에 있다.

from math import sqrt, sin, cos, pi

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def lemke(M, q, maxit=5000):
    """LCP: w = Mz + q >= 0, z >= 0, w'z = 0 (레믹 피벗)"""
    n = len(q)
    if min(q) >= 0:
        return [0.0] * n
    T = [[float(i == j) for j in range(n)]
         + [-M[i][j] for j in range(n)] + [-1.0, q[i]]
         for i in range(n)]
    basis, art = list(range(n)), 2 * n
    r = min(range(n), key=lambda i: q[i])
    col = art
    for _ in range(maxit):
        p = T[r][col]
        T[r] = [a / p for a in T[r]]
        for i in range(n):
            f = T[i][col]
            if i != r and f != 0.0:
                T[i] = [a - f * b for a, b in zip(T[i], T[r])]
        out, basis[r] = basis[r], col
        if out == art:
            z = [0.0] * n
            for i, b in enumerate(basis):
                if n <= b < art:
                    z[b - n] = T[i][-1]
            return z
        col = out + n if out < n else out - n
        cand = [i for i in range(n) if T[i][col] > 1e-11]
        if not cand:
            return None                      # 광선 종료
        key = lambda i: [T[i][j] / T[i][col]
                         for j in [-1] + list(range(n))]
        best = min(key(i) for i in cand)
        tol = 1e-9 * (1 + abs(best[0]))
        tie = [i for i in cand
               if key(i)[0] - best[0] < tol]
        r = next((i for i in tie if basis[i] == art), None)
        if r is None:
            r = min(tie, key=key)
    return None

def pgs(A, b, it=200):
    """투영 가우스-자이델: w = Ax + b >= 0, x >= 0 (대칭 A)"""
    n, x = len(b), [0.0] * len(b)
    for _ in range(it):
        for i in range(n):
            s = b[i] + sum(A[i][j] * x[j] for j in range(n))
            x[i] = max(0.0, x[i] - s / A[i][i])
    return x

def contact_lcp(Minv, N, D, mu, h, phi, vs):
    """스튜어트-트링클 꼴: z = (c_n, beta, lambda)"""
    nc, md = len(N), len(D) // len(N)
    G = N + D
    W = matmul(matmul(G, Minv), transpose(G))
    ng, nz = nc * (1 + md), nc * (2 + md)
    A = [[0.0] * nz for _ in range(nz)]
    for i in range(ng):
        A[i][:ng] = W[i]
    for k in range(nc):
        for i in range(md):
            A[nc + k * md + i][ng + k] = 1.0
            A[ng + k][nc + k * md + i] = -1.0
        A[ng + k][k] = mu[k]
    b = [dot(g, vs) for g in G]
    b = [b[k] + phi[k] / h for k in range(nc)] + b[nc:]
    z = lemke(A, b + [0.0] * nc)
    if z is None:
        return None
    imp = [dot(r, z) for r in transpose(G)]
    vp = [a + dot(r, imp) for a, r in zip(vs, Minv)]
    return z, vp

def impact_lcp(Minv, N, v, e):
    """뉴턴 반발: 0 <= c, N v+ + e N v- >= 0 (c 상보)"""
    A = matmul(matmul(N, Minv), transpose(N))
    b = [(1 + e) * dot(n, v) for n in N]
    c = lemke(A, b)
    imp = [dot(r, c) for r in transpose(N)]
    return c, [a + dot(r, imp) for a, r in zip(v, Minv)]

def f_hc(k, a, x, v):
    """헌트-크로슬리: x 침투, v 접근 속도 (n = 3/2)"""
    if x <= 0:
        return 0.0
    return max(0.0, k * x ** 1.5 * (1 + 1.5 * a * v))

def f_kv(k, c, x, v):
    """켈빈-포크트: 선형 스프링-댐퍼 (인장 허용)"""
    return k * x + c * v if x > 0 else 0.0

def collide(m, F, v0, dt):
    """RK4로 x'' = -F(x, x')/m, x(0)=0, x'(0)=v0 -> 이력, e"""
    s, hist = (0.0, v0), [(0.0, v0, F(0.0, v0))]
    f = lambda s: (s[1], -F(s[0], s[1]) / m)
    while True:
        k1 = f(s)
        k2 = f((s[0] + dt / 2 * k1[0], s[1] + dt / 2 * k1[1]))
        k3 = f((s[0] + dt / 2 * k2[0], s[1] + dt / 2 * k2[1]))
        k4 = f((s[0] + dt * k3[0], s[1] + dt * k3[1]))
        t = tuple(s[i] + dt / 6 * (k1[i] + 2 * (k2[i] + k3[i])
                                   + k4[i]) for i in (0, 1))
        if t[0] <= 0.0:
            w = s[0] / (s[0] - t[0])
            vend = s[1] + w * (t[1] - s[1])
            return hist, -vend / v0, (len(hist) + w) * dt
        s = t
        hist.append((s[0], s[1], F(s[0], s[1])))

def hydro_sphere(E, R, d):
    """압력장 p=E(1-r/R) 구와 강체 평면: 힘, 패치 반지름"""
    F = pi * E * d * d * (1 - 2 * d / (3 * R))
    return F, sqrt(2 * R * d - d * d)

def hydro_quad(E, R, d, n=4000):
    """패치 위 압력 적분 F = int p dA (중점 구적)"""
    dr = sqrt(2 * R * d - d * d) / n
    p = lambda r: E * (1 - sqrt((R - d) ** 2 + r * r) / R)
    return sum(p((i + .5) * dr) * 2 * pi * (i + .5) * dr * dr
               for i in range(n))

def pyramid_cap(mu, md, psi):
    """내접 md각 피라미드의 방향 psi 마찰 용량 / (mu N)"""
    w = 2 * pi / md
    off = (psi % w) - w / 2
    return cos(w / 2) / cos(off)

① LCP 풀이기.

레믹 풀이기는 계수가 부족한 마찰 없는 LCP를 포함해 800개 모두 광선 종료 없이 풀었다. 상보성 잔차가 $10^{-12}$ 이하이고 마찰 접촉 LCP에서 원뿔 조건과 관통 조건 위반이 없다. 광선 종료가 처음에는 500개 중 2개에서 나타났다. 원인은 사전식 비율 검정의 동률 허용오차가 너무 작아 거의 같은 비율 두 행 중 잘못된 행을 고른 것이었고 상대 허용오차 $10^{-9}$로 바꾸자 사라졌다. 따라서 이 풀이기는 수치 허용오차에 민감하며 해의 존재 증명은 아니다.

② 다면체 마찰 원뿔의 슬립 임계.

슬립 임계는 질점에 수평 힘을 늘려 가며 이분법으로 구했다. 꼭짓점 방향에서는 $\mu mg$이고 모서리 중점 방향에서는 닫힌 식 $\cos(\pi/m_d)\mu mg$와 $7\times10^{-8}$ 안에서 같다(4변 $0.7071$, 8변 $0.9239$, 16변 $0.9808$, 32변 $0.9952$). 슬립이 시작된 직후 미끄러짐 방향은 힘의 방향과 최대 $\pi/m_d$까지 어긋난다(4변 $42.75^\circ$, 상한 $45^\circ$). 미끄러짐 속도가 근처 모서리의 법선 방향으로 정렬되기 때문이며 정확한 쿨롱 법칙에서는 속도가 알짜 힘과 같은 방향이다. 다면체가 방향에 따라 마찰 용량이 다른 이방 마찰이 되는 것이다.

③ 강체 접촉력의 비유일성.

네 모서리의 법선 힘은 합 $mg$와 모멘트 평형을 만족하는 해가 무수히 많다. 접촉 행렬 $A$의 특이값(numpy svd)은 $5.48,\ 4.00,\ 2.92,\ 3.9\times10^{-16}$이어서 계수 3이다(바닥 평면의 접촉이 세 자유도를 제어하는데 접촉점은 4개). 세 해는 모두 상보성 잔차 $10^{-17}$ 이하의 유효한 LCP 해이며 사후 속도는 같다. 레믹 해는 한 모서리의 힘이 0인 꼭짓점 해이고 PGS 해는 시작값과 갱신 순서에 의존하며 최소 노름 해는 모서리 4개 모두에 양의 힘을 준다. 강체 접촉은 이 선택을 모형이 정하지 못하므로 풀이기가 정한다.

④ 2R 팔의 충격과 유효 질량.

벽의 법선 방향에서 말단의 유효 질량 $\Lambda_n=1/(N M^{-1}N^{T})$이 충격량을 정한다. 팔이 거의 펴진 $q_2=10^\circ$에서 반경 방향 벽은 $\Lambda_n=20.06$ kg으로 팔 전체 질량 4 kg의 5.0배이고 충격량이 $26.07$ N·s로 커진다. 같은 자세에서 접선 방향 벽은 $0.597$ kg이다. 말단이 벽에 닫힘 속도 1 m/s로 닿아도 방향에 따라 손실 에너지가 $0.27$ J에서 $9.13$ J까지 달라진다. 충돌 직후 법선 속도는 모든 경우 $e=0.3$ m/s이다. 이 자세 의존성이 충돌 안전 해석에서 유효 질량을 따로 보는 이유이다(26.9, 27.10의 ISO/TS 15066 본문은 읽지 못했다 (확인 필요)). 이 표에서 반경 방향의 $q_2$와 접선 방향의 $180^\circ-q_2$ 값이 같은 것은 링크 길이·질량이 같은 이 팔의 대칭 때문으로 보며 일반 성질로 확인하지 않았다.

⑤ 헤르츠·헌트-크로슬리 충돌.

탄성($\alpha=0$) 충돌의 $e$는 $1$이고 최대 침투와 접촉 시간이 닫힌 식과 같다(접촉 시간 상수 $2.9433$). 접촉 시간은 $77\sim169$ µs이다. 명시적 적분기가 안정하려면 $\Delta t$가 $38\sim84$ µs 이하여야 하므로 1 ms 간격의 명시적 시뮬레이션으로는 이 강철 충돌의 접촉 구간을 해상할 수 없다. 헌트-크로슬리 반발계수는 $\alpha v_0\le0.04$에서 $1-\alpha v_0$와 $1.6\times10^{-3}$ 이내이고 $\alpha v_0=0.4$에서는 $0.712$로 $0.6$과 $0.11$ 차이가 난다(1차 근사의 범위가 좁다). 반발계수가 속도에 의존하는 것이 헌트-크로슬리의 특징이며 켈빈-포크트의 $e$는 속도와 무관하다.

세 모형의 힘-침투 궤적에서 헌트-크로슬리는 접촉 시작과 끝에서 힘이 0이고 폐곡선이 닫힌다. 켈빈-포크트는 접촉 직후 $1563.8$ N(자기 최대 $6333.7$ N의 25 %)로 힘이 불연속하고 접촉이 끝날 때 $-1108.2$ N의 인장력을 낸다. 인장력을 0으로 제한하면(인장 제거) 접촉 종료 시점이 앞당겨져 $e$가 $0.7120$에서 $0.7288$로 달라진다. 헌트-크로슬리의 $e$가 켈빈-포크트와 같도록 $\zeta$를 맞췄으므로($0.7120$) 같은 $e$를 내는 두 모형의 힘 이력이 이처럼 다르다. 헤르츠 탄성의 최대 힘은 $9060.5$ N이다.

⑥ 하이드로엘라스틱 구-평면.

패치 위 압력의 수치 적분이 닫힌 식과 $10^{-7}$ 이내로 같다. 같은 계수 $E_h=E^*$로 두면 하이드로엘라스틱 힘이 헤르츠의 $0.23$배($\delta=0.5$ mm)에서 $0.84$배($8$ mm)로 변한다. 지수가 다르므로 두 모형을 한 침투 깊이에서 맞추면 다른 깊이에서 어긋난다.

⑦ 라이브러리 대조(실행 2026-11-05). SciPy 1.18.1의 nnls를 마찰 없는 LCP(제곱 노름 최소화 동치)와, solve_ivp(DOP853)를 헌트-크로슬리·켈빈-포크트 충돌과, quad를 접촉 시간 상수와 하이드로엘라스틱 패치 적분과, Pinocchio 4.1.0의 crba·computeFrameJacobian·impulseDynamics를 2R 팔의 질량행렬·야코비안·충격과, MuJoCo 3.17.0을 피라미드·타원 원뿔의 이탈 힘·상자 모서리 힘·미끄럼 감속과, Drake 1.58.0(kHydroelastic, 구는 순응, 평면은 강체)을 하이드로엘라스틱 구-평면 힘·패치 면적과 같은 입력에서 비교했다. 스크립트는 robotics/fig-src/3.12-libcheck.py(Drake는 별도 가상환경에서 drake 인자로 먼저 실행), 출력은 robotics/fig-src/data/3.12-libcheck-output.txt에 있다.

SciPy·Pinocchio 대조는 모두 $2\times10^{-10}$ 이하에서 같다. MuJoCo의 피라미드 원뿔은 이 장의 이방성을 재현하지만 이탈 힘의 비가 닫힌 식보다 $0.9\,\%$ 크다($0.7723$ 대 $0.7654$, $0.7136$ 대 $0.7071$). 이탈 힘을 속도 $0.05$ m/s로 판정했으므로 램프 오버슈트가 포함되며(꼭짓점 방향에서 $1.0194$) 이 오차의 일부다. 램프 중 고착 구간에서도 mm/s 수준의 느린 미끄러짐(크리프)이 관측되었으며 이는 MuJoCo의 소프트 구속의 결과로 해석한다(이 장은 원인을 분리하지 않았다). 타원 원뿔은 방향에 무관하게 $0.9987\sim0.9993$으로 정확한 원과 같다. 미끄럼 감속도는 피라미드 $4.875$·타원 $4.889$ m/s²로 $\mu g=4.905$ m/s²보다 $0.6\,\%$·$0.3\,\%$ 작다. 모서리 힘 4개에서 MuJoCo의 피라미드는 최소 노름 해와 $2\times10^{-4}$ N, 타원은 $6.6\times10^{-3}$ N 안에서 같다. MuJoCo 문서에 따르면 정규화 $R_{ii}$의 대각 근사는 초기 자세의 말단 관성에서 만들므로 모서리마다 같은 값이 되어 균일한 정규화의 극한(최소 노름 해)이 나온다고 설명할 수 있으나 이 해석은 이 장의 추정이다. Drake의 하이드로엘라스틱 힘은 메시 해상도를 $10\to5\to2.5\to1.25$ mm로 줄일 때 해석식의 $0.8420\to0.9595\to0.9898\to0.9974$로 수렴하며 오차가 해상도마다 약 $1/4$로 줄어 메시 크기의 제곱에 비례하는 것과 일치한다. 이 결과로 구의 압력장이 $p=E_h(1-r/R)$이고 힘이 $\delta^2$ 법칙을 따른다는 이 장의 식이 Drake의 구현과 같다고 본다.

한계. 질점·상자·2R 팔·구-평면의 단일 또는 소수 접촉이며 접촉 검출·다물체 충돌 순서·구름 마찰·비틀림 마찰·마찰 충격(뉴턴 대 푸아송)은 다루지 않았다. 시간 간격 $h$에 대한 수렴과 접촉 파라미터의 실물 식별도 시험하지 않았다. 레믹 풀이기는 변수 40개까지만 시험했고 대규모 접촉에는 PGS·ADMM 등이 필요하다(16.1).

견해: 접촉이 많은 시뮬레이션과 모델 기반 제어에는 MuJoCo형 볼록 소프트 접촉(또는 Drake의 볼록 근사)을 기본으로 두고, 엄격한 LCP는 충격량이 중요한 소수 접촉의 해석·검증에 사용한다고 본다. 근거는 이 장에서 중복 접촉의 강체 법선 힘이 풀이기마다 다른 반면 소프트 모형은 하나의 분포(최소 노름 해와 $2\times10^{-4}$ N 이내)를 정한다는 점, 다면체 근사가 마찰 이방성($0.7071$)을 만드는 반면 타원 원뿔은 방향에 무관하다는 점이다. 강체 접촉의 사실성(마찰 최대 소산·이탈 거동)이 중요한 과제에서는 달라질 수 있으며 이 장은 실물 접촉 시험을 하지 않았다.

3.12.6 설계 트레이드오프 — 접촉 모형의 선택

접촉을 모형화하는 방식은 다음 중에서 고른다.

• ① 강체 상보성(LCP·시간 간격). 침투가 없고 충격이 정의된다. 중복 접촉에서 힘이 비유일하고 다면체 근사가 이방성을 만들며 마찰 최대 소산이 비볼록이다.
• ② 연성 점 접촉(스프링-댐퍼, 헤르츠·헌트-크로슬리). 침투의 함수로 힘을 정하므로 풀이가 상미분방정식이고 반발계수를 감쇠로 조정한다. 강성이 크면 적분 간격이 작아지고 면 접촉의 압력 분포를 알 수 없다.
• ③ 볼록 소프트(MuJoCo)·하이드로엘라스틱(Drake). 정규화로 해가 유일하며 후자는 접촉 면적과 압력 분포·모멘트를 준다. 소프트 파라미터가 물리량이 아니며 실물 보정이 필요하다.

견해: 로봇 학습용 시뮬레이션처럼 수천 개 접촉을 빠르게 풀어야 하면 ③ 중 볼록 소프트를, 평평한 접촉면(발바닥·상자 바닥·평면 손가락)의 모멘트가 중요한 해석에서는 하이드로엘라스틱을 우선 검토한다고 본다. 근거는 이 장에서 중복 접촉의 힘 분포를 정하는 것이 모형의 정규화이고(③) 구-평면에서 하이드로엘라스틱이 패치 면적과 힘을 닫힌 식으로 주며 Drake와 $0.3\,\%$ 이내(메시 1.25 mm, $\delta=4$ mm)로 같다는 점이다. 하이드로엘라스틱 계수의 실물 보정은 이 장이 다루지 않았고 소프트 접촉의 크리프가 정밀 조작에 미치는 영향도 시험하지 않았다.

3.12.7 응용 축 비교

이 표는 이 책의 정리이며 응용별 수치 요구는 이 장이 확인하지 않았다.

3.12.8 관련 규격·프로토콜

ISO/TS 15066의 본문과 부위별 한계값은 읽지 못했다 (확인 필요). 접촉 파라미터의 단위(소프트 구속의 시간 상수·임피던스, 하이드로엘라스틱 계수 Pa, 소산 s/m)는 엔진마다 다르므로 모델 이식 때 대응을 확인해야 한다.

3.12.9 관련 규제·인증 요구와의 연결

접촉·충돌 동역학 자체를 규정하는 법령·인증 요건은 확인하지 못했다. 협동 운전에서 인체와 로봇의 충돌은 접촉 힘·압력의 한계로 평가되며 유효 질량과 접근 속도가 그 입력이다(이 장은 $\Lambda_n$과 충격량의 식을 제공한다). 한계값과 측정법은 26.9, 인증 현장 시험은 41.9, 협동 운전 방식은 28.4과 이어진다. 이 장은 규정의 조항을 확인하지 않았다 (확인 필요).

3.12.10 한계·유의점·미해결

• 모형의 범위. 구름·비틀림 마찰(MuJoCo condim 4·6), 마찰 충격, 다접촉 동시 충격, 점탄성·접착은 다루지 않았다.
• 원문 미확인. LCP의 원문 식·해 존재 증명, 헌트-크로슬리·헤르츠의 원문 형식, 하이드로엘라스틱 원논문의 압력장 정의는 읽지 못했다.
• 존재와 유일성. 이 장의 레믹 풀이기가 광선 종료 없이 푼 것은 시험한 표본에서의 사실이다. 마찰 있는 강체 쿨롱 모형은 연속 시간에서 해가 없거나 여럿일 수 있다(MR).
• 다면체 근사. 변 수를 늘리면 이방성이 줄지만 LCP 크기가 커진다. MuJoCo 문서는 소프트 접촉에서는 변 수를 늘려도 타원 원뿔 해로 수렴하지 않는다고 설명한다.
• 파라미터의 실물 대응. 소프트 계수(강성·감쇠·소프트 구속 시간 상수·하이드로엘라스틱 계수)는 재료 물성과 같지 않아 실험 보정이 필요하다(16.12).
• 시간 간격. $h$에 대한 LCP 해의 수렴, 접촉 순서 변화의 영향은 시험하지 않았다.
• 충돌 검사. 접촉 목록과 거리·법선의 계산은 11.2의 주제이다.

3.12.11 체크포인트·연구 질문·출처

체크포인트:

1. 단일 법선 접촉에서 뉴턴 반발 법칙으로 충격량 $c=(1+e)\Lambda_ns$와 에너지 손실 $\tfrac12(1-e^2)\Lambda_ns^2$를 유도하라.
2. 내접 $m_d$각 피라미드의 방향별 마찰 용량 $\cos(\pi/m_d)/\cos(\psi-\psi_{\text{mid}})$를 유도하고 슬립 직후 방향 오차가 $\pi/m_d$ 이하인 이유를 설명하라.
3. 마찰 없는 LCP가 QP의 최적 조건과 같음을 보이고 $A$가 특이할 때 $c$는 비유일하고 $v^+$는 유일한 이유를 설명하라.
4. 헤르츠 접촉에서 $\delta_{\max}$와 접촉 시간 $T_c=c_T\delta_{\max}/v_0$를 에너지 보존에서 유도하라.
5. 켈빈-포크트가 접촉 시작에서 힘이 불연속하고 종료에서 인장력을 내는 이유와 헌트-크로슬리가 이를 없애는 이유를 설명하라.
6. 압력장 $p=E_h(1-r/R)$인 구-평면 접촉력 $\pi E_h\delta^2(1-2\delta/3R)$를 $\rho\,d\rho=r\,dr$ 치환으로 유도하고 패치 반지름을 구하라.

연구 질문:

• 마찰 충격(뉴턴·푸아송·에너지 반발계수)을 다접촉 LCP에서 에너지 일관성을 보장하며 정식화하는 방법은 무엇인가(이 장은 시험하지 않았다)?
• 하이드로엘라스틱 계수와 소산을 실물 접촉 시험(발바닥·손가락)에서 식별하는 절차와 그 불확실성은 어떻게 정량화하는가?
• 소프트 접촉의 크리프(고착 구간의 느린 미끄러짐)를 모형 파라미터에서 어떻게 분리·제한하는가?

1차 출처(조회 2026-11-05).

• K. M. Lynch, F. C. Park, Modern Robotics, Cambridge University Press, 2017(사전 인쇄본 2017-05-03) — 12.2.1·12.3절을 직접 읽었다.
• MuJoCo 공식 문서 "Computation"(mujoco.readthedocs.io, stable) — 직접 읽었다. E. Todorov, "Convex and analytically-invertible dynamics with contacts and constraints: Theory and implementation in MuJoCo," ICRA 2014, 6054–6061, doi:10.1109/icra.2014.6907751 — 서지만 확인. E. Todorov, T. Erez, Y. Tassa, "MuJoCo: A physics engine for model-based control," IROS 2012, 5026–5033, doi:10.1109/iros.2012.6386109 — 서지만 확인.
• Drake 문서 "Hydroelastic Contact User Guide"·"Modeling Compliant Contact"(drake.mit.edu, 조회 2026-11-05) — 직접 읽었다. R. Elandt, E. Drumwright, M. Sherman, A. Ruina, "A pressure field model for fast, robust approximation of net contact force and moment between nominally rigid objects," IROS 2019, 8238–8245, doi:10.1109/iros40897.2019.8968548 — 서지만 확인. J. Masterjohn 외, "Velocity level approximation of pressure field contact patches," IEEE RA-L 7(4), 11593–11600, 2022, doi:10.1109/lra.2022.3203845 — 서지만 확인. A. Castro, F. Permenter, X. Han, "An unconstrained convex formulation of compliant contact," IEEE T-RO 39(2), 1301–1320, 2023, doi:10.1109/tro.2022.3209077 — 서지만 확인.
• D. E. Stewart, J. C. Trinkle, "An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction," IJNME 39(15), 2673–2691, 1996, doi:10.1002/(sici)1097-0207(19960815)39:15<2673::aid-nme972>3.0.co;2-i — 서지만 확인. M. Anitescu, F. A. Potra, "Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems," Nonlinear Dynamics 14(3), 231–247, 1997, doi:10.1023/a:1008292328909 — 서지만 확인. J. S. Pang, J. C. Trinkle, Mathematical Programming 73(2), 199–226, 1996, doi:10.1007/bf02592103 — 서지만 확인. D. E. Stewart, "Rigid-body dynamics with friction and impact," SIAM Review 42(1), 3–39, 2000, doi:10.1137/s0036144599360110 — 서지만 확인.
• H. Hertz, "Ueber die Berührung fester elastischer Körper," J. reine angew. Math. 92, 156–171, 1882, doi:10.1515/crll.1882.92.156 — 서지만 확인. K. H. Hunt, F. R. E. Crossley, "Coefficient of restitution interpreted as damping in vibroimpact," J. Appl. Mech. 42(2), 440–445, 1975, doi:10.1115/1.3423596 — 서지만 확인. R. Marhefka, D. Orin, IEEE Trans. SMC-A 29(6), 566–572, 1999, doi:10.1109/3468.798060 — 서지만 확인. W. J. Stronge, "Rigid body collisions with friction," Proc. R. Soc. Lond. A 431, 169–181, 1990, doi:10.1098/rspa.1990.0125 — 서지만 확인. K. L. Johnson, Contact Mechanics, Cambridge University Press, 1985, doi:10.1017/cbo9781139171731 — 서지만 확인.
• Pinocchio 4.1.0(impulseDynamics 설명, crba, computeFrameJacobian), SciPy 1.18.1(optimize.nnls, integrate.solve_ivp·quad), numpy 2.5.4, MuJoCo 3.17.0, Drake 1.58.0(pydrake)의 계산 대조.

다음 장. 3.13는 부유기저 다물체 동역학과 중심 동역학(centroidal dynamics)을 다룬다. 이 장의 접촉력과 마찰 원뿔이 부유기저의 운동량 변화율을 결정하는 외력이 된다.