제3부 동역학

3.9 유연 관절·유연 링크 동역학: 스퐁 유연 관절 모델·가정 모드법

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

이 장은 감속기·케이블·링크의 탄성이 로봇 동역학에 상태를 더하는 경우를 다룬다. 관절의 탄성을 회전 스프링으로 모델링한 스퐁 유연 관절 모델(링크 각 $q$와 모터 각 $\theta$를 함께 두는 $2n$ 자유도), 감속비가 모터 회전자의 관성을 $G^2$배로 반영하는 반사 관성, 1R 유연 관절의 공진·반공진 닫힌 식과 주파수 응답, 강성이 커질수록 강체 모델로 환원되는 특이 섭동의 정도와 그 오차의 $1/K$ 비례, 유연 관절의 역동역학이 링크 각의 4계 미분을 요구한다는 점과 이를 사용한 피드포워드, 연속체 링크를 유한 개의 모드로 줄이는 가정 모드법(Rayleigh–Ritz)의 수렴, 회전하는 유연 링크에서 허브(관절 구동기)의 관성이 링크의 진동수를 바꾸는 경계 조건 효과를 계산으로 확인한다. 강체 운동방정식은 3.4, 라그랑주 유도는 3.3, 마찰은 3.8, 동역학 파라미터 식별은 3.7, 진동과 모드 해석은 4.8, 직렬 탄성 구동기는 5.18, 유연 관절의 진동 억제 제어는 12.11, 계산 토크는 12.5, 연속체·소프트 로봇은 3.10와 이어진다.

배우는 것:

• 스퐁 모델 — 운동에너지 $T=\tfrac12\dot q^\top M(q)\dot q+\tfrac12\dot\theta^\top J\dot\theta$, 위치에너지에 스프링 $\tfrac12(q-\theta)^\top K(q-\theta)$, 방정식 두 개.
• 반사 관성 — 감속비 $G$에서 $J=G^2J_{\text{rotor}}$, 링크 관성 $M$과의 비가 공진비 $\omega_r/\omega_a=\sqrt{1+M/J}$를 정한다.
• 주파수 응답 — 모터 각 전달함수의 반공진, 링크 각 전달함수의 공진, 병치(collocated)와 비병치.
• 특이 섭동 — $K\to\infty$의 강체 환원, 강체 토크의 추종 오차가 $1/K$로 줄어드는 것, 4계 미분 피드포워드.
• 가정 모드법 — 연속체 링크의 모드 근사, 외팔보 진동수 방정식, 다항식 허용 함수의 수렴.
• 회전 링크의 경계 조건 — 허브 관성에 따라 첫 진동수가 핀-자유에서 고정-자유로 이동하는 것.

선수 지식: 3.1(강체 동역학), 3.3(라그랑주 정식화), 3.4(운동방정식), 1.2(일반화 고유값 문제의 선형대수).

특히 중요한 독자: 파동기어·케이블 구동 관절의 진동과 정착 시간을 다루는 제어 엔지니어, 로봇 팔의 모드 해석과 경량화 설계를 하는 기구 엔지니어, 시뮬레이터에서 관절 탄성을 모델링하는 개발자.

출처 방침. 기준 시점은 2026-11-05이다. 직접 읽은 자료: Modern Robotics(Lynch·Park, 사전 인쇄본 2017-05) 8.9.2절(겉보기 관성 $G^2I_{\text{rotor}}$)과 8.9.5절(파동기어의 플렉스플라인이 모터 회전자와 링크 사이의 비교적 강한 비틀림 스프링으로 모델링될 수 있다는 서술, 링크의 유한한 강성이 링크를 따른 진동으로 나타난다는 서술, 유연 관절·링크가 동역학의 상태를 늘려 동역학과 제어를 크게 복잡하게 한다는 서술). 스퐁 모델의 식, 1R의 전달함수와 공진·반공진, 특이 섭동 환원과 4계 미분 토크, 외팔보 진동수 방정식, 가정 모드법의 행렬, 회전 링크의 질량 행렬은 이 장이 에너지와 보 이론에서 직접 유도하고 numpy·SciPy·MuJoCo와 같은 입력에서 비교해 확인했다. 서지만 확인한 것(본문은 읽지 못했다): Spong(1987, 탄성 관절 로봇의 모델링과 제어, J. Dynamic Systems, Measurement, and Control), Book(1984, 유연 매니퓰레이터 팔의 재귀 라그랑주 동역학, IJRR), De Luca·Siciliano(1991, 평면 다링크 경량 로봇의 닫힌 형식 동역학), Cannon·Schmitz(1984, 유연 1링크 로봇의 끝점 제어 초기 실험), Tomei(1991, 탄성 관절 로봇의 단순 PD 제어기), Ozgoli·Taghirad(2006, 유연 관절 로봇 제어 개관), Albu-Schäffer 외(2007, DLR 경량 로봇의 설계와 제어), Pratt·Williamson(1995, 직렬 탄성 구동기), Ott 외(2008, 유연 관절 로봇의 수동성 기반 임피던스 제어), De Luca·Book(Springer Handbook of Robotics, 유연 요소를 가진 로봇 장). 이들의 내용은 제목과 서지에서 확인한 범위로만 서술한다. 스퐁 모델의 가정(회전자의 운동에너지가 자신의 회전만으로 이루어진다는 것)은 널리 알려진 모델의 가정으로 서술하며 Spong의 원문에서 확인하지 못했다 (확인 필요). 질량·강성·감속비는 모두 이 장이 정한 예시이며 어떤 제품의 제원이 아니다. 수치는 모두 이 장의 순수 파이썬 합성 계산이다.

3.9.2 개념 — 관절과 링크의 탄성은 자유도를 늘린다

결론부터 적으면 다음과 같다. 감속기(특히 파동기어)와 케이블·벨트는 완전히 강하지 않아서 모터 쪽 각 $\theta$와 링크 쪽 각 $q$가 서로 다르고, 두 각의 차 $q-\theta$가 스프링의 변형이다. 이 차를 무시하지 않는 모델이 유연 관절 모델이고 관절 하나당 자유도가 2개가 되어 상태가 두 배로 늘어난다. 모터 쪽에서 본 전달함수는 링크의 고유 진동수(회전자를 고정한 상태의 진동수) $\omega_a$에서 0이 되는 반공진과 그보다 높은 공진 $\omega_r$을 가지며, 두 진동수의 비는 모터 반사 관성 $J$와 링크 관성 $M$의 비만으로 $\omega_r/\omega_a=\sqrt{1+M/J}$로 정해진다. 강성 $K$를 키우면 강체 모델(두 관성의 합 $J+M$)로 환원되고 환원에 따른 링크 각의 오차가 $1/K$에 비례한다. 링크 자체의 탄성은 무한 자유도이고 가정 모드법은 링크의 휨 변위를 몇 개의 모양 함수의 합으로 근사해 유한 자유도로 줄인다. 이때 링크가 매달린 관절의 경계 조건(구동기를 고정할지, 토크 0으로 자유롭게 둘지)이 진동수를 바꾼다. 감속비가 큰 관절에서 회전자의 겉보기 관성 $G^2J_{\text{rotor}}$는 회전자 자체보다 크게 나타나 링크 관성과 같은 크기가 될 수 있다(Lynch·Park 8.9.2절). 관절 변형의 토크-비틀림 곡선이 보이는 히스테리시스는 이 장의 선형 모델에 들어 있지 않다.

이 장의 계산 결과를 미리 적으면 다음과 같다. ① 예시 1R 유연 관절($M=0.5$, $J=0.2$ kg·m², $K=2000$ N·m/rad)에서 반공진은 $63.246$ rad/s, 공진은 $118.322$ rad/s, 비는 $\sqrt{1+M/J}=1.8708$이며 numpy 다항식 근과 $10^{-16}$ 안에서 같고, MuJoCo 두 힌지·힘줄 스프링의 자유 진동수가 닫힌 식과 $1.8\times10^{-6}$ 안에서 같다. ② 같은 팔을 강체 모델 토크 $(J+M)\ddot q+mgl\sin q$로 구동하면 1 Hz 기준 궤적(진폭 1 rad)의 링크 각 추종 오차가 $K=10^3\to10^6$ N·m/rad에서 $4.4\times10^{-3}\to4.4\times10^{-6}$ rad으로 $1005$배 줄고($1/K$ 비례), 링크 각의 4계 미분을 포함한 유연 관절 역동역학 토크는 모든 $K$에서 $10^{-11}$ rad 수준(수치 적분 오차)이다. ③ 외팔보의 다항식 가정 모드법은 $N=3$에서 첫 진동수 $\mu_1=1.875387$(정확한 값 $1.875104$), $N=8$에서 처음 세 진동수가 정확한 값과 $5\times10^{-6}$ 안에서 같다. ④ 회전 링크의 첫 진동수 $\mu_1$는 허브 관성 비 $r=0$에서 핀-자유 값 $3.9266$이고 $r\to\infty$에서 고정-자유 값 $1.8751$로 내려간다.

3.9.3 원리

3.9.3.1 스퐁 유연 관절 모델

각 관절의 모터(와 감속기 입력 쪽)의 각 $\theta\in\mathbb R^n$을 링크 쪽으로 환산한 값(감속비로 나눈 값)으로 두고, 링크 각 $q\in\mathbb R^n$과 대각 강성 행렬 $K$의 선형 스프링이 둘을 잇는다고 하자. 모터 회전자의 운동에너지가 자신의 축 둘레 회전만으로 이루어진다고 가정하면(스퐁 모델의 가정으로 알려져 있다 (확인 필요)) 에너지는 다음이다.

$$T=\tfrac12\,\dot q^\top M(q)\,\dot q+\tfrac12\,\dot\theta^\top J\,\dot\theta,\qquad V=V_g(q)+\tfrac12\,(q-\theta)^\top K\,(q-\theta).$$

여기서 $M(q)$는 링크 질량행렬(3.4), $J$는 반사 회전자 관성의 대각 행렬이다. 오일러-라그랑주 식(3.3)에서 토크 $\tau$가 $\theta$에 작용하는 입력일 때 다음이 나온다.

$$M(q)\ddot q+C(q,\dot q)\dot q+g(q)+K(q-\theta)=0,\qquad J\ddot\theta+K(\theta-q)=\tau.$$

첫 식은 링크 동역학(스프링 토크가 입력), 둘째 식은 모터 동역학(스프링 토크가 부하)이다. 관절 변형 $q-\theta$의 감쇠 $D(\dot q-\dot\theta)$는 두 식에 부호를 달리해 더할 수 있고(이 장의 계산은 무감쇠), 마찰(3.8)은 둘째 식에 더한다. 두 식은 강체 모델 $M\ddot q+C\dot q+g=\tau_{\text{joint}}$(3.4)에서 관절 토크 $\tau_{\text{joint}}=K(\theta-q)$를 스프링으로 주는 구조이며, $K\to\infty$이면 $\theta\to q$이고 두 식의 합이 $(M+J)\ddot q+C\dot q+g=\tau$가 된다.

3.9.3.2 반사 관성

감속비 $G$의 기어헤드에서 모터 회전자의 각속도는 링크 쪽 각속도의 $G$배이므로 회전자 운동에너지가 $\tfrac12I_{\text{rotor}}(G\dot\theta)^2=\tfrac12\,(G^2I_{\text{rotor}})\,\dot\theta^2$이다. $G^2I_{\text{rotor}}$를 겉보기 관성(반사 관성)이라 하며(Lynch·Park 8.9.2절) 스퐁 모델의 $J$가 이 값이다. 회전자 관성 $2\times10^{-5}$ kg·m², $G=100$이면 $J=0.2$ kg·m²로 회전자의 $10^4$배이고 링크 관성 $M=0.5$ kg·m²와 같은 크기가 된다. 감속비가 클수록 질량행렬이 대각에 가까워지는 효과(Lynch·Park)는 이 $J$의 증가에서 나온다. 이상 기어(효율 1)에서는 링크 쪽으로 환산한 입력 토크가 $\tau=G\,\tau_{\text{motor}}$이고 모터 쪽에서 측정한 비틀림 변위는 $G$로 나누어 $\theta$에 넣는다.

3.9.3.3 1R 유연 관절의 주파수 응답

중력을 무시하고 $M\ddot q+K(q-\theta)=0$, $J\ddot\theta+K(\theta-q)=\tau$를 라플라스 변환하면 다음이다.

$$\frac{\Theta(s)}{T(s)}=\frac{Ms^2+K}{s^2\bigl[JMs^2+K(J+M)\bigr]},\qquad \frac{Q(s)}{T(s)}=\frac{K}{s^2\bigl[JMs^2+K(J+M)\bigr]}.$$

극은 강체 모드의 이중 영점과 $\pm j\omega_r$, $\omega_r=\sqrt{K(J+M)/(JM)}$이다. 모터 각 전달함수는 $\pm j\omega_a$, $\omega_a=\sqrt{K/M}$에 영점(반공진)을 가지며 이는 모터를 고정했을 때의 링크 진동수이다. 링크 각 전달함수는 영점이 없고 상대 차수가 4이다. 입력 $\tau$와 출력 $\theta$가 같은 축에 있는 병치(collocated) 구성은 영점이 극보다 먼저 나타나 위상이 $-180^\circ$ 아래로 내려가지 않고, 링크 각을 출력으로 보는 비병치 구성은 공진 뒤에 위상이 $-360^\circ$로 내려가 제어가 어렵다. $\omega_r^2/\omega_a^2=1+M/J$이므로 $J$가 $M$보다 작을수록 공진이 반공진보다 높고, $J\gg M$이면 두 진동수가 같아져 구분되지 않는다.

3.9.3.4 강체 환원과 유연 관절의 역동역학

첫 식에서 $\theta=q+K^{-1}\bigl(M\ddot q+C\dot q+g\bigr)$이므로 $\theta-q=O(1/K)$이다(느린 다양체). 이를 둘째 식에 넣으면 토크가 링크 각만의 식이 된다.

$$\tau=J\ddot\theta+\bigl(M\ddot q+C\dot q+g\bigr),\qquad \ddot\theta=\ddot q+\frac{d^2}{dt^2}\Bigl[K^{-1}\bigl(M\ddot q+C\dot q+g\bigr)\Bigr].$$

$\ddot\theta$가 $\tfrac{d^2}{dt^2}(M\ddot q)$를 포함하므로 $q$의 4계 도함수가 필요하다. 따라서 기준 궤적이 4번 미분 가능해야 하고(스냅까지 연속) 강체 모델의 토크 $(J+M)\ddot q+C\dot q+g$는 $K^{-1}$ 항을 버린 것이다. 링크 각이 이 모델의 모든 상태와 입력을 결정하는 출력이라는 사실은 이 장의 식에서 직접 확인된다(평탄 출력과의 관계는 문헌에서 확인하지 못했다 (확인 필요)). Spong(1987)이 이 모델의 특이 섭동 환원을 다룬 것으로 알려져 있다(본문은 읽지 못했다 (확인 필요)). 강체 토크로 구동한 링크 각 오차 $e=q-q_d$는 $O(1/K)$로 줄어든다.

3.9.3.5 유연 링크와 가정 모드법

연속체 링크는 무한 자유도이다. 균일 보(Euler–Bernoulli, 질량 선밀도 $\rho A$, 굽힘 강성 $EI$, 길이 $L$)의 횡변위 $w(x,t)$는 $\rho A\,w_{tt}+EI\,w_{xxxx}=0$을 만족하고, 고정-자유 경계($w(0)=w'(0)=0$, $w''(L)=w'''(L)=0$)에서 변수 분리를 하면 진동수 방정식이 다음이다.

$$\cos\mu\cosh\mu+1=0,\qquad \omega_i=\mu_i^2\sqrt{\frac{EI}{\rho AL^4}},\quad \mu_i=\beta_iL=1.8751,\ 4.6941,\ 7.8548,\ \ldots$$

가정 모드법은 $w(x,t)=\sum_{k=1}^{N}\phi_k(x)\,q_k(t)$로 두고 에너지에서 방정식을 만든다. 운동에너지 $\tfrac12\int\rho A\,w_t^2dx$에서 질량행렬 $M_{kl}=\int\rho A\,\phi_k\phi_l\,dx$, 변형에너지 $\tfrac12\int EI\,(w'')^2dx$에서 강성행렬 $K_{kl}=\int EI\,\phi_k''\phi_l''\,dx$가 나온다. 허용 함수 $\phi_k$는 기하 경계 조건($\phi_k(0)=\phi_k'(0)=0$)만 만족하면 되고(Rayleigh–Ritz), 이 장은 $\phi_k=(x/L)^{k+1}$ ($k=1\ldots N$)을 사용했다. $\rho A=L=EI=1$로 정규화하면 $M_{kl}=1/(k+l+3)$, $K_{kl}=(k+1)k(l+1)l/(k+l-1)$이다. 일반화 고유값 문제 $Kv=\lambda Mv$의 $\mu_i=\lambda_i^{1/4}$는 정확한 $\mu_i$의 상한으로 수렴한다(허용 함수 공간이 커질수록 진동수가 줄어든다). 다링크 평면 경량 로봇의 닫힌 형식 동역학 모델(De Luca·Siciliano 1991)과 유연 매니퓰레이터 팔의 재귀 라그랑주 동역학(Book 1984)이 제목에서 확인되며, 두 논문이 사용한 방법(가정 모드법 여부 포함)은 확인하지 못했다.

3.9.3.6 회전하는 유연 링크와 경계 조건

관절에서 허브 각 $\theta$로 회전하는 보의 점 $x$의 속도 성분은 $x\dot\theta+\dot w$(변형이 작고 원심 강성화를 무시하는 선형화)이다. 허브 관성 $I_h$를 더한 운동에너지에서 좌표 $(\theta,q_1,\ldots,q_N)$의 질량행렬은 다음이다.

$$\mathcal M=\begin{bmatrix}I_h+\rho AL^3/3 & m_k^\top\\ m_k & M_{kl}\end{bmatrix},\qquad m_k=\rho A\!\int_0^L x\,\phi_k\,dx=\frac{\rho A L^2}{k+3},$$

강성행렬은 $\theta$ 성분이 0인 블록 대각이다. 이 시스템의 고유값은 강체 모드의 0 하나와 휨 진동수이다. 허브가 고정($I_h\to\infty$, $\dot\theta=0$)이면 외팔보의 $\mu_i$(고정-자유)이고, 허브가 자유($I_h=0$, 토크 0)이면 핀-자유 보의 근 $\tan\mu=\tanh\mu$, $\mu=3.9266,\ 7.0686,\ 10.2102,\ldots$이다. 허브의 관성이 중간이면 진동수가 그 사이에 놓인다. 경계 조건이 구동기의 임피던스로 정해지므로 같은 링크라도 토크 제어(관절이 자유)와 위치 제어(관절이 고정에 가까움)에서 첫 진동수가 두 배 이상 다르다. 이 장의 선형화는 회전 속도가 커질 때 나타나는 원심 강성화와 기하학적 비선형을 포함하지 않는다.

3.9.3.7 감쇠·히스테리시스·백래시

실제 전달계는 감쇠와 비선형 강성을 가진다. 선형 감쇠 $D$는 공진 봉우리의 높이를 정하고, 토크-비틀림 곡선의 히스테리시스는 진폭에 따라 달라지는 감쇠와 강성으로 나타난다. 파동기어의 토크-비틀림 곡선이 이런 비선형성과 히스테리시스를 보이는 것으로 알려져 있으나 이 장은 확인하지 못했다 (확인 필요). 감속기의 간극(백래시)은 4.3·4.6에서, 마찰의 동적 모형(Dahl·LuGre의 강모 모형이 히스테리시스 루프를 만든다)은 3.8에서 다루며 이 장은 이 항들을 계산하지 않았다.

3.9.4 알고리즘

유연 관절 시뮬레이션(1R, 입력 tau(t)):
    상태 s = (q, qd, th, thd);  RK4로 적분;  dt ≤ 0.05 / ω_r (ω_r = sqrt(K (J+M) / (J M)))
유연 관절 피드포워드(기준 q_d(t)가 C⁴):
    th_d = q_d + (M q_d'' + g(q_d)) / K
    tau = J th_d'' + M q_d'' + g(q_d),  th_d'' = q_d'' + (M q_d'''' + d²g/dt²) / K
    초기 상태를 (q_d, q_d', th_d, th_d')(0)로 맞춘다
가정 모드법(외팔보):
    M_kl ← 1/(k+l+3);  K_kl ← (k+1)k(l+1)l/(k+l−1)  (k, l = 1..N)
    M = L Lᵀ (촐레스키);  S ← L⁻¹ K L⁻ᵀ;  S의 고유값(야코비) λ_i;  μ_i ← λ_i^{1/4}
회전 링크: 좌표 (θ, q);  M[0][0] ← I_h + 1/3;  M[0][k] ← 1/(k+3);  K의 θ 행·열 ← 0

3.9.5 풀이 예제 — 유연 관절의 공진, 강체 환원, 가정 모드법, 회전 링크

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

구현. flex_rhs·rk4(스퐁 1R 모의), ref·tau_rigid·tau_flex·init_state(기준 궤적과 강체·유연 관절 토크), jacobi_eigs·gen_eigs·beam_mats(가정 모드법의 행렬과 일반화 고유값)가 아래 코드이다. 주파수 응답, 강성 소인, 회전 링크의 허브 관성 소인, 그림은 스크립트에 있다. 외팔보와 핀-자유 보의 근은 이분법으로 구했다.

from math import sin, cos, sqrt, pi

def flex_rhs(s, tau, P):
    """스퐁 1R: M q'' + mgl sin q + K (q - th) = 0,
    J th'' + K (th - q) = tau.  s = (q, qd, th, thd)"""
    M, J, K, mgl = P
    q, qd, th, thd = s
    return [qd, -(mgl * sin(q) + K * (q - th)) / M,
            thd, (tau - K * (th - q)) / J]

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

def ref(t, A, W):
    """q_d = A (1 - cos W t) / 2 와 1~4계 도함수"""
    c, s = cos(W * t), sin(W * t)
    return (A * (1 - c) / 2, A * W * s / 2,
            A * W ** 2 * c / 2, -A * W ** 3 * s / 2,
            -A * W ** 4 * c / 2)

def tau_rigid(t, A, W, P):
    """강체 모델 토크: (J + M) q'' + mgl sin q"""
    M, J, K, mgl = P
    q, qd, qdd, _, _ = ref(t, A, W)
    return (J + M) * qdd + mgl * sin(q)

def tau_flex(t, A, W, P):
    """유연 관절 역동역학. th = q + (M q'' + g)/K 이므로
    tau = J th'' + M q'' + g, th''에 q의 4계 미분이 든다"""
    M, J, K, mgl = P
    q, qd, qdd, q3, q4 = ref(t, A, W)
    g = mgl * sin(q)
    gdd = mgl * (cos(q) * qdd - sin(q) * qd ** 2)
    return J * (qdd + (M * q4 + gdd) / K) + M * qdd + g

def init_state(A, W, P):
    """기준 궤적과 일치하는 초기 상태 (q, qd, th, thd)"""
    M, J, K, mgl = P
    q, qd, qdd, q3, q4 = ref(0.0, A, W)
    th = q + (M * qdd + mgl * sin(q)) / K
    thd = qd + (M * q3 + mgl * cos(q) * qd) / K
    return [q, qd, th, thd]

def jacobi_eigs(A, sweeps=60):
    """대칭 행렬의 고유값(순환 야코비), 오름차순"""
    n, A = len(A), [r[:] for r in A]
    for _ in range(sweeps):
        if sum(A[i][j] ** 2 for i in range(n)
               for j in range(i)) < 1e-30:
            break
        for p in range(n - 1):
            for q in range(p + 1, n):
                if abs(A[p][q]) < 1e-300:
                    continue
                x = (A[q][q] - A[p][p]) / (2 * A[p][q])
                sg = 1 if x >= 0 else -1
                t = sg / (abs(x) + sqrt(x * x + 1))
                c = 1 / sqrt(t * t + 1)
                s = t * c
                for k in range(n):
                    u, v = A[k][p], A[k][q]
                    A[k][p] = c * u - s * v
                    A[k][q] = s * u + c * v
                for k in range(n):
                    u, v = A[p][k], A[q][k]
                    A[p][k] = c * u - s * v
                    A[q][k] = s * u + c * v
    return sorted(A[i][i] for i in range(n))

def gen_eigs(K, M):
    """K v = lam M v: M = L L^T 로 L^-1 K L^-T 를 대칭화"""
    n = len(M)
    L = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(i + 1):
            s = M[i][j] - sum(L[i][k] * L[j][k]
                              for k in range(j))
            L[i][j] = sqrt(s) if i == j else s / L[j][j]
    def fwd(b):
        x = [0.0] * n
        for i in range(n):
            x[i] = (b[i] - sum(L[i][k] * x[k]
                               for k in range(i))) / L[i][i]
        return x
    X = [fwd([K[i][j] for i in range(n)]) for j in range(n)]
    S = [fwd([X[i][j] for i in range(n)]) for j in range(n)]
    return jacobi_eigs([[(S[i][j] + S[j][i]) / 2
                         for j in range(n)]
                        for i in range(n)])

def beam_mats(N, Ih=None):
    """외팔보(rho A = L = EI = 1), 허용 함수 phi_k = x^(k+1),
    k = 1..N. Ih를 주면 허브 회전 좌표 th를 맨 앞에 더한
    회전 링크의 (M, K)이다(허브 관성 Ih)"""
    rg = range(1, N + 1)
    M = [[1 / (k + l + 3) for l in rg] for k in rg]
    K = [[(k + 1) * k * (l + 1) * l / (k + l - 1) for l in rg]
         for k in rg]
    if Ih is None:
        return M, K
    row = [1 / (k + 3) for k in rg]
    M = [[Ih + 1 / 3] + row] + [[row[k - 1]] + M[k - 1]
                                for k in rg]
    K = [[0.0] * (N + 1)] + [[0.0] + K[k - 1] for k in rg]
    return M, K

주파수 응답 격자 탐색(간격 $0.1\,\%$)의 극값이 닫힌 식과 $10^{-4}$ 안에서 같다. $J/M=0.4$이므로 공진이 반공진의 $1.871$배이다. 회전자의 반사 관성이 링크 관성과 같은 크기($0.2$ 대 $0.5$)여서 공진과 반공진이 가까운 값에 놓인다.

강체 모델 토크로 구동했을 때 링크 각 오차는 $K$에 반비례하고($K=10^3$과 $10^6$의 오차 비 $1005$, 같은 궤적에서 반공진 $\omega_a$는 $45$ rad/s에서 $1414$ rad/s로 31.6배) 4계 미분을 포함한 토크는 모든 $K$에서 $10^{-11}$ rad으로 수치 적분 오차 수준이다. 두 토크의 RMS 차는 $K=10^4$에서 $4.8\times10^{-3}$ N·m로 강체 토크 RMS $8.97$ N·m의 $0.05\,\%$이다. 이 작은 토크 차이가 링크 각 오차 $4.4\times10^{-4}$ rad(진폭의 $0.044\,\%$)에 해당한다.

다항식 가정 모드법은 상한으로 수렴한다($\mu_1$: $N=1$ 2.1147, $N=3$ 1.8754, $N=8$ 1.875104). 고차 진동수는 더 많은 허용 함수가 필요하다($\mu_3$는 $N=3$에서 $10.87$, $N=8$에서 $7.854762$로 정확한 값 $7.854757$과 같은 자리까지). 다항식 허용 함수의 질량행렬은 힐베르트 행렬이어서 $N$이 커지면 조건수가 커지므로(이 장은 $N=8$까지) SciPy의 일반화 고유값과 $2\times10^{-10}$ 안에서 같았다.

회전 링크의 첫 진동수 $\mu_1$은 허브 관성 비 $r=0$(토크 0, 핀-자유)에서 $3.9266$이고 $r=1000$에서 $1.8756$으로 고정-자유 값 $1.8751$에 가까워진다. 진동수 비율로는 $(3.9266/1.8751)^2=4.39$배이다. 높은 모드($\mu_3$)는 허브 관성에 덜 민감하고 $r=0$의 $\mu_3=10.2132$는 $N=8$ 근사이며 정확한 근은 $10.2102$이다.

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

모든 항목이 일치한다. MuJoCo는 모터와 링크를 두 힌지로, 스프링을 두 관절 좌표의 차를 변형량으로 하는 힘줄 강성으로 구성했고 중력 없이 자유 진동수만 비교했다. 따라서 이 대조는 선형 모델의 진동수를 확인한 것이며 중력과 비선형 항이 있는 응답은 SciPy Radau와의 비교(두 적분기 차이 $1.0\times10^{-11}$ rad)로 확인했다.

한계. 선형 스프링, 무감쇠, 1R, 균일 보, 작은 변형이다. 관절 감쇠·히스테리시스·백래시·비선형 강성(소프트 윈드업)은 계산하지 않았다. 다관절 유연 팔의 질량행렬과 코리올리 항, 링크 변형이 큰 비선형 동역학, 유한요소 모델(FEM)은 다루지 않았다. 4계 미분 피드포워드는 기준 궤적이 정확히 알려진 이상적인 경우이고 모형 파라미터 오차와 토크 한계를 고려하지 않았다. 다항식 허용 함수는 $N$이 크면 조건이 나빠지므로 실제 모드 해석에서는 직교 모드 함수나 FEM을 사용한다.

견해: 관절 강성이 기준 궤적의 지배 주파수에 비해 충분히 크면 강체 모델로 충분하고, 그렇지 않으면 링크 각의 4계 미분을 포함한 피드포워드나 병치 구성의 제어(모터 각 피드백, 입력 성형 등; 12.11)를 더하는 것이 낫다고 본다. 근거는 이 장의 한 예시에서 반공진이 기준 궤적 주파수($\Omega=6.28$ rad/s)의 $22.5$배($K=10^4$)일 때 강체 토크의 링크 각 오차가 진폭의 $0.044\,\%$, $7.1$배($K=10^3$)일 때 $0.44\,\%$로 $1/K$에 비례했다는 점이다. 허용 오차와 기준 궤적의 스펙트럼, 관성비 $J/M$, 감쇠에 따라 달라지고 이 장은 실물 관절을 시험하지 않았다. 유연 링크는 첫 진동수가 구동기의 경계 조건에 따라 두 배 이상 달라지므로 모드 해석(4.8)에서 허브 경계를 구동기의 제어 모드(토크·위치)와 함께 정해야 한다.

3.9.6 설계 트레이드오프 — 모델의 자유도 선택

관절·링크의 탄성을 모델에 넣는 방식은 다음 중에서 고른다.

• ① 강체 모델. 탄성을 무시한다. $n$ 자유도, 계산 토크와 식별(3.7)이 단순하다. 강성이 충분하지 않으면 진동과 위치 오차가 남는다.
• ② 유연 관절 모델(스퐁). $2n$ 자유도. 감속기·케이블·토크 센서의 탄성을 포함하고 모터 쪽 제어와 링크 쪽 출력을 구분한다. 역동역학에 4계 미분이 필요하고 상태가 늘어난다.
• ③ 유연 링크 포함(가정 모드법·FEM). 링크 휨의 모드를 추가한다. 긴 링크·경량 링크에 필요하고 모드 수와 경계 조건을 정해야 한다.

견해: 파동기어나 케이블을 사용하는 관절은 ②를, 링크 길이 대 단면 치수가 큰 장대 팔은 ③을 기본으로 두고, 강성이 큰 짧은 팔은 ①로 두는 것이 낫다고 본다. 근거는 이 장에서 ②의 모터 쪽과 링크 쪽 전달함수가 서로 다른 구조(반공진 대 공진)를 가진다는 점과, ①의 오차가 $1/K$로 줄어 강성이 충분하면 무시할 수 있다는 점이다. 이 장은 세 방식의 계산 비용과 정확도를 실물에서 비교하지 않았다.

3.9.7 응용 축 비교

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

3.9.8 관련 규격·프로토콜

유연 관절·링크 모델 자체를 규정한 국제 표준은 확인하지 못했다 (확인 필요).

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

관절·링크 탄성의 모델링을 규정한 법령·인증은 확인하지 못했다. 직렬 탄성 구동기와 토크 센서를 사용하는 관절은 협동 운전에서 힘 제한 기능(26.7)과 충돌 안전(26.9)의 구현 방식과 연관되나 어느 표준이 구성 요소의 탄성 특성을 요구하는지는 확인하지 못했다 (확인 필요).

3.9.10 한계·유의점·미해결

• 선형 스프링, 무감쇠. 감쇠·히스테리시스·백래시·소프트 윈드업(비선형 강성)은 계산하지 않았다.
• 1R과 외팔보. 다관절 유연 팔의 질량행렬·코리올리 항과 링크 변형이 큰 비선형 동역학은 다루지 않았다.
• 회전 링크의 선형화. 원심 강성화와 기하학적 비선형을 무시했다.
• 허용 함수. 다항식은 $N$이 크면 조건이 나쁘다. 직교 모드 함수·FEM은 다루지 않았다.
• 4계 미분 피드포워드. 이상적인 모델·궤적에서의 결과이며 파라미터 오차와 토크 한계를 고려하지 않았다.
• 스퐁 모델의 가정. 회전자 운동에너지가 자신의 회전만으로 이루어진다는 가정과 그 위반(자이로 결합)의 영향은 계산하지 않았다.

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

체크포인트:

1. 스퐁 모델의 에너지에서 오일러-라그랑주 식을 유도하고 $K\to\infty$에서 $(M+J)\ddot q+C\dot q+g=\tau$가 됨을 보여라.
2. 1R 모델의 전달함수 $\Theta/T$와 $Q/T$를 구하고 반공진 $\omega_a=\sqrt{K/M}$, 공진 $\omega_r=\sqrt{K(J+M)/(JM)}$, 비 $\sqrt{1+M/J}$를 확인하라.
3. 유연 관절의 역동역학 토크에 링크 각의 4계 도함수가 들어가는 이유를 $\theta=q+K^{-1}(M\ddot q+C\dot q+g)$에서 설명하라.
4. 외팔보의 진동수 방정식 $\cos\mu\cosh\mu=-1$을 경계 조건에서 유도하고 가정 모드법의 $N=1$ 근사 $\mu_1=2.1147$을 직접 계산하라.
5. 회전 링크에서 허브가 자유일 때와 고정일 때의 경계 조건과 첫 진동수 $3.9266$, $1.8751$의 차이를 설명하라.
6. 반사 관성 $G^2I_{\text{rotor}}$가 감속비에 따라 링크 관성과 같은 크기가 되는 조건을 구하라.

연구 질문:

• 파동기어의 비선형 강성과 히스테리시스가 있을 때 유연 관절 모델의 파라미터를 식별하고 제어하는 방법은 무엇인가(이 장은 시험하지 않았다)?
• 유연 링크의 모드 수를 기준 궤적의 주파수 성분과 허용 오차에서 어떻게 정하는가?
• 직렬 탄성 구동기(5.18)의 스프링 강성을 과제의 충격·에너지 요구에서 어떻게 설계하는가?

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

• K. M. Lynch, F. C. Park, Modern Robotics: Mechanics, Planning, and Control, Cambridge University Press, 2017, 사전 인쇄본(2017-05) 8.9.2·8.9.5절 — 직접 읽음.
• M. W. Spong, "Modeling and control of elastic joint robots," Journal of Dynamic Systems, Measurement, and Control 109(4), 310–318, 1987, doi:10.1115/1.3143860 — 서지만 확인.
• W. J. Book, "Recursive Lagrangian dynamics of flexible manipulator arms," International Journal of Robotics Research 3(3), 87–101, 1984, doi:10.1177/027836498400300305 — 서지만 확인.
• A. De Luca, B. Siciliano, "Closed-form dynamic model of planar multilink lightweight robots," IEEE Transactions on Systems, Man, and Cybernetics 21(4), 826–839, 1991, doi:10.1109/21.108300 — 서지만 확인.
• R. H. Cannon, E. Schmitz, "Initial experiments on the end-point control of a flexible one-link robot," International Journal of Robotics Research 3(3), 62–75, 1984, doi:10.1177/027836498400300303 — 서지만 확인.
• P. Tomei, "A simple PD controller for robots with elastic joints," IEEE Transactions on Automatic Control 36(10), 1208–1213, 1991, doi:10.1109/9.90238 — 서지만 확인.
• S. Ozgoli, H. D. Taghirad, "A survey on the control of flexible joint robots," Asian Journal of Control 8(4), 332–344, 2006, doi:10.1111/j.1934-6093.2006.tb00285.x — 서지만 확인.
• A. Albu-Schäffer, S. Haddadin, C. Ott, A. Stemmer, T. Wimböck, G. Hirzinger, "The DLR lightweight robot: design and control concepts for robots in human environments," Industrial Robot 34(5), 376–385, 2007, doi:10.1108/01439910710774386 — 서지만 확인.
• G. A. Pratt, M. M. Williamson, "Series elastic actuators," Proc. IEEE/RSJ IROS 1995, 399–406, doi:10.1109/iros.1995.525827 — 서지만 확인.
• C. Ott, A. Albu-Schäffer, A. Kugi, G. Hirzinger, "On the passivity-based impedance control of flexible joint robots," IEEE Transactions on Robotics 24(2), 416–429, 2008, doi:10.1109/tro.2008.915438 — 서지만 확인.
• A. De Luca, W. J. Book, "Robots with flexible elements," in Springer Handbook of Robotics, Springer, 2008, 287–319, doi:10.1007/978-3-540-30301-5_14 — 서지만 확인.
• numpy 2.5.4, SciPy 1.18.1(linalg.eigh, optimize.brentq, integrate.solve_ivp), MuJoCo 3.17.0.

다음 장. 3.10는 연속체·소프트 로봇의 동역학(코세라 막대·FEM 소프트 모델)을 다룬다. 이 장의 보 모델(작은 변형, 선형)에서 큰 변형의 비선형 보로 확장한다.