AEGLOSLAB 로봇공학 전서(全書)

제3부 동역학

3.3 라그랑주 정식화와 에너지 기반 유도

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

이 장은 로봇의 운동방정식을 운동 에너지와 위치 에너지에서 유도하는 라그랑주 정식화를 다룬다. 일반화 좌표·일반화 힘과 오일러-라그랑주 방정식, 직렬 개연쇄의 운동 에너지 $K=\frac12\dot q^{T}M(q)\dot q$와 질량행렬 $M=\sum_i J_{ib}^{T}G_iJ_{ib}$, 크리스토펠 기호로 나타나는 코리올리·원심 항, 에너지 균형(수동성)과 레일리 감쇠, 단일 강체의 오일러 방정식을 에너지에서 얻는 오일러-푸앙카레 변분, 그리고 이 정식화가 3.2의 RNEA와 같은 토크를 내는지와 계산량이 어떻게 다른지를 구현하고 확인한다. 링크 하나의 동역학은 3.1, 힘에서 출발하는 재귀 알고리즘은 3.2, 질량행렬·코리올리·중력항의 성질은 3.4, 효율적 알고리즘은 3.5, 라이브러리는 3.6, 유연 관절·링크는 3.9, 연속체는 3.10, 구속 동역학과 라그랑주 승수는 3.11, 부유 기저는 3.13, 에너지 기반 제어는 12.5와 이어진다.

배우는 것:

  • 오일러-라그랑주 방정식 — 일반화 좌표·일반화 힘, 가상일 원리에서의 유도.
  • 직렬 개연쇄의 에너지 — $K=\frac12\sum V_i^{T}G_iV_i$, 질량행렬 $\sum J_{ib}^{T}G_iJ_{ib}$, 균일 중력의 위치 에너지.
  • 크리스토펠 기호 — $\tau=M\ddot q+C\dot q+g(q)$의 전개와 $\dot M-2C$의 반대칭.
  • 에너지 균형 — $\dot E=\dot q^{T}\tau$, 수동성, 레일리 감쇠.
  • 오일러-푸앙카레 — 오일러 방정식의 변분 유도.
  • 정식화의 비교 — RNEA와 같은 토크, 곱셈 횟수와 기호식 크기의 증가.

선수 지식: 1.7(리군·리대수 I), 2.8(물체 야코비안), 3.1(공간 관성·오일러 방정식), 3.2(RNEA).

특히 중요한 독자: 로봇의 운동방정식을 기호로 유도하거나 제어기의 안정성을 에너지로 증명하는 연구자, 동역학 코드의 정확성을 독립적으로 검증하려는 개발자, 라그랑주와 뉴턴-오일러의 관계를 정리하려는 학생.

출처 방침. 기준 시점은 2026-11-05이다. 직접 읽은 자료: Lynch·Park, Modern Robotics(사전 인쇄본 2017-05판) 8.1절(라그랑주 정식화: 기본 개념과 예제·일반 정식화·질량행렬의 이해·라그랑주 대 뉴턴-오일러), 8.4절 앞부분(운동 에너지의 링크 합과 $M=\sum J_{ib}^{T}G_iJ_{ib}$)을 텍스트로 읽었다. 오일러-라그랑주 방정식의 유도는 이 자료가 "동역학 교재에 있다"고 적은 부분이며 이 장은 변분 계산을 직접 전개했다. 서지만 확인한 것(Crossref로 실재를 확인했고 본문은 읽지 못했다): Marsden·Ratiu(1994, Introduction to Mechanics and Symmetry), Uicker(1967, 1969, 공간 링크의 동적 해석), Hollerbach(1980, 재귀 라그랑주 정식화와 연산량 비교), Featherstone(2008), Carpentier 외(2019, Pinocchio), Todorov 외(2012, MuJoCo). 오일러-푸앙카레 방정식의 일반 이론과 문헌의 연산량 수치는 읽지 않았다. 확인하지 못한 것: 케인(Kane) 방법·기브스-아펠 방정식 등 다른 정식화의 문헌. 링크 치수·질량은 모두 이 장이 정한 예시이며 어떤 제품의 제원이 아니다. 수치는 모두 이 장의 순수 파이썬 계산이다(SymPy 기호식 크기는 라이브러리 대조에서 계산했다).

3.3.2 개념 — 두 스칼라 함수에서 운동방정식을 얻는다

결론부터 적으면 다음과 같다. 라그랑주 정식화는 계의 운동 에너지 $K(q,\dot q)$와 위치 에너지 $P(q)$라는 두 스칼라 함수에서 운동방정식을 얻는다. 뉴턴-오일러 방법이 링크마다 힘과 모멘트의 평형을 세우고 관절 구속력을 소거해야 하는 것과 달리, 일반화 좌표(관절 변수)를 그대로 사용하므로 구속력이 처음부터 방정식에 나타나지 않는다. 직렬 개연쇄에서는 $K=\frac12\dot q^{T}M(q)\dot q$이므로 방정식이 $M(q)\ddot q+C(q,\dot q)\dot q+g(q)=\tau$ 꼴로 나오며 코리올리·원심 항은 질량행렬의 편미분에서 나온다. 같은 로봇에 두 정식화를 적용하면 같은 관절 토크가 얻어진다(이 장이 수치로 확인한다). 차이는 용도와 계산량이다. 라그랑주 정식화는 에너지 균형·수동성처럼 제어기의 안정성을 증명하는 구조와 기호식 유도, 구속·연속체·유연체로의 확장에 적합하고, RNEA는 수치 계산이 관절 수에 선형이다. 직접 구현한 라그랑주 계산의 곱셈 횟수는 관절 수에 대해 지수가 $2.5$에서 $3.2$로 늘며 $n=24$에서 RNEA의 약 $940$배이다.

링크 Gi, 관절 축 Ai 에너지 경로 (3.3) 힘 경로 (3.2) K = q̇TM(q)q̇/2 Vi, V̇i (전방) P(q) = −∑mi g·pci GiV̇i − adTGiVi L = K − P Fi (후방 합) 오일러-라그랑주 방정식 τi = FiTAi 같은 관절 토크 τ
같은 관절 토크에 이르는 두 경로. 왼쪽 파란 경로는 운동 에너지 $K$와 위치 에너지 $P$에서 라그랑주 함수 $L=K-P$를 만들고 오일러-라그랑주 방정식을 전개한다. 오른쪽 주황 경로는 링크의 속도·가속도와 렌치를 전달하는 RNEA(3.2)이다. 두 경로는 링크 질량 특성 $G_i$와 관절 축 $A_i$에서 출발해 같은 관절 토크 $\tau$로 모인다(이 장의 계산에서 $3.7\times10^{-11}$ N·m 이내).

이 장의 계산 결과를 미리 적으면 다음과 같다. ① 오일러-라그랑주(수치 도함수)가 RNEA와 $n=2\sim6$의 임의 사슬에서 $3.7\times10^{-11}$ N·m 이하로 같고, $\sum J^{T}GJ$로 만든 질량행렬이 RNEA 열 구성과 $2.7\times10^{-15}$ 이하로 같으며 $\dot M-2C$가 $1.4\times10^{-12}$ 이하로 반대칭이다. ② 평면 2R의 닫힌 식(Modern Robotics 식 (8.9))과 $3.0\times10^{-12}$ N·m 안에서 같다. ③ 곱셈 횟수 비(라그랑주 대 RNEA)가 $n=2$의 $9.3$에서 $n=24$의 $943.5$로 커진다. ④ 수동 낙하 1 s의 에너지가 RK4 간격 $4\to2\to1$ ms에서 $15.0$·$15.5$분의 1로 줄어 4차 수렴을 하고($9.3\times10^{-9}$ 상대 변화 @ 1 ms), 레일리 감쇠의 에너지 손실이 $-\int\dot q^{T}B\dot q\,dt$와 $5.8\times10^{-8}$ J 안에서 같다. ⑤ 오일러 방정식의 해에서 오일러-푸앙카레 변분이 $-3.6\times10^{-14}$이고 해가 아닌 궤적에서는 $-2.0\times10^{-2}$이다. ⑥ Pinocchio·MuJoCo·modern_robotics·SymPy·SciPy와 $3.8\times10^{-11}$ 이하(RK4와 DOP853의 궤적 차이 제외)에서 같다.

3.3.3 원리

3.3.3.1 일반화 좌표와 오일러-라그랑주 방정식

계의 형상을 독립인 $n$개 변수 $q\in\mathbb R^{n}$(일반화 좌표)로 나타내고, 일률이 $\tau^{T}\dot q$가 되도록 일반화 힘 $\tau$를 정의한다. 개연쇄의 모든 관절이 구동되면 $q$를 관절 변수로 둘 수 있고 회전 관절의 $\tau$는 토크, 직동 관절의 $\tau$는 힘이다. 라그랑주 함수를 $L(q,\dot q)=K(q,\dot q)-P(q)$로 두면 운동방정식은 다음이다.

$$\frac{d}{dt}\frac{\partial L}{\partial\dot q}-\frac{\partial L}{\partial q}=\tau.$$

유도는 가상일 원리(라그랑주-달랑베르 원리)에서 나온다. 끝점을 고정한 임의의 변분 $\delta q(t_0)=\delta q(t_1)=0$에 대해 $\int_{t_0}^{t_1}(\delta L+\tau^{T}\delta q)\,dt=0$이고, $\delta L=\frac{\partial L}{\partial q}\delta q+\frac{\partial L}{\partial\dot q}\delta\dot q$에서 둘째 항을 부분적분하면 $\int\left(\frac{\partial L}{\partial q}-\frac{d}{dt}\frac{\partial L}{\partial\dot q}+\tau\right)^{T}\delta q\,dt=0$이 된다. $\delta q$가 임의이므로 괄호가 0이며 위 식이 얻어진다. 가상일을 하지 않는 구속력(이상 구속의 관절 반력)은 $\delta q$ 방향에서 일이 0이므로 식에 나타나지 않는다. 이 유도는 이 장이 직접 전개한 것이며 일반 이론은 Marsden·Ratiu(1994)에 있는 것으로 확인했으나 읽지 못했다.

3.3.3.2 직렬 개연쇄의 운동 에너지와 위치 에너지

링크 $i$의 운동 에너지는 3.1의 $\frac12V_i^{T}G_iV_i$이고 $V_i=J_{ib}(q)\dot q$($J_{ib}$는 프레임 $\{i\}$ 성분의 $6\times n$ 물체 야코비안으로 $i$열 뒤는 0)이므로 다음이 성립한다.

$$K=\frac12\dot q^{T}M(q)\dot q,\qquad M(q)=\sum_{i=1}^{n}J_{ib}^{T}(q)\,G_i\,J_{ib}(q),\qquad P(q)=-\sum_{i=1}^{n}m_i\,g\cdot p_{ci}(q).$$

$g$는 기저 프레임 성분의 중력 벡터, $p_{ci}(q)$는 링크 $i$의 질량 중심의 기저 프레임 위치이다. $G_i$가 양의 정부호이고 $J_{ib}$의 열이 선형 독립인 범위에서 $M$은 양의 정부호이다(Modern Robotics 8.1.3, 8.4). 질량행렬은 자세에 따라 달라지므로 같은 크기의 관절 가속도도 방향에 따라 다른 토크를 요구한다.

3.3.3.3 방정식의 전개와 크리스토펠 기호

$\partial L/\partial\dot q=M\dot q$이고 $\frac{d}{dt}(M\dot q)=M\ddot q+\dot M\dot q$이며 $\partial L/\partial q_k=\frac12\dot q^{T}\frac{\partial M}{\partial q_k}\dot q-\frac{\partial P}{\partial q_k}$이므로 $k$번째 식은 다음이다.

$$\tau_k=\sum_jM_{kj}\ddot q_j+\sum_{i,j}\Gamma_{kij}\dot q_i\dot q_j+\frac{\partial P}{\partial q_k},\qquad \Gamma_{kij}=\frac12\left(\frac{\partial M_{kj}}{\partial q_i}+\frac{\partial M_{ki}}{\partial q_j}-\frac{\partial M_{ij}}{\partial q_k}\right).$$

$\Gamma_{kij}$를 제1종 크리스토펠 기호라 하며 코리올리 행렬은 $C_{kj}=\sum_i\Gamma_{kij}\dot q_i$, 중력항은 $g_k(q)=\partial P/\partial q_k$이다. 따라서 $\tau=M\ddot q+C\dot q+g(q)$이다. $i=j$인 항이 원심 항, $i\ne j$인 항이 코리올리 항이다(Modern Robotics 8.1.2). 이 분해에서 $\dot M-2C$의 $(k,j)$ 성분은 $\sum_i\left(\frac{\partial M_{ij}}{\partial q_k}-\frac{\partial M_{ik}}{\partial q_j}\right)\dot q_i$이고 $k$와 $j$를 바꾸면 부호가 바뀌므로 $\dot M-2C$는 반대칭이다. 이 성질의 의미와 제어에서의 사용은 3.4의 주제이며 이 장은 다음 소절의 에너지 균형에 사용한다.

3.3.3.4 에너지 균형·수동성·감쇠

총 에너지 $E=K+P$를 시간으로 미분하면 $\dot E=\dot q^{T}M\ddot q+\frac12\dot q^{T}\dot M\dot q+\dot q^{T}g(q)$이고 운동방정식 $M\ddot q=\tau-C\dot q-g$를 대입하면 다음이다.

$$\dot E=\dot q^{T}\tau+\frac12\dot q^{T}\left(\dot M-2C\right)\dot q=\dot q^{T}\tau.$$

계에 들어가는 일률이 곧 에너지의 증가율이다. 입력 $\tau$, 출력 $\dot q$인 계가 저장함수 $E$에 대해 수동적이며(위치 에너지가 아래로 유계일 때), 이 구조가 에너지 기반 제어기의 안정성 증명에 사용된다(3.4, 12.5). 마찰·감쇠처럼 에너지를 소산하는 힘은 레일리 소산 함수 $\mathcal R=\frac12\dot q^{T}B\dot q$($B\succeq0$)로 넣으며 $\tau_{\mathrm{total}}=\tau-B\dot q$이므로 $\dot E=\dot q^{T}\tau-\dot q^{T}B\dot q$이다. 수동 낙하에서 $E$가 보존되는지가 모델과 적분기의 정확성을 점검하는 지표이다(풀이 예제 ④).

3.3.3.5 오일러-푸앙카레 — 오일러 방정식을 에너지에서 얻는다

단일 강체의 라그랑주 함수는 회전 운동 에너지 $l(\omega)=\frac12\omega^{T}I\omega$이며 변수는 자세 $R\in SO(3)$이다. 물체 각속도는 $\hat\omega=R^{T}\dot R$이고 변분 $\delta R=R\hat\eta$($\eta$는 물체 프레임 성분의 변분, 끝점에서 0)에 대해 $\delta\hat\omega=\dot{\hat\eta}+[\hat\omega,\hat\eta]$, 곧 $\delta\omega=\dot\eta+\omega\times\eta$이다. 작용 $A=\int l\,dt$의 변분은 다음이다.

$$\delta A=\int(I\omega)\cdot(\dot\eta+\omega\times\eta)\,dt=\int\left[-\frac{d}{dt}(I\omega)+(I\omega)\times\omega\right]\cdot\eta\,dt.$$

$(I\omega)\cdot(\omega\times\eta)=\eta\cdot((I\omega)\times\omega)$를 사용했다. $\eta$가 임의이므로 정상 조건은 $I\dot\omega=I\omega\times\omega=-\omega\times I\omega$, 곧 3.1의 오일러 방정식이다. 토크가 있으면 우변에 일반화 힘 $\tau$가 더해진다. 이 장은 보조정리 $\delta\omega=\dot\eta+\omega\times\eta$를 유한 차분으로, 오일러 방정식의 해에서 $\delta A=0$임을 구적으로 확인했다(풀이 예제 ⑤). 일반 리군에서의 오일러-푸앙카레 방정식은 Marsden·Ratiu(1994)의 주제이며 읽지 못했다.

3.3.3.6 라그랑주와 뉴턴-오일러의 관계, 확장

두 정식화는 같은 운동방정식을 낸다(풀이 예제 ①). 뉴턴-오일러는 관절 반력 $F_i$를 부산물로 주고(3.2), 라그랑주는 일반화 좌표 방향의 토크 $\tau$만 주며 반력은 나타나지 않는다. 일반화 좌표가 독립이 아닌 폐연쇄·접촉 구속은 구속을 라그랑주 승수로 방정식에 넣어 처리하며 3.11에서 다룬다. 유연 관절·유연 링크·연속체(3.9, 3.10)와 부유 기저(3.13)로의 확장도 해당 장의 주제이고 이 장은 구현하지 않았다.

3.3.3.7 계산량과 기호식의 크기

수치 도함수로 $\Gamma$를 만드는 구현은 $\partial M/\partial q_k$를 얻으려고 $k$마다 $M$을 여러 번 평가한다. $M$ 하나의 조립이 $O(n^{3})$이므로 전체가 $O(n^{4})$에 접근한다. 이 장의 구현은 5점 중앙 차분이라 $4n+1$번의 $M$ 평가를 하고, 측정한 지수는 $n$ 사이에서 $2.5$에서 $3.2$로 늘어난다(풀이 예제 ③). 기호식은 항 수가 급증한다. 평면 $n$링크의 오일러-라그랑주 식(미정리)의 연산 수는 $n=2,3,4$에서 $492$·$2{,}449$·$8{,}097$이다(SymPy count_ops, 라이브러리 대조에서 계산). 재귀적으로 정식화해 연산량을 줄인 라그랑주 방법은 Hollerbach(1980)가 다룬 것으로 서지를 확인했으나 이 장은 그 내용을 읽지 못했고 구현하지 않았다 (확인 필요).

3.3.4 알고리즘

입력: 링크 {M_{i,i−1}, A_i, G_i}, 중력 g, q, q̇, q̈
링크 자세:     T_{i,i−1} ← exp(−[A_i] q_i) M_{i,i−1};   T_{0,i} ← Π T_{k−1,k}
물체 야코비안: J_i[:, i] ← A_i;   J_i[:, j] ← Ad_{T_{i,j}} A_j   (j < i)
질량행렬:      M(q) ← Σ_i J_iᵀ G_i J_i
위치 에너지:   P(q) ← −Σ_i m_i g · p_ci(q)
도함수:        ∂M/∂q_k, ∂P/∂q_k ← 5점 중앙 차분
토크:          τ_k ← Σ_j M_kj q̈_j + Σ_ij Γ_kij q̇_i q̇_j + ∂P/∂q_k
에너지:        E ← ½ q̇ᵀ M q̇ + P

3.3.5 풀이 예제 — 오일러-라그랑주와 RNEA의 일치, 에너지 균형, 오일러-푸앙카레

결론부터 적으면 다음과 같다. ① 오일러-라그랑주와 RNEA가 같은 토크를 낸다. ② 평면 2R 닫힌 식과 일치한다. ③ 라그랑주 계산의 곱셈 횟수가 RNEA의 $9.3\sim943.5$배이다. ④ 수동 낙하에서 에너지가 4차 수렴으로 보존되고 감쇠 손실이 소산 일률의 적분과 같다. ⑤ 오일러 방정식의 해가 작용의 정상점이다. ⑥ 라이브러리 5종과 일치한다.

구현. link_poses·body_jacobians·mass_matrix_el·potential·grad·euler_lagrange·energy. 사슬 생성·RNEA·질량행렬은 3.2의 구현을, 행렬 연산은 3.1의 구현을 사용한다. 크리스토펠 $C$·수동 낙하 적분(simulate)·곱셈 횟수 계수·오일러-푸앙카레 구적·도해는 스크립트에 있다.

항목값
입력 1임의 사슬: 3.2 입력 2와 같은 생성(수정 DH, 점질량 12개로 만든 관성), $n=2\sim6$, 각 30개, 상태 $q\sim U(-\pi,\pi)$·$\dot q\sim U(-2,2)$·$\ddot q\sim U(-4,4)$(시드 65)
입력 2평면 2R(끝점 질량 $m_1=1.5$, $m_2=0.8$ kg, $L_1=0.5$, $L_2=0.4$ m, 중력 $-y$), 상태 100개
입력 3곱셈 횟수: $n=2,3,4,6,8,12,16,24$ 임의 사슬(dot·mv·mm·cross의 스칼라 곱셈)
입력 4수동 낙하: 3.2의 예제 6R 팔(입력 5)의 앞 세 관절(질량 4·4·3 kg), $q(0)=(0.3,-0.2,0.9)$ rad, 정지에서 출발, 중력 $-z$, RK4 $h=4,2,1$ ms·1 s. 감쇠 $B=\mathrm{diag}(0.8,0.8,0.3)$ N·m·s/rad
입력 5오일러-푸앙카레: 직육면체($0.30\times0.20\times0.10$ m, 2 kg) $I=\mathrm{diag}(0.00833,0.01667,0.02167)$ kg·m², $\omega(0)=(1,6,2)$ rad/s, $\eta(t)=\sin(\pi t/T)\,(1,0.5,-0.3)$, $T=1$ s, 심프슨 $h=1$ ms. 보조정리: $R(t)=\exp[(0.7\sin t+0.2,\,-0.4\sin t+0.5,\,0.9\sin t-0.3)]$, $\eta(t)=(0.8\sin t,\,0.3\cos2t,\,-0.5t)$, $\varepsilon=10^{-4}$
입력 6라이브러리 대조: 임의 6R 20개×상태 3개, SymPy 평면 $n=2,3,4$($L_i=0.4+0.05(i-1)$ m, $m_i=2.0-0.3(i-1)$ kg 균일 막대), DOP853 rtol $10^{-13}$ (시드 651)
def link_poses(ch, q):
    """T_{k,k-1}(q_k) 목록과 T_{0,k}(q) 목록"""
    E3 = [[1.0 if i == j else 0.0 for j in range(3)]
          for i in range(3)]
    Ts, T0, T = [], [], (E3, [0.0] * 3)
    for k in range(len(q)):
        E = screw_exp(ch['A'][k], -q[k])
        Tk = pose_mul(E, ch['M'][k])
        Ts.append(Tk)
        T = pose_mul(T, inv_pose(Tk))
        T0.append(T)
    return Ts, T0

def body_jacobians(ch, q):
    """J[i][j]: 링크 i 물체 야코비안의 j열, V_i = J[i] qd"""
    n = len(q)
    Ts, _ = link_poses(ch, q)
    J = []
    for i in range(n):
        cols = [[0.0] * 6 for _ in range(n)]
        cols[i] = list(ch['A'][i])
        T = Ts[i]
        for j in range(i - 1, -1, -1):
            cols[j] = mv(adjoint(*T), ch['A'][j])
            T = pose_mul(T, Ts[j])
        J.append(cols)
    return J

def mass_matrix_el(ch, q):
    """M(q) = sum_i J_i^T G_i J_i"""
    n = len(q)
    M = [[0.0] * n for _ in range(n)]
    for i, cols in enumerate(body_jacobians(ch, q)):
        GJ = [mv(ch['G'][i], c) for c in cols]
        for a in range(i + 1):
            for b in range(i + 1):
                M[a][b] += dot(cols[a], GJ[b])
    return M

def potential(ch, q, g):
    """P = -sum m_i g . p_ci (g: 기저 프레임 중력 벡터)"""
    _, T0 = link_poses(ch, q)
    P = 0.0
    for T, (m, c, _) in zip(T0, ch['links']):
        p = [a + b for a, b in zip(mv(T[0], c), T[1])]
        P -= m * dot(g, p)
    return P

def grad(f, x, h=1e-3):
    """5점 중앙 차분"""
    out = []
    for k in range(len(x)):
        v = []
        for s in (2, 1, -1, -2):
            y = list(x)
            y[k] += s * h
            v.append(f(y))
        out.append((-v[0] + 8 * v[1] - 8 * v[2]
                    + v[3]) / (12 * h))
    return out

def euler_lagrange(ch, q, qd, qdd, g):
    """tau = M qdd + Gam qd qd + dP/dq (크리스토펠)"""
    n = len(q)
    M = mass_matrix_el(ch, q)
    dM = []                       # dM[k][a][b] = dM_ab/dq_k
    for k in range(n):
        v = []
        for s in (2, 1, -1, -2):
            y = list(q)
            y[k] += s * 1e-3
            v.append(mass_matrix_el(ch, y))
        dM.append([[(-v[0][a][b] + 8 * v[1][a][b]
                     - 8 * v[2][a][b] + v[3][a][b]) / 0.012
                    for b in range(n)] for a in range(n)])
    dP = grad(lambda y: potential(ch, y, g), q)
    tau = []
    for k in range(n):
        c = 0.0
        for i in range(n):
            for j in range(n):
                gam = 0.5 * (dM[i][k][j] + dM[j][k][i]
                             - dM[k][i][j])
                c += gam * qd[i] * qd[j]
        tau.append(dot(M[k], qdd) + c + dP[k])
    return tau, M, dM

def energy(ch, q, qd, g):
    M = mass_matrix_el(ch, q)
    return 0.5 * dot(qd, mv(M, qd)) + potential(ch, q, g)

① 오일러-라그랑주 대 RNEA.

관절 수 $n$오일러-라그랑주 $\tau$ 대 RNEA$M=\sum J_i^{T}G_iJ_i$ 대 RNEA 열크리스토펠 $C\dot q$ 대 RNEA ($g=0$, $\ddot q=0$)$\dot M-2C$의 비반대칭 $\lvert S+S^{T}\rvert$
25.5e-120.0e+003.6e-131.1e-16
38.0e-128.9e-161.0e-121.8e-13
41.1e-111.8e-152.7e-124.4e-13
51.7e-112.7e-153.8e-121.2e-12
63.7e-112.7e-157.6e-121.4e-12

세 방법이 같은 값을 낸다. 오일러-라그랑주의 오차 $10^{-11}$은 5점 차분(간격 $10^{-3}$)의 절단·반올림 오차로 추정하며(간격을 바꿔 확인하지 않았다) 차분을 사용하지 않는 질량행렬(둘째 열)은 $10^{-15}$이다. $\dot M-2C$의 반대칭은 크리스토펠 기호에서 나오는 구조이며 $10^{-12}$ 수준으로 확인된다.

② 평면 2R 닫힌 식.

평면 2R(끝점 질량) 닫힌 식 Modern Robotics (8.9) 대 이 장의 오일러-라그랑주, 임의 상태 100개최대 오차
$m_1=1.5$, $m_2=0.8$ kg, $L_1=0.5$, $L_2=0.4$ m, 중력 $-y$3.0e-12 N·m

Modern Robotics 식 (8.9)와 3차원 구현이 일치하므로 일반 구현의 부호·프레임 규약이 맞는다. 같은 식을 SymPy의 기호 오일러-라그랑주로 유도한 결과와도 $5.3\times10^{-15}$로 일치한다(⑥).

③ 곱셈 횟수.

관절 수 $n$RNEA 곱셈 횟수오일러-라그랑주(수치 도함수) 곱셈 횟수비율인접 $n$ 사이 지수 $\log(\text{횟수비})/\log(n\text{비})$
27146,6319.3—
31,07118,38117.22.51
41,42839,17827.42.63
62,142119,33755.72.75
82,856271,87695.22.86
124,284908,526212.12.98
165,7122,208,712386.73.09
248,5688,083,788943.53.20
2 3 4 6 8 12 24 103 104 105 106 107 108 관절 수 n 곱셈 횟수 RNEA 오일러-라그랑주(수치 도함수)
관절 수 $n$에 대한 곱셈 횟수(양대수 축). RNEA(파랑)는 $n$에 정비례하고, 수치 도함수로 구현한 오일러-라그랑주(주황)는 $n$이 클수록 기울기가 커져 $n=24$에서 RNEA의 약 940배이다.

RNEA의 곱셈은 $n$에 정비례하고($357n$) 오일러-라그랑주는 지수가 $n$과 함께 $2.51$에서 $3.20$으로 커진다. 점근 지수는 $4$로 예상되나 $n=24$까지만 측정했다. 이 값은 이 장의 구현(5점 차분, $\sum J^{T}GJ$를 전체 행렬 곱으로 조립)에서 센 것이며 재귀적 라그랑주나 해석적 도함수를 사용하면 줄어들 것으로 예상한다(측정하지 않았다). 실행 시간은 측정하지 않았다.

④ 에너지 균형.

수동 낙하 1 s (앞 세 관절, $\tau=0$), RK4 간격 $h$ (s)$\max\lvert E(t)-E(0)\rvert$ (J)상대 ($/\lvert E(0)\rvert$)
0.0045.56e-052.14e-06
0.0023.71e-061.43e-07
0.0012.40e-079.26e-09
0.00 0.25 0.50 0.75 1.00 0 10 20 30 시간 (s) 에너지 (J) 운동 K 위치 P 합 K + P
예제 팔의 앞 세 관절이 정지에서 중력으로 낙하하는 1 s 동안의 에너지(RK4, $h=10$ ms). 운동 에너지(파랑)와 위치 에너지(주황 파선)가 서로 바뀌고 합 $K+P$(초록 점선)가 약 $26.0$ J로 유지된다.

수동 낙하의 에너지 변화가 RK4 간격 $4\to2\to1$ ms에서 $15.0$·$15.5$분의 1로 줄어 4차 방법의 이론 비 $16$에 가깝다. 보존이 정확하다는 것이 아니라 오차가 적분 오차로 설명된다는 뜻이며, 운동방정식(RNEA와 질량행렬 풀이)과 에너지(라그랑주의 $\frac12\dot q^{T}M\dot q+P$)를 서로 다른 경로로 계산했으므로 두 정식화의 일치도 같이 확인한다. 도해에서 운동 에너지(최대 $23.3$ J)와 위치 에너지(최대 $26.0$ J)가 서로 바뀌고 합이 $26.0$ J로 유지된다.

레일리 감쇠 $B=\mathrm{diag}(0.8,0.8,0.3)$ N·m·s/rad, 1 s, $h=1$ ms값
에너지 변화 $E(1)-E(0)$ (J)-22.358752
$-\int\dot q^{T}B\dot q\,dt$ (심프슨, J)-22.358752
차이 (J)5.8e-08

감쇠가 있는 경우 에너지 변화 $-22.358752$ J이 $-\int\dot q^{T}B\dot q\,dt=-22.358752$ J과 $5.8\times10^{-8}$ J 안에서 같다. $\dot E=\dot q^{T}\tau-\dot q^{T}B\dot q$가 맞음을 보인다.

⑤ 오일러-푸앙카레 변분.

단일 강체의 오일러-푸앙카레 변분 ($I=\mathrm{diag}(0.00833,0.01667,0.02167)$, $\omega(0)=(1,6,2)$ rad/s, $T=1$ s)값
보조정리 $\partial_\varepsilon\omega_\varepsilon$ 대 $\dot\eta+\omega\times\eta$ (세 시각, 유한 차분) 최대 차이7.2e-10
오일러 방정식 해에서 $\delta A=\int I\omega\cdot(\dot\eta+\omega\times\eta)\,dt$-3.6e-14 (적분 규모 9.5e-02)
같은 $\eta$에서 오일러 방정식을 만족하지 않는 $\omega(t)=\omega(0)$ (일정)-2.0e-02

보조정리 $\delta\omega=\dot\eta+\omega\times\eta$가 $7.2\times10^{-10}$ 안에서 성립하고, 오일러 방정식의 해에서는 작용의 변분이 적분 규모($9.5\times10^{-2}$)의 $10^{-12}$ 배 이하($-3.6\times10^{-14}$)인 반면 오일러 방정식을 만족하지 않는 일정한 $\omega$에서는 $-2.0\times10^{-2}$이다. 오일러 방정식이 작용의 정상 조건임을 수치로 보이는 것이다. 이는 하나의 $\eta$에 대한 확인이며 모든 $\eta$에 대한 증명은 위 변분 계산이 한다.

⑥ 라이브러리 대조(실행 2026-11-05). Pinocchio 4.1.0의 computeKineticEnergy·computePotentialEnergy·computeMechanicalEnergy·computeCoriolisMatrix·computeGeneralizedGravity·aba, MuJoCo 3.17.0의 mj_forward와 에너지 플래그(energy), modern_robotics의 MassMatrix·VelQuadraticForces·GravityForces, SymPy 1.14.0의 기호 오일러-라그랑주와 count_ops, SciPy 1.18.1의 solve_ivp(DOP853)를 같은 입력에서 비교했다. 스크립트는 robotics/fig-src/3.3-libcheck.py, 출력은 robotics/fig-src/data/3.3-libcheck-output.txt에 있다.

비교결과
임의 6R 사슬 20개×상태 3개: 운동 에너지 $\tfrac12\dot q^{T}M\dot q$ 대 Pinocchio computeKineticEnergy: 최대 차이 (J)8.9e-15
같은 입력: 위치 에너지 $P$ 대 Pinocchio computePotentialEnergy · 총 에너지 대 computeMechanicalEnergy: 최대 차이 (J)4.3e-14 · 4.3e-14
같은 입력: 크리스토펠 코리올리 행렬 $C$ 대 Pinocchio computeCoriolisMatrix: 최대 차이4.8e-12
같은 입력: $\partial P/\partial q$ 대 Pinocchio computeGeneralizedGravity: 최대 차이 (N·m)3.7e-11
같은 입력: 운동·위치 에너지 대 MuJoCo energy(mj_forward, 에너지 플래그): 최대 차이 (J)2.1e-14 · 5.7e-14
같은 입력: $M$ · $C\dot q$ · $\partial P/\partial q$ 대 modern_robotics MassMatrix · VelQuadraticForces · GravityForces: 최대 차이2.8e-14 · 1.0e-11 · 3.8e-11
SymPy 기호 오일러-라그랑주(평면 $n=2,3,4$링크 균일 막대, 상태 20개)를 대입한 $\tau$ 대 이 장의 RNEA: 최대 차이 (N·m)1.4e-14
SymPy 기호 오일러-라그랑주(평면 2R 끝점 질량, 상태 100개) 대 Modern Robotics 식 (8.9): 최대 차이 (N·m)5.3e-15
앞 세 관절 수동 낙하 1 s: SciPy DOP853(rtol $10^{-13}$, Pinocchio aba) 대 이 장의 RK4(RNEA·질량행렬 풀이)의 상태 최대 차이, $h=4,2,1$ ms1.1e-04 · 8.2e-06 · 5.6e-07
같은 DOP853 해: Pinocchio 총 에너지의 상대 변화 · 이 장의 에너지 식 대 Pinocchio 에너지1.0e-14 · 2.1e-14 J

기호식 크기는 평면 $n=2,3,4$링크의 미정리 오일러-라그랑주 식에서 492,·2,449,·8,097(SymPy count_ops)이며 같은 $n$의 RNEA 곱셈은 $714$·$1{,}071$·$1{,}428$이다. 에너지·코리올리 행렬·중력 벡터·질량행렬·SymPy 식이 모두 $3.8\times10^{-11}$ 이하에서 같다. 궤적의 상태 차이($1.1\times10^{-4}$·$8.2\times10^{-6}$·$5.6\times10^{-7}$)는 RK4 간격에 따른 4차 수렴(감소비 $13.4$·$14.6$)으로 DOP853과 같은 해에 수렴하며, DOP853 해에서 Pinocchio의 총 에너지가 $1.0\times10^{-14}$ 상대 변화로 보존되고 이 장의 에너지 식이 Pinocchio 값과 $2.1\times10^{-14}$ J 안에서 같다. 이상 징후는 없었다. Pinocchio의 코리올리 행렬은 이 장의 크리스토펠 형태와 같은 행렬로 확인되었다($C\dot q$가 같은 다른 $C$ 정의가 있을 수 있으나 $\dot M-2C$의 반대칭을 이 장이 확인한 것은 크리스토펠 형태이므로 행렬 자체를 비교했다).

한계. 이상 구속·균일 중력·1자유도 관절의 직렬 개연쇄이다. 마찰·회전자 관성·탄성은 다루지 않았고(3.8, 3.9) 구속 승수(3.11)와 부유 기저(3.13)는 구현하지 않았다. 도함수는 수치 차분이며 해석적 도함수·자동 미분은 3.6이 다룬다. 오일러-푸앙카레는 SO(3)의 단일 강체만 확인했고 일반 리군 이론은 읽지 못했다. 에너지 균형은 앞 세 관절 하나의 궤적에서 확인했다.

견해: 운동방정식을 기호로 유도하거나 제어기의 안정성을 에너지로 증명할 때는 라그랑주 정식화를, 실시간 토크 계산에는 RNEA(3.2)를 사용하는 것이 맞다고 본다. 근거는 이 장에서 두 방법이 같은 토크를 내므로(최대 $3.7\times10^{-11}$ N·m) 선택은 용도와 비용의 문제이고, 라그랑주 구현이 RNEA보다 곱셈이 $n=24$에서 약 $940$배 많으며 기호식이 $n=2\to4$에서 $492\to8{,}097$로 커진다는 점이다. 반면 에너지 균형·수동성은 라그랑주 형태에서 한 줄로 얻어진다. 이 장은 재귀 라그랑주나 해석적 도함수 같은 효율적 라그랑주 방법은 구현하지 않았으므로 "라그랑주는 느리다"를 일반 결론으로 삼지 않는다.

3.3.6 설계 트레이드오프 — 운동방정식을 얻는 방법

로봇의 운동방정식을 얻는 방법은 다음 중에서 고른다.

  • ① 라그랑주 기호 유도. $K$·$P$를 기호로 적고 오일러-라그랑주를 전개한다. 저자유도에서 식이 짧고 구조(수동성·대칭)가 드러난다. 자유도가 늘면 항이 급증한다.
  • ② 라그랑주 수치 구현. $M=\sum J^{T}GJ$와 도함수를 수치로 평가한다. 모델 변경이 쉽고 에너지를 같이 얻는다. 연산량이 크다.
  • ③ 뉴턴-오일러 재귀(RNEA). 링크당 일정한 연산이며 관절 반력을 부산물로 준다(3.2). 에너지 구조는 따로 구성해야 한다.
방법연산량부산물적합한 경우
① 기호 유도$n$에 따라 식 급증구조·안정성 증명3자유도 이하의 해석
② 수치 라그랑주$n=24$에서 RNEA의 약 940배(이 구현)에너지·$C$검증·교육
③ RNEA$357n$관절 반력실시간 계산

견해: 검증 목적에는 ③의 결과를 ②나 SymPy ①로 독립 확인하는 것이 유용하다고 본다. 근거는 이 장에서 두 경로(힘·에너지)를 독립으로 구현해 $3.7\times10^{-11}$ N·m 안에서 일치했으므로 한쪽의 부호·프레임 오류가 다른 쪽과의 비교에서 드러난다는 점이다. 방법 사이의 실행 시간 비교는 이 장이 수행하지 않았다.

3.3.7 응용 축 비교

응용 축라그랑주 정식화가 맡는 일요구관련 장
산업·제조저자유도 기구의 기호 모델·에너지 기반 제어 안정성 증명모델의 에너지 구조12.5
물류·이동해당하는 고유 요구 없음——
도로 주행해당하는 고유 요구 없음——
서비스해당하는 고유 요구 없음——
의료·헬스케어원리 수준에서 해당하는 고유 요구를 확인하지 못했다——
휴머노이드·범용접촉 구속과 부유 기저를 포함한 승수 정식화구속 조건 처리3.11, 3.12, 3.13
필드유연 링크·대형 팔 모델로의 확장(해당 장 참조, 이 장은 구현하지 않았다)탄성 모델3.9, 7.6
국방·안보원리 수준에서 해당하는 고유 요구 없음——

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

3.3.8 관련 규격·프로토콜

문서발행무엇을 정하나상세 장
computeKineticEnergy·computePotentialEnergy·computeCoriolisMatrixPinocchio 4.1.0에너지와 코리올리 행렬3.6
MuJoCo 에너지 플래그MuJoCo 3.17.0energy에 퍼텐셜·운동 에너지3.6
MassMatrix·VelQuadraticForces·GravityForcesmodern_robotics$M$·$C\dot q$·$g$ 분해3.6
SymPy count_ops·diffSymPy 1.14.0기호 오일러-라그랑주와 식 크기1.18
수학 표기·기호표이 책$M\ddot q+C\dot q+g=\tau$A.2

에너지 정식화를 규정한 국제 표준은 확인하지 못했다 (확인 필요).

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

라그랑주 정식화 자체를 규정하는 법령·인증 요건은 확인하지 못했다. 협동 운전에서 충돌 시 전달되는 에너지와 유효 질량이 인체 접촉 한계의 입력이 될 수 있고(3.4, 26.9), 이 값을 계산하는 질량행렬이 이 장의 $M(q)$이다. 이를 규정하는 규격 조항은 확인하지 않았다 (확인 필요).

3.3.10 한계·유의점·미해결

  • 모델 범위. 이상 구속·균일 중력·직렬 개연쇄이다. 마찰·회전자 관성·감속기 탄성은 포함하지 않았다.
  • 구현의 효율. 수치 차분·전체 행렬 곱 구현이며 재귀 라그랑주·해석적 도함수의 연산량은 측정하지 않았다. 문헌의 비교는 읽지 못했다.
  • 수치 도함수의 정확도. 5점 차분으로 $10^{-11}$ 수준이며 극단적 조건(특이 자세·큰 질량비)에서의 정확도는 시험하지 않았다.
  • 구속과 확장. 라그랑주 승수, 접촉, 부유 기저, 연속체는 구현하지 않았다(3.11~3.13, 3.10).
  • 일반 이론. 오일러-푸앙카레의 일반 리군 이론과 다른 정식화(Kane·Gibbs-Appell)는 확인하지 못했다 (확인 필요).

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

체크포인트:

  1. 가상일 원리에서 오일러-라그랑주 방정식을 유도하라.
  2. $K=\frac12\dot q^{T}M\dot q$에서 $\tau_k=\sum M_{kj}\ddot q_j+\sum\Gamma_{kij}\dot q_i\dot q_j+\partial P/\partial q_k$를 전개하라.
  3. $\dot M-2C$가 반대칭임을 보이고 이로부터 $\dot E=\dot q^{T}\tau$를 유도하라.
  4. 레일리 소산 함수가 있을 때 에너지 변화율을 구하라.
  5. $\delta\omega=\dot\eta+\omega\times\eta$를 이용해 오일러 방정식을 유도하라.
  6. 라그랑주 정식화에서 관절 반력이 나타나지 않는 이유와 폐연쇄에서 승수가 필요한 이유를 설명하라.
  7. 평면 2R에서 $M$·$C$·$g$를 손으로 전개해 Modern Robotics 식 (8.9)를 얻어라.

연구 질문:

  • 재귀 라그랑주 방법과 RNEA의 곱셈 횟수·실행 시간을 같은 사슬에서 비교하면 어떤 $n$에서 교차하는가(이 장은 구현하지 않았다)?
  • 해석적 도함수나 자동 미분(3.6)으로 $\partial M/\partial q$를 구하면 이 장의 수치 차분 구현의 연산량과 정확도는 어떻게 달라지는가?
  • 에너지 기반 제어기의 수동성 조건이 관절 마찰과 모델 오차에서 어느 정도까지 유지되는가?

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

  • K. M. Lynch, F. C. Park, Modern Robotics: Mechanics, Planning, and Control, Cambridge University Press, 2017, doi:10.1017/9781316661239 — 8.1절·8.4절 앞부분을 사전 인쇄본(2017-05판) 텍스트로 읽었다.
  • J. E. Marsden, T. S. Ratiu, Introduction to Mechanics and Symmetry, Springer, 1994, doi:10.1007/978-1-4612-2682-6 — 서지만 확인.
  • J. J. Uicker, "Dynamic force analysis of spatial linkages," J. Applied Mechanics 34(2), 418–424, 1967, doi:10.1115/1.3607699; "Dynamic behavior of spatial linkages: Part 1—Exact equations of motion; Part 2—Small oscillations about equilibrium," J. Engineering for Industry 91(1), 251–265, 1969, doi:10.1115/1.3591539 — 서지만 확인.
  • J. M. Hollerbach, "A recursive Lagrangian formulation of manipulator dynamics and a comparative study of dynamics formulation complexity," IEEE Trans. Systems, Man, and Cybernetics 10(11), 730–736, 1980, doi:10.1109/tsmc.1980.4308393 — 서지만 확인.
  • R. Featherstone, Rigid Body Dynamics Algorithms, Springer, 2008, doi:10.1007/978-1-4899-7560-7; J. Carpentier 외, "The Pinocchio C++ library," 2019, doi:10.1109/sii.2019.8700380; E. Todorov, T. Erez, Y. Tassa, "MuJoCo," IROS, 2012, doi:10.1109/iros.2012.6386109 — 서지만 확인.
  • Pinocchio 4.1.0, MuJoCo 3.17.0, modern_robotics 1.1.1, SymPy 1.14.0, SciPy 1.18.1(integrate.solve_ivp), numpy 2.5.4.

다음 장. 3.4은 매니퓰레이터 운동방정식의 질량행렬·코리올리·중력항이 갖는 성질을 다룬다. 이 장의 $M(q)$·$C(q,\dot q)$·$g(q)$가 대칭·양의 정부호, 반대칭 성질, 매개변수 선형성 같은 구조를 가지며 유효 질량·동적 조작성과 어떻게 이어지는지를 정리한다.