11장. 힘이 움직임을 만든다: 동역학 첫걸음
모터가 토크를 내면 팔은 어떻게 움직일까? 거꾸로, 이렇게 움직이려면 토크가 얼마나 들까? 이 장은 강체 하나의 뉴턴·오일러 방정식과 관성 텐서를 익히고, 링크를 따라 바깥으로 갔다가 안으로 돌아오며 관절 토크를 구한다. 강의 11강(74분)에 해당한다.
11.1 동역학은 무엇을 묻나
지금까지는 팔의 모양과 움직임만 다뤘다. 관절 각도로 손끝 위치를 구하고(순기구학), 관절 속도로 손끝 속도를 구했다(자코비안). 어디에도 질량이나 힘은 나오지 않았다. 이것이 기구학(kinematics)이다.
이제 팔에 무게가 있다고 본다. 모터가 토크를 내면 링크는 얼마나 빨리 가속할까? 이 물음을 다루는 것이 동역학(dynamics)이다.
1.4절에서 "힘을 주면 팔이 어떻게 움직이나"라고 이름만 붙인 네 번째 질문이 바로 이것이다.
로봇 팔의 동역학은 관절 변수로 적는다. 관절 토크 τ와 관절 각도 q, 각속도 q̇, 각가속도 q̈ 사이의 관계다. 이것을 관절 공간 동역학(joint space dynamics)이라 한다.
두 방향의 물음
- 순동역학(forward dynamics): 토크 τ를 주면 가속도 q̈는 얼마인가? q̈를 조금씩 쌓으면 속도와 각도가 나온다. 컴퓨터 안에서 로봇을 움직여 보는 시뮬레이션이 이것이다.
- 역동역학(inverse dynamics): 이런 q, q̇, q̈로 움직이려면 관절마다 토크가 얼마 들까? 10장에서 짠 궤적을 따라가게 하는 제어와, 충분히 센 모터 고르기에 쓴다.
잠깐, 생각해 보기. 서보를 고른다. 사양표의 "최대 토크"와 견줄 숫자는 순동역학과 역동역학 가운데 어느 쪽에서 나올까?
역동역학이다. 팔이 해야 할 가장 힘든 움직임(가장 빠른 동작, 가장 무거운 짐)을 넣고 토크를 구한다. 그 값보다 넉넉히 센 서보를 고른다.
왜 어려운가
링크 하나를 돌리는 것은 쉽다. 그런데 팔은 링크가 사슬처럼 이어져 있어서, 한 관절이 움직이면 다른 링크를 밀고 당긴다. 강의 첫머리에 나오는 휴머노이드가 손을 조금만 움직여도 몸 전체에 힘이 퍼지는 것도 같은 까닭이다.
그래서 두 단계로 나눈다. 먼저 강체(rigid body) 하나의 동역학을 익히고, 그다음 강체들을 관절로 이어 붙인다. 이어 붙이는 방법은 여러 가지인데, 강의는 두 가지를 본다.
- 뉴턴-오일러(Newton–Euler) 방법: 링크마다 힘과 모멘트를 따진다. 8장의 정적 힘 전파에 관성력을 더한 것이다(11.8~11.9절).
- 라그랑주(Lagrange) 방법: 팔 전체의 에너지로 한 번에 식을 얻는다(11.10절, 12장).
무거운 관절은 흔들림을 걸러 낸다
강의는 80 kg을 드는 6관절 팔의 실험을 보여 준다. 다른 관절을 빠르게 움직이는 동안 한 관절은 제자리를 지키게 했다. 받침 쪽의 무거운 관절은 오차가 작았고, 손목 쪽의 가벼운 관절은 오차가 컸다.
교수는 이것을 트럭과 파리에 빗댄다. 달리는 트럭에 파리가 부딪혀도 트럭은 꿈쩍하지 않지만, 파리에게는 큰일이다. 관성이 큰 관절은 다른 링크가 미는 힘을 걸러 내고, 관성이 작은 관절은 그 힘에 흔들린다.
강의는 정찰 로봇 영상(2000)으로 시작한다. 본론은 4분 18초부터다.
영상에서 보기: 04:18 이제 동역학 · 05:20 두 가지 방법 · 10:32 자세마다 바뀌는 관성 · 14:08 트럭과 파리
11.2 운동 방정식의 모양
어떤 방법으로 구하든, 관절이 n개인 팔의 동역학은 늘 같은 모양으로 나온다. 강의는 이 식을 먼저 보여 주고 시작한다.
q, q̇, q̈, τ는 숫자 n개짜리 벡터다. 항을 하나씩 보자.
- τ: 관절마다 모터가 내는 힘이다. 회전 관절이면 토크(N·m), 직동 관절이면 힘(N)이다. 둘을 묶어 일반화 힘(generalized force)이라 하고, 강의는 Γ(감마)로도 쓴다.
- M(q) q̈: 가속도에 맞서는 관성의 몫이다. M은 n×n 질량 행렬(mass matrix)이다. 관절이 하나면 "질량 × 가속도"의 질량 자리다.
- V(q, q̇): 속도 때문에 생기는 힘이다. 원심력(centrifugal force)과 코리올리 힘(Coriolis force)이 여기 든다. 팔이 멈춰 있으면(q̇ = 0) 0이다.
- G(q): 중력을 버티는 몫이다. 자세 q에만 달렸다.
질량 행렬: 자세마다 바뀌고, 관절끼리 얽힌다
1.4절의 그림 1-5에서 팔꿈치를 접으면 어깨가 가벼워졌다. 그 "어깨가 느끼는 무거움"이 M의 첫 칸 M₁₁이다. 그래서 M은 자세 q에 따라 바뀐다.
M에는 대각선 밖의 칸도 있다. M₂₁은 "어깨만 가속할 때 팔꿈치가 버텨야 하는 몫"이다. 한 관절의 움직임이 다른 관절을 민다는 뜻이고, 이것을 얽힘(coupling)이라 한다(11.9절).
직동 관절 세 개가 x, y, z를 따로 미는 팔(문형 로봇)은 다르다. 관절 1은 팔 전체를, 관절 2는 그 뒤의 두 링크를, 관절 3은 마지막 링크만 민다. 그래서 M은 대각선에 질량의 합이 놓인 상수 행렬이다.
속도 항: 마이크 돌리기
교수는 마이크 줄을 잡고 빙빙 돌린다. 도는 동안 마이크는 바깥으로 당겨져 떠 있다가, 멈추면 떨어진다. 바깥으로 당기는 이 힘이 원심력이다.
링크가 여럿이면 두 관절 속도의 곱에서 나오는 코리올리 힘도 생긴다. 그래서 V의 모든 항은 q̇₁², q̇₁q̇₂처럼 속도 두 개의 곱이다.
11.10절에서 보듯 V는 M을 미분한 것에서 나온다. 그래서 M이 상수인 문형 로봇은 V = 0이다.
운동 에너지 행렬
질량 m인 점이 속도 v로 가면 운동 에너지는 ½mv²다. 팔 전체의 운동 에너지도 같은 모양이다.
그래서 M을 운동 에너지 행렬(kinetic energy matrix)이라고도 부른다.
M은 대칭 행렬이다. 움직이는 팔의 운동 에너지는 언제나 0보다 크므로, 어떤 q̇ ≠ 0을 넣어도 q̇ᵀMq̇ > 0이다. 이런 행렬을 양의 정부호(positive definite)라 한다.
수학 더하기: q̇ᵀ M q̇. 관절이 2개면 q̇ᵀMq̇ = M₁₁q̇₁² + 2M₁₂q̇₁q̇₂ + M₂₂q̇₂²로 풀린다. 행렬과 벡터를 곱했지만 결과는 숫자 하나다. M이 대칭(M₁₂ = M₂₁)이라 가운데 항이 두 번 나온다.
중력 항: 세워 든 링크
종이 팔은 책상 면을 따라 수평으로 움직인다. 관절 축이 위를 향하고 링크의 높이는 바뀌지 않는다. 중력은 일을 하지 않고 책상이 무게를 받치므로, 종이 팔은 G = 0이다.
중력을 보려면 다른 예가 필요해서, 이 장은 세워 든 링크를 쓴다. 길이 10 cm, 질량 100 g인 고른 막대를 아래 끝의 모터로 세로 평면에서 돌린다. 각도 θ는 수평에서 위로 잰다.
관절이 하나라 V = 0이고, M은 숫자 하나 I₀다(11.5절에서 구한다). 운동 방정식은 다음과 같다.
I₀ = 3.33 × 10⁻⁴ kg·m², m = 0.1 kg, g = 9.8 m/s², lc = 0.05 m(축에서 질량 중심까지)다. 중력 항은 무게 mg에 수평 거리 lc cos θ를 곱한 것이다. 링크가 수평이면 가장 커서 0.049 N·m이고, 곧게 서면 0이다.
해 볼 것 세 가지:
- θ = 0°(수평)에서 τ를 0.049로 맞춘다. θ̈ = 0이다. 수평으로 버티는 데 0.049 N·m가 든다.
- τ = 0으로 두고 θ를 바꿔 본다. 수평에서 θ̈ = −147 rad/s²로 가장 빨리 떨어지고, 곧게 선 90°에서는 0이다.
- θ를 90°보다 크게 넘긴다. cos θ < 0이라 중력 항의 부호가 바뀐다. 링크가 왼쪽으로 넘어가니 모터도 반대쪽으로 버텨야 한다.
알아 두기. 멈춰 있어도 중력 토크는 계속 든다. 서보가 링크를 수평으로 든 채 가만히 있으면 0.049 N·m(약 0.5 kg·cm)를 쉬지 않고 낸다. 무거운 팔을 오래 들고 있으면 서보가 뜨거워지는 까닭이다.
영상에서 보기: 14:40 운동 방정식 · 16:12 M은 질량 행렬 · 17:45 마이크 돌리기 · 20:21 운동 에너지 ½ q̇ᵀMq̇
11.3 점 하나: 운동량과 각운동량
동역학의 출발점은 뉴턴 법칙이다. 질량 m인 점에 힘 F를 주면 a = F/m으로 가속한다. 질량이 클수록 가속하기 어렵다.
강의는 같은 법칙을 조금 다르게 적는다. mv를 운동량(linear momentum)이라 하면, 뉴턴 법칙은 "운동량의 변화율 = 힘"이다.
이렇게 적는 까닭은 회전에도 똑같은 모양의 법칙이 있기 때문이다. 회전에서는 운동량 대신 각운동량을 쓴다.
모멘트를 잡으면 각운동량
움직이지 않는 점 O를 하나 잡고, O에서 점까지 화살표 p를 긋는다. 뉴턴 법칙 m v̇ = F의 양변 앞에 p ×를 붙인다. 오른쪽은 O 둘레의 모멘트(moment) N = p × F로, 8장의 토크와 같은 것이다.
왼쪽은 p × mv의 변화율과 같다. 곱의 미분으로 확인한다.
ṗ = v이고, 나란한 두 벡터의 외적은 0이다(6.6절). 그래서 첫 항이 사라진다.
수학 더하기: 곱의 미분. 두 양의 곱이 시간에 따라 바뀌면, 하나씩 번갈아 미분해서 더한다. (ab)′ = a′b + ab′다.
외적과 행렬 곱에도 그대로 통한다. 다만 외적은 순서를 지켜서 (a × b)′ = a′ × b + a × b′다.
p × mv를 각운동량(angular momentum)이라 하고 φ(파이)로 쓴다. 강의는 운동량에도 φ를 쓰지만, 이 장은 헷갈리지 않게 각운동량만 φ로 쓴다. 이제 두 법칙이 나란히 선다.
운동량의 변화율은 힘이고, 각운동량의 변화율은 모멘트다. 강체로 넘어가도 이 두 줄이 뼈대다.
영상에서 보기: 31:22 점 하나에 힘을 주면 · 32:28 운동량 · 34:03 각운동량 만들기 · 36:09 변화율은 모멘트
11.4 강체의 각운동량과 관성 텐서
강체는 점이 아주 많이 모인 것이다. 점 i의 질량을 mi, O에서의 위치를 pi라 하자. 강체가 O를 지나는 축 둘레로 ω로 돌면, 점마다의 속도는 6.6절의 vi = ω × pi다(그림 11-3 (나)).
각운동량을 점마다 구해 모두 더한다.
ω를 밖으로 꺼낸다
이 합에서 점마다 다른 것은 mi와 pi뿐이다. ω는 모든 점이 같다. ω를 합 밖으로 꺼낼 수 있으면, 남는 합은 강체의 모양과 질량만으로 정해진다.
요령은 둘이다. 먼저 순서를 바꾸면 부호가 바뀐다: ω × p = −p × ω. 그다음 6.7절의 외적 연산자로 p × 를 p̂로 쓴다.
이제 ω가 오른쪽 끝에 홀로 섰다. 합 밖으로 꺼낸다.
괄호 안의 3×3 행렬 I가 관성 텐서(inertia tensor)다. 운동량 = 질량 × 속도처럼, 각운동량 = 관성 텐서 × 각속도다.
점이 셀 수 없이 많으면 합 대신 적분을 쓴다. 작은 부피 dV의 질량은 밀도 ρ를 곱한 dm = ρ dV다.
수학 더하기: 적분은 아주 잘게 나눈 합. ∫는 "잘게 나누어 모두 더하기"다. 길이 10 cm, 100 g인 막대를 조각 1000개로 나눈다고 하자. 조각 하나는 0.1 g이고, 가운데에서 x만큼 떨어져 있다.
조각마다 0.1 × x²을 구해 모두 더하면 833.33 g·cm²다. 더 잘게 나누면 mL²/12 = 833.33…에 끝없이 다가간다. 이 끝값이 적분이다.
칸을 풀어 쓰기
−p̂p̂는 계산하기 쉬운 모양으로 바꿀 수 있다. p = (x, y, z)이면 pᵀp = x² + y² + z²는 숫자 하나이고, ppᵀ는 3×3 행렬이다. I₃는 3×3 단위행렬이다.
이것을 적분하면 관성 텐서의 칸이 나온다.
- 대각선의 Ixx, Iyy, Izz를 관성 모멘트(moment of inertia)라 한다. y² + z²은 그 조각에서 x축까지 거리의 제곱이다. 그래서 Ixx는 "x축 둘레로 돌리기 어려운 정도"로, 1.4절의 "거리의 제곱"이 여기서 나온다.
- 대각선 밖의 Ixy, Ixz, Iyz를 관성 곱(product of inertia)이라 한다. 두 축이 얽힌 정도다. 텐서에는 빼기 부호를 붙여 넣는다.
잠깐, 생각해 보기. 연필을 손가락 사이에서 길이 방향 축으로 굴리기는 쉽다. 가운데를 잡고 바람개비처럼 돌리기는 훨씬 힘들다. 관성 텐서의 어느 칸이 이 차이를 말해 주나?
대각선 칸이다. 연필이 x 방향으로 길면 조각들이 x축에 바짝 붙어 있어 Ixx가 작다. y축이나 z축에서는 멀리 떨어진 조각이 많아 Iyy, Izz가 크다.
그림 11-4 (나)의 막대는 16.7 g·cm²와 842 g·cm²로 약 50배 차이 난다.
성질 두 가지
- 대칭이다. i행 j열과 j행 i열이 같다. 그래서 독립인 숫자는 9개가 아니라 6개다.
- 양의 정부호다. 회전 운동 에너지는 ½ ωᵀ I ω이고, 질량이 있는 물체가 돌면 이 값은 0보다 크다. 교수는 0이 나오는 경우는 블랙홀뿐이라고 농담한다.
상자와 막대
x, y, z 방향으로 a, b, c인 고른 상자를 생각하자. 질량 중심에 원점을 두고 축을 모서리와 나란히 잡으면 적분이 쉽게 끝난다.
관성 곱이 모두 0인 까닭은 대칭이다. (x, y, z)에 조각이 있으면 (−x, y, z)에도 같은 조각이 있다. 두 조각의 xy는 크기가 같고 부호가 반대라서 더하면 지워진다.
가는 막대는 굵기 b, c가 0에 가까운 상자다. 길이 방향 축 둘레는 거의 0이고, 가로 축 둘레는 mL²/12다. 세워 든 링크(10 cm, 100 g)라면 이렇다.
알아 두기. 단위를 조심한다. 1 g·cm² = 10⁻³ kg × 10⁻⁴ m² = 10⁻⁷ kg·m²다. 작은 팔은 g·cm²로 적어야 숫자를 읽기 좋지만, 토크(N·m)를 계산할 때는 kg·m²로 바꿔 넣는다.
영상에서 보기: 37:13 점을 모두 더한다 · 43:02 φ = Iω · 47:17 −p̂p̂ 풀어 쓰기 · 51:00 관성 모멘트와 관성 곱
11.5 평행축 정리: 다른 점에서 본 관성
관성 텐서는 어느 점에서 쟀는지에 따라 다르다. 같은 막대도 가운데 축보다 끝 축으로 돌리기가 더 어렵다. 그런데 관절 축은 보통 링크의 끝에 있다.
매번 새로 적분할 필요는 없다. 질량 중심 C에서 한 번 구해 두면, 다른 점의 관성은 더하기 한 번으로 나온다.
축 하나라면: I = Ic + m d²
질량 중심을 지나는 축 둘레의 관성이 Ic다. 그 축과 나란하고 d만큼 떨어진 축 둘레의 관성은 다음과 같다. 이것이 평행축 정리(parallel axis theorem)다.
막대로 까닭을 보자. 질량 중심에서 x만큼 떨어진 조각은 새 축에서 x + d만큼 떨어져 있다. 거리의 제곱을 질량으로 더하면 세 덩어리가 나온다.
가운데 덩어리가 0인 까닭은 C가 질량 중심이기 때문이다. 질량 중심은 막대가 균형을 잡는 점이라, 양쪽의 x dm이 서로 지워진다.
해 볼 것 세 가지:
- d = 0: 질량 중심 축이다. I = Ic = 833 g·cm²다.
- d = 5 cm: 막대 끝, 곧 세워 든 링크의 모터 축이다. I₀ = 833 + 2500 = 3333 g·cm² = mL²/3이다. 11.2절의 I₀ = 3.33 × 10⁻⁴ kg·m²가 이것이다.
- d = 14 cm: 축이 막대에서 멀리 떨어진다. m d²가 거의 전부다. 멀리 있는 물체는 점 하나로 봐도 된다는 뜻이다.
잠깐, 생각해 보기. 그림 1-5는 링크를 질량 중심에 모인 점으로 보았다. 무엇을 빠뜨렸고, 얼마나 큰가?
링크마다의 Ic를 빠뜨렸다. 링크 1(10 cm, 20 g)은 167 g·cm², 링크 2(8 cm, 15 g)는 80 g·cm²다.
팔을 다 편 자세에서 어깨가 느끼는 관성은 19640에서 19887 g·cm²로 1.3% 커질 뿐이다. 손에 든 50 g이 어깨에서 18 cm나 떨어져 있어 전체의 80%를 넘게 차지하기 때문이다.
3차원: 관성 곱도 함께
3차원에서는 A에서 C까지의 화살표 pc = (xc, yc, zc) 전체를 쓴다. 두 좌표계의 축은 나란하다고 둔다.
칸으로 풀면 Ic + m d²와 같은 모양이다. 관성 곱에는 좌표 두 개의 곱이 더해진다.
강의 예제: 정육면체의 꼭짓점
한 변 a, 질량 m인 고른 정육면체다. 질량 중심에서는 상자 공식에 a를 넣어 세 축 모두 m(a² + a²)/12 = ma²/6이고, 관성 곱은 0이다.
이제 원점을 한 꼭짓점 A로 옮기고, 축은 세 모서리를 따라 잡는다. 그러면 pc = (a/2, a/2, a/2)다.
- 대각선: ma²/6 + m(a²/4 + a²/4) = ma²/6 + ma²/2 = 2ma²/3
- 관성 곱: 0 + m(a/2)(a/2) = ma²/4. 텐서에는 −ma²/4로 들어간다.
질량 중심에서는 대각 행렬이던 것이 꼭짓점에서는 관성 곱을 갖는다. 꼭짓점에서 보면 모든 조각이 x > 0, y > 0 쪽에 있어서 xy가 지워지지 않는다.
영상에서 보기: 51:31 대칭이면 관성 곱이 0 · 52:33 평행축 정리 · 53:37 정육면체 예제 · 54:40 다른 점에서는 관성 곱이 생긴다
11.6 좌표계를 돌리면: R I Rᵀ와 주축
평행축 정리는 원점만 옮겼다. 축의 방향이 다르면 어떻게 될까?
링크의 관성은 보통 링크에 붙은 좌표계로 한 번 적어 둔다. 그런데 링크가 돌면, 받침에서 본 관성 텐서의 숫자는 바뀐다.
각운동량 식에서 바로 나온다
각운동량과 각속도는 둘 다 화살표다. {C}의 숫자를 {A}의 숫자로 바꾸려면 회전행렬 R = ACR을 곱한다(2장).
마지막 단계에서 Cω = Rᵀ Aω를 썼다(2.7절). 그러니 {A}의 숫자로 쓴 관성 텐서는 다음과 같다.
원점도 다르면 먼저 돌리고(R I Rᵀ), 그다음 평행축 정리로 옮긴다. pc는 A에서 C까지의 화살표를 {A}의 숫자로 적은 것이다.
판 하나로 보기
가로 6 cm, 세로 2 cm, 질량 12 g인 얇은 판이 책상 위에 누워 있다. 판에 붙인 좌표계 {C}의 x′축은 긴 쪽, y′축은 짧은 쪽이다. 상자 공식에 두께 0을 넣는다.
판을 책상 위에서 α만큼 돌린다. R은 z축 회전 RZ(α)다. 슬라이더로 α를 바꿔 R I Rᵀ의 숫자를 보자.
해 볼 것 네 가지:
- α = 0°: 대각 행렬 그대로다. φ = (4, 0, 0)으로 ω와 나란하다.
- α = 30°: Ixx = 12, Iyy = 28이고 관성 곱 Ixy = 13.9가 생긴다. φ = (12, −13.9, 0)은 ω에서 49° 비켜 있다.
- α = 90°: 다시 대각 행렬이다. 긴 쪽이 y축을 따라 놓여 Ixx = 36이다.
- φ와 ω가 나란해지는 α를 찾아본다. −90°, 0°, 90°뿐이다. ω가 판에 붙은 축 x′나 y′와 겹칠 때다.
어느 α에서든 Izz = 40이다. 판을 책상 위에서 돌려도 책상에 수직인 축 둘레 관성은 그대로다.
주축
관성 곱이 모두 0이 되는 세 축을 주축(principal axes)이라 한다. 주축으로 잰 관성 텐서는 대각 행렬이고, 대각 칸을 주관성 모멘트(principal moments of inertia)라 한다. 판에서는 x′, y′, z가 주축이다.
어떤 모양의 강체든 서로 수직인 주축 세 개가 반드시 있다. 모양이 대칭이면 대칭축이 주축이다. 그래서 링크의 관성은 주축에 맞춘 좌표계로 대각 행렬 하나를 적어 두고, 필요할 때 R I Rᵀ로 돌린다.
수학 더하기: 고유벡터. 행렬 A에 벡터 u를 곱했더니 방향은 그대로이고 길이만 λ배가 되면(Au = λu), u를 A의 고유벡터(eigenvector), λ를 고유값(eigenvalue)이라 한다. 주축은 관성 텐서의 고유벡터이고, 주관성 모멘트는 고유값이다. 그림 11-6에서 φ = Iω가 ω와 나란해지는 ω가 바로 고유벡터다.
영상에서 보기: 46:15 링크마다 I를 구할 줄 알아야 한다
11.7 강체 하나의 두 방정식: 뉴턴과 오일러
강체는 옮겨 가면서 돌기도 한다. 이 움직임을 두 식으로 나누어 적는다. 옮겨 가는 쪽은 뉴턴 방정식(Newton's equation), 도는 쪽은 오일러 방정식(Euler's equation)이다.
- F는 강체에 걸린 바깥 힘의 합, v̇c는 질량 중심의 가속도다.
- N은 질량 중심 둘레의 모멘트 합, Ic는 질량 중심에서 잰 관성 텐서다.
뉴턴 방정식은 질량 중심에
강체를 이루는 점마다 F = ma가 성립한다. 모두 더하면 점끼리 밀고 당기는 힘은 짝으로 지워진다(작용·반작용). 바깥 힘의 합만 남고, 오른쪽 Σ miai는 m v̇c가 된다.
그래서 강체는 질량 중심 한 점에 질량이 모두 모인 것처럼 옮겨 간다. "짝으로 지워진다"는 생각은 11.8절에서 링크 사이의 힘에 다시 쓴다.
오일러 방정식: 회전의 F = ma
11.3절에서 각운동량의 변화율은 모멘트였다. 각운동량은 φ = Iω이니 미분하면 된다. 그런데 미분하면 항이 둘 나온다.
- Ic ω̇: 각가속도의 몫이다. 질량 × 가속도와 같은 모양이다.
- ω × Icω: 각속도 두 개의 곱이다. 각가속도가 0이어도 생긴다. 회전에서 생기는 원심력·코리올리 항이고, 11.2절의 V로 이어진다.
둘째 항이 생기는 까닭은 I가 강체에 붙어 함께 돌기 때문이다. 그림 11-6에서 φ는 ω에서 비켜 있었다.
강체가 돌면 비켜 있던 φ도 같이 돌아서, ω가 일정해도 φ의 방향은 계속 바뀐다. 방향을 바꾸는 데 모멘트가 든다.
수학 더하기: 왜 ω × Iω인가. 받침에서 본 관성은 I₀ = R Ic Rᵀ다(11.6절). R의 세 열은 강체에 붙은 축이고, 축 끝은 ω × (축)의 속도로 움직인다(6.6절). 그래서 Ṙ = ω̂R이다.
곱의 미분으로 I₀를 미분하면 İ₀ = ω̂I₀ − I₀ω̂다(ω̂ᵀ = −ω̂). 그러면 d(I₀ω)/dt = I₀ω̇ + ω̂I₀ω − I₀ω̂ω이고, ω̂ω = ω × ω = 0이다. I₀ω̇ + ω × I₀ω가 남는다.
예: 판을 기울여 돌리면
그림 11-6의 판을 α = 30°로 두고, 책상의 x축 둘레로 ω = (10, 0, 0) rad/s로 일정하게 돌린다. ω̇ = 0이라 첫 항은 0이다. 둘째 항을 계산한다.
z 성분만 10 × (−138.6) = −1386으로 남는다. 단위는 g·cm²/s²이고, 10⁻⁷을 곱하면 −1.39 × 10⁻⁴ N·m다.
각속도는 그대로인데도 축이 z 방향 모멘트를 계속 내야 한다. 판의 주축이 회전축과 어긋나 있기 때문이다.
알아 두기. 자동차 바퀴의 동적 밸런스가 이 이야기다. 바퀴의 주축이 차축과 조금만 어긋나도 ω × Iω가 생기고, 이 모멘트가 바퀴와 함께 돌며 차축을 흔든다. 그래서 정비소에서 작은 추를 붙여 주축을 차축에 맞춘다.
평면 팔에서는 간단해진다
종이 팔처럼 평면에서 도는 링크는 ω가 늘 책상에 수직(z)이다. 링크가 책상 면을 기준으로 대칭이면 z가 주축이라 Iω = Izzω로 ω와 나란하다. 그러면 ω × Iω = 0이고 오일러 방정식은 N = Izz ω̇ 하나로 줄어든다(세워 든 링크도 같다).
영상에서 보기: 43:34 뉴턴: 운동량의 변화율 = 힘 · 44:06 오일러: 각운동량의 변화율 = 모멘트 · 44:40 Iω̇ 말고 한 항이 더 · 45:12 강체 하나의 두 식
11.8 링크 하나의 자유물체도와 내부 힘
이제 강체를 관절로 잇는다. 8장에서는 멈춰 있는 팔의 관절을 모두 끊고 링크마다 힘의 평형을 따졌다. 동역학도 똑같이 하되, 링크가 가속하므로 힘의 합이 0이 아니라 관성력과 같아야 한다.
링크 i에 걸리는 힘
링크 i는 관절 i에서 링크 i−1과, 관절 i+1에서 링크 i+1과 닿는다. 다음처럼 이름을 붙인다.
- fi, ni: 링크 i−1이 관절 i에서 링크 i에 주는 힘과 모멘트다.
- −fi+1, −ni+1: 링크 i+1이 관절 i+1에서 링크 i에 되돌려 주는 힘과 모멘트다. 작용·반작용이라 부호만 반대다.
이렇게 몸 하나만 떼어 놓고 걸리는 힘을 모두 그린 그림을 자유물체도(free-body diagram)라 한다. 그림의 무게 mig는 11.9절의 요령으로 따로 넣는다. 그래서 아래 식에는 관절의 힘만 남는다.
pc i는 관절 i에서 링크 i의 질량 중심까지, pi+1은 관절 i에서 관절 i+1까지의 화살표다. 11.7절의 Fi = miv̇c i, Ni = Ic iω̇i + ωi × Ic iωi를 쓰면 링크 i의 두 방정식은 다음과 같다.
모멘트 식은 질량 중심 둘레로 잡았다. 관절 i는 질량 중심에서 −pc i, 관절 i+1은 pi+1 − pc i만큼 떨어져 있다. 그래서 두 힘이 각각 그 화살표와의 외적만큼 모멘트를 보탠다.
왜 내부 힘을 없애야 하나
fi, ni는 3 + 3 = 6개의 숫자인데, 회전 관절의 모터가 내는 것은 그중 하나인 관절 축 둘레의 토크뿐이다. 나머지 다섯은 관절의 베어링과 링크의 뼈대가 버틴다. 관절이 그 방향의 움직임을 막고 있기 때문이다(2.1절의 "관절은 5개를 막는다").
이 다섯 개는 모터가 고를 수 없고, 움직임도 정하지 않는다. 링크를 붙잡아 사슬 모양을 지키는 구속력(constraint force)일 뿐이다. 그래서 운동 방정식에 τ와 q만 남기려면 이것들을 지워야 한다.
지우는 방법은 8장과 같다. ni를 관절 축 Zi 위에 내린다(내적).
직동 관절이면 모터가 축 방향으로 미는 힘을 내므로, ni 대신 fi를 내린다: τi = fiᵀZi.
잠깐, 생각해 보기. 세워 든 링크를 수평으로 든 채 멈춰 있다. 관절이 링크에 주는 힘 f₁은 얼마인가? 모터가 이 힘을 내나?
f₁은 위로 0.98 N으로, 링크의 무게와 같다. 이 힘은 모터가 아니라 베어링이 받치고, 모터는 축 둘레의 토크 0.049 N·m만 낸다. 무게는 받침으로, 토크만 모터로 간다.
알아 두기. 지운다고 쓸모없는 값은 아니다. 뉴턴-오일러 방법은 계산하는 동안 fi, ni를 모두 구해 준다. 베어링과 링크가 얼마나 튼튼해야 하는지 정할 때 이 값을 쓴다.
영상에서 보기: 21:57 정적 평형 다시 보기 · 22:30 움직이면 관성력이 더해진다 · 24:34 축에 내려 내부 힘을 지운다 · 45:43 관절은 구속이다
11.9 뉴턴-오일러 순환 알고리즘
링크마다 식을 세웠으니 차례로 풀면 된다. 순서는 사슬을 두 번 지나가는 것이다. 이것을 뉴턴-오일러 순환 알고리즘(recursive Newton–Euler algorithm)이라 한다.
① 바깥으로: 속도와 가속도
6.8절의 속도 전파를 그대로 쓰고, 한 번 더 미분해 가속도도 넘긴다. 회전 관절만 적는다. 모든 화살표는 받침 {0}의 숫자이고, Zi+1은 관절 i+1의 축이다.
v̇i+1은 관절 i+1의 가속도다. ω × (ω × p)는 원을 따라 돌 때 중심 쪽을 향하는 구심 가속도(centripetal acceleration)다. 마지막으로 질량 중심의 가속도를 구해 관성력을 계산한다.
받침 {0}의 숫자로 계산하므로 Ic도 11.6절의 R I Rᵀ로 돌려 넣는다. 직동 관절이면 ω는 그대로이고, 가속도에 ḋ, d̈가 든 항이 더해진다(강의 유인물).
수학 더하기: 가속도는 속도를 한 번 더 미분한 것. ωi+1 식을 시간으로 미분하면 ω̇i+1 식이 된다. 축 Zi+1은 링크 i에 붙어 ωi로 돌기 때문에 Żi+1 = ωi × Zi+1이다. 그래서 곱의 미분으로 θ̇Z를 미분하면 θ̈Z와 θ̇ ωi × Z가 나온다.
② 안으로: 힘과 모멘트
11.8절의 두 식을 fi, ni에 대해 정리하면, 다음 링크의 값으로 이번 링크의 값을 구하는 식이 된다. 끝에서 시작해 받침 쪽으로 간다.
출발점은 손끝이다. 손이 아무것도 잡지 않았으면 fn+1 = 0, nn+1 = 0이다. 무언가를 밀고 있으면 그 반작용을 넣는다.
강의 유인물은 6.8절처럼 링크마다 자기 좌표계의 숫자로 적고 회전행렬로 넘긴다. 숫자를 적는 방법만 다르고 내용은 같다.
중력은 받침을 위로 가속해서 넣는다
중력을 링크마다 따로 더하지 않는 요령이 있다. 받침이 위로 g만큼 가속하고 있다고 놓고 시작한다. z가 위라면 v̇₀ = (0, 0, g)다.
엘리베이터가 위로 가속하면 몸이 무거워지는 것과 같은 이치다. 모든 링크가 받침의 가속도를 물려받으므로, 링크마다 무게 mg가 저절로 F에 들어간다. 강의에서 학생이 바로 이 답을 말한다.
손으로 풀기: 세워 든 링크
관절이 하나라 n = 1이다. x는 수평, y는 위, z는 종이 밖이다. 링크의 방향은 u = (cos θ, sin θ, 0), 질량 중심은 pc = lc u다.
① 바깥으로 간다. 중력 때문에 v̇₀ = (0, g, 0)으로 둔다. 관절 1은 받침 원점에 있어 v̇₁ = v̇₀다.
- ω₁ = θ̇ ẑ, ω̇₁ = θ̈ ẑ
- 질량 중심의 가속도는 아래 식이다. ω₁ × (ω₁ × pc) = −θ̇² pc를 썼다.
- F₁ = m v̇c, N₁ = Ic θ̈ ẑ (평면이라 ω × Iω = 0)
② 안으로 온다. 끝에는 아무것도 없어 f₂ = 0, n₂ = 0이다. 그래서 f₁ = F₁, n₁ = N₁ + pc × F₁이다.
pc × F₁의 z 성분을 계산하면 θ̇²가 든 항은 서로 지워진다. 구심 가속도는 pc와 나란해서 모멘트를 만들지 못하기 때문이다.
11.2절의 식이 나왔다. Ic + m lc²에 평행축 정리가 저절로 들어 있는 것도 눈여겨보자.
종이 팔에 돌려 보기
종이 팔(링크 1 = 10 cm, 20 g, 링크 2 = 8 cm, 15 g, 손에 든 물건 50 g)에 같은 알고리즘을 돌린다. 수평으로 움직이니 중력은 넣지 않는다(v̇₀ = 0). 링크 2와 물건은 한 덩어리로 묶는다. 질량 65 g, 질량 중심은 팔꿈치에서 7.08 cm이고, 그 둘레 관성은 평행축 정리로 265 g·cm²다.
멈춘 팔(q̇ = 0)에 가속도 q̈ = (1, 0)을 넣으면 τ = M q̈, 곧 M의 첫 열이 나온다. q̈ = (0, 1)이면 둘째 열이다. 뉴턴-오일러로 질량 행렬을 한 열씩 채울 수 있다.
팔을 쭉 편 θ₂ = 0°에서는 이렇게 나온다(파이썬 확인 코드).
어깨를 0.25초 동안 90° 돌린다고 하자. 앞 절반은 가속, 뒤 절반은 감속한다. 그러면 1.57 rad = θ̈ × 0.25²/4에서 θ̈₁ ≈ 100 rad/s²가 든다. 팔꿈치 각은 그대로 둔다(θ̈₂ = 0).
해 볼 것 세 가지:
- θ₂ = 0°: τ₁ = 0.199 N·m(약 2 kg·cm)로, 작은 서보의 최대 토크와 비슷하다. 팔꿈치는 각도가 바뀌지 않는데도 τ₂ = 0.081 N·m를 내야 한다.
- θ₂를 150°까지 접는다. τ₁은 0.027 N·m로 약 7분의 1이 된다. 그림 1-5와 같은 이야기다.
- θ₂ = 90°에서는 M₂₁ = M₂₂다. 팔꿈치 모터를 끄면 θ̈₂ = −θ̈₁이 되어, 링크 2는 방향을 지킨 채 옮겨 가기만 한다.
잠깐, 생각해 보기. 팔을 쭉 편 채 팔꿈치 모터를 끄고 어깨를 빨리 돌리면, 팔꿈치는 어느 쪽으로 꺾일까?
θ₂ = 0°에서 θ̈₂ = −(8120/3520) × 100 ≈ −231 rad/s²로, 뒤로 처진다. 링크 2와 물건이 제자리에 남으려 해서 어깨가 도는 반대쪽으로 꺾인다. 팔을 휘두르면 힘 뺀 손목이 늦게 따라오는 것과 같다.
계산량과 한계
관절마다 바깥으로 한 번, 안으로 한 번 지나가니 계산량은 관절 수에 비례한다. 강의는 6관절 팔이면 바깥으로 0~5, 안으로 6~1을 돌고, 휴머노이드 ASIMO는 25까지 간다고 말한다. 식은 복잡해 보여도 컴퓨터에게는 빠르다.
대신 결과가 숫자로만 나온다. 교수도 슬라이드의 식을 따라가다 "길을 잃었다"며 웃는다. 어느 항이 원심력이고 어느 항이 얽힘인지 한눈에 보려면 라그랑주 방법이 낫다.
영상에서 보기: 29:13 속도는 바깥으로, 힘은 안으로 · 55:42 속도·가속도 순환식 · 58:20 되넘기며 축에 내리기 · 59:54 중력은 어떻게 넣나
11.10 라그랑주 방정식 첫걸음
뉴턴-오일러 방법은 링크마다 힘을 따지고 내부 힘을 지웠다. 라그랑주 방법은 처음부터 내부 힘을 보지 않고, 팔 전체의 에너지 두 가지만 쓴다.
- 운동 에너지(kinetic energy) K: 움직임의 에너지다. 링크 i의 몫은 ½ mivc iᵀvc i + ½ ωiᵀIc iωi이고, 모두 더한 것이 K다. 회전 몫에 관성 텐서가 다시 나온다.
- 위치 에너지(potential energy) U: 높이의 에너지다. 링크 i의 몫은 mighi(hi는 질량 중심의 높이)이고, 모두 더한 것이 U다. 자세 q에만 달렸다.
둘의 차이를 라그랑지안(Lagrangian) L이라 한다.
라그랑주 방정식(Lagrange equations)은 관절마다 하나씩, 다음과 같이 생겼다.
q와 q̇를 서로 다른 변수로 보고 L을 편미분한다(6.2절). 그다음 첫 항은 d/dt로 시간에 대해 한 번 더 미분한다.
세워 든 링크로 해 보기
강의처럼 관절 하나로 먼저 해 본다. 세워 든 링크는 축에 고정되어 돌기만 한다.
운동 에너지는 ½ m vc² + ½ Icθ̇² = ½(m lc² + Ic)θ̇²라서 축 둘레 관성 I₀ 하나로 적힌다. 질량 중심의 높이는 lc sin θ다.
한 조각씩 계산한다.
- ∂L/∂θ̇ = I₀θ̇. θ̇²를 θ̇로 미분하면 2θ̇이고, ½과 만나 θ̇가 된다. ½mv²를 v로 미분하면 운동량 mv가 나오는 것과 같다.
- d/dt(I₀θ̇) = I₀θ̈. I₀는 상수다.
- ∂L/∂θ = −∂U/∂θ = −m g lc cos θ. K에는 θ가 없다.
넣으면 11.2절, 11.9절과 같은 식이 나온다.
내부 힘 f₁은 한 번도 나오지 않았다. 뉴턴-오일러보다 훨씬 짧다.
잠깐, 생각해 보기. 중력 항 m g lc cos θ는 위치 에너지 U를 θ로 미분한 것이다. 링크가 곧게 선 θ = 90°에서 중력 항이 0인 까닭을 U로 설명해 보자.
θ = 90°에서 높이 lc sin θ가 가장 크다. 꼭대기에서는 θ를 조금 바꿔도 높이가 거의 그대로라 기울기 ∂U/∂θ가 0이다. 그래서 모터가 버틸 토크도 0이다.
관절이 여럿이면: M, V, G가 보인다
이번에는 벡터로 한꺼번에 쓰고 K = ½ q̇ᵀM(q)q̇를 넣는다. U는 q̇에 달리지 않으므로 K와 U를 나누어 적는다.
세 조각을 하나씩 본다.
- ∂K/∂q̇ = M q̇. 숫자 하나일 때 ½mv²를 v로 미분하면 mv인 것과 같다. 관절이 2개면 K = ½(M₁₁q̇₁² + 2M₁₂q̇₁q̇₂ + M₂₂q̇₂²)다. q̇₁로 미분하면 M₁₁q̇₁ + M₁₂q̇₂로, M q̇의 첫 줄이다.
- d/dt(M q̇) = M q̈ + Ṁ q̇. 곱의 미분이다. 자세가 바뀌면 M도 바뀌므로 Ṁ이 생긴다.
- ∂K/∂q. K에서 q에 달린 것은 M뿐이다. i번째 칸은 ½ q̇ᵀ(∂M/∂qi)q̇다.
모아서 M q̈ + V + G = τ의 모양과 맞춰 보면 두 항의 정체가 드러난다.
- V의 모든 항은 q̇ 두 개의 곱이다(Ṁ = Σ (∂M/∂qj) q̇j도 q̇를 하나 품는다). 멈추거나(q̇ = 0) M이 상수이면 V는 사라진다. 11.2절에서 말한 성질이 여기서 나온다.
- G는 위치 에너지의 기울기다. 팔이 아래로 내려가려는 몫이고, 모터는 그만큼 버텨야 한다.
결국 동역학을 얻는 데 필요한 것은 둘뿐이다. 운동 에너지 속의 M과 위치 에너지 U다. M을 알면 V는 미분으로 나온다.
다음 장: M을 직접 구하기
K = ½ q̇ᵀMq̇의 M은 어떻게 구할까? 링크마다 운동 에너지 ½ m vcᵀvc + ½ ωᵀIcω를 구해 모두 더한다. 그 합을 q̇ᵀ(…)q̇ 모양으로 정리하면 괄호 안이 M이다.
vc와 ω는 자코비안으로 q̇에 이어진다. 1장에서 교수가 말한 대로 동역학도 자코비안으로 만든다. 이 계산(명시적 형태)과 V를 칸마다 적는 법, 전체 운동 방정식은 12장에서 다룬다.
알아 두기. 뉴턴-오일러와 라그랑주는 같은 답을 준다. 세워 든 링크에서 두 방법이 같은 식을 냈고, 종이 팔에서도 V가 맞아떨어진다(파이썬 확인 코드). 컴퓨터로 빨리 계산할 때는 뉴턴-오일러, 식의 구조를 보고 제어를 설계할 때는 라그랑주를 주로 쓴다.
이것으로 11강 내용이 끝난다.
영상에서 보기: 1:01:59 회전 운동 에너지 ½ ωᵀIω · 1:04:38 L = K − U · 1:07:45 ∂K/∂q̇ = Mq̇ · 1:12:01 운동 방정식의 모양
11.11 연습문제
중력 가속도는 g = 9.8 m/s²로 둔다. 답은 문제 8개 아래에 있다. 답을 보기 전에 손으로 끝까지 풀어 본다.
- x, y, z 방향으로 4 cm, 2 cm, 2 cm인 60 g 상자의 질량 중심 관성 텐서를 구하라. 세 축 가운데 어느 축으로 돌리기가 가장 쉬운가?
- 종이 팔의 링크 2(8 cm, 15 g, 가는 막대)를 생각하자. 축은 모두 책상에 수직이다.
- 질량 중심을 지나는 축 둘레의 관성 Ic
- 팔꿈치 축 둘레의 관성
- 팔을 쭉 폈을 때(θ₂ = 0°) 어깨 축 둘레의 관성
- 한 변 2 cm, 질량 24 g인 정육면체(강의 예제)의 관성 텐서를 질량 중심과 한 꼭짓점에서 구하라. 꼭짓점에서는 축을 세 모서리를 따라 잡는다.
- 11.6절의 판(6 cm × 2 cm, 12 g)을 α = 30° 돌렸다. R I Rᵀ를 계산하라. ω = (1, 0, 0) rad/s일 때 각운동량 φ를 구하라.
- 문제 4의 판을 ω = (10, 0, 0) rad/s로 일정하게 돌린다. 필요한 모멘트 N을 구하라. α = 0°이면 어떻게 되나?
- 세워 든 링크(10 cm, 100 g)가 θ = 60°에 있다.
- 멈춰 버틸 토크는 얼마인가?
- 모터가 0.05 N·m를 내면 θ̈는?
- τ = 0으로 놓으면 θ̈는? 어느 쪽으로 도나?
- 세워 든 링크가 수평(θ = 0°)을 θ̇ = 2 rad/s로 지나는 순간 θ̈ = 0이다. 뉴턴-오일러로 v̇c, 관절이 주는 힘 f₁, 토크 τ를 구하라. 중력은 받침을 위로 가속해서 넣는다.
- 질량 0.2 kg인 상자를 위아래로 미는 직동 관절 하나가 있다. 관절 변수는 높이 d다. 라그랑주 방정식으로 관절 힘 f의 식을 구하고, 위로 2 m/s²로 가속하는 데 드는 f를 구하라.
답
답 1. m/12 = 5 g이다. Ixx = 5(2² + 2²) = 40, Iyy = 5(4² + 2²) = 100, Izz = 5(4² + 2²) = 100 g·cm²다. 관성 곱은 대칭이라 0이다.
긴 쪽인 x축 둘레가 가장 쉽다. 조각들이 x축에 가깝기 때문이다.
답 2.
- Ic = mL²/12 = 15 × 64/12 = 80 g·cm²
- 팔꿈치는 질량 중심에서 4 cm다. 80 + 15 × 4² = 320 g·cm²로, mL²/3과 같다.
- 팔을 펴면 질량 중심은 어깨에서 10 + 4 = 14 cm다. 80 + 15 × 14² = 80 + 2940 = 3020 g·cm²다. 그림 1-5는 80을 빼고 2940으로 보았다.
답 3. 질량 중심에서는 ma²/6 = 24 × 4/6 = 16 g·cm²인 대각 행렬이다. 꼭짓점에서는 pc = (1, 1, 1) cm다. 대각선은 16 + 24 × 2 = 64(= 2ma²/3), 관성 곱은 24 × 1 × 1 = 24(= ma²/4)다.
답 4. c = cos 30° ≈ 0.866, s = sin 30° = 0.5로 R diag(4, 36) Rᵀ를 칸마다 쓴다.
Ixx = 4c² + 36s² = 3 + 9 = 12, Iyy = 4s² + 36c² = 1 + 27 = 28이다. 대각선 밖은 (4 − 36)cs = −32 × 0.433 = −13.9다.
관성 곱은 Ixy = 13.9다. φ = Iω는 첫 열 그대로 (12, −13.9, 0) g·cm²/s이고, ω와 49° 어긋난다.
답 5. ω̇ = 0이라 N = ω × Iω다. Iω = (120, −138.6, 0)이고, ω × Iω = (0, 0, 10 × (−138.6)) = (0, 0, −1386) g·cm²/s², 약 −1.39 × 10⁻⁴ N·m다.
α = 0°이면 Iω = (40, 0, 0)으로 ω와 나란해서 N = 0이다. 주축 둘레로 돌리면 모멘트가 필요 없다.
답 6. m g lc = 0.049 N·m, cos 60° = 0.5, I₀ = 3.33 × 10⁻⁴ kg·m²다.
- 버틸 토크: 0.049 × 0.5 = 0.0245 N·m
- θ̈ = (0.05 − 0.0245)/(3.33 × 10⁻⁴) = 76.5 rad/s². 위로 가속한다.
- θ̈ = −0.0245/(3.33 × 10⁻⁴) = −73.5 rad/s². θ가 줄어드는 쪽, 곧 수평 쪽으로 넘어진다.
답 7. θ = 0이라 pc = (0.05, 0, 0) m, θ̈ = 0이다.
- v̇c = (0, 9.8, 0) + 0 − 2² × (0.05, 0, 0) = (−0.2, 9.8, 0) m/s²
- f₁ = F₁ = 0.1 × v̇c = (−0.02, 0.98, 0) N. 무게를 받치는 0.98 N과, 링크를 축 쪽으로 당겨 원을 돌게 하는 0.02 N이다.
- τ = (pc × F₁)z = 0.05 × 0.98 − 0 × (−0.02) = 0.049 N·m. 구심 몫은 축을 지나서 토크에 보태지 않는다.
답 8. K = ½ m ḋ², U = m g d다. ∂L/∂ḋ = m ḋ이고, 시간으로 미분하면 m d̈다. ∂L/∂d = −m g다.
M = m(상수), V = 0, G = mg다. 위로 2 m/s²라면 f = 0.2 × (2 + 9.8) = 2.36 N이다.
파이썬으로 확인
파이썬(numpy)이 있으면 아래를 붙여 넣어 돌려 본다. 주석의 값이 그대로 나온다. -0.은 아주 작은 계산 오차로 생긴 0이다.
newton_euler는 11.9절의 두 번 지나가기를 평면 팔에 맞게 줄인 것이다. 관절 축이 모두 z라서 각속도는 숫자 하나이고, ω × Iω = 0이다. mass_matrix는 멈춘 팔에 단위 가속도를 하나씩 넣어 M을 한 열씩 채운다.
import numpy as np
np.set_printoptions(precision=4, suppress=True)
D = np.pi / 180 # 도 → 라디안
g = 9.8 # m/s²
def box_inertia(m, a, b, c):
"""고른 상자(x 방향 a, y 방향 b, z 방향 c)의 질량 중심 관성 텐서."""
return m / 12 * np.diag([b * b + c * c, a * a + c * c, a * a + b * b])
def parallel_axis(I_c, m, p):
"""평행축 정리: 질량 중심에서 p만큼 옮긴 점의 관성 텐서 (축은 나란하다)."""
p = np.asarray(p, float)
return I_c + m * ((p @ p) * np.eye(3) - np.outer(p, p))
def rot_z(deg):
t = np.radians(deg)
return np.array([[np.cos(t), -np.sin(t), 0], [np.sin(t), np.cos(t), 0], [0, 0, 1]])
def inertia_by_pieces(P, m):
"""조각마다 (pᵀp)I − ppᵀ 에 조각 질량을 곱해 더한다 (11.4절의 적분을 합으로)."""
dm = m / len(P)
return dm * (np.sum(P ** 2) * np.eye(3) - P.T @ P)
# 11.4 막대(10 cm, 100 g)를 조각 1000개로 나눠 더하기
x = (np.arange(1000) + 0.5) / 100 - 5 # 조각 가운데의 위치 (cm)
print(round(np.sum(0.1 * x ** 2), 4)) # 833.3325
print(box_inertia(100, 10, 1, 1)) # 대각 16.6667 841.6667 841.6667 (그림 11-4)
# 11.5 강의 예제와 문제 3: 한 변 2 cm, 24 g인 정육면체
I_C = box_inertia(24, 2, 2, 2) # 대각 16
print(parallel_axis(I_C, 24, [1, 1, 1])) # [[ 64. -24. -24.] [-24. 64. -24.] [-24. -24. 64.]]
c = (np.arange(40) + 0.5) * 2 / 40 # 꼭짓점에서 잰 작은 조각의 위치
P = np.array(np.meshgrid(c, c, c)).reshape(3, -1).T
print(inertia_by_pieces(P, 24)) # 대각 63.99, 나머지 -24. (조각 64000개)
# 11.6 문제 4: 판 6 cm × 2 cm, 12 g을 30° 돌리기
I30 = rot_z(30) @ box_inertia(12, 6, 2, 0) @ rot_z(30).T
print(I30) # [[ 12. -13.8564 0.] [-13.8564 28. 0.] [0. 0. 40.]]
print(I30 @ [1, 0, 0]) # [ 12. -13.8564 0. ]
print(np.linalg.eigh(I30)[0]) # [ 4. 36. 40.] 주관성 모멘트
# 11.7 문제 5: ω = (10, 0, 0)으로 일정하게 돌릴 때의 ω × Iω
w = np.array([10, 0, 0])
print(np.cross(w, I30 @ w)) # [ 0. 0. -1385.6406]
def newton_euler(q, qd, qdd, links, grav=(0, 0)):
"""평면 팔의 뉴턴-오일러 (관절 축은 모두 z).
links: (질량, 길이, 관절에서 질량 중심까지, I_c) 목록. 돌려주는 값: τ와 관절 1의 힘 f₁."""
th = w = wd = 0.0
acc = -np.array(grav, float) # 받침 가속도 = 중력의 반대 (위로 g)
F, N, pc, pn = [], [], [], []
for (m, L, lc, Ic), q_, qd_, qdd_ in zip(links, q, qd, qdd):
th, w, wd = th + q_, w + qd_, wd + qdd_ # ① 바깥으로: ω, ω̇
u = np.array([np.cos(th), np.sin(th)])
z_u = np.array([-u[1], u[0]]) # z × u
acc_at = lambda k: acc + wd * k * z_u - w ** 2 * k * u
F.append(m * acc_at(lc)) # F = m v̇_c
N.append(Ic * wd) # 평면이라 ω × Iω = 0
pc.append(lc * u)
pn.append(L * u)
acc = acc_at(L) # 다음 관절의 가속도 v̇
cross = lambda a, b: a[0] * b[1] - a[1] * b[0]
f, n, tau = np.zeros(2), 0.0, np.zeros(len(links))
for i in reversed(range(len(links))): # ② 안으로: f, n, τ = n·z
n = N[i] + n + cross(pc[i], F[i]) + cross(pn[i], f)
f = F[i] + f
tau[i] = n
return tau, f
# 11.2, 11.9 세워 든 링크: 10 cm, 100 g, 세로 평면 (중력은 −y)
m, L, lc = 0.1, 0.1, 0.05
I0 = m * L ** 2 / 3
rod = [(m, L, lc, m * L ** 2 / 12)]
up = (0, -g)
print(round(I0, 7), round(m * g * lc, 4)) # 0.0003333 0.049
print(newton_euler([30 * D], [0], [0], rod, up)[0]) # [0.0424] 30°에서 버틸 토크
for th, thd, thdd in [(0.3, 1.5, -20.0), (2.0, -4.0, 100.0)]:
ne = newton_euler([th], [thd], [thdd], rod, up)[0][0]
print(round(ne, 6), round(I0 * thdd + m * g * lc * np.cos(th), 6))
# 0.040145 0.040145 / 0.012942 0.012942 (뉴턴-오일러 = 손으로 구한 식)
G60 = m * g * lc * np.cos(60 * D) # 문제 6
print(round(G60, 4), round((0.05 - G60) / I0, 1), round(-G60 / I0, 1)) # 0.0245 76.5 -73.5
print(*newton_euler([0], [2], [0], rod, up)) # [0.049] [-0.02 0.98] 문제 7: τ, f₁
# 순동역학(시뮬레이션): 수평에서 τ = 0으로 놓으면 언제, 얼마나 빨리 맨 아래를 지나나
th, thd, t, dt = 0.0, 0.0, 0.0, 1e-5
while th > -np.pi / 2:
thdd = (0 - m * g * lc * np.cos(th)) / I0 # θ̈ = (τ − G) / I₀
thd += thdd * dt
th += thd * dt
t += dt
print(round(t, 3), round(thd, 2), round(-np.sqrt(2 * m * g * lc / I0), 2))
# 0.153 -17.15 -17.15 (0.153초 뒤, 에너지 보존으로 구한 속도와 같다)
# 11.9 종이 팔: 링크 2와 손에 든 50 g을 한 덩어리로
def combine(parts):
"""한 링크 위의 (질량, 관절에서 거리, 자기 I_c) 조각들을 한 덩어리로 (평행축 정리)."""
mt = sum(p[0] for p in parts)
xc = sum(p[0] * p[1] for p in parts) / mt
return mt, xc, sum(p[2] + p[0] * (p[1] - xc) ** 2 for p in parts)
m2, x2, I2 = combine([(0.015, 0.04, 0.015 * 0.08 ** 2 / 12), (0.050, 0.08, 0.0)])
print(m2, round(x2 * 100, 2), round(I2 * 1e7, 1)) # 0.065 7.08 264.6
paper = [(0.020, 0.10, 0.05, 0.020 * 0.10 ** 2 / 12), (m2, 0.08, x2, I2)]
def mass_matrix(q):
"""멈춘 팔(q̇ = 0)에 단위 가속도를 하나씩 넣어 M을 한 열씩 채운다."""
return np.column_stack([newton_euler(q, [0, 0], e, paper)[0] for e in np.eye(2)])
print(round(4 * (np.pi / 2) / 0.25 ** 2, 1)) # 100.5 0.25초에 90°: θ̈ ≈ 100 rad/s²
for t2 in (0, 90, 150):
M = mass_matrix([20 * D, t2 * D])
print(t2, M * 1e7, M[:, 0] * 100) # M (g·cm²), θ̈₁ = 100일 때 (τ₁, τ₂) N·m
# 0 [[19886.6667 8120. ] [ 8120. 3520. ]] [0.1989 0.0812]
# 90 [[10686.6667 3520. ] [ 3520. 3520. ]] [0.1069 0.0352]
# 150 [[2719.233 -463.7169] [-463.7169 3520. ]] [ 0.0272 -0.0046]
M0 = mass_matrix([20 * D, 0])
print(round(-M0[1, 0] / M0[1, 1] * 100, 1)) # -230.7 팔꿈치 모터를 끄면 θ̈₂
# 11.10 라그랑주와 맞춰 보기: V = Ṁq̇ − ½[q̇ᵀ(∂M/∂qᵢ)q̇]
q, qd = np.array([20 * D, 90 * D]), np.array([10.0, 0.0])
V_ne = newton_euler(q, qd, [0, 0], paper)[0] # q̈ = 0이면 τ = V
h = 1e-6
dM = [(mass_matrix(q + h * e) - mass_matrix(q - h * e)) / (2 * h) for e in np.eye(2)]
Mdot = dM[0] * qd[0] + dM[1] * qd[1]
V_lag = Mdot @ qd - 0.5 * np.array([qd @ dM[i] @ qd for i in range(2)])
print(V_ne, V_lag) # [0. 0.046] [-0. 0.046]
# 운동 에너지: ½ q̇ᵀMq̇ = 링크마다 ½ m v_c² + ½ I_c ω² 의 합
q, qd = np.array([20 * D, 60 * D]), np.array([1.0, 2.0])
u1 = np.array([-np.sin(q[0]), np.cos(q[0])])
u12 = np.array([-np.sin(q[0] + q[1]), np.cos(q[0] + q[1])])
v1 = qd[0] * 0.05 * u1
v2 = qd[0] * 0.10 * u1 + (qd[0] + qd[1]) * x2 * u12
K = 0.5 * 0.020 * v1 @ v1 + 0.5 * paper[0][3] * qd[0] ** 2 + 0.5 * m2 * v2 @ v2 + 0.5 * I2 * (qd[0] + qd[1]) ** 2
print(round(0.5 * qd @ mass_matrix(q) @ qd * 1e7, 2), round(K * 1e7, 2)) # 26323.33 26323.33
스탠퍼드 숙제
스탠퍼드 강의 페이지의 Homework 5에서 1번(관성 텐서를 다른 좌표계로 옮기기)은 이 장으로 풀 수 있다. (a)는 11.6절의 식을 이끌어 내는 문제다. (b)는 상자 공식, (c)는 R I Rᵀ와 평행축 정리를 차례로 쓴다.
2번(RP 팔의 라그랑주 운동 방정식)은 12장을 배운 뒤에 푼다. 풀이는 같은 페이지의 Solution 5에 있다.
다 풀었으면 11강 영상을 한 번 본다.