제3부 동역학

3.8 마찰: 쿨롱·스트리벡·LuGre 모델과 관절 마찰

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

이 장은 로봇 관절의 기어헤드와 베어링에서 생기는 마찰 토크를 모델링하고 식별하며 보상하는 방법을 다룬다. 정상 상태 마찰 곡선(쿨롱·점성·스트리벡 효과)의 음의 기울기 구간, 변수 투영(VP)으로 스트리벡 파라미터를 식별할 때 소인 속도 설계가 식별 가능성과 정밀도를 어떻게 정하는지, 마찰이 상태 하나를 가지는 LuGre 동적 모형이 프리슬라이딩 변위·이탈 마찰의 속도 의존·마찰 지연 히스테리시스를 만드는 과정, 마찰이 있는 1관성 서보에서 적분 동작이 일으키는 헌팅(한계 순환)과 LuGre 논문의 관측기 기반 보상, 쿨롱·점성 모형으로 동역학을 식별할 때 스트리벡 잔차가 관성 파라미터에 주는 편의를 계산으로 확인한다. 동역학 파라미터 식별은 3.7, 운동방정식의 구조는 3.4, 유연 관절은 3.9, 접촉 마찰(쿨롱 마찰 원뿔)은 3.12, 관절·모션 제어의 마찰 보상은 12.3, 서보 루프는 5.9, 파동기어와 감속기는 5.10·5.11, 베어링·윤활은 4.2과 이어진다.

배우는 것:

• 정상 상태 마찰 곡선 — 쿨롱·점성·스트리벡, 음의 기울기 구간 $(\omega_{-},\omega_{+})$의 닫힌 식(Lambert $W$).
• 정지 마찰과 영속도 — 정적 사상이 영속도에서 다가(多價)인 이유, 쿨롱 모형의 불연속.
• LuGre 모형 — 상태 $z$ 하나의 강모 모형, 유계성·소산성, 프리슬라이딩·이탈 마찰·마찰 지연.
• 식별 — 부호 대칭 평균으로 중력·오프셋 제거, VP 최소제곱, 소인 설계, 점근 공분산.
• 보상과 헌팅 — 적분 동작이 일으키는 한계 순환, 관측기 기반 마찰 보상과 양실수(SPR) 조건.
• 동역학 식별과의 결합 — 마찰 모형 오차가 관성·중력 파라미터에 만드는 편의.

선수 지식: 3.4(운동방정식), 1.12(비선형 최소제곱), 3.7(회귀 행렬·기저 파라미터).

특히 중요한 독자: 관절 토크 지령에 마찰 보상을 넣는 제어 엔지니어, 고감속비 관절의 저속 정밀 제어를 다루는 설계자, 마찰 모형을 시뮬레이터 모델에 반영하는 개발자(16.1, 16.12).

출처 방침. 기준 시점은 2026-11-05이다. 직접 읽은 자료: Modern Robotics(Lynch·Park, 사전 인쇄본 2017-05) 8.9.4절(마찰 모형, 도 8.13의 여섯 가지 속도 의존 모형, 기어헤드 마찰이 감속비와 함께 커지는 경향)과 8.9.5절(관절·링크의 유연성); Canudas de Wit·Olsson·Åström·Lischinsky의 "A new model for control of systems with friction"(IEEE TAC 40(3), 1995, HAL 공개본 전문 — 모형식, 유계성·소산성, 선형화, 프리슬라이딩·마찰 지연·이탈 마찰·스틱-슬립·PID 헌팅 모의, 관측기 기반 보상과 정리 1); Swevers·Al-Bender·Ganseman·Prajogo의 "An integrated friction model structure with improved presliding behavior for accurate friction compensation"(IEEE TAC 45(4), 2000, KU Leuven 공개본의 초록·서론·2절 앞부분). 서지만 확인한 것(본문은 읽지 못했다): Armstrong-Hélouvry·Dupont·Canudas de Wit(1994, Automatica, 마찰 모형·분석·보상 개관), Olsson 외(1998, European J. Control, 마찰 모형과 보상), Dupont·Hayward·Armstrong·Altpeter(2002, 단일 상태 탄소성 마찰 모형), Karnopp(1985, 스틱-슬립 모의), Haessig·Friedland(1991, 마찰의 모델링과 모의), Bona·Indri(2005, 로봇 마찰 보상 개관). 스트리벡의 원 논문은 서지를 확인하지 못했다 (확인 필요). 마찰 파라미터의 값(쿨롱·정지·스트리벡 속도·점성·LuGre $\sigma_0$·$\sigma_1$)은 모두 이 장이 정한 예시이며 어떤 제품의 값이 아니다. 수치는 모두 이 장의 순수 파이썬 합성 계산이다.

3.8.2 개념 — 마찰은 속도의 함수이면서 변위의 함수이다

결론부터 적으면 다음과 같다. 관절 마찰의 가장 단순한 모형은 쿨롱(속도 부호에 비례하는 일정 크기)과 점성(속도에 비례) 마찰의 합이고 파라미터에 대해 선형이어서 동역학 회귀 행렬(3.7)의 열로 들어간다. 그러나 저속에서는 속도가 커질수록 마찰이 줄어드는 스트리벡 효과가 있어 속도 사상의 기울기가 음이 되는 구간이 생기고, 이 구간은 스틱-슬립과 헌팅의 원인이 된다. 정적 속도 사상은 영속도에서 값이 하나로 정해지지 않으며(Lynch·Park 도 8.13: 영속도에서 마찰이 $[-b_{\text{static}},b_{\text{static}}]$의 임의 값) 접촉면이 탄성 강모처럼 변형하는 프리슬라이딩(presliding) 영역에서는 마찰이 속도가 아니라 변위의 함수이다. LuGre 모형은 평균 강모 변형 $z$ 하나의 미분 방정식으로 두 영역을 하나로 묶고, 정상 상태에서는 임의의 정적 곡선(스트리벡 포함)으로 환원된다. 제어의 관점에서는 속도가 0 근처를 지나는 모든 과제(방향 전환, 저속 추종, 정지 유지)에서 정적 모형의 불연속이 오차를 만들고, 적분 동작이 있는 서보는 마찰이 정지 마찰 이상이 될 때까지 적분이 쌓였다가 미끄러지는 한계 순환(헌팅)을 일으킬 수 있다.

이 장의 계산 결과를 미리 적으면 다음과 같다. ① 예시 파라미터($F_c=1.0$, $F_s=1.5$ N·m, $\omega_s=0.01$ rad/s, $\sigma_2=0.5$ N·m·s/rad)의 정상 상태 곡선은 $\omega\in(5.0\times10^{-5},\ 0.0249)$ rad/s에서 기울기가 음이고 최솟값 $1.0135$ N·m이며, 두 근이 Lambert $W$ 닫힌 식·brentq와 $10^{-16}$ 안에서 같다. ② 스트리벡 4파라미터를 VP로 식별한 결과가 SciPy least_squares와 상대 $7\times10^{-8}$ 안에서 같고, 넓은 소인에서 점근 표준편차가 $F_c$ $0.0112$·$F_s$ $0.0227$ N·m이며 몬테카를로 비는 $0.87\sim1.05$이다. 속도를 $0.2$ rad/s 이상으로만 소인하면 $F_s$·$\omega_s$가 식별되지 않는다. ③ LuGre 프리슬라이딩에서 강모 변형이 $F/\sigma_0$와 $0.1\,\%$ 안에서 같고(변위는 3.1배), 활주 직전 마찰 최댓값이 토크 증가율 $0.1\to100$ N·m/s에서 $1.489\to1.173$ N·m로 줄며 마찰 지연 루프 면적이 주파수에 따라 늘어난다. ④ PID 서보의 헌팅이 진폭 $0.0685$ rad·주기 $23.9$ s(Kv=6)이고, SPR 조건을 만족하는 $K_v=150$에서 관측기 보상은 오차를 $10^{-19}$ rad로 줄인다(조건을 만족하지 못하는 $K_v=6$에서도 $2\times10^{-4}$ rad 진동으로 줄인다). ⑤ 쿨롱·점성 모형으로 식별한 1링크 관성이 최대 속도 $0.2$ rad/s의 궤적에서 $+91.8\,\%$ 어긋난다.

3.8.3 원리

3.8.3.1 정상 상태 마찰 곡선

일정 속도 $\omega$로 움직일 때의 마찰 토크 $F_{ss}(\omega)$를 정상 상태 마찰이라 한다. 점성 마찰 $F_v\omega$, 쿨롱 마찰 $F_c\operatorname{sgn}\omega$, 정지·운동 마찰의 차이, 스트리벡 효과는 Lynch·Park 도 8.13(a)~(f)가 여섯 모형으로 정리한다. 이 장은 LuGre 논문이 스트리벡 효과의 매개화로 소개한 지수 형태(논문의 식 (4), $\sigma_0g(v)=F_c+(F_s-F_c)e^{-(v/v_s)^2}$)를 사용한다.

$$F_{ss}(\omega)=\bigl[F_c+(F_s-F_c)\,e^{-(\omega/\omega_s)^2}\bigr]\operatorname{sgn}\omega+\sigma_2\,\omega,$$

여기서 $F_c$는 쿨롱 마찰 수준, $F_s$는 정지 마찰 수준($F_s\ge F_c$), $\omega_s$는 스트리벡 속도, $\sigma_2$는 점성 계수이다. 기울기는 $f_{ss}^\prime(\omega)=-2(F_s-F_c)\,\omega\,\omega_s^{-2}e^{-(\omega/\omega_s)^2}+\sigma_2$이고 $r=\omega/\omega_s$로 두면 $f_{ss}^\prime=0$은 $r\,e^{-r^2}=\kappa$, $\kappa=\sigma_2\omega_s/\bigl(2(F_s-F_c)\bigr)$와 같다. $r\,e^{-r^2}$의 최댓값이 $1/\sqrt{2e}\approx0.429$이므로 음의 기울기 구간은 $\kappa<0.429$일 때 존재하고, 두 근은 Lambert $W$ 함수로 다음처럼 나타낸다.

$$\omega_\pm=\omega_s\sqrt{-\tfrac12\,W_k(-2\kappa^2)},\qquad k=0\ (\omega_-),\ k=-1\ (\omega_+).$$

$\omega_-<\omega<\omega_+$에서 마찰이 속도와 함께 줄어드는 음의 감쇠가 나타나며 이 구간이 저속에서의 불안정(스틱-슬립)의 기원이다(LuGre 논문의 서론이 스트리벡 효과의 저속 불안정화 작용을 서술한다). 이 장의 예시에서는 $\kappa=0.005$이므로 구간이 $\omega_s$보다 훨씬 넓게 나타나 $\omega_+\approx2.49\,\omega_s$이다.

3.8.3.2 영속도와 정적 모형의 한계

$\operatorname{sgn}$ 함수는 $\omega=0$에서 불연속이고 정적 사상은 영속도에서 값이 하나가 아니다. 시뮬레이션에서는 영속도 근방의 작은 대역 $\lvert\omega\rvert<\epsilon$을 정지로 취급하는 방법(Karnopp 1985의 제목이 스틱-슬립 모의를 다룬다; 본문은 읽지 못했다)이나 $\tanh(\omega/\epsilon)$ 같은 평활화를 사용하는데, 전자는 정지·미끄럼 전환 논리를, 후자는 $\epsilon$에 따라 저속 오차와 수치 강성을 바꾼다. MuJoCo의 frictionloss는 일정 토크에서 종단 속도가 $(\tau-F_c)/\sigma_2$와 같고 $F_c$ 아래의 토크에서는 거의 정지하는 연성 구속처럼 거동했다(이 장은 동작만 확인했고 문서는 읽지 못했다). 정적 모형은 방향 전환 때의 지연, 정지 상태에서 힘이 작을 때의 미소 변위, 이탈 마찰의 힘 증가율 의존을 설명하지 못한다(LuGre 논문 서론).

3.8.3.3 LuGre 동적 마찰 모형

Canudas de Wit 외(1995)는 접촉면을 탄성 강모(bristle)로 연결된 두 강체로 보고 강모의 평균 변형 $z$의 미분 방정식을 제안했다. 이 장의 기호(회전 관절, 각속도 $\omega$)로 다음이다.

$$\dot z=\omega-\frac{\sigma_0\lvert\omega\rvert}{h(\omega)}\,z,\qquad h(\omega)=F_c+(F_s-F_c)e^{-(\omega/\omega_s)^2},\qquad F=\sigma_0z+\sigma_1\dot z+\sigma_2\omega.$$

$\sigma_0$는 강모 강성(N·m/rad), $\sigma_1$은 강모 감쇠, $\sigma_2$는 점성이다. 정상 상태($\omega$ 일정)에서 $z_{ss}=h(\omega)\operatorname{sgn}\omega/\sigma_0$이므로 $F=F_{ss}(\omega)$로 앞 절의 곡선과 같다. 논문이 보인 성질은 다음이다. 속성 1: $0<h/\sigma_0\le a$이고 $\lvert z(0)\rvert\le a$이면 $\lvert z(t)\rvert\le a$(강모 변형이 유계). 속성 2: 속도에서 $z$로의 사상이 $V=z^2/2$에 대해 소산적. 선형화: $z=0$·$\omega=0$ 근방에서 $\dot z\approx\omega$이므로 $F\approx\sigma_0\theta+(\sigma_1+\sigma_2)\dot\theta$이고 정지 영역은 강성 $\sigma_0$·감쇠 $\sigma_1+\sigma_2$의 스프링-댐퍼로 거동한다. 관성 $J$와 연결하면 고유 각진동수 $\sqrt{\sigma_0/J}$, 감쇠비 $(\sigma_1+\sigma_2)/(2\sqrt{\sigma_0J})$이다. 논문은 프리슬라이딩 변위, 마찰 지연(Hess·Soom의 실험에 대응하는 히스테리시스), 변하는 이탈 마찰, 스틱-슬립을 모의로 보였다. $g=F_c/\sigma_0$이고 $\sigma_1=\sigma_2=0$이면 Dahl 모형이 된다고 논문이 서술한다. Swevers 외(2000)는 LuGre 모형이 정상 상태 거동과 속도 반전의 매끄러운 전이를 잘 기술하지만 프리슬라이딩에서 소산이 지나치게 크고 변위-힘 전이 곡선의 모양이 모형 구조로 고정되며 비국소 기억 히스테리시스를 포함하지 않는다고 지적한다(같은 논문은 KUKA IR 361 로봇 1축의 프리슬라이딩 토크 측정으로 이를 보였다고 서술하며, Armstrong이 LuGre 모형의 정확한 예측이 어떤 장치에서도 입증되지 않았다고 말했다고 인용한다).

3.8.3.4 관절 마찰의 성질

관절 마찰은 기어헤드와 베어링에서 나오며 Lynch·Park(8.9.4절)은 베어링 하중, 정지 시간, 온도가 마찰에 영향을 주고 기어헤드 마찰이 감속비 $G$와 함께 커지는 경향이 있다고 서술한다. 이상 기어(효율 1)를 가정하면 모터 쪽 마찰 토크 $F_m$은 동력 균형 $F_m\omega_m=F_j\omega_j$, $\omega_m=G\omega_j$에서 관절 쪽 $F_j=GF_m$이므로 감속비가 클수록 관절 쪽의 쿨롱 수준이 커진다. 모터 전류에서 토크를 구하는 구성(6.3)에서는 이 마찰이 측정 토크에 포함되어 식별 대상이 된다. 파동기어 관절은 모터 회전자와 링크 사이에 비교적 강한 비틀림 스프링을 두어 모델링할 수 있다(Lynch·Park 8.9.5절; 3.9). 이 장은 마찰의 위치 의존(기어 치 리플)·온도·하중 의존의 함수 형태를 정하지 않았다 (확인 필요).

3.8.3.5 정상 상태 마찰의 식별

속도를 일정하게 유지하며 토크를 측정하면 $\tau(\omega)=F_{ss}(\omega)+g(q)+o$(중력과 오프셋)이다. 같은 자세에서 $+\omega$와 $-\omega$로 측정해 $\bigl(\tau(+\omega)-\tau(-\omega)\bigr)/2=F_{ss}(\omega)$를 구하면 중력·오프셋이 소거된다(자세가 같을 때). $F_{ss}$는 $\omega_s$에 대해 비선형이고 나머지 세 계수는 선형이므로 변수 투영(VP)을 사용한다. $\omega_s$를 고정하면 $F_{ss}=a_1\operatorname{sgn}\omega+a_2e^{-(\omega/\omega_s)^2}\operatorname{sgn}\omega+\sigma_2\omega$가 $(a_1,a_2,\sigma_2)=(F_c,F_s-F_c,\sigma_2)$에 대해 선형이어서 정규방정식 하나로 풀고, 잔차 제곱합을 $\log\omega_s$의 1변수 함수로 보아 격자와 황금분할로 최소화한다. 점근 공분산은 $\sigma^2(J^\top J)^{-1}$, $J$는 $F_{ss}$의 파라미터 미분이다($\partial F/\partial F_c=\operatorname{sgn}\omega\,(1-e^{-r^2})$, $\partial F/\partial F_s=\operatorname{sgn}\omega\,e^{-r^2}$, $\partial F/\partial\omega_s=(F_s-F_c)\operatorname{sgn}\omega\,e^{-r^2}\,2\omega^2/\omega_s^3$, $\partial F/\partial\sigma_2=\omega$). 소인 속도가 $\omega_s$의 수 배보다 큰 영역만 있으면 지수 항이 0이어서 $J$의 두 열이 0이 되고 $F_s$·$\omega_s$를 식별할 수 없다. 동적 파라미터 $\sigma_0$·$\sigma_1$은 정지 영역의 소신호 응답(고유 진동수와 감쇠비)이나 시간 영역 모의 적합으로 구하며 이 장은 식별을 구현하지 않았다.

3.8.3.6 마찰 보상과 헌팅

마찰이 있는 1관성 서보 $J\ddot\theta=u-F$에서 PID $u=-K_pe-K_v\omega-K_i\!\int e$는 정지 마찰 이상의 적분이 쌓일 때까지 정지했다가 미끄러져 목표를 지나치는 한계 순환(헌팅)을 일으킬 수 있고 LuGre 논문이 이를 모의로 보였다. 논문은 마찰 상태 $z$를 관측하는 보상을 제안한다.

$$\dot{\hat z}=\omega-\frac{\sigma_0\lvert\omega\rvert}{h(\omega)}\hat z-k\,e,\quad \hat F=\sigma_0\hat z+\sigma_1\dot{\hat z}+\sigma_2\omega,\quad u=-H(s)e+\hat F+J\ddot x_d,\ \ k>0.$$

정리 1(논문): $G(s)=(\sigma_1s+\sigma_0)/(Js^2+H(s))$가 엄밀 양실수(SPR)이면 관측기 오차 $F-\hat F$와 위치 오차 $e$가 0으로 수렴한다(모형과 파라미터가 정확히 알려진 경우; 민감도는 논문이 다루지 않았다고 서술한다). $H(s)=K_vs+K_p$이면 $\operatorname{Re}G(j\omega)$의 분자가 $\sigma_0(K_p-J\omega^2)+\sigma_1K_v\omega^2$이므로 SPR 조건은 $\sigma_1K_v>\sigma_0J$이다(이 장이 논문의 정리에 적용해 유도했다). 논문은 모형 파라미터를 정확히 알 때 이 관측기 보상으로 위치·속도 추종 오차가 0으로 수렴함을 보였다.

3.8.3.7 마찰 모형 오차와 동역학 식별

3.7의 회귀 행렬에 마찰 열 $[\operatorname{sgn}\dot q,\ \dot q,\ 1]$을 더하면 쿨롱·점성·오프셋이 선형 파라미터로 식별된다. 참 마찰이 스트리벡 항을 가지면 모형 오차 $\delta(\dot q)=F_{ss}-F_{c}'\operatorname{sgn}\dot q-F_v'\dot q-o'$가 남는다. $\delta$는 $\dot q$에 대해 홀함수이고 속도가 0을 지나는 방향 전환 근처에서 크며 그 지점에서 $\ddot q$도 크므로 $\delta$가 $\ddot q$ 열과 상관되어 관성 파라미터에 편의를 만든다. 편의의 크기와 부호는 궤적에 따라 다르다.

3.8.4 알고리즘

부호 대칭 평균:  같은 자세·속도 ±ω에서 τ를 측정 → F_ss(ω) = (τ(+ω) − τ(−ω)) / 2
스트리벡 VP 적합(속도 v_i, 마찰 F_i):
    for ω_s in 로그 격자(60):  A ← [sgn v, e^{−(v/ω_s)²} sgn v, v];  c ← 정규방정식;  잔차 제곱합 기록
    최솟값 둘레 구간에서 황금분할로 log ω_s 정밀화
    F_c ← c₁,  F_s ← c₁ + c₂,  σ₂ ← c₃
LuGre 단계(RK4):  ż = ω − σ0|ω| z / h(ω);  F = σ0 z + σ1 ż + σ2 ω
관측기 기반 보상 서보(오차 e, 속도 ω):
    ẑ' = ω − σ0|ω| ẑ / h(ω) − k e;  F̂ = σ0 ẑ + σ1 ẑ' + σ2 ω;  u = −Kp e − Kv ω + F̂
    SPR 조건(PD): σ1 Kv > σ0 J

3.8.5 풀이 예제 — 정상 상태 곡선, 식별, LuGre, 헌팅과 보상, 마찰 모형 오차

결론부터 적으면 앞 절의 ①~⑤이다. 입력은 모두 예시 가정이다.

구현. f_ss·dfss(정상 상태 곡선과 기울기), lugre·rk4·servo(LuGre와 PID·관측기 서보), ss_basis·lin_fit·stribeck_fit(VP 적합)이 아래 코드이다. 소인 설계, 점근 공분산, 프리슬라이딩·이탈·마찰 지연 모의, 헌팅 모의, 1링크 편의 계산, 그림은 스크립트에 있다. 아래 코드는 1.1의 solve_lu를 사용한다.

from math import exp, log, sqrt, sin, cos, pi

def sgn(x):
    return (x > 0) - (x < 0)

def hfun(w, Fc, Fs, ws):
    """h(w) = Fc + (Fs - Fc) exp(-(w/ws)^2) = s0 g(w)"""
    return Fc + (Fs - Fc) * exp(-(w / ws) ** 2)

def f_ss(w, Fc, Fs, ws, s2):
    """정상 상태 마찰: 쿨롱 + 스트리벡 + 점성"""
    return hfun(w, Fc, Fs, ws) * sgn(w) + s2 * w

def dfss(w, Fc, Fs, ws, s2):
    """정상 상태 곡선의 기울기 (w != 0)"""
    e = exp(-(w / ws) ** 2)
    return -2 * (Fs - Fc) * w / ws ** 2 * e + s2

def lugre(w, z, P):
    """LuGre: z' = w - s0|w|z/h(w), F = s0 z + s1 z' + s2 w"""
    s0, s1, s2, Fc, Fs, ws = P
    zd = w - s0 * abs(w) * z / hfun(w, Fc, Fs, ws)
    return zd, s0 * z + s1 * zd + s2 * w

def rk4(f, st, dt):
    k1 = f(st)
    k2 = f([a + dt / 2 * b for a, b in zip(st, k1)])
    k3 = f([a + dt / 2 * b for a, b in zip(st, k2)])
    k4 = f([a + dt * b for a, b in zip(st, k3)])
    return [a + dt / 6 * (b + 2 * c + 2 * d + e)
            for a, b, c, d, e in zip(st, k1, k2, k3, k4)]

def servo(st, P, J, g, mode):
    """상태 (e, w, z, ie, zh), 목표 0, g = (Kp, Kv, Ki, k).
    'pid': u = -Kp e - Kv w - Ki ie.  'obs': u = -Kp e - Kv w
    + F^, F^는 zh' = w - s0|w|zh/h(w) - k e 의 관측기 마찰"""
    e, w, z, ie, zh = st
    Kp, Kv, Ki, k = g
    zd, F = lugre(w, z, P)
    if mode == 'pid':
        u, zhd = -Kp * e - Kv * w - Ki * ie, 0.0
    else:
        s0, s1, s2, Fc, Fs, ws = P
        h = hfun(w, Fc, Fs, ws)
        zhd = w - s0 * abs(w) * zh / h - k * e
        u = -Kp * e - Kv * w + s0 * zh + s1 * zhd + s2 * w
    return [w, (u - F) / J, zd, e, zhd]

def ss_basis(v, ws):
    """ws 고정 시 선형 열: sgn, exp(-(v/ws)^2) sgn, v"""
    return [[sgn(x), exp(-(x / ws) ** 2) * sgn(x), x]
            for x in v]

def lin_fit(A, y):
    """정규방정식 (A^T A) c = A^T y. (해, 잔차 제곱합)"""
    n = len(A[0])
    G = [[sum(r[i] * r[j] for r in A) for j in range(n)]
         for i in range(n)]
    b = [sum(r[i] * t for r, t in zip(A, y))
         for i in range(n)]
    c = solve_lu(G, b)[0]
    return c, sum((sum(a * x for a, x in zip(r, c)) - t) ** 2
                  for r, t in zip(A, y))

def stribeck_fit(v, F):
    """변수 투영: ws는 로그 격자 + 황금분할, (Fc, Fs-Fc, s2)는
    선형 최소제곱. (Fc, Fs, ws, s2, 잔차 제곱합)"""
    cost = lambda lw: lin_fit(ss_basis(v, exp(lw)), F)[1]
    lo, hi = log(1e-3), log(1.0)
    grid = [lo + i * (hi - lo) / 60 for i in range(61)]
    k = min(range(61), key=lambda i: cost(grid[i]))
    a, b = grid[max(k - 1, 0)], grid[min(k + 1, 60)]
    r = (sqrt(5) - 1) / 2
    for _ in range(60):
        c1, c2 = b - r * (b - a), a + r * (b - a)
        if cost(c1) < cost(c2):
            b = c2
        else:
            a = c1
    ws = exp((a + b) / 2)
    c, rss = lin_fit(ss_basis(v, ws), F)
    return c[0], c[0] + c[1], ws, c[2], rss

정상 상태 곡선의 최솟값은 $\omega=0.02492$ rad/s에서 $1.0135$ N·m이고 이탈 마찰 $F_s=1.5$ N·m보다 $0.49$ N·m 낮다. 쿨롱·점성만으로 이 소인 자료를 적합하면 잔차 RMS가 $0.159$ N·m, 최대 $0.339$ N·m로 마찰 수준 $\sim1$ N·m의 $16\,\%$이다.

넓은 소인의 점근 표준편차는 $F_c$ $0.011$, $F_s$ $0.023$ N·m, $\omega_s$ $7\times10^{-4}$ rad/s(참값의 $7\,\%$), $\sigma_2$ $0.015$ N·m·s/rad이고 몬테카를로 표준편차가 $10\,\%$ 안팎으로 같다($F_c$의 비 $0.87$은 200회 표본 표준편차의 상대 오차 $5\,\%$보다 크며 $\omega_s$ 추정의 비선형성이 한 원인일 수 있다). 고속만 소인하면 두 열이 정확히 0이어서 $F_s$·$\omega_s$를 식별할 수 없다. 저속 위주 소인은 $\sigma_2$의 표준편차가 $1.23$으로 넓은 소인의 80배이다. 소인은 $\omega_s$의 수 분의 1에서 수십 배까지 포함해야 한다.

프리슬라이딩에서 강모 변형 $z$는 선형 예측 $F/\sigma_0$와 $0.1\,\%$ 안에서 같다. 변위 $\theta$는 $z$의 3.1배인데 $z$가 $F_s/\sigma_0=1.5\times10^{-4}$ rad에 가까워져 포화하고 이후에는 미소 속도로 크리프하기 때문이다. 활주 직전 마찰은 토크 증가율이 클수록 작아진다($1.489\to1.173$ N·m, 증가율 $0.1\to100$ N·m/s). 이 경향의 방향이 논문의 모의나 실험과 같은지는 논문 그림의 수치를 읽지 못해 확인하지 못했다. 마찰 지연 루프의 면적은 $\omega_v=1\to100$ rad/s에서 $9.3\times10^{-6}\to1.1\times10^{-3}$으로 커지고 정상 상태 곡선과의 RMS 차도 $2.3\times10^{-4}\to2.0\times10^{-2}$ N·m로 커진다(논문이 서술한 "루프 폭이 주파수와 함께 증가"와 같은 방향). 강모의 지연 시간 상수가 $h/(\sigma_0\lvert\omega\rvert)\approx3\text{–}5$ ms이므로 예시 강성에서는 $\omega_v\sim100$ rad/s에서 루프가 뚜렷하다. 루프의 방향은 $\sigma_1$에 의존한다. $\sigma_1=0$이면 속도 $0.02$ rad/s(음의 기울기 구간)를 지날 때 증가 구간의 마찰($1.0764$ N·m)이 감소 구간($1.0120$)보다 높아 순수 지연의 기대와 같고, LuGre 논문이 Hess·Soom의 실험으로 서술한 "감소 속도의 마찰이 더 낮다"와 같은 방향이다. 감쇠비 $0.5$를 주는 $\sigma_1=99.5$에서는 $\sigma_1\dot z$ 항이 지연 효과를 넘어 방향이 반대($0.9635$ 대 $1.0261$ N·m)이다. 논문 그림의 모의 조건과 $\sigma_1$ 값을 읽어내지 못해 논문의 루프 방향과 직접 비교하지 못했다 (확인 필요).

PID 서보는 목표 $1$ rad로의 계단에서 $\pm0.0685$ rad 진폭, 주기 $23.9$ s의 한계 순환(헌팅)에 들어간다($K_v=6$; 정상 구간 오차가 거의 $\pm0.067$ rad에 머물며 부호만 바뀐다). $K_v=150$에서는 진폭이 $0.68$ rad로 커진다(100 s 동안 주기 몇 개만 관측). 관측기 보상은 SPR 조건($\sigma_1K_v>\sigma_0J$, $K_v>100.5$)을 만족하는 $K_v=150$에서 오차가 $10^{-19}$ rad로 수렴하고, 조건을 만족하지 못하는 $K_v=6$에서도 헌팅은 사라지지만 진폭 $2.05\times10^{-4}$ rad, 주기 $0.109$ s의 진동이 남는다(논문 정리 1의 보장 밖). 관측기는 모형과 파라미터를 정확히 안다고 가정한다.

쿨롱·점성·오프셋으로 적합한 관성은 최대 속도 $0.2$ rad/s의 궤적에서 $0.959$ kg·m²($+91.8\,\%$), $0.05$ rad/s에서 $+19.5\,\%$이고 $\ge1$ rad/s에서는 $+0.8\,\%$ 이하이다. 중력 파라미터 $k_g$는 $+1.3\,\%$ 이하로 덜 어긋난다. 스트리벡 항을 포함하고 $\omega_s$를 아는 적합은 $10^{-13}$ 상대 오차로 참값을 복원하므로 편의는 최소제곱의 결함이 아니라 모형 오차에서 나온다. 편의의 크기는 이 궤적(기본파와 2차 고조파)에서의 값이다.

라이브러리 대조. 같은 입력에서 SciPy·numpy·MuJoCo와 비교한 결과는 다음이다.

모든 항목이 일치한다. RK4와 DOP853의 차이($10^{-5}$~$10^{-6}$)는 $\lvert\omega\rvert$의 비평활성과 $\Delta t=10^{-4}$ s에서 오는 것이고 헌팅 진폭·주기의 결론을 바꾸지 않는다. MuJoCo의 frictionloss는 구속 풀이기에서 연성으로 처리되어 고착 토크(0.5·0.9 N·m)에서도 속도가 $10^{-3}$ rad/s 수준으로 남는다. 이 장에서 쿨롱 마찰의 이론값과 MuJoCo의 차이는 이 연성 때문이다.

한계. 합성 자료이고 마찰 파라미터는 모두 가정이다. 실제 관절은 위치(치 리플)·온도·하중·정지 시간에 따라 마찰이 변하며(Lynch·Park) 이 장은 이를 모델링하지 않았다. 헌팅 모의는 한 관성(강체)이어서 관절 유연성(3.9)이 있으면 진동수와 안정성이 달라진다. 식별은 정상 상태 곡선만 구현했고 $\sigma_0$·$\sigma_1$의 식별과 비국소 기억 히스테리시스(Leuven 모형)는 다루지 않았다. 관측기 보상은 파라미터를 정확히 아는 가정에서만 시험했다. 이탈 마찰의 방향과 크기는 이 장의 $\sigma_0$·$\sigma_1$·$J$ 값에 의존한다.

견해: 관절 토크 제어에서 마찰 모형은 정상 상태 스트리벡 곡선(변수 투영 식별)을 기본으로 두고, 저속 정밀 추종·정지 유지·헌팅이 문제인 과제에서만 LuGre 같은 동적 모형과 관측기를 더하는 것이 낫다고 본다. 근거는 이 장에서 ① 쿨롱·점성 모형이 스트리벡 곡선과 $0.159$ N·m RMS(마찰 수준의 $16\,\%$)만큼 어긋났고 저속 궤적의 관성 식별을 $+91.8\,\%$까지 편의시켰으며, ② 스트리벡 4파라미터가 넓은 소인 50점에서 상대 표준편차 $1\text{–}7\,\%$($F_c$ $1.1\,\%$, $F_s$ $1.5\,\%$, $\omega_s$ $7\,\%$, $\sigma_2$ $3\,\%$)로 식별되고 SciPy와 일치했다는 점, ③ 동적 보상은 모형이 정확하고 SPR 조건을 만족할 때에만 오차를 $10^{-19}$ rad로 줄였고 파라미터 오차를 시험하지 않았다는 점이다. 실물에서는 마찰의 하중·온도 의존으로 모형이 달라지므로 이 장의 값을 설계에 사용할 수 없다.

3.8.6 설계 트레이드오프 — 마찰 모형의 선택

마찰 모형은 다음 중에서 고른다.

• ① 쿨롱 + 점성(+ 오프셋). 파라미터 2~3개, 파라미터에 대해 선형이어서 동역학 회귀 행렬(3.7)의 열로 들어간다. 저속과 방향 전환에서 오차가 크다.
• ② 정상 상태 스트리벡. 파라미터 4개, $\omega_s$만 비선형(VP). 저속 속도 사상을 설명하지만 영속도에서 값이 하나가 아니다.
• ③ LuGre 동적 모형. 파라미터 6개와 상태 1개. 프리슬라이딩·마찰 지연·이탈 마찰의 변화를 모사하고 관측기 설계에 사용된다. 상태를 측정할 수 없어 관측기가 필요하고 구현의 수치 강성이 있다.
• ④ Leuven 통합 구조. 비국소 기억 히스테리시스를 가진 프리슬라이딩과 활주를 전환 함수 없이 잇는 모형(Swevers 외 2000, 초록·서론·2절 앞부분만 읽었다). 파라미터가 많고 이 장은 구현하지 않았다.

견해: 시뮬레이터(16.1)와 계산 토크 피드포워드에는 ①을 넣고, 저속·정지 정밀도가 요구되는 관절만 ②로 올리며, ③은 헌팅·스틱-슬립이 실측으로 확인된 경우에만 사용하는 것이 낫다고 본다. 근거는 이 장에서 ①이 고속에서는 편의가 $0.8\,\%$ 이하였지만 저속에서 $+91.8\,\%$까지 커졌다는 점과, ③의 보상이 정확한 모형을 전제로 했다는 점이다. 이 장은 실물에서 모형 간 성능을 비교하지 않았다.

3.8.7 응용 축 비교

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

3.8.8 관련 규격·프로토콜

마찰 모형이나 마찰 식별 절차를 규정한 국제 표준은 확인하지 못했다 (확인 필요).

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

관절 마찰 모형 자체를 규정한 법령·인증은 확인하지 못했다. 토크 센서 없이 모터 전류와 모델로 접촉을 검출하는 제품에서는 마찰 모델 오차가 검출 임계값의 오경보·미검출에 영향을 줄 수 있으나(9.27) 어느 표준이 이를 요구하는지는 확인하지 못했다 (확인 필요). 충돌 안전은 26.9, 산업 로봇 안전 표준은 27.10에서 다룬다.

3.8.10 한계·유의점·미해결

• 파라미터가 모두 가정이다. 실제 마찰은 위치·온도·하중·정지 시간에 따라 변한다.
• 한 관성. 헌팅·보상 모의는 강체 1관성이다. 유연 관절(3.9)·다관절 결합은 다루지 않았다.
• 관측기의 민감도. 모형 파라미터 오차에 대한 관측기 보상의 거동은 논문도 다루지 않았고 이 장도 시험하지 않았다.
• 동적 파라미터. $\sigma_0$·$\sigma_1$의 식별을 구현하지 않았다.
• 이탈 마찰의 경향. 증가율에 따른 감소가 문헌 그림과 같은 방향인지 확인하지 못했다.
• Leuven 모형. 비국소 기억 히스테리시스 모형은 초록·서론·2절 앞부분만 읽었고 구현하지 않았다.
• 연성 구속 마찰. 시뮬레이터에서 쿨롱 마찰의 연성 처리(MuJoCo)는 고착 토크 근처에서 미소 속도를 남긴다.

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

체크포인트:

1. 스트리벡 정상 상태 곡선의 기울기 $f_{ss}^\prime=0$이 $r\,e^{-r^2}=\kappa$로 정리됨을 보이고 음의 기울기 구간이 존재할 조건 $\kappa<1/\sqrt{2e}$를 유도하라.
2. LuGre 모형에서 $\omega$가 일정할 때 $z_{ss}=h(\omega)\operatorname{sgn}\omega/\sigma_0$임을 보이고 $F=F_{ss}(\omega)$를 확인하라.
3. 정지 영역의 선형화 $F\approx\sigma_0\theta+(\sigma_1+\sigma_2)\dot\theta$에서 고유 진동수와 감쇠비를 구하라.
4. 변수 투영에서 $\omega_s$를 고정할 때 나머지 세 계수가 선형인 이유와 소인이 $\omega_s$보다 훨씬 클 때 자코비안 두 열이 0이 되는 이유를 설명하라.
5. $H(s)=K_vs+K_p$일 때 $G(s)=(\sigma_1s+\sigma_0)/(Js^2+H(s))$의 SPR 조건이 $\sigma_1K_v>\sigma_0J$임을 유도하라.
6. 마찰 모형 오차가 관성 파라미터에 편의를 만드는 조건을 $\ddot q$ 열과 모형 오차의 상관으로 설명하라.

연구 질문:

• 파라미터 오차와 하중·온도 의존이 있을 때 관측기 기반 마찰 보상의 안정 영역은 어떻게 바뀌는가(이 장은 시험하지 않았다)?
• 고감속비 파동기어 관절에서 마찰의 위치 의존(치 리플)과 유연성(3.9)을 함께 고려한 모형을 실측으로 어떻게 식별하는가?
• 헌팅이 일어나는 적분 이득과 마찰 파라미터의 경계를 닫힌 식으로 구할 수 있는가?

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

• K. M. Lynch, F. C. Park, Modern Robotics: Mechanics, Planning, and Control, Cambridge University Press, 2017, 사전 인쇄본(2017-05) 8.9.4·8.9.5절 — 직접 읽음.
• C. Canudas de Wit, H. Olsson, K. J. Åström, P. Lischinsky, "A new model for control of systems with friction," IEEE Transactions on Automatic Control 40(3), 419–425, 1995, doi:10.1109/9.376053 — 직접 읽음(HAL 공개본 hal-05367408).
• J. Swevers, F. Al-Bender, C. G. Ganseman, T. Prajogo, "An integrated friction model structure with improved presliding behavior for accurate friction compensation," IEEE Transactions on Automatic Control 45(4), 675–686, 2000, doi:10.1109/9.847103 — 초록·서론·2절 앞부분만 읽음.
• B. Armstrong-Hélouvry, P. Dupont, C. Canudas de Wit, "A survey of models, analysis tools and compensation methods for the control of machines with friction," Automatica 30(7), 1083–1138, 1994, doi:10.1016/0005-1098(94)90209-7 — 서지만 확인.
• H. Olsson, K. J. Åström, C. Canudas de Wit, M. Gäfvert, P. Lischinsky, "Friction models and friction compensation," European Journal of Control 4(3), 176–195, 1998, doi:10.1016/s0947-3580(98)70113-x — 서지만 확인.
• P. Dupont, V. Hayward, B. Armstrong, F. Altpeter, "Single state elastoplastic friction models," IEEE Transactions on Automatic Control 47(5), 787–792, 2002, doi:10.1109/tac.2002.1000274 — 서지만 확인.
• D. Karnopp, "Computer simulation of stick-slip friction in mechanical dynamic systems," Journal of Dynamic Systems, Measurement, and Control 107(1), 100–103, 1985, doi:10.1115/1.3140698 — 서지만 확인.
• D. A. Haessig, B. Friedland, "On the modeling and simulation of friction," Journal of Dynamic Systems, Measurement, and Control 113(3), 354–362, 1991, doi:10.1115/1.2896418 — 서지만 확인.
• B. Bona, M. Indri, "Friction compensation in robotics: an overview," Proc. 44th IEEE CDC, 4360–4367, 2005, doi:10.1109/cdc.2005.1582848 — 서지만 확인.
• SciPy 1.18.1(optimize.least_squares·brentq, integrate.solve_ivp, special.lambertw), numpy 2.5.4, MuJoCo 3.17.0.

다음 장. 3.9은 유연 관절과 유연 링크의 동역학(스퐁 유연 관절 모델·가정 모드법)을 다룬다. 이 장의 마찰 모형이 강체 관절을 가정한 것과 달리 감속기의 탄성이 상태를 늘리는 경우를 다룬다.