PIDゲインの決め方とC言語実装｜Arduinoで10パターン実測した

式は読んだ、コードも写した。それでも、目の前の装置でKp・Ki・Kdをいくつにすればいいのかが決まらない。

ゲインは、先に制御対象を測ってから計算で出します。勘で決めるものではありません。この記事では、LEDの明るさをフォトレジスタで測って合わせる装置をArduinoで組み、Kp・Ki・Kdを10パターン振って実測した波形を並べます。部品は4点、測定は全部で約2分です。

＜自己紹介＞筆者は現役の車載組み込みソフトエンジニアです。仕事ではRH系・SH系のマイコンで制御ソフトを書いています。この記事の数値は、すべて手元のArduino Uno R3で測ったものです。

まず、PID制御がどういう形をしているかを1枚で見ておきます。

やっていることは「測る → くらべる → 動かす → また測る」の繰り返しだけです。輪の中でPIDが担当するのは制御器の箱ひとつで、残りは装置の側にあります。

今回作るもの：部品4点・配線10分

LEDを光らせて、その明るさをフォトレジスタで測り、目標の明るさになるようPWMを調節します。使うのはArduino Uno R3のほか、LED・220Ω・フォトレジスタ・10kΩの4点だけです。トランジスタも電池も要りません。配線は10分ほどで終わります。

この装置を選んだ理由は2つあります。
1つは速さ。1パターンが5〜9秒で終わるので、9本を約2分で取り切れます。ゲインを振って比べる作業は回数が要るので、1分待つ相手だと同じことが1日仕事になります。
もう1つは明るさが対数で返ってくること。同じ装置のまま、どの明るさで使うかによって必要なゲインが変わるところまで見えます。「ゲインの決め方」を扱う記事なら、ここまで見せられる相手のほうが役に立ちます（動作点が変わると、ゲインも変わる）。

つなぎ方は2つに分かれます。出す側はD3から220Ωを通してLEDへ。(5V−2V)÷220Ω＝13.6mAです。ATmega328Pデータシート の Absolute Maximum Ratings で1ピンあたりの上限は40mAです。ただしこれは「これを超えたら壊れる」値で、常用の目安は20mA以下。この電流ならトランジスタは要りません。測る側は5V─フォトレジスタ─A0─10kΩ─GNDの分圧で、明るいほどフォトレジスタの抵抗が下がってA0の電圧が上がります。

LEDとフォトレジスタは向かい合わせにして、上から黒いテープか紙コップで覆ってください。部屋の照明が入ると、それがそのまま外乱になって毎回ちがう波形になります。

教科書の式をC言語にする

PIDの式は3つの項の足し算。連続系の式を、そのままC言語に置き換えます。

積分は「足していく変数」、微分は「前回との差を時間で割る」だけです。∫と d/dt という記号が難しく見えるだけで、コードにすると2行です。

float pid_update(pid_t *p, float setpoint, float meas)
{
    float error = setpoint - meas;

    /* P項：偏差にそのまま比例 */
    float p_term = p->kp * error;

    /* I項：偏差を足していく。dtを掛けるので単位は「秒あたり」 */
    p->integral += error * p->dt;
    float i_term = p->ki * p->integral;

    /* D項：前回からの変化。dtで割るので単位は「1秒あたり」 */
    float d_term = p->kd * (error - p->prev_error) / p->dt;
    p->prev_error = error;

    float u = p_term + i_term + d_term;

    /* 出せる範囲に丸める */
    if (u > p->out_max) { u = p->out_max; }
    if (u < p->out_min) { u = p->out_min; }
    return u;
}

これが動く最小形。ただしこのままだと2つの問題が起きます。積分ワインドアップと微分キックで、どちらもあとの節で実測を出します。

制御周期をタイマ割り込みで作る

PID制御は「一定の間隔で回す」ことが前提です。式のdtが一定でないと、積分も微分も値が狂います。

loop()の中でdelay(10)を使うと、そのあいだの処理時間が加わって周期がばらつきます。タイマ割り込みで一定の周期を作り、割り込みではフラグを立てるだけにします。

const float DT_S = 0.01f;          /* 制御周期 10ms */
volatile bool g_tick = false;

void timer1_begin_10ms(void)
{
    noInterrupts();
    TCCR1A = 0;
    TCCR1B = 0;
    TCNT1  = 0;
    /* 16MHz / 256 = 62500Hz。10ms = 625カウント */
    OCR1A  = 625 - 1;
    TCCR1B |= (1 << WGM12);        /* CTCモード */
    TCCR1B |= (1 << CS12);         /* 分周256 */
    TIMSK1 |= (1 << OCIE1A);
    interrupts();
}

ISR(TIMER1_COMPA_vect)
{
    g_tick = true;                 /* フラグを立てるだけ */
}

void loop(void)
{
    if (!g_tick) { return; }
    g_tick = false;
    /* ここでPIDを1回まわす */
}

タイマ割り込みで周期を作る手順そのものは、タイマ割り込みで一定周期を作る手順で分周比の計算から扱っています。

間に合っているかは、自分で数えて確かめます。1周期ぶんの処理を終えた時点で次の割り込みフラグが立っていたら、間に合っていません。

if (g_tick) { g_overrun++; }   /* 次のtickが来てしまった＝間に合っていない */

ゲインを決める前に、制御対象を測る

ここが「ゲインの決め方」の入口です。Kpをいくつにするかは、制御対象がどれだけ反応するかで決まります。反応の大きさを知らずに数字を決めることはできません。

測るのは3つです。

1. 静ゲインPWMを1上げると、測定値が何％上がるか
2. むだ時間と時定数出力を変えてから、測定値が動きはじめるまでと、落ち着くまで
3. 安定限界Kpをどこまで上げると、振動が止まらなくなるか

3つ目だけは計算で出せません。上げていって、振動した値を見つけるしかありません。そのかわり、これがいちばん確かな数字になります。

まず静ゲインです。LEDを消したときと全開にしたときの明るさを測れば出ます。

この数字があると、目標値に必要な操作量が計算で出ます。目標25％なら (25−0.48) ÷ 0.173 ＝ 142カウント。上限255に収まるので届きます。目標60％なら344カウント必要で、上限を超えるので絶対に届きません。ワインドアップの節ではこれを利用します。

2点で出した傾きは、その場の傾きではない

ところが、実際に操作量を142にしても25％にはなりません。本当は95カウントで25％に届きます。計算が1.5倍ずれています。

理由は、消灯と全開の2点を結んだ直線が平均の傾きでしかないことです。確かめるために、操作量を5カウントずつ上げながら明るさを測りました（各点で0.3秒待ってから16回平均、52点）。

一番急なところと一番なだらかなところで3.5倍違います。暗いほど1カウントの効きが大きく、明るくなるほど鈍ります。フォトレジスタの抵抗が光に対して対数的に変わるためで、明るさの読み方をどう作ってもこの性質は残ります。

「振動していない」の基準線を先に取る

このあと波形を「振れ幅」で比べます。そのとき、何％なら振動していないと言えるのかが決まっていないと比べられません。

そこでPIDを止めて、操作量を固定したまま、制御と同じ10ms周期で2秒ぶん記録しました。

操作量を動かしていなくても、明るさは1.9％ぶん揺れています。ADCのばらつきとLEDの明滅です。消灯すると0.09％まで下がるので、揺れの出どころはLEDを点けていること自体だと分かります。
つまり振れ幅2％前後は「振動していない」と読みます。以降の波形はこの線と比べてください。

速いのは取り柄だが、そのままではPIDの見せ場が無い

操作量を0から128へ一段だけ上げて、明るさの立ち上がりを1ms刻みで350点記録しました。制御周期の10msでは2点しか入らないので、この測定だけ別のルーチンで撮っています。3回繰り返した平均が次の値です。

/* むだ時間：測定値をN サンプル遅らせてからPIDへ渡す */
float dead_apply(float v, int n)
{
    if (n <= 0) {
        return v;              /* 遅らせない */
    }
    if (n >= DEAD_MAX) {
        n = DEAD_MAX - 1;      /* バッファからはみ出さない */
    }
    g_dead_buf[g_dead_idx] = v;
    g_dead_idx = (g_dead_idx + 1) % DEAD_MAX;
    int k = (g_dead_idx - 1 - n + DEAD_MAX * 2) % DEAD_MAX;
    return g_dead_buf[k];
}

むだ時間L＝50ms、時定数T＝24.7ms、静ゲインK＝0.211が分かれば、Ziegler-Nicholsのステップ応答法でゲインの出発点が計算で出ます。

安定限界を実測して、計算と突き合わせる

計算で出た2.8が妥当かどうかは、どこまで上げたら振動するかを知らないと判断できません。そこでKiとKdを0にしたP制御だけで、Kpを0.25から8まで12段階に振りました。

境界がここまで鋭いのは、むだ時間が支配的な系の特徴です。時定数が支配的な系なら、振れ幅はもっとなだらかに増えます。

この記事で使う6つの指標

以降の波形は、次の6つの数字で比べます。定義を先に決めておかないと、同じ波形から違う結論が出ます。

Kpを振る：小さいと届かない、大きいと振動する

まずP制御だけで動かします。変えるのはKpだけで、他は何も変えていません。

float p_term = 1.0f * error;

float p_term = 8.0f * error;

Kiを足す：残った偏差が消える

Kiを足すと、P制御で17.67％残っていた偏差が−0.05％まで消えます。積み上げた差が、Pだけでは届かなかった分を埋めるからです。

p->integral += error * p->dt;
float i_term = 28.0f * p->integral;

Kdを足す：収まっていた制御が振動に変わった

教科書では、D項は行き過ぎを抑える働きをすると説明されます。同じKp・Kiのまま、Kdだけを0から0.070にして測りました。この0.070は、さきほどZiegler-Nicholsの式が出した値そのものです。

理由はむだ時間にあります。D項は50ms前の変化に反応するので、押すタイミングがずれます。今回の装置はむだ時間50ms・時定数24.7msで、L/T＝2.0の「むだ時間が支配的」な系です。ずれた向きに押す力を足したので、ループが自分で振動を保つようになりました。

ただし、ここまではすべて1回ずつの測定です。1回の測定でゲインを決めないで同じ条件を4回ずつ繰り返し、この4つの差が測定のばらつきでは説明できないことを確かめています。

積分ワインドアップ：届かない目標を出したあと

ここがPID制御でいちばん実害の出るところです。目標値が装置の能力を超えていると、偏差がいつまでも残るため積分が際限なく育ちます。

先ほど計算したとおり、この装置で目標60％を出すには344カウント必要で、上限255では届きません。LEDを全開にしても44.67％までしか出ないので、どうやっても届かない目標です。2.0秒でそこへ上げ、5.0秒で元へ戻します。

p->integral += error * p->dt;   /* 無条件に足す */

float i_cand = p->ki * (p->integral + error * p->dt);
float u_cand = p_term + i_cand + d_term;

/* 出せる範囲を超えていて、しかもさらに押し込む向きのときだけ、足さない */
bool pushing_further = (u_cand > p->out_max && error > 0.0f) ||
                       (u_cand < p->out_min && error < 0.0f);
if (!pushing_further) {
    p->integral += error * p->dt;
}

D項を使えるようにする：ローパスと微分キック対策

ここまでD項は悪者でしたが、直し方は2つあります。ローパスフィルタと、微分の元を偏差から測定値へ変えることです。

この2つは直す相手が違います。まとめて入れると、どちらが直したのか分からなくなります。1つずつ入れて測りました。目標値は2.5秒で15％から32％へ変えています。

まずローパスを入れる

D項は測定値の変化を見るので、センサのノイズをそのまま増幅します。それなら、微分した値を平らにしてから使えば済むはずです。

/* 1次のローパス。alphaが小さいほど平らになる（この記事は0.1） */
p->d_filtered += p->d_lpf_alpha * (d_raw - p->d_filtered);
float d_term = p->d_filtered;

残るのは、目標値を変えた瞬間だけ

ローパスで振動は止まりました。それでも1か所だけ跳ねます。D項を偏差から計算していると、目標値を変えた瞬間に偏差が段差になり、微分が跳ね上がります。

float d_raw = p->kd * (error - p->prev_error) / p->dt;
p->d_filtered += p->d_lpf_alpha * (d_raw - p->d_filtered);

/* 目標値の変化に反応しないよう、測定値の変化だけを見る（符号は反転） */
float d_raw = -p->kd * (meas - p->prev_meas) / p->dt;
p->d_filtered += p->d_lpf_alpha * (d_raw - p->d_filtered);

Kdを上げすぎると、微分項がノイズを増幅する

Kdを10倍にすると、微分項もちょうど10倍になります。D項は測定値の変化を見ているので、ノイズも同じ倍率で増えます。

float d_raw = -0.70f * (meas - p->prev_meas) / p->dt;

1回の測定でゲインを決めない

ここまでの数値は、すべて1回ずつの測定です。装置に触っていなくても、測るたびに値は動きます。その動きより小さい差は、差ではありません。

PIとPIDを4回ずつ、交互に測りました。順番の影響が出ないよう、PI→PID を1組として4回繰り返しました。PIを4回続けてからPIDを4回にすると、その間に装置がずれたぶんが差に化けます。

上の4行は、Kdを足すとどれも悪くなることを別の角度から示しています。いちばん下の行だけが、同じ「行き過ぎ」を測り方だけ変えたものです。ここがこの節でいちばん大事なところなので、順に見ていきます。

ゲインを変えて1回測り、良くなったから採用する。この手順では、実際にはばらつきを見ているだけということが起こります。

この結論が当てはまる範囲

この記事の数字がそのまま使えるのは、L/T＝2.0、つまりむだ時間が時定数の2倍ある相手までです。手順は他の対象へ持っていけますが、数字は持っていけません。

あなたの制御対象に読み替える

この記事は光でやりましたが、測る量と動かす量が決まっていれば、中身は同じです。自分の対象に置き換えるときは、次の6つが何になるかを先に決めてください。

ゲインを決める手順：どこから手を付けるか

ここまでの実測を、手を動かす順に並べ直します。この順でやると、うまくいかなかったときにどこへ戻ればいいかが分かります。

1. 静ゲインを測る操作量を0と最大にして、測定値がどこからどこまで動くか。この記事なら 0.48％〜44.67％。ただしゲインの計算には、使う動作点でステップを入れて測った傾きを使う（この記事は0.211）
2. むだ時間と時定数を測る操作量にステップを入れて、動きはじめるまで（L）と落ち着くまで（T）。L/Tが1より大きければ、Dは生のままでは使えないと思っておく
3. 計算で出発点を出すa＝K・L÷T から Kp＝1.2÷a、Ti＝2L、Td＝0.5L。ここはまだ出発点で、答えではありません
4. Ki＝0・Kd＝0にしてKpだけ振る振動しはじめる値を見つけ、その半分から3分の1に戻します。偏差は残りますが、それはIで埋めます
5. Kiを足すTi＝Kp÷Ki。むだ時間の3〜4倍から始め、短くしていって振動が出たら戻します
6. Kdは生のまま試さないTd＝Kd÷Kpがむだ時間の半分あたりを試します。ただし最初からローパスと測定値微分を入れてください。この記事では、生の微分だと振れ幅が30.58％まで荒れたのに、Kdを変えずにローパスを入れるだけで2.30％まで収まりました
7. 各段階で複数回測る差がばらつきより小さいなら、その変更に意味があるかは分かりません。指標も1つに絞らないでください

Ti・Td・制御周期を「時間」で決める

ここまでの実測から、実際に手を動かすときの勘所をまとめます。

制御周期は時定数の何分の1にするか

目安は時定数の5分の1から20分の1です。細かくすれば良くなるというものではありません。粗すぎると制御が追いつかず、細かすぎると1周期あたりの変化がノイズに埋もれます。

この装置は時定数が24.7ms、制御周期が10msなので、時定数のなかに2.5サンプルしか入っていません。かなり粗い部類です。
ただしこの記事の波形を振動させているのは、足した50msのむだ時間のほう。遅れが50msでサンプリングが10msなら、支配しているのはむだ時間です。

KiとKdは「時間」で考えると決めやすい

Ki・Kdの数字は、そのままでは大きいのか小さいのか分かりません。時間に直すと、むだ時間や時定数と直接くらべられます。

Ziegler-Nicholsのステップ応答法は、この2つをむだ時間Lから直接出します。PI制御ならTi＝3.3L、PID制御ならTi＝2L・Td＝0.5L。時間のほうはLだけで決まります。ただしKi＝Kp÷Ti、Kd＝Kp×Td なので、ゲインに直すにはKpが要り、Kpを出すには静ゲインKと時定数Tも要ります。3つとも測らないと数字にはなりません。
そしてこの記事で見たとおり、出た値が使えるかどうかは別の話です。Ti＝0.1秒はそのまま当たりましたが、Td＝25msは微分の作り方を直すまで使えませんでした。

動作点が変わると、ゲインも変わる

静ゲインを52点で測った表のとおり、1カウントの効きは動作点で3.5倍変わります。変わるのは静ゲインだけではありません。装置には一切触らず、ステップを入れる場所だけを変えて確かめました。

明るい側では、静ゲインが半分になり、時定数も半分になりました。フォトレジスタ（CdS）は明るいほど反応が速く、そのぶん1カウントの効きが鈍ります。

安定限界の探し方

Kpを上げていって振動しはじめる値を見つけ、そこから戻します。このとき、KiとKdは0のままにしてください。Iが入っていると、振動しているのかIが遅れて追いついているのかが見分けられません。

上の実測では、Kp＝5.5までは振れ幅3.03％で静かなのに、Kp＝6.0で21.74％へ跳びました。0.5上げただけです。実用のKpは、この値の半分から3分の1に取ります。この記事の2.8は安定限界の49％にあたります。

全体のコード

実測に使ったPID本体です。10パターンはすべて、このコードのゲインとフラグを変えただけで、本体は1行も書き換えていません。

typedef struct {
    float kp, ki, kd;
    float dt;              /* 制御周期[s]。一定であることが前提 */
    float integral;
    float prev_error;
    float prev_meas;
    float d_filtered;      /* 微分項のローパス後の値 */
    float d_lpf_alpha;     /* 1.0でフィルタなし */
    float out_min, out_max;
    bool  anti_windup;     /* 飽和するなら積分しない */
    bool  deriv_on_meas;   /* 測定値を微分する */
    bool  started;         /* 初回の微分を0にするため */
} pid_t;

float pid_update(pid_t *p, float setpoint, float meas)
{
    float error = setpoint - meas;

    /* --- P項 --- */
    float p_term = p->kp * error;

    /* --- D項：偏差から取るか、測定値から取るか --- */
    float d_raw;
    if (!p->started) {
        d_raw = 0.0f;              /* 初回。前回値が無いので微分しない */
    } else if (p->deriv_on_meas) {
        d_raw = -p->kd * (meas - p->prev_meas) / p->dt;
    } else {
        d_raw = p->kd * (error - p->prev_error) / p->dt;
    }
    p->d_filtered += p->d_lpf_alpha * (d_raw - p->d_filtered);
    float d_term = p->d_filtered;

ここまでがPとDです。Dはderiv_on_measで「偏差を微分」と「測定値を微分」を切り替え、そのあとローパスを通します。微分キックの2つの対策が、この数行に入っています。

残るI項は、足す前に「足したらどうなるか」を計算します。これが積分ワインドアップの対策です。

    /* --- I項：足す前に、足した結果が出せる範囲か確かめる --- */
    float i_cand = p->ki * (p->integral + error * p->dt);
    float u_cand = p_term + i_cand + d_term;

    bool saturating = (u_cand > p->out_max) || (u_cand < p->out_min);
    bool pushing_further = (u_cand > p->out_max && error > 0.0f) ||
                           (u_cand < p->out_min && error < 0.0f);

    /* 飽和していても、戻る向きの偏差なら積分してよい */
    if (!(p->anti_windup && saturating && pushing_further)) {
        p->integral += error * p->dt;
    }
    float i_term = p->ki * p->integral;

    p->prev_error = error;
    p->prev_meas  = meas;
    p->started    = true;

    float u = p_term + i_term + d_term;
    if (u > p->out_max) { u = p->out_max; }
    if (u < p->out_min) { u = p->out_min; }
    return u;
}

うまく動かないとき

練習問題

まとめ

この記事で書いたPIDのコードは、そのまま他の制御対象にも使えます。変えるのはゲインと制御周期だけ。ただし制御対象が変われば静ゲインもむだ時間も変わるので、ゲインは測り直してください。