PIDチューニング・チートシート

このシートは、関節、車輪、ジンバルの制御ループを立ち上げる人向けです。その上に開発を積み重ねられる程度のチューニングを得ることを目的としています。コードで書くときに当サイトが推奨するコントローラの形、段階的な手順、Ziegler-Nicholsの出発点とその限界をまとめました。各対策の効果を示すシミュレーションも1つ載せています。概念については、PID制御のレッスンと、用語集のPIDコントローラ・ステップ応答の項目から始めてください。

コードに書くべきコントローラ

形式	式	備考
並列形（ゲイン）	u = Kp e + Ki integral(e) + Kd de/dt	多くのコードやドライブが採用する形
標準形（時定数）	u = Kp (e + integral(e) / Ti + Td de/dt)	Ki = Kp / Ti、Kd = Kp Td。チューニング則はこの形で書かれている
目標値重み付き	u = Kp (beta r - y) + Ki integral(r - y) - Kd dy/dt	0から1の間のbetaで、外乱抑制を変えずに目標値応答を調整する。積分項は誤差全体に作用させる必要がある

下のコードで使っている、サンプリング周期Tsでの離散形は次のとおりです。

e_k   = r_k - y_k
P_k   = Kp (beta r_k - y_k)
D_k   = a D_(k-1) - (1 - a) Kd (y_k - y_(k-1)) / Ts      a = Tf / (Tf + Ts),  Tf = Td / N
u_raw = P_k + I_k + D_k
u_k   = clamp(u_raw, -u_max, u_max)
I_(k+1) = I_k + Ki Ts e_k     unless u_k is saturated and e_k pushes further into saturation

ゲインの値よりも重要な点が3つあります。

微分は誤差ではなく測定値に対してかけます。こうすると、rのステップ変化で微分キックが生じません。
微分項にはローパスフィルタをかけます。ÅströmとMurrayは、Tf = Td / N、Nは5から20の間としています。コードではN = 10を使います。
アンチワインドアップ。ここではクランプ（条件付き積分）を示します。バックカリキュレーションは、飽和による誤差をゲインを介して積分器にフィードバックする方式です。ÅströmとMurrayは、そのゲインをKiの倍数に選ぶことを提案しています。

ロボット関節のチューニング手順

符号を確認します。小さな正の出力を指令し、関節がエンコーダの正方向に動くことを確かめてください。符号が逆だと、どのゲインでも不安定に見えます。
先に制限を設定します。出力飽和、位置ウィンドウ、ドライブの電流制限です。飽和は応答を変えるので、実際の制限を入れた状態でチューニングしてください。
まずPのみで始めます。応答が速く、ただし明らかにオーバーシュートが大きすぎる状態になるまでKpを上げ、その値を半分にします。出発点となる数値が必要なら、後述の限界ゲイン法を使ってください。
Dを加えてオーバーシュートを減衰させます。オーバーシュートが減るか、モータからノイズが出始めるまでTdを大きくします。ノイズが出る場合は、Nが大きすぎるか、この大きさのDに対してエンコーダの分解能が粗すぎます。
Iは最後に、負荷外乱の抑制のために加えます。長いTi（遅い積分）から始め、負荷ステップからの回復が十分速くなるまで短くします。アンチワインドアップは最初のテストから有効にしておいてください。
目標値応答はオーバーシュートするが負荷からの回復は適切、という状態になったら、ゲインではなくbetaを下げます。ロボットは通常なめらかな軌道に追従するので、ランプ状の目標値でも同じ効果が得られます。
極端な条件で再調整します。最大可搬重量、腕を伸ばした姿勢、低温時と高温時です。軽い負荷で腕を畳んだ状態で調整したゲインは楽観的です。

Ziegler-Nichols

1942年に発表された古典的な2つの規則を、ÅströmとMurrayがまとめた表（Table 11.1）のとおりに示します。

制御器	ステップ応答法（Kp、Ti、Td）	限界ゲイン法（Kp、Ti、Td）
P	1/a	0.5 Ku
PI	0.9/a, tau/0.3	0.45 Ku, Tu/1.2
PID	1.2/a, tau/0.5, 0.5 tau	0.6 Ku, Tu/2, Tu/8

ステップ応答法では、aとtauを開ループステップ応答の最も急な接線から求めます。tauは見かけのむだ時間、a/tauは最大傾きです。限界ゲイン法では、Pのみのゲインをループが持続振動するまで上げます。そのときのゲインがKu、振動の周期がTuです。

ロボットに使う前に、次の点に注意してください。

ÅströmとMurrayは2つの欠点を挙げています。規則が使うプロセス情報が少なすぎること、そして得られるループのロバスト性が不足することです。下のシミュレーションのように、数十%のオーバーシュートを想定してください。
実機のアームを持続振動させると、ギアを傷めたり機械的ストッパに当たったりするおそれがあります。シミュレーション上で行うか、関節を固定してごく小さな振幅で行うか、振幅を制限するリレー試験を使ってください。
モータを囲む位置ループは、すでに積分器を含んでいます。そのため積分動作は、目標値ステップのたびにオーバーシュートを加えます。Iを弱めるか、目標値重み付けを使うことを想定してください。
サンプリングとフィルタリングは遅延を加えます。そのため、ディジタルループで測定したKuは連続時間の値より小さくなります。

シミュレーション

モータ遅れ（0.1 s）とアンプ遅れ（0.02 s）を介して駆動される関節角を考えます。伝達関数はG(s) = 1 / (s (0.1 s + 1)(0.02 s + 1))です。ループは1 kHzで動作し、出力制限は20です。t = 1.5 sから一定の負荷がかかります。このコードは、Pのみのループで二分法によりKuとTuを求め、Ziegler-Nicholsを適用します。その上で、アンチワインドアップ、目標値重み付け、手動でデチューンしたゲインの組を比較します。オーバーシュートと整定時間（2%バンド）は、負荷がかかる前の目標値ステップで測定します。

import numpy as np

DT = 0.001                      # controller period (s), 1 kHz
T_M, T_A, K = 0.10, 0.02, 1.0   # motor and amplifier time constants (s), plant gain
U_MAX = 20.0                    # actuator limit
LOAD = 2.0                      # constant load disturbance applied from t = 1.5 s


def plant_step(x: np.ndarray, u: float, load: float, sub: int = 10) -> np.ndarray:
    """Joint angle: K / (s (T_M s + 1)(T_A s + 1)), load enters at the motor."""
    pos, vel, drive = x
    h = DT / sub
    for _ in range(sub):
        drive += h * (K * u - drive) / T_A
        vel += h * (drive - load - vel) / T_M
        pos += h * vel
    return np.array([pos, vel, drive])


def make_pid(kp: float, ki: float, kd: float, tf: float, u_max: float, anti_windup: bool = True, beta: float = 1.0):
    s = {"i": 0.0, "d": 0.0, "y_prev": None}
    alpha = tf / (tf + DT)                        # first-order low-pass on the derivative

    def step(r: float, y: float) -> float:
        e = r - y
        y_prev = y if s["y_prev"] is None else s["y_prev"]
        s["d"] = alpha * s["d"] - (1 - alpha) * kd * (y - y_prev) / DT   # derivative on measurement
        s["y_prev"] = y
        u_unsat = kp * (beta * r - y) + s["i"] + s["d"]   # setpoint weight beta on P
        u = float(np.clip(u_unsat, -u_max, u_max))
        if not anti_windup or u == u_unsat or np.sign(e) != np.sign(u_unsat):
            s["i"] += ki * e * DT                 # clamping: freeze I while saturated and pushing further
        return u

    return step


def simulate(pid, t_end: float = 3.0, r: float = 1.0, load_on: bool = True) -> np.ndarray:
    x, ys = np.zeros(3), []
    for k in range(int(t_end / DT)):
        x = plant_step(x, pid(r, x[0]), LOAD if (load_on and k * DT >= 1.5) else 0.0)
        ys.append(x[0])
    return np.array(ys)


def growing(kp: float) -> bool:
    y = simulate(make_pid(kp, 0, 0, 1e-3, 1e9), t_end=4.0, r=0.01, load_on=False) - 0.01
    return np.abs(y[2000:]).max() >= np.abs(y[1000:2000]).max()


lo, hi = 1.0, 200.0
for _ in range(30):                               # bisection on P-only gain
    mid = 0.5 * (lo + hi)
    lo, hi = (lo, mid) if growing(mid) else (mid, hi)
ku = hi
y = simulate(make_pid(ku, 0, 0, 1e-3, 1e9), t_end=4.0, r=0.01, load_on=False) - 0.01
tu = np.diff(np.where(np.diff(np.sign(y[1000:])) > 0)[0]).mean() * DT
print(f"Ku = {ku:.1f} (continuous-time analytic {(T_M + T_A) / (T_M * T_A * K):.1f}), Tu = {tu:.3f} s")


def report(name: str, y: np.ndarray) -> None:
    step = y[:1500]
    outside = np.where(np.abs(step - 1.0) > 0.02)[0]
    settle = (outside[-1] + 1) * DT if len(outside) else 0.0
    print(f"{name:30s} overshoot {max(0.0, 100 * (step.max() - 1)):5.1f} %  settling {settle:.2f} s  "
          f"error at 3 s {1 - y[-1]:+.4f} rad")



kp, ti, td = 0.6 * ku, 0.5 * tu, 0.125 * tu       # classic Ziegler-Nichols PID
report("Z-N, no anti-windup", simulate(make_pid(kp, kp / ti, kp * td, td / 10, U_MAX, False)))
y_zn = simulate(make_pid(kp, kp / ti, kp * td, td / 10, U_MAX))
report("Z-N, clamping", y_zn)
report("Z-N, clamping, beta = 0", simulate(make_pid(kp, kp / ti, kp * td, td / 10, U_MAX, beta=0.0)))
kp2, ti2, td2 = 0.5 * ku, 1.0 * tu, 0.2 * tu      # detuned by hand from the Z-N start
y_tuned = simulate(make_pid(kp2, kp2 / ti2, kp2 * td2, td2 / 10, U_MAX, beta=0.3))
report("tuned, clamping, beta = 0.3", y_tuned)
for name, y in (("Z-N", y_zn), ("tuned", y_tuned)):
    print(name, [[round(k * DT, 2), round(float(y[k - 1]), 3)] for k in range(60, 3001, 60)])

出力の最初の数行は次のとおりです。

Ku = 58.7 (continuous-time analytic 60.0), Tu = 0.284 s
Z-N, no anti-windup            overshoot  83.8 %  settling 1.07 s  error at 3 s -0.0002 rad
Z-N, clamping                  overshoot  41.0 %  settling 0.87 s  error at 3 s -0.0002 rad
Z-N, clamping, beta = 0        overshoot  16.0 %  settling 0.86 s  error at 3 s -0.0002 rad
tuned, clamping, beta = 0.3    overshoot   0.0 %  settling 0.64 s  error at 3 s +0.0000 rad

測定したKuは58.7で、連続系プラントの60.0より小さくなっています。これは1 msのサンプル遅延によるものです。出力制限20は初期の比例動作（0.6 Kuに1 radを掛けて約35）より小さいので、各ステップの開始時にアクチュエータが飽和します。アンチワインドアップがないと、飽和中も積分が増え続け、オーバーシュートは2倍になります。クランプでその大部分が解消します。目標値重み付けは、負荷に対するループの応答を変えずに、オーバーシュートをさらに減らします。デチューンした組（0.5 Ku、Ti = Tu、Td = 0.2 Tu、beta = 0.3）は、オーバーシュートなしで0.64 sで2%バンドに入り、負荷による誤差も除去します。

t = 0での1 radのステップと、t = 1.5 sでの負荷ステップに対する関節角。クランプ付きのZiegler-Nicholsゲインと、目標値重み0.3の手動調整ゲインを比較する。点は60 msごとで、上のコードが出力した値。

症状と対処

症状	考えられる原因	対処
大きな移動の後に増大する遅い振動	積分ワインドアップ	アンチワインドアップ、または積分状態の制限
速いうなり、モータの発熱、耳で聞こえるノイズ	Dが大きすぎる、Nが大きすぎる、またはエンコーダの量子化	TdまたはNを下げる、測定値をフィルタする、エンコーダの分解能を上げる
いつまでも消えない定常偏差	Iがない、クランプでIが凍結している、または静止摩擦	Iを追加・確認する、摩擦フィードフォワードを加える
目標値付近の小さな振動	積分動作と静止摩擦によるスティックスリップ、またはバックラッシ	Iに不感帯を設ける、摩擦補償、Kiを下げる
目標値ステップでのみオーバーシュートし、負荷ステップでは問題ない	ステップ状の目標値に積分が作用している	目標値重み付け、または目標値をランプ状にする
先週は問題なかったゲインで振動する	可搬重量、温度、ファームウェアのループ周期の違い	極端な条件で再試験し、必要ならゲインスケジューリング

さらに読む

K. J. Åström and R. M. Murray、「Feedback Systems: An Introduction for Scientists and Engineers」第2版、Princeton University Press、無料オンライン版。第11章でPIDチューニング、ワインドアップ、実装を扱っています。
J. G. Ziegler and N. B. Nichols, “Optimum settings for automatic controllers”, 1942, doi:10.1115/1.4019264。
制御ハブとアクチュエータ・駆動系ハブ。
測定値をDに渡す前にノイズを取り除くためのカルマンフィルタ・チートシート。