Skip to content

Split View: 로봇 팔 제어 루프 — 궤적 생성부터 PID와 중력 보상까지

✨ Learn with Quiz
|

로봇 팔 제어 루프 — 궤적 생성부터 PID와 중력 보상까지

들어가며 — 목표 각도를 그대로 써 넣으면 팔이 튑니다

역기구학 계산이 끝났다고 하겠습니다. 손끝을 어디에 놓고 싶은지 정했고, 각 관절이 몇 도가 되어야 하는지도 나왔습니다. 어깨 0도, 팔꿈치 90도입니다.

그래서 코드에 이렇게 씁니다.

shoulder.write(0);
elbow.write(90);

팔이 튀어 나갑니다. 목표 각도 근처에서 두어 번 흔들리고, 팔 전체가 부르르 떨다가 멈춥니다. 물건을 들고 있었다면 놓쳤을 겁니다.

각도는 정확히 맞았는데 움직임이 엉망입니다. 이 글은 그 간극에 관한 것입니다.

문제의 정체는 이렇습니다. 위 두 줄이 관절에 요구한 것은 "0.02초 안에 90도만큼 이동하라"입니다. 서보 입장에서는 평균 각속도 4500도/초, 그리고 출발과 도착 순간의 가속도가 무한대입니다. 실제로 낼 수 있는 토크에는 상한이 있으므로 관절은 자기가 낼 수 있는 최대 토크로 밀어붙이다가 목표를 지나치고, 되돌아오다 또 지나칩니다.

제어는 두 층으로 나뉩니다. 위층은 시간에 따라 목표를 만들어 주는 궤적 생성이고, 아래층은 그 목표를 실제 각도가 따라가게 하는 피드백 루프입니다. 위층이 없으면 아래층이 아무리 좋아도 팔은 튑니다. 아래층이 없으면 위층이 아무리 매끄러워도 중력이 팔을 끌어내립니다.

사다리꼴 속도 프로파일 — 가장 단순한 궤적

궤적 생성의 가장 오래되고 가장 널리 쓰이는 형태입니다. 이름 그대로 속도를 사다리꼴로 만듭니다. 일정한 가속도로 올리고, 최대 속도로 순항하고, 일정한 가속도로 내립니다.

숫자를 넣어 보겠습니다. 어깨를 90도, 즉 1.5708 rad 움직입니다. 이 관절이 낼 수 있는 최대 속도를 1.0 rad/s, 최대 가속도를 2.0 rad/s²로 잡습니다.

가속 구간에서 최대 속도까지 걸리는 시간은,

t_a = v_max / a_max = 1.0 / 2.0 = 0.5초

그동안 이동한 각도는,

d_a = 0.5 × a_max × t_a² = 0.5 × 2.0 × 0.25 = 0.25 rad

감속 구간도 대칭이므로 같은 0.25 rad입니다. 둘을 합쳐 0.5 rad이 가감속에 쓰이고, 남은 각도는,

1.5708 - 0.5 = 1.0708 rad

이 구간을 최대 속도로 지나가므로,

t_c = 1.0708 / 1.0 = 1.0708초

전체 이동 시간은,

T = 0.5 + 1.0708 + 0.5 = 2.0708초

앞의 0.02초짜리 지령과 비교하면 100배 느립니다. 그리고 이 100배가 정확히 팔이 튀지 않는 이유입니다.

여기서 중요한 예외가 하나 있습니다. 이동 거리가 짧으면 최대 속도에 도달하기 전에 감속을 시작해야 합니다. 가감속에 필요한 최소 거리가 0.5 rad이었으므로, 0.2 rad을 움직인다면 순항 구간이 없는 삼각형 프로파일이 됩니다. 이때 도달하는 최고 속도는,

v_peak = √(a_max × d) = √(2.0 × 0.2) = 0.6325 rad/s
T = 2 × v_peak / a_max = 2 × 0.6325 / 2.0 = 0.6325초

이 분기를 코드에서 빼먹으면 짧은 이동에서 프로파일이 목표를 지나쳐 버립니다. 아래 구현에 그 판정이 들어 있습니다.

import numpy as np


def trapezoid(dq, v_max, a_max):
    """사다리꼴(또는 삼각형) 속도 프로파일의 구간 시간을 계산합니다."""
    dq = abs(dq)
    d_a = v_max ** 2 / (2 * a_max)      # 가속 구간에서 이동하는 거리
    if 2 * d_a >= dq:                    # 최대 속도에 도달하지 못합니다
        v_peak = np.sqrt(a_max * dq)
        t_a = v_peak / a_max
        return t_a, 0.0, t_a, v_peak
    t_a = v_max / a_max
    t_c = (dq - 2 * d_a) / v_max
    return t_a, t_c, t_a, v_max


def sample(t, dq, v_max, a_max):
    """시각 t에서의 위치와 속도를 돌려줍니다. 부호는 마지막에 붙입니다."""
    sign = 1.0 if dq >= 0 else -1.0
    t_a, t_c, t_d, v_p = trapezoid(dq, v_max, a_max)
    T = t_a + t_c + t_d
    t = min(max(t, 0.0), T)
    if t < t_a:
        q, v = 0.5 * a_max * t * t, a_max * t
    elif t < t_a + t_c:
        q, v = 0.5 * a_max * t_a ** 2 + v_p * (t - t_a), v_p
    else:
        td = t - t_a - t_c
        q = 0.5 * a_max * t_a ** 2 + v_p * t_c + v_p * td - 0.5 * a_max * td * td
        v = v_p - a_max * td
    return sign * q, sign * v


DQ, VMAX, AMAX = np.radians(90), 1.0, 2.0
ta, tc, td, vp = trapezoid(DQ, VMAX, AMAX)
print(f"가속 {ta:.4f}s  순항 {tc:.4f}s  감속 {td:.4f}s  전체 {ta+tc+td:.4f}s")
for t in (0.0, 0.25, 0.5, 1.0, 1.5708, 2.0708):
    q, v = sample(t, DQ, VMAX, AMAX)
    print(f"  t={t:.4f}s  q={np.degrees(q):8.4f}도  v={v:.5f} rad/s")

실행 결과입니다.

가속 0.5000s  순항 1.0708s  감속 0.5000s  전체 2.0708s
  t=0.0000s  q=  0.0000도  v=0.00000 rad/s
  t=0.2500s  q=  3.5810도  v=0.50000 rad/s
  t=0.5000s  q= 14.3239도  v=1.00000 rad/s
  t=1.0000s  q= 42.9718도  v=1.00000 rad/s
  t=1.5708s  q= 75.6763도  v=0.99999 rad/s
  t=2.0708s  q= 90.0000도  v=0.00000 rad/s

사다리꼴 프로파일에는 남는 문제가 하나 있습니다. t=0t=0.5에서 가속도가 0에서 2.0으로, 2.0에서 0으로 계단처럼 뜁니다. 가속도의 미분인 저크가 그 순간 무한대입니다. 가벼운 팔에서는 잘 안 보이지만 링크가 길고 가늘면 이 충격이 구조를 때려서 눈에 보이는 잔진동이 남습니다.

5차 다항식 궤적 — 가속도까지 이어 붙이기

저크 문제를 없애는 표준적인 방법은 위치를 시간의 다항식으로 두고, 양 끝에서 위치·속도·가속도를 모두 지정하는 것입니다. 조건이 여섯 개이므로 미지수도 여섯 개, 곧 5차 다항식이 됩니다.

q(t) = a0 + a1·t + a2·t² + a3·t³ + a4·t⁴ + a5·t⁵

양 끝에서 속도와 가속도를 0으로 두면 계수가 닫힌 형태로 나옵니다. 이동량을 Δq, 전체 시간을 T라고 하면,

a0 = q0,  a1 = 0,  a2 = 0
a3 =  10·Δq / T³
a4 = -15·Δq / T⁴
a5 =   6·Δq / T⁵

여기서 실무적으로 중요한 두 숫자가 나옵니다. 정규화한 시간 τ = t/T로 보면 속도는 τ = 0.5에서, 가속도는 τ = 0.2113τ = 0.7887에서 최대가 되고, 그 값은 각각 이렇습니다.

v_max = 1.875 × Δq / T
a_max = (10/√3) × Δq / T² = 5.7735 × Δq / T²

이제 앞의 사다리꼴과 같은 조건으로 비교해 보겠습니다. Δq = 1.5708 radT = 2.0708초에 움직이면,

v_max = 1.875 × 1.5708 / 2.0708 = 1.4223 rad/s
a_max = 5.7735 × 1.5708 / 2.0708² = 2.1149 rad/s²

같은 시간을 쓰는데 최대 속도가 1.0이 아니라 1.4223입니다. 42퍼센트 더 빠른 속도를 요구합니다. 반대로 최대 속도를 1.0으로 묶으면,

T = 1.875 × 1.5708 / 1.0 = 2.9452초

2.0708초짜리 이동이 2.9452초가 됩니다. 역시 42퍼센트 차이입니다.

항목사다리꼴5차 다항식
같은 시간(2.0708초)일 때 최대 속도1.0000 rad/s1.4223 rad/s
같은 시간일 때 최대 가속도2.0000 rad/s²2.1149 rad/s²
최대 속도를 1.0으로 묶었을 때 소요 시간2.0708초2.9452초
저크시작·끝에서 무한대어디서나 유한
모터 성능 활용률높음낮음
계산량조건 분기 3개계수 6개, 분기 없음

이 표가 말하는 것은 어느 쪽이 우월하다는 것이 아닙니다. 산업용 로봇이 사다리꼴 계열을 쓰는 이유는 사이클 타임이 곧 생산성이기 때문이고, 정밀 측정 장비나 카메라 짐벌이 다항식을 쓰는 이유는 잔진동이 곧 품질이기 때문입니다. 두 성질을 섞은 S-커브 프로파일이 실제로 가장 많이 쓰이는데, 사다리꼴의 가감속 모서리만 다항식으로 둥글린 형태라고 보면 됩니다.

import numpy as np


def quintic(q0, qf, T):
    """양 끝의 속도와 가속도가 0인 5차 다항식 계수를 돌려줍니다."""
    d = qf - q0
    return np.array([q0, 0.0, 0.0, 10 * d / T ** 3, -15 * d / T ** 4, 6 * d / T ** 5])


def evaluate(coef, t):
    powers = np.array([t ** i for i in range(6)])
    dpow = np.array([i * t ** (i - 1) if i >= 1 else 0.0 for i in range(6)])
    ddpow = np.array([i * (i - 1) * t ** (i - 2) if i >= 2 else 0.0 for i in range(6)])
    return coef @ powers, coef @ dpow, coef @ ddpow


T = 2.0
c = quintic(0.0, np.radians(90), T)
print("계수:", np.round(c, 6))
for t in (0.0, 0.5, 1.0, 1.5, 2.0):
    q, v, a = evaluate(c, t)
    print(f"  t={t}s  q={np.degrees(q):8.4f}도  v={v:8.5f} rad/s  a={a:8.5f} rad/s²")
print(f"  예측 최대 속도 {1.875*np.radians(90)/T:.6f}, 최대 가속도 {10/np.sqrt(3)*np.radians(90)/T**2:.6f}")

실행 결과입니다.

계수: [ 0.        0.        0.        1.963495 -1.472622  0.294524]
  t=0.0s  q=  0.0000도  v= 0.00000 rad/s  a= 0.00000 rad/s²
  t=0.5s  q=  9.3164도  v= 0.82835 rad/s  a= 2.20893 rad/s²
  t=1.0s  q= 45.0000도  v= 1.47262 rad/s  a= 0.00000 rad/s²
  t=1.5s  q= 80.6836도  v= 0.82835 rad/s  a=-2.20893 rad/s²
  t=2.0s  q= 90.0000도  v= 0.00000 rad/s  a= 0.00000 rad/s²
  예측 최대 속도 1.472622, 최대 가속도 2.267249

t=1.0에서 정확히 절반인 45도를 지나고 속도가 1.4726으로 최대이며, 양 끝에서 속도와 가속도가 모두 0입니다. 공식이 예측한 최대 속도와 표본 값이 일치합니다.

관절공간 보간과 작업공간 보간이 만드는 서로 다른 경로

여기서 갈림길이 하나 있습니다. 궤적을 관절 각도에 대해 만들 것인가, 손끝 위치에 대해 만들 것인가.

관절공간 보간은 시작 각도와 끝 각도 사이를 위 프로파일로 채웁니다. 각 관절이 독립적으로 자기 궤적을 따라가면 됩니다. 계산이 싸고, 관절 속도·가속도 한계를 직접 지킬 수 있고, 특이점을 신경 쓸 필요가 없습니다.

작업공간 보간은 시작 위치와 끝 위치 사이를 직선으로 잇고, 그 직선 위의 점마다 역기구학을 풀어 관절 각도를 얻습니다. 손끝이 실제로 직선을 그립니다.

두 방법이 만드는 경로가 얼마나 다른지 계산해 보겠습니다. 순기구학 편에서 쓴 것과 같은 2링크 팔입니다. 위팔 L1 = 0.20 m, 아래팔 L2 = 0.15 m.

출발 자세를 (30도, 45도), 도착 자세를 (-30도, 45도)로 두면 손끝은,

출발 (0.212028, 0.244889) m
도착 (0.318094, -0.061177) m

관절공간에서 정확히 절반인 (0도, 45도)의 손끝 위치는,

x = 0.20·cos(0) + 0.15·cos(45도) = 0.20 + 0.106066 = 0.306066
y = 0.20·sin(0) + 0.15·sin(45도) = 0 + 0.106066 = 0.106066

작업공간에서 직선의 중점은 두 끝점의 평균이므로,

x = (0.212028 + 0.318094) / 2 = 0.265061
y = (0.244889 - 0.061177) / 2 = 0.091856

두 점 사이의 거리는,

√((0.306066 - 0.265061)² + (0.106066 - 0.091856)²)
= √(0.041005² + 0.014210²)
= √(0.00168141 + 0.00020193)
= 0.043398 m = 43.4 mm

같은 두 점을 잇는데 중간에서 43밀리미터가 벌어집니다. 팔의 전체 길이가 350밀리미터이니 12퍼센트가 넘습니다. 컵에 물을 따르는 동작이라면 이 차이가 성공과 실패를 가릅니다. 반대로 상자 위 A 지점에서 B 지점으로 옮기기만 하는 동작이라면 43밀리미터 부풀어 오르는 경로가 오히려 장애물을 피해 주기도 합니다.

기준관절공간 보간작업공간 보간
손끝 경로예측하기 어려운 곡선직선(또는 지정한 곡선)
관절 속도 한계직접 지킬 수 있음간접적, 위반하기 쉬움
특이점통과해도 문제없음통과하면 관절 속도가 발산
매 주기 계산다항식 평가 한 번역기구학 풀이 한 번
자세 해가 바뀌는 문제없음elbow-up과 elbow-down 사이를 넘나들 수 있음
쓰는 곳이동, 대기 자세 복귀용접, 도포, 삽입, 따르기

마지막 줄이 실무에서 가장 아픈 항목입니다. 직선 위의 점마다 역기구학을 독립적으로 풀면, 어떤 점에서 갑자기 다른 해가 선택되어 팔꿈치가 반대로 뒤집힐 수 있습니다. 손끝은 여전히 직선 위에 있지만 팔 전체가 한 주기 만에 반대편으로 넘어갑니다. 해결책은 매번 새로 푸는 대신 직전 해에서 가장 가까운 해를 고르도록 강제하는 것입니다.

PID — 세 항이 각각 고치는 문제

궤적이 매 주기 목표 각도를 내놓으면, 아래층이 그 목표를 따라가야 합니다.

관절 하나를 이렇게 모형화하겠습니다. 회전 관성 J = 0.02 kg·m², 점성 마찰 b = 0.05 N·m·s/rad, 그리고 순기구학과 정적 토크 편에서 구한 어깨 중력 토크 1.79 N·m가 아래로 걸립니다. 목표는 90도입니다.

J·q̈ = τ - b·q̇ - τ_g

세 항을 하나씩 켜 보겠습니다.

비례항은 오차에 비례해 밀어 줍니다. Kp = 20으로 두고 이것만 켜면 관절이 목표를 크게 지나치고 한참 진동합니다. 마찰만으로는 에너지를 못 빼기 때문입니다.

미분항은 속도에 비례해 반대로 당깁니다. 물리적으로는 인공 마찰을 만드는 것입니다. Kd = 1.0을 더하면 오버슈트가 사라집니다. 그런데 목표에 도달하지 못하고 5.13도 아래에서 멈춥니다.

이 5.13도는 우연한 값이 아닙니다. 정지 상태에서는 속도가 0이므로 미분항이 0이고, 비례항 혼자 중력과 균형을 이룹니다.

Kp × e = τ_g
e = 1.79 / 20 = 0.0895 rad = 5.129도

손끝이 어깨에서 0.35미터 떨어져 있으므로,

0.0895 rad × 0.35 m = 0.0313 m = 31.3 mm

31밀리미터 아래로 처집니다. 팔은 조용히 멈춰 있고 오차도 안정적인데, 그냥 틀린 자리에 있습니다.

적분항은 이 잔류 오차를 시간에 걸쳐 누적해서 없앱니다. Ki = 40을 더하면 정상상태 오차가 0.07도로 떨어집니다. 대신 오버슈트가 4.75도 생깁니다. 적분항은 과거를 기억하므로 반응이 한 박자 늦습니다.

import numpy as np

J_INERTIA = 0.02     # kg·m²
B_VISCOUS = 0.05     # N·m·s/rad
TAU_GRAVITY = 1.79   # N·m, 어깨가 수평일 때 중력이 만드는 토크
DT, T_END = 0.001, 2.0
REF = np.radians(90)


def simulate(kp, ki, kd, feedforward=0.0):
    q, dq, integral = 0.0, 0.0, 0.0
    log = []
    for _ in range(int(T_END / DT)):
        e = REF - q
        integral += e * DT
        tau = kp * e + ki * integral + kd * (0.0 - dq) + feedforward
        ddq = (tau - B_VISCOUS * dq - TAU_GRAVITY) / J_INERTIA
        dq += ddq * DT
        q += dq * DT
        log.append(q)
    return np.array(log)


for name, (kp, ki, kd, ff) in {
    "P only":      (20, 0, 0.0, 0.0),
    "PD":          (20, 0, 1.0, 0.0),
    "PID":         (20, 40, 1.0, 0.0),
    "PD+gravity":  (20, 0, 1.0, TAU_GRAVITY),
}.items():
    log = simulate(kp, ki, kd, ff)
    err = np.degrees(REF - log[-1])
    over = max(0.0, np.degrees(log.max() - REF))
    print(f"{name:<11} 최종 {np.degrees(log[-1]):8.4f}도  정상상태오차 {err:8.4f}도  오버슈트 {over:7.4f}도")

실행 결과입니다.

P only      최종  78.4275도  정상상태오차  11.5725도  오버슈트 69.8285도
PD          최종  84.8720도  정상상태오차   5.1280도  오버슈트  0.0000도
PID         최종  90.0718도  정상상태오차  -0.0718도  오버슈트  4.7486도
PD+gravity  최종  90.0000도  정상상태오차  -0.0000도  오버슈트  0.6777도

PD의 정상상태 오차 5.1280도가 손으로 구한 5.129도와 일치합니다. 이론과 시뮬레이션이 같은 값을 내면 모형이 맞다는 뜻입니다.

마지막 줄에 0.6777도의 오버슈트가 남아 있는 것도 의미가 있습니다. 중력을 완전히 지워 버리면 남는 것은 순수한 2차 시스템이고, 이 게인에서의 감쇠비를 계산하면,

ζ = (Kd + b) / (2·√(Kp·J)) = 1.05 / (2·√(20 × 0.02)) = 1.05 / 1.2649 = 0.830

0.83이므로 이론적인 오버슈트는 약 0.9퍼센트, 90도에 대해 0.84도입니다. 시뮬레이션의 0.68도와 같은 자릿수입니다. PD만 켰을 때 오버슈트가 0이었던 것은 제어기가 잘해서가 아니라 중력이 팔을 계속 아래로 당겨 브레이크 역할을 했기 때문입니다. 중력 보상은 그 공짜 브레이크를 없애기도 하므로, 보상을 켠 뒤에는 Kd를 다시 봐야 합니다.

적분항을 쓸 때 반드시 함께 넣어야 하는 장치가 하나 있습니다. 팔이 물리적으로 막혀 있거나 모터가 토크 한계에 걸리면 오차가 줄지 않는데, 적분항은 그동안에도 계속 누적됩니다. 나중에 장애물이 치워지는 순간 어마어마하게 커진 적분항이 팔을 날려 버립니다. 이것이 적분 포화이고, 대책은 출력이 포화되면 적분을 멈추거나 적분값 자체에 상한을 두는 것입니다.

// 아두이노에서의 관절 하나 PID. 고정 주기로 도는 것이 전제입니다.
const float KP = 20.0f, KI = 40.0f, KD = 1.0f;
const float DT = 0.005f;             // 200Hz. 이 값이 실제 주기와 같아야 합니다
const float I_LIMIT = 3.0f;          // 적분 포화 방지 상한 (N·m 환산)
const float TAU_LIMIT = 8.0f;        // 모터가 낼 수 있는 토크 상한

float integral = 0.0f;
float prevMeasured = 0.0f;

float pidStep(float target, float measured, float gravityFeedforward) {
  float error = target - measured;

  // 미분은 오차가 아니라 측정값으로 계산합니다. 목표가 계단처럼 바뀌면
  // 오차의 미분이 순간적으로 폭발해 출력이 튀기 때문입니다.
  float derivative = -(measured - prevMeasured) / DT;
  prevMeasured = measured;

  float unsaturated = KP * error + KI * integral + KD * derivative + gravityFeedforward;

  // 출력이 포화된 방향으로는 적분을 더 쌓지 않습니다.
  bool pushingIntoLimit =
      (unsaturated > TAU_LIMIT && error > 0) || (unsaturated < -TAU_LIMIT && error < 0);
  if (!pushingIntoLimit) {
    integral += error * DT;
    if (integral > I_LIMIT) integral = I_LIMIT;
    if (integral < -I_LIMIT) integral = -I_LIMIT;
  }

  float tau = KP * error + KI * integral + KD * derivative + gravityFeedforward;
  if (tau > TAU_LIMIT) tau = TAU_LIMIT;
  if (tau < -TAU_LIMIT) tau = -TAU_LIMIT;
  return tau;
}

미분을 오차가 아니라 측정값으로 계산한 부분을 눈여겨봐 주십시오. 목표가 계단처럼 바뀌는 순간 오차의 미분은 이론상 무한대이고, 실제로는 한 주기 만에 출력이 상한까지 튑니다. 측정값의 미분을 쓰면 목표가 어떻게 바뀌든 이 항은 실제 관절 속도만 봅니다.

피드포워드와 중력 보상

위 시뮬레이션의 마지막 줄로 돌아가겠습니다. PD에 중력 토크를 그대로 더해 주었더니 정상상태 오차가 0.0000도가 되었습니다. 적분항을 하나도 쓰지 않고서입니다.

이것이 피드포워드입니다. 피드백은 오차가 생긴 뒤에야 반응합니다. 그런데 중력은 예측 가능합니다. 관절 각도만 알면 중력이 얼마나 걸리는지 지금 계산할 수 있습니다. 굳이 처지기를 기다렸다가 고칠 이유가 없습니다.

어깨의 중력 토크를 각도의 함수로 써 보겠습니다. 각 질량의 수평 거리에 무게를 곱해 더하면 됩니다.

τ_g(θ1, θ2) = g × [ m1·(L1/2)·cos(θ1)
                  + m2·(L1·cos(θ1) + (L2/2)·cos(θ1+θ2))
                  + (m_grip + m_pay)·(L1·cos(θ1) + L2·cos(θ1+θ2)) ]

앞 편의 값(m1 = 0.15, m2 = 0.10, 그리퍼 0.15, 페이로드 0.25 kg)을 넣으면,

어깨 각도팔꿈치 각도중력 토크
0도 (수평)0도1.7903 N·m
30도0도1.5505 N·m
60도0도0.8952 N·m
90도 (수직)0도0.0000 N·m
0도45도1.5964 N·m
0도90도1.1282 N·m

수직으로 세우면 정확히 0이고, 팔꿈치를 접으면 무게중심이 안쪽으로 오면서 줄어듭니다. 이 표가 곧 피드포워드 항입니다. 매 주기 현재 각도로 이 식을 계산해서 제어 출력에 더하면, 피드백은 모형이 놓친 나머지만 담당하면 됩니다.

여기서 정직해야 할 부분이 있습니다. 이 식은 모형입니다. 실제 팔의 질량 분포를 정확히 모르면 계산값이 어긋나고, 어긋난 만큼은 여전히 피드백이 처리해야 합니다. 그래도 1.79 N·m 전부를 피드백이 감당하는 것과 0.2 N·m만 감당하는 것은 완전히 다른 문제입니다. 필요한 Kp가 그만큼 작아지고, Kp가 작아지면 노이즈 증폭과 진동도 함께 줄어듭니다.

피드포워드로 넣을 수 있는 항은 중력만이 아닙니다. 궤적에서 이미 목표 속도와 목표 가속도를 알고 있으므로,

τ_ff = J·q̈_desired + b·q̇_desired + τ_g(q_desired)

이렇게 세 항을 다 넣으면 모형이 완벽할 때 피드백 출력이 0이 됩니다. 이것을 계산 토크 제어라고 부르고, 산업용 로봇 제어기가 하는 일이 대체로 이것입니다. 피드백은 모형이 틀린 만큼만 일하게 만드는 것이 제어 설계의 큰 방향입니다.

제어 주기와 지연이 안정성을 갉아먹는 방식

여기까지의 이야기는 전부 연속 시간이었습니다. 실제 제어기는 이산적으로 돌고, 그 사실이 안정성에 직접 영향을 줍니다.

디지털 제어 루프에는 최소 두 가지 지연이 있습니다. 하나는 영차 홀드에서 오는 평균 반 주기, 다른 하나는 센서를 읽고 계산해서 출력하기까지의 시간입니다. 둘을 합쳐 대략 1.5 주기로 잡는 것이 실무적인 어림입니다.

순수 지연이 위상에 미치는 영향은 주파수에 비례합니다.

위상 지연(라디안) = ω × T_delay

루프 이득이 1이 되는 주파수, 즉 교차 주파수를 20 rad/s(약 3.2Hz)로 잡고 계산해 보겠습니다.

제어 주기유효 지연(1.5주기)20 rad/s에서의 위상 지연
1000Hz1.50 ms1.72도
200Hz7.50 ms8.59도
100Hz15.00 ms17.19도
50Hz30.00 ms34.38도
50Hz:   20 × 0.030 = 0.600 rad = 34.38도
1000Hz: 20 × 0.0015 = 0.030 rad = 1.72도
차이:   32.66도

제어 주기를 1킬로헤르츠에서 50헤르츠로 낮추면 위상 여유 32.7도가 그냥 사라집니다. 통상적인 설계 목표가 45도에서 60도 사이인데, 32.7도를 잃으면 남는 것이 거의 없습니다. 게인은 하나도 안 건드렸는데 시스템이 발진합니다.

증상의 형태가 특징적입니다. 팔이 목표 근처에서 일정한 주파수로 떨고, 게인을 낮추면 멎고, 다시 올리면 같은 주파수로 떱니다. 그 주파수가 제어 주기와 관계있다면 원인은 게인이 아니라 지연입니다.

주기가 흔들리는 것도 같은 문제입니다. 아두이노에서 loop() 안에 PID를 넣고 delay(5)로 주기를 맞추는 코드를 흔히 보는데, 실제 주기는 5밀리초에 계산 시간이 더해진 값이고 시리얼 출력이라도 하면 훌쩍 늘어납니다. 미분항과 적분항이 모두 DT로 나누고 곱하므로, DT가 실제와 다르면 게인이 그 비율만큼 틀어집니다. 하드웨어 타이머 인터럽트로 고정 주기를 만들거나, 최소한 실제 경과 시간을 측정해 쓰는 편이 안전합니다.

지연을 만드는 것은 소프트웨어만이 아닙니다. 통신도 지연입니다. 하나의 버스에 데이지 체인으로 물리는 스마트 서보라면 관절 여섯 개의 상태를 한 바퀴 읽는 데 걸리는 시간이 곧 최소 제어 주기가 됩니다. I2C로 엔코더를 읽는다면 I2C 버스의 특성 때문에 클럭 스트레칭과 슬레이브 응답 시간이 주기에 그대로 들어옵니다.

시뮬레이터에서 맞춘 게인이 실기에서 깨지는 이유

이 글의 시뮬레이션에서 Kp = 20, Kd = 1.0은 아주 잘 동작했습니다. 같은 값을 실제 팔에 넣으면 대체로 깨집니다. 이유가 몇 가지 있고, 전부 모형에 없던 것들입니다.

첫째, 백래시입니다. 시뮬레이션에서 모터 축과 관절 축은 같은 각도입니다. 실물에서는 감속기 유격만큼 어긋납니다. 이 구간에서는 모터가 돌아도 관절이 안 움직이므로 제어기 입장에서는 이득이 0인 구간입니다. 유격을 지나 이가 맞물리는 순간 이득이 갑자기 정상으로 돌아옵니다. 이득이 구간마다 다른 시스템에서 하나의 게인 세트가 모든 구간에서 잘 동작하기는 어렵습니다. 증상은 목표 근처에서의 저주파 리밋 사이클, 즉 작은 폭으로 계속 왔다 갔다 하는 것입니다.

둘째, 관절 유연성입니다. 링크를 강체로 두었지만 실물은 휩니다. 하모닉 드라이브는 원리상 탄성체를 쓰고, 3D 프린트 링크는 눈에 보이게 휘고, 벨트는 늘어납니다. 이 탄성이 모터와 링크 사이에 공진을 만들고, 그 공진 주파수 위쪽에서는 위상이 180도 뒤집힙니다. 시뮬레이션의 강체 모형에는 이 극점이 아예 없습니다.

셋째, 마찰이 점성이 아닙니다. 모형의 b·q̇는 속도에 비례하는 점성 마찰인데, 실제 관절에서 지배적인 것은 속도와 무관한 쿨롱 마찰과, 정지 상태에서 더 큰 정지 마찰입니다. 저속에서 스틱슬립, 즉 붙었다가 미끄러지기를 반복하는 현상이 나옵니다. 아주 느리게 움직이라는 명령이 가장 어려운 명령인 이유입니다.

넷째, 엔코더 양자화와 미분 노이즈입니다. 12비트 엔코더는 한 바퀴를 4096으로 나누므로 분해능이 0.0879도입니다. 200Hz에서 이 한 눈금 차이를 미분하면,

0.0879도 / 0.005초 = 17.6도/초

관절이 완전히 정지해 있어도 측정값이 한 눈금 흔들리면 미분항은 17.6도/초에 해당하는 속도를 봅니다. Kd가 크면 이 잡음이 그대로 출력으로 나가서 모터가 웁니다. 그래서 미분항에는 거의 항상 저역통과 필터가 붙습니다.

다섯째, 전원입니다. 시뮬레이션의 모터는 요구한 토크를 즉시 냅니다. 실물에서는 모터를 구동하는 전원이 전류를 못 대면 토크가 안 나옵니다. 여러 관절이 동시에 가속하는 순간 전압이 무너지고, 제어기는 자기가 명령한 토크가 나갔다고 믿은 채로 다음 계산을 합니다. 팔이 무거운 자세로 갈 때만 이상해진다면 게인이 아니라 전류를 먼저 재 보십시오.

실무적인 순서는 이렇습니다. 시뮬레이터에서 얻는 것은 최종 게인이 아니라 게인의 자릿수와 구조입니다. 어느 항이 필요한지, 대략 어느 크기인지, 궤적이 관절 한계 안에 들어오는지를 확인하는 도구입니다. 실기에서는 중력 보상을 먼저 켜고, KiKd를 0으로 둔 채 Kp를 진동이 시작될 때까지 올린 뒤 절반으로 내리고, 그다음 Kd를 올려 진동을 잡고, 정상상태 오차가 남으면 마지막에 Ki를 아주 작게 넣습니다. 순서를 바꾸면 무엇이 무엇을 고쳤는지 알 수 없게 됩니다.

마치며 — 좋은 제어기는 대체로 좋은 궤적입니다

이 글에서 다룬 대책들은 층이 다릅니다.

궤적 생성은 애초에 불가능한 것을 요구하지 않는 일입니다. 90도를 0.02초에 가라는 명령은 어떤 제어기로도 잘 수행할 수 없습니다. 2.07초에 가라고 하면 평범한 제어기도 잘 따라갑니다.

피드포워드는 알고 있는 것을 미리 넣어 주는 일입니다. 중력은 매 순간 계산할 수 있는데 굳이 오차가 생기기를 기다릴 이유가 없습니다.

피드백은 나머지를 담당합니다. 모형이 놓친 마찰, 예상 못 한 외력, 부품 편차. 이 나머지가 작을수록 게인이 작아도 되고, 게인이 작을수록 시스템이 관대해집니다.

튀는 팔을 만났을 때 먼저 봐야 할 것은 게인이 아니라 지령입니다. 게인을 아무리 만져도 고쳐지지 않던 문제가 궤적 한 줄로 사라지는 경우가 실제로 많습니다. 그다음이 지연이고, 게인은 대체로 마지막입니다.

이 순서를 스스로 확인할 수 있는 실험이 하나 있습니다. 제어 주기를 절반으로 줄여 보십시오. 증상이 나아지면 지연 문제이고, 그대로면 모형 문제입니다. 게인은 그 둘을 구분한 다음에 만지는 것입니다. 나머지 배경이 되는 수학은 로봇공학에 필요한 수학 편에 순서대로 정리해 두었습니다.

The Robot Arm Control Loop: From Trajectory Generation to PID and Gravity Compensation

Introduction — Write the Target Angle Straight In and the Arm Jerks

Say inverse kinematics has finished its work. You've decided where to put the end effector, and you now know what angle each joint needs. Shoulder 0 degrees, elbow 90 degrees.

So you write this into your code:

shoulder.write(0);
elbow.write(90);

The arm lurches. It shakes back and forth a couple times near the target angle, trembles all over, and stops. If it had been holding something, it would have dropped it.

The angle was exactly right, but the motion is a mess. This post is about that gap.

Here's what's actually going on. Those two lines demanded of the joint: "move 90 degrees within 0.02 seconds." From the servo's point of view, that's an average angular velocity of 4500 degrees per second, and infinite acceleration at the instant of starting and stopping. Real torque has a ceiling, so the joint pushes with the maximum torque it can produce, overshoots the target, comes back, and overshoots again.

Control splits into two layers. The upper layer generates a target that changes over time — trajectory generation — and the lower layer is a feedback loop that makes the actual angle follow that target. Without the upper layer, no matter how good the lower layer is, the arm jerks. Without the lower layer, no matter how smooth the upper layer is, gravity drags the arm down.

Trapezoidal Velocity Profile — The Simplest Trajectory

This is the oldest and most widely used shape of trajectory generation. As the name says, it shapes velocity into a trapezoid: ramp up at constant acceleration, cruise at max velocity, ramp down at constant acceleration.

Let's plug in numbers. Move the shoulder 90 degrees, or 1.5708 rad. Take this joint's max velocity as 1.0 rad/s and max acceleration as 2.0 rad/s².

The time to reach max velocity during the acceleration phase is:

t_a = v_max / a_max = 1.0 / 2.0 = 0.5s

The angle traveled during that time is:

d_a = 0.5 × a_max × t_a² = 0.5 × 2.0 × 0.25 = 0.25 rad

The deceleration phase is symmetric, so it's also 0.25 rad. Together they use up 0.5 rad, leaving:

1.5708 - 0.5 = 1.0708 rad

This stretch is covered at max velocity, so:

t_c = 1.0708 / 1.0 = 1.0708s

The total move time is:

T = 0.5 + 1.0708 + 0.5 = 2.0708s

Compare this to the 0.02-second command from earlier and it's 100 times slower. And that 100x factor is exactly why the arm doesn't jerk.

There's an important exception here. If the move distance is short, deceleration has to start before max velocity is even reached. Since the minimum distance needed for accel-plus-decel was 0.5 rad, moving only 0.2 rad produces a triangular profile with no cruise phase at all. The peak velocity reached is:

v_peak = √(a_max × d) = √(2.0 × 0.2) = 0.6325 rad/s
T = 2 × v_peak / a_max = 2 × 0.6325 / 2.0 = 0.6325s

Miss this branch in your code and the profile overshoots the target on short moves. The implementation below includes this check.

import numpy as np


def trapezoid(dq, v_max, a_max):
    """Computes the segment times of a trapezoidal (or triangular) velocity profile."""
    dq = abs(dq)
    d_a = v_max ** 2 / (2 * a_max)      # distance covered during the acceleration phase
    if 2 * d_a >= dq:                    # max velocity never gets reached
        v_peak = np.sqrt(a_max * dq)
        t_a = v_peak / a_max
        return t_a, 0.0, t_a, v_peak
    t_a = v_max / a_max
    t_c = (dq - 2 * d_a) / v_max
    return t_a, t_c, t_a, v_max


def sample(t, dq, v_max, a_max):
    """Returns the position and velocity at time t. The sign is applied at the end."""
    sign = 1.0 if dq >= 0 else -1.0
    t_a, t_c, t_d, v_p = trapezoid(dq, v_max, a_max)
    T = t_a + t_c + t_d
    t = min(max(t, 0.0), T)
    if t < t_a:
        q, v = 0.5 * a_max * t * t, a_max * t
    elif t < t_a + t_c:
        q, v = 0.5 * a_max * t_a ** 2 + v_p * (t - t_a), v_p
    else:
        td = t - t_a - t_c
        q = 0.5 * a_max * t_a ** 2 + v_p * t_c + v_p * td - 0.5 * a_max * td * td
        v = v_p - a_max * td
    return sign * q, sign * v


DQ, VMAX, AMAX = np.radians(90), 1.0, 2.0
ta, tc, td, vp = trapezoid(DQ, VMAX, AMAX)
print(f"accel {ta:.4f}s  cruise {tc:.4f}s  decel {td:.4f}s  total {ta+tc+td:.4f}s")
for t in (0.0, 0.25, 0.5, 1.0, 1.5708, 2.0708):
    q, v = sample(t, DQ, VMAX, AMAX)
    print(f"  t={t:.4f}s  q={np.degrees(q):8.4f}°  v={v:.5f} rad/s")

Here's the output.

accel 0.5000s  cruise 1.0708s  decel 0.5000s  total 2.0708s
  t=0.0000s  q=  0.0000°  v=0.00000 rad/s
  t=0.2500s  q=  3.5810°  v=0.50000 rad/s
  t=0.5000s  q= 14.3239°  v=1.00000 rad/s
  t=1.0000s  q= 42.9718°  v=1.00000 rad/s
  t=1.5708s  q= 75.6763°  v=0.99999 rad/s
  t=2.0708s  q= 90.0000°  v=0.00000 rad/s

A trapezoidal profile still has one problem left over. At t=0 and t=0.5, acceleration jumps like a step from 0 to 2.0, then from 2.0 back to 0. Jerk — the derivative of acceleration — is infinite at those instants. A light arm barely notices, but a long, slender link takes that shock and shows visible residual vibration.

Quintic Polynomial Trajectory — Chaining Acceleration Together Too

The standard way to eliminate the jerk problem is to write position as a polynomial in time, and specify position, velocity, and acceleration at both ends. That's six conditions, so six unknowns, which means a 5th-order (quintic) polynomial.

q(t) = a0 + a1·t + a2·t² + a3·t³ + a4·t⁴ + a5·t⁵

With velocity and acceleration set to zero at both ends, the coefficients come out in closed form. With Δq as the total move and T as the total time:

a0 = q0,  a1 = 0,  a2 = 0
a3 =  10·Δq / T³
a4 = -15·Δq / T⁴
a5 =   6·Δq / T⁵

Two practically important numbers fall out here. In normalized time τ = t/T, velocity peaks at τ = 0.5, and acceleration peaks at τ = 0.2113 and τ = 0.7887, with these values:

v_max = 1.875 × Δq / T
a_max = (10/√3) × Δq / T² = 5.7735 × Δq / T²

Now let's compare against the trapezoidal profile under the same conditions. Move Δq = 1.5708 rad in T = 2.0708s:

v_max = 1.875 × 1.5708 / 2.0708 = 1.4223 rad/s
a_max = 5.7735 × 1.5708 / 2.0708² = 2.1149 rad/s²

Same total time, but the max velocity required is 1.4223, not 1.0 — 42 percent faster. Flip it around: cap max velocity at 1.0 instead, and:

T = 1.875 × 1.5708 / 1.0 = 2.9452s

The 2.0708-second move becomes 2.9452 seconds. Same 42 percent difference again.

ItemTrapezoidalQuintic polynomial
Max velocity at the same total time (2.0708s)1.0000 rad/s1.4223 rad/s
Max acceleration at the same total time2.0000 rad/s²2.1149 rad/s²
Time needed with max velocity capped at 1.02.0708s2.9452s
JerkInfinite at start and endFinite everywhere
Motor performance utilizationHighLow
Computation cost3 conditional branches6 coefficients, no branches

What this table says isn't that one is superior. Industrial robots use the trapezoidal family because cycle time directly is throughput, and precision measurement equipment or camera gimbals use polynomials because residual vibration directly is quality. The S-curve profile, a blend of both, is actually the most widely used in practice — think of it as a trapezoidal profile with its accel/decel corners rounded off by a polynomial.

import numpy as np


def quintic(q0, qf, T):
    """Returns the coefficients of a 5th-order polynomial whose velocity and acceleration are both zero at both ends."""
    d = qf - q0
    return np.array([q0, 0.0, 0.0, 10 * d / T ** 3, -15 * d / T ** 4, 6 * d / T ** 5])


def evaluate(coef, t):
    powers = np.array([t ** i for i in range(6)])
    dpow = np.array([i * t ** (i - 1) if i >= 1 else 0.0 for i in range(6)])
    ddpow = np.array([i * (i - 1) * t ** (i - 2) if i >= 2 else 0.0 for i in range(6)])
    return coef @ powers, coef @ dpow, coef @ ddpow


T = 2.0
c = quintic(0.0, np.radians(90), T)
print("Coefficients:", np.round(c, 6))
for t in (0.0, 0.5, 1.0, 1.5, 2.0):
    q, v, a = evaluate(c, t)
    print(f"  t={t}s  q={np.degrees(q):8.4f}°  v={v:8.5f} rad/s  a={a:8.5f} rad/s²")
print(f"  predicted max velocity {1.875*np.radians(90)/T:.6f}, max acceleration {10/np.sqrt(3)*np.radians(90)/T**2:.6f}")

Here's the output.

Coefficients: [ 0.        0.        0.        1.963495 -1.472622  0.294524]
  t=0.0s  q=  0.0000°  v= 0.00000 rad/s  a= 0.00000 rad/s²
  t=0.5s  q=  9.3164°  v= 0.82835 rad/s  a= 2.20893 rad/s²
  t=1.0s  q= 45.0000°  v= 1.47262 rad/s  a= 0.00000 rad/s²
  t=1.5s  q= 80.6836°  v= 0.82835 rad/s  a=-2.20893 rad/s²
  t=2.0s  q= 90.0000°  v= 0.00000 rad/s  a= 0.00000 rad/s²
  predicted max velocity 1.472622, max acceleration 2.267249

At t=1.0, it's exactly halfway through at 45 degrees, with velocity peaking at 1.4726, and velocity and acceleration are both zero at both ends. The predicted max velocity matches the sampled value.

The Different Paths Joint-Space and Task-Space Interpolation Produce

Here's a fork in the road. Should you generate the trajectory in terms of joint angles, or in terms of end-effector position?

Joint-space interpolation fills in the profile above between the start and end angle. Each joint follows its own trajectory independently. It's cheap to compute, joint velocity and acceleration limits are directly respected, and you don't need to worry about singularities.

Task-space interpolation connects the start and end position with a straight line, and solves inverse kinematics at every point along that line to get the joint angles. The end effector actually traces a straight line.

Let's calculate exactly how different the paths these two methods produce can be. Same 2-link arm as in the forward kinematics post: upper arm L1 = 0.20 m, forearm L2 = 0.15 m.

With the start pose at (30°, 45°) and the end pose at (-30°, 45°), the end effector is at:

Start (0.212028, 0.244889) m
End   (0.318094, -0.061177) m

In joint space, the exact midpoint (0°, 45°) puts the end effector at:

x = 0.20·cos(0) + 0.15·cos(45°) = 0.20 + 0.106066 = 0.306066
y = 0.20·sin(0) + 0.15·sin(45°) = 0 + 0.106066 = 0.106066

In task space, the midpoint of the straight line is the average of the two endpoints:

x = (0.212028 + 0.318094) / 2 = 0.265061
y = (0.244889 - 0.061177) / 2 = 0.091856

The distance between the two points is:

√((0.306066 - 0.265061)² + (0.106066 - 0.091856)²)
= √(0.041005² + 0.014210²)
= √(0.00168141 + 0.00020193)
= 0.043398 m = 43.4 mm

The two methods connect the exact same two points, and diverge by 43 millimeters in the middle. The arm's total length is 350 millimeters, so that's over 12 percent. If the motion is pouring water into a cup, this difference is the line between success and failure. On the other hand, if the motion is just moving from point A to point B above a box, a path that bulges out by 43 millimeters can actually help clear an obstacle.

CriterionJoint-space interpolationTask-space interpolation
End-effector pathA curve that's hard to predictA straight line (or whatever curve you specify)
Joint velocity limitsDirectly respectedIndirect, easy to violate
SingularitiesFine passing throughJoint velocity diverges passing through
Per-cycle computationOne polynomial evaluationOne inverse-kinematics solve
Solution-flip problemNoneCan jump between elbow-up and elbow-down
Used forPoint-to-point moves, returning to a rest poseWelding, coating, insertion, pouring

The last row is where this hurts most in practice. Solve inverse kinematics independently at every point along the line, and at some point a different solution can suddenly get selected, flipping the elbow to the opposite side. The end effector is still on the straight line, but the whole arm has jumped to the other side within a single cycle. The fix is to force each solve to pick the solution closest to the previous one, instead of solving fresh every time.

PID — What Each of the Three Terms Actually Fixes

Once the trajectory produces a target angle every cycle, the lower layer needs to make the actual angle follow it.

Let's model a single joint like this. Rotational inertia J = 0.02 kg·m², viscous friction b = 0.05 N·m·s/rad, and the shoulder's gravity torque of 1.79 N·m from the forward kinematics and static torque post pulling downward. The target is 90 degrees.

J·q̈ = τ - b·q̇ - τ_g

Let's switch on the three terms one at a time.

The proportional term pushes proportional to the error. Set Kp = 20 and turn on only this term, and the joint blows well past the target and oscillates for a long while. Friction alone can't remove enough energy.

The derivative term pulls back proportional to velocity. Physically, it's creating artificial friction. Add Kd = 1.0 and the overshoot disappears. But it never reaches the target — it stops 5.13 degrees short.

This 5.13 degrees isn't a random value. At rest, velocity is zero, so the derivative term is zero, and the proportional term alone balances against gravity.

Kp × e = τ_g
e = 1.79 / 20 = 0.0895 rad = 5.129°

Since the end effector is 0.35 meters from the shoulder:

0.0895 rad × 0.35 m = 0.0313 m = 31.3 mm

It sags 31 millimeters. The arm sits quietly, its error stable, but it's simply in the wrong place.

The integral term eliminates this leftover error by accumulating it over time. Add Ki = 40 and steady-state error drops to 0.07 degrees. In exchange, it produces 4.75 degrees of overshoot. The integral term remembers the past, so its response lags a beat.

import numpy as np

J_INERTIA = 0.02     # kg·m²
B_VISCOUS = 0.05     # N·m·s/rad
TAU_GRAVITY = 1.79   # N·m, torque gravity produces when the shoulder is horizontal
DT, T_END = 0.001, 2.0
REF = np.radians(90)


def simulate(kp, ki, kd, feedforward=0.0):
    q, dq, integral = 0.0, 0.0, 0.0
    log = []
    for _ in range(int(T_END / DT)):
        e = REF - q
        integral += e * DT
        tau = kp * e + ki * integral + kd * (0.0 - dq) + feedforward
        ddq = (tau - B_VISCOUS * dq - TAU_GRAVITY) / J_INERTIA
        dq += ddq * DT
        q += dq * DT
        log.append(q)
    return np.array(log)


for name, (kp, ki, kd, ff) in {
    "P only":      (20, 0, 0.0, 0.0),
    "PD":          (20, 0, 1.0, 0.0),
    "PID":         (20, 40, 1.0, 0.0),
    "PD+gravity":  (20, 0, 1.0, TAU_GRAVITY),
}.items():
    log = simulate(kp, ki, kd, ff)
    err = np.degrees(REF - log[-1])
    over = max(0.0, np.degrees(log.max() - REF))
    print(f"{name:<11} final {np.degrees(log[-1]):8.4f}°  steady-state error {err:8.4f}°  overshoot {over:7.4f}°")

Here's the output.

P only      final  78.4275°  steady-state error  11.5725°  overshoot 69.8285°
PD          final  84.8720°  steady-state error   5.1280°  overshoot  0.0000°
PID         final  90.0718°  steady-state error  -0.0718°  overshoot  4.7486°
PD+gravity  final  90.0000°  steady-state error  -0.0000°  overshoot  0.6777°

PD's steady-state error of 5.1280 degrees matches the 5.129 degrees we worked out by hand. When theory and simulation agree, that means the model is right.

The 0.6777 degrees of overshoot left over on the last line means something too. Strip gravity out entirely and what's left is a pure second-order system, and computing the damping ratio at these gains:

ζ = (Kd + b) / (2·√(Kp·J)) = 1.05 / (2·√(20 × 0.02)) = 1.05 / 1.2649 = 0.830

At 0.83, the theoretical overshoot is about 0.9 percent, or 0.84 degrees on a 90-degree move. The same order of magnitude as the simulation's 0.68 degrees. The reason PD alone had zero overshoot wasn't that the controller was doing a great job — it was that gravity kept pulling the arm downward, acting as a free brake. Gravity compensation removes that free brake too, so once you turn it on, you need to look at Kd again.

There's one device you absolutely need to add whenever you use the integral term. If the arm is physically blocked, or the motor hits its torque limit, the error stops shrinking — but the integral term keeps accumulating the whole time regardless. The moment the obstacle clears later, that hugely inflated integral term flings the arm violently. This is integral windup, and the fix is to freeze the integral once the output saturates, or cap the integral value itself.

// PID for a single joint on an Arduino. Assumes it runs on a fixed period.
const float KP = 20.0f, KI = 40.0f, KD = 1.0f;
const float DT = 0.005f;             // 200Hz. This value must match the actual period.
const float I_LIMIT = 3.0f;          // Integral windup guard, in N·m terms
const float TAU_LIMIT = 8.0f;        // Torque limit the motor can produce

float integral = 0.0f;
float prevMeasured = 0.0f;

float pidStep(float target, float measured, float gravityFeedforward) {
  float error = target - measured;

  // The derivative is computed from the measurement, not the error. If the target
  // jumps like a step, the derivative of the error would momentarily explode
  // and the output would spike.
  float derivative = -(measured - prevMeasured) / DT;
  prevMeasured = measured;

  float unsaturated = KP * error + KI * integral + KD * derivative + gravityFeedforward;

  // Don't keep accumulating the integral in the direction that's pushing into the limit.
  bool pushingIntoLimit =
      (unsaturated > TAU_LIMIT && error > 0) || (unsaturated < -TAU_LIMIT && error < 0);
  if (!pushingIntoLimit) {
    integral += error * DT;
    if (integral > I_LIMIT) integral = I_LIMIT;
    if (integral < -I_LIMIT) integral = -I_LIMIT;
  }

  float tau = KP * error + KI * integral + KD * derivative + gravityFeedforward;
  if (tau > TAU_LIMIT) tau = TAU_LIMIT;
  if (tau < -TAU_LIMIT) tau = -TAU_LIMIT;
  return tau;
}

Pay attention to computing the derivative from the measurement rather than the error. The instant a step-changing target arrives, the derivative of the error is theoretically infinite, and in practice the output spikes to its limit within a single cycle. Using the derivative of the measurement means this term only ever looks at the joint's actual velocity, no matter how the target changes.

Feedforward and Gravity Compensation

Let's go back to the last line of the simulation above. Adding the gravity torque straight into PD brought steady-state error to 0.0000 degrees — without using the integral term at all.

This is feedforward. Feedback only reacts after an error has already appeared. But gravity is predictable. Knowing only the joint angle, you can compute how much gravity is pulling right now. There's no reason to wait for the arm to sag and then correct it.

Let's write the shoulder's gravity torque as a function of angle. Multiply each mass's horizontal distance by its weight and add them up.

τ_g(θ1, θ2) = g × [ m1·(L1/2)·cos(θ1)
                  + m2·(L1·cos(θ1) + (L2/2)·cos(θ1+θ2))
                  + (m_grip + m_pay)·(L1·cos(θ1) + L2·cos(θ1+θ2)) ]

Plug in the values from the earlier post (m1 = 0.15, m2 = 0.10, gripper 0.15, payload 0.25 kg):

Shoulder angleElbow angleGravity torque
0° (horizontal)1.7903 N·m
30°1.5505 N·m
60°0.8952 N·m
90° (vertical)0.0000 N·m
45°1.5964 N·m
90°1.1282 N·m

Standing it up vertically, it's exactly zero, and folding the elbow shrinks it as the center of mass moves inward. This table is exactly the feedforward term. Compute this formula every cycle with the current angle and add it to the control output, and feedback only has to handle whatever the model missed.

Something worth being honest about here: this equation is a model. If it doesn't precisely match the arm's real mass distribution, the computed value will be off, and feedback still has to handle whatever's off by. Even so, having feedback carry the full 1.79 N·m versus only 0.2 N·m is a completely different problem. The Kp you need shrinks accordingly, and a smaller Kp means less noise amplification and less vibration too.

Gravity isn't the only thing you can feed forward. Since the trajectory already knows the target velocity and target acceleration:

τ_ff = J·q̈_desired + b·q̇_desired + τ_g(q_desired)

Put all three terms in and, when the model is perfect, feedback output goes to zero. This is called computed torque control, and it's largely what an industrial robot's controller does. Making feedback work only as hard as the model is wrong is the general direction of control design.

How Control Period and Latency Eat Away at Stability

Everything up to this point was continuous time. A real controller runs discretely, and that fact directly affects stability.

A digital control loop has at least two sources of delay: the zero-order hold's average half-period, and the time it takes to read the sensor, compute, and produce the output. A common practical rule of thumb is to lump these together as roughly 1.5 periods.

Pure delay's effect on phase is proportional to frequency.

Phase lag (radians) = ω × T_delay

Take the crossover frequency — where the loop gain equals 1 — as 20 rad/s (about 3.2Hz) and calculate:

Control rateEffective delay (1.5 periods)Phase lag at 20 rad/s
1000Hz1.50 ms1.72°
200Hz7.50 ms8.59°
100Hz15.00 ms17.19°
50Hz30.00 ms34.38°
50Hz:   20 × 0.030 = 0.600 rad = 34.38°
1000Hz: 20 × 0.0015 = 0.030 rad = 1.72°
Difference: 32.66°

Drop the control rate from 1 kilohertz to 50 hertz and 32.7 degrees of phase margin just vanishes. A typical design target is 45 to 60 degrees, so losing 32.7 degrees leaves almost nothing. Not a single gain got touched, and yet the system oscillates.

The symptom has a distinctive shape. The arm trembles near the target at a fixed frequency, quiets down if you lower the gain, and trembles at the same frequency again if you raise it back up. If that frequency has anything to do with the control period, the cause is delay, not gain.

A wobbling period is the same kind of problem. It's common to see code that puts PID inside loop() on an Arduino and matches the period with delay(5) — but the actual period is 5 milliseconds plus whatever computation time got added on top, and it stretches noticeably further if there's even a serial print in there. Since the derivative and integral terms both divide and multiply by DT, if DT doesn't match reality, the gains get skewed by that same ratio. Using a hardware timer interrupt for a fixed period, or at minimum measuring actual elapsed time and using that, is the safer approach.

Software isn't the only source of delay. Communication is delay too. If you're daisy-chaining smart servos on a single bus, the time it takes to read the state of all six joints in one pass becomes your minimum control period. If you're reading an encoder over I2C, I2C bus characteristics mean clock stretching and slave response time land directly in your period.

Why Gains Tuned in a Simulator Break on Real Hardware

In this post's simulation, Kp = 20 and Kd = 1.0 worked beautifully. Put those same values on a real arm and they usually break. There are several reasons, and every one of them is something the model didn't have.

First, backlash. In the simulation, the motor shaft and the joint shaft are at the same angle. In reality, they're off by however much slop the gearbox has. In that region, the motor turns but the joint doesn't move, so as far as the controller's concerned, gain is zero in that zone. The instant the teeth mesh again, gain snaps back to normal. A system whose gain differs region to region is hard for a single set of gains to handle well everywhere. The symptom is a low-frequency limit cycle near the target — small oscillations back and forth that never settle.

Second, joint flexibility. Links were treated as rigid, but real ones flex. A harmonic drive uses an elastic element by design, a 3D-printed link visibly bends, and a belt stretches. This flexibility creates resonance between the motor and the link, and above that resonant frequency, phase flips 180 degrees. The simulation's rigid-body model doesn't have this pole at all.

Third, friction isn't viscous. The model's b·q̇ is viscous friction proportional to velocity, but what actually dominates in a real joint is Coulomb friction, independent of velocity, and static friction, which is even larger at rest. At low speed you get stick-slip — sticking, then slipping, repeatedly. This is exactly why commanding very slow motion is the hardest command to execute well.

Fourth, encoder quantization and derivative noise. A 12-bit encoder divides one revolution into 4096 steps, giving a resolution of 0.0879 degrees. Differentiate a single tick of noise at 200Hz:

0.0879° / 0.005 second = 17.6°/second

Even with the joint completely still, if the reading wobbles by one tick, the derivative term sees a velocity equivalent to 17.6 degrees per second. With a large Kd, this noise goes straight through to the output and the motor whines. That's why the derivative term almost always gets a low-pass filter.

Fifth, the power supply. In the simulation, the motor produces the requested torque instantly. In reality, if the power driving the motor can't deliver the current, torque doesn't show up. The instant several joints accelerate at once, voltage collapses, and the controller keeps computing its next step believing it commanded that torque and got it. If things only go strange when the arm moves into a heavy pose, measure current before you touch the gains.

The practical order looks like this. What you get from a simulator isn't your final gains — it's the order of magnitude and structure of the gains. It's a tool for confirming which terms you need, roughly what size they should be, and whether the trajectory stays within joint limits. On real hardware, turn on gravity compensation first, hold Ki and Kd at zero and raise Kp until it starts oscillating, then halve it; next raise Kd to damp the oscillation; and only at the very end, if steady-state error remains, add in a very small Ki. Change the order and you lose track of what fixed what.

Conclusion — A Good Controller Is Usually a Good Trajectory

The countermeasures covered in this post operate at different layers.

Trajectory generation is about never demanding the impossible in the first place. No controller can execute "get to 90 degrees in 0.02 seconds" well. Tell it to get there in 2.07 seconds, and even an ordinary controller follows along just fine.

Feedforward is about feeding in what you already know ahead of time. Gravity can be computed at every instant — there's no reason to wait for an error to appear first.

Feedback handles everything else: friction the model missed, unexpected external forces, part-to-part variation. The smaller that remainder is, the smaller the gains can be, and the smaller the gains are, the more forgiving the whole system becomes.

When you run into a jerky arm, the first thing to check isn't the gains — it's the command. Problems that no amount of gain-tweaking ever fixed disappear all the time with a single line of trajectory. Delay comes next, and gains are usually last.

There's one experiment you can run yourself to confirm this order: cut the control period in half. If the symptom improves, it's a delay problem; if it doesn't change, it's a model problem. Gains are what you touch after you've told those two apart. The math underlying all of this is laid out in order in the math you need for robotics post.