Chỉnh PID: Ziegler–Nichols, thử rơ le tìm Ku Tu, công thức SIMC và chỉnh tay theo đáp ứng

Cách chỉnh thông số PID: ảnh hưởng của Kp, Ki, Kd; Ziegler–Nichols bằng thử rơ le tìm Ku, Tu; công thức SIMC; so sánh trên mô hình lò nhiệt; chỉnh tay; ESP32 tự dò PID cho lò nhiệt dùng SSR.

Viết bộ điều khiển PID thì dễ, chỉnh PID mới khó: Kp, Ki, Kd chọn sai thì lò nhiệt vọt quá nhiệt độ, động cơ rung, mực nước dao động mãi. Bài này giải thích từng hệ số ảnh hưởng thế nào, rồi đi qua ba cách chỉnh: phương pháp Ziegler–Nichols dùng thử rơ le để tìm hệ số tới hạn Ku và chu kỳ Tu, công thức SIMC dựa trên mô hình quá trình, và chỉnh tay theo đáp ứng. Các cách được so sánh trên mô hình lò nhiệt bằng script Python. Bài có chương trình ESP32 tự dò thông số bằng thử rơ le cho lò nhiệt dùng SSR.

1. Bộ điều chỉnh PID là gì, ba hệ số làm gì

Bộ điều chỉnh PID tính tín hiệu điều khiển (công suất, duty) từ sai lệch e = giá trị đặt − giá trị đo:

u = Kp × e + Ki × (tổng e theo thời gian) + Kd × (tốc độ thay đổi)

• P (tỉ lệ): sai lệch càng lớn, tác động càng mạnh. Chỉ có P thường còn sai lệch tĩnh: cần một sai lệch nhỏ để giữ được công suất.
• I (tích phân): cộng dồn sai lệch theo thời gian, đẩy sai lệch tĩnh về 0; nhiều quá gây vọt lố, dao động chậm.
• D (vi phân): phản ứng với tốc độ thay đổi, giống phanh, giảm vọt lố; nhưng khuếch đại nhiễu cảm biến.

Cách viết PID trên ESP32 (vi phân theo giá trị đo, chống bão hòa tích phân) có ở bài giữ tốc độ động cơ DC ổn định bằng PID trên ESP32. Bài này tập trung vào chọn số.

2. Ziegler–Nichols và thử rơ le

Phương pháp Ziegler–Nichols vòng kín: chỉ dùng P, tăng Kp tới khi hệ dao động đều không tắt dần. Kp lúc đó là hệ số tới hạn Ku, chu kỳ dao động là Tu; rồi tính Kp, Ti, Td theo bảng. Đẩy hệ thật tới bờ mất ổn định thì nguy hiểm (lò nhiệt vọt nhiệt, máy rung mạnh), nên ngày nay thường dùng thử rơ le (Åström–Hägglund) để tìm Ku, Tu an toàn hơn:

1. Thay bộ điều khiển bằng một rơ le: đo cao hơn điểm đặt thì giảm công suất xuống mức thấp, thấp hơn thì tăng lên mức cao (mức nền ± d).
2. Hệ dao động đều quanh điểm đặt với biên độ a và chu kỳ Tu, biên độ được giới hạn bởi d.
3. Tính Ku = 4d / (π a); Tu chính là chu kỳ dao động.

Hai lưu ý:

• Ku từ thử rơ le là giá trị gần đúng. Với lò nhiệt mô phỏng trong bài, nó thấp hơn giá trị tính đúng khoảng 18 %, còn Tu gần như chính xác.
• Bộ số Ziegler–Nichols được thiết kế để dao động tắt dần nhanh nên khá "gắt": vọt lố 30–40 % là bình thường. Coi đó là điểm xuất phát, sau đó chỉnh tay (mục 4).

3. So sánh trên mô hình lò nhiệt

Nhiều quá trình (lò nhiệt, bể nước, nhiệt độ phòng) mô tả gần đúng bằng mô hình bậc nhất có trễ (FOPDT):

• độ lợi K: nhiệt độ tăng bao nhiêu mỗi % công suất;
• hằng số thời gian T;
• thời gian trễ L.

Ba thông số này đo được bằng một thử nghiệm nhảy bậc: đổi công suất một bước, ghi lại nhiệt độ theo thời gian. Khi đã có K, T, L thì dùng được công thức SIMC (Skogestad). Chọn τc = L cho đáp ứng nhanh vừa phải, τc lớn hơn cho êm hơn.

Script dưới đây làm cả hai việc trên mô hình lò nhiệt K = 2 °C/%, T = 60 s, L = 10 s:

• thử rơ le để lấy Ku, Tu, so với giá trị tính đúng;
• tính ba bộ hệ số, mô phỏng điểm đặt tăng 30 °C rồi thêm tải;
• đo vọt lố, thời gian ổn định, mức tụt nhiệt khi có tải.

# PYWORLD.VN - So sánh các cách chỉnh PID trên mô hình lò nhiệt bậc nhất có trễ (FOPDT): Ziegler–Nichols, SIMC
# PYWORLD.VN - Thử nghiệm rơ le (Åström–Hägglund) để dò Ku, Tu giống chương trình ESP32 trong bài, so với giá trị tính đúng
# PYWORLD.VN - Chạy: python chinh_pid.py --K 2 --T 60 --L 10   (K °C mỗi %, T hằng số thời gian, L thời gian trễ, giây)
"""Chỉnh PID cho quá trình K·e^(−Ls)/(Ts+1): thử rơ le, Ziegler–Nichols, SIMC; mô phỏng nhảy bậc và nhiễu tải."""
from __future__ import annotations

import argparse
import math

DT_SIM, DT_DK = 0.1, 1.0                       # bước mô phỏng, chu kỳ bộ điều khiển (giây)


class QuaTrinh:
    """Lò nhiệt: y là độ tăng nhiệt so với môi trường; u là % công suất, tác động sau L giây."""

    def __init__(self, K, T, L):
        self.K, self.T = K, T
        self.hang = [0.0] * max(1, int(round(L / DT_SIM)))
        self.y = 0.0

    def buoc(self, u, nhieu=0.0):
        self.hang.append(u)
        u_tre = self.hang.pop(0)
        self.y += (self.K * u_tre - nhieu - self.y) / self.T * DT_SIM
        return self.y


def ku_tu_dung(K, T, L):
    """Ku, Tu tính đúng: tìm w để pha của K e^(-jwL)/(jwT+1) bằng -180 độ."""
    lo, hi = 1e-6, math.pi / L
    for _ in range(200):
        w = (lo + hi) / 2
        if w * L + math.atan(w * T) < math.pi:
            lo = w
        else:
            hi = w
    w = (lo + hi) / 2
    return math.sqrt(1 + (w * T) ** 2) / K, 2 * math.pi / w


def thu_role(K, T, L, d=10.0, tre_tren=0.0, t_max=3000.0):
    """Thử rơ le quanh điểm đặt 30 (độ tăng): u = nen ± d, nen là công suất giữ được điểm đặt. Trả về (Ku, Tu) từ 4 chu kỳ cuối."""
    qt = QuaTrinh(K, T, L)
    dat, t = 30.0, 0.0
    nen = dat / K
    u = nen + d
    # dua len gan diem dat bang cong suat nen
    while qt.y < dat:
        qt.buoc(100.0)
        t += DT_SIM
    dinh, day, chuyen = [], [], []
    y_max, y_min = -1e9, 1e9
    while t < t_max:
        y = qt.buoc(u)
        t += DT_SIM
        y_max, y_min = max(y_max, y), min(y_min, y)
        if u > nen and y > dat + tre_tren:
            u = nen - d
            chuyen.append(t)
            day.append(y_min)
            y_max = -1e9
        elif u < nen and y < dat - tre_tren:
            u = nen + d
            dinh.append(y_max)
            y_min = 1e9
    a = (sum(dinh[-4:]) / 4 - sum(day[-4:]) / 4) / 2
    tu = (chuyen[-1] - chuyen[-5]) / 4
    return 4 * d / (math.pi * a), tu


def mo_phong(K, T, L, Kp, Ki, Kd, t_max=1500.0, t_nhieu=900.0, nhieu=10.0):
    """PID rời rạc giống chương trình ESP32: D theo giá trị đo, chống bão hòa tích phân, ra 0..100 %."""
    qt = QuaTrinh(K, T, L)
    dat, tp, y_truoc, u = 30.0, 0.0, 0.0, 0.0
    kq = []
    buoc = int(round(DT_DK / DT_SIM))
    for n in range(int(t_max / DT_SIM)):
        t = n * DT_SIM
        if n % buoc == 0:
            e = dat - qt.y
            uu = Kp * e + Ki * tp - Kd * (qt.y - y_truoc) / DT_DK
            y_truoc = qt.y
            if not ((uu >= 100 and e > 0) or (uu <= 0 and e < 0)):
                tp += e * DT_DK
            u = min(100.0, max(0.0, uu))
            kq.append((t, qt.y, u))
        qt.buoc(u, nhieu if t >= t_nhieu else 0.0)
    return kq


def danh_gia(kq, dat=30.0, t_nhieu=900.0):
    tren = [p for p in kq if p[0] < t_nhieu]
    vot = (max(p[1] for p in tren) - dat) / dat * 100
    on_dinh = None
    for i in range(len(tren) - 1, -1, -1):
        if abs(tren[i][1] - dat) > 0.02 * dat:
            on_dinh = tren[i + 1][0] if i + 1 < len(tren) else None
            break
    sau = [p for p in kq if p[0] >= t_nhieu]
    tut = dat - min(p[1] for p in sau)
    return vot, on_dinh, tut


def bo_chinh(K, T, L, ku, tu):
    tau_c = L                                        # SIMC: chọn hằng số thời gian vòng kín bằng thời gian trễ
    kc = T / (K * (tau_c + L))
    ti = min(T, 4 * (tau_c + L))
    return {"ZN PID (thu role)": (0.6 * ku, 0.6 * ku / (tu / 2), 0.6 * ku * tu / 8),
            "ZN PI (thu role)": (0.45 * ku, 0.45 * ku / (tu / 1.2), 0.0),
            "SIMC PI (tau_c = L)": (kc, kc / ti, 0.0)}


def main() -> None:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--K", type=float, default=2.0, help="độ lợi quá trình, °C mỗi %% công suất")
    ap.add_argument("--T", type=float, default=60.0, help="hằng số thời gian, s")
    ap.add_argument("--L", type=float, default=10.0, help="thời gian trễ, s")
    a = ap.parse_args()
    ku_d, tu_d = ku_tu_dung(a.K, a.T, a.L)
    ku_r, tu_r = thu_role(a.K, a.T, a.L)
    print(f"Qua trinh K={a.K:g} C/%, T={a.T:g} s, L={a.L:g} s")
    print(f"Ku, Tu tinh dung   : Ku = {ku_d:.3f} %/C, Tu = {tu_d:.1f} s")
    print(f"Ku, Tu thu ro le   : Ku = {ku_r:.3f} %/C, Tu = {tu_r:.1f} s  (lech {100 * (ku_r / ku_d - 1):+.0f} %, {100 * (tu_r / tu_d - 1):+.0f} %)")
    print("Nhay bac dat +30 C luc 0 s; tai them (mat 10 C) luc 900 s:")
    for ten, (kp, ki, kd) in bo_chinh(a.K, a.T, a.L, ku_r, tu_r).items():
        vot, od, tut = danh_gia(mo_phong(a.K, a.T, a.L, kp, ki, kd))
        print(f"  {ten:20s} Kp {kp:6.3f} Ki {ki:7.4f} Kd {kd:6.3f} | vot lo {vot:5.1f} %, on dinh 2 % sau {od:5.0f} s, tai lam tut {tut:4.1f} C")


if __name__ == "__main__":
    main()

Kết quả trên mô hình này:

• SIMC PI vọt lố chỉ 6 % và ổn định nhanh nhất.
• Bộ Ziegler–Nichols PID kéo nhiệt độ về nhanh nhất khi có tải (tụt ít nhất), đổi lại vọt lố lớn khi đổi điểm đặt.

Không có bộ số "đúng" cho mọi trường hợp: chọn theo việc quan trọng hơn là không vượt nhiệt độ hay là giữ nhiệt độ khi có tải.

4. Chỉnh tay theo đáp ứng

Khi không đo được mô hình, hoặc để tinh chỉnh sau Ziegler–Nichols, làm theo trình tự:

1. Đặt Ki = 0, Kd = 0. Tăng dần Kp tới khi đáp ứng bắt đầu dao động, rồi lùi về khoảng một nửa.
2. Tăng dần Ki tới khi sai lệch tĩnh hết trong thời gian chấp nhận được; vọt lố tăng quá mức thì giảm Ki.
3. Nếu còn vọt lố, thêm Kd từ nhỏ. Thấy tín hiệu điều khiển nhảy loạn (do nhiễu cảm biến) thì giảm Kd hoặc lọc giá trị đo trước khi tính D.
4. Thử cả đổi điểm đặt và thêm tải (mở cửa lò, thêm vật nung), vì hai việc này có thể cần bộ số khác nhau.

Mỗi lần chỉ đổi một hệ số, khoảng 20–50 %, và ghi lại đồ thị (Serial Plotter) để so sánh. Với động cơ, hệ thống nhanh, trình tự giống hệt nhưng thời gian mỗi lần thử chỉ vài giây; ví dụ chỉnh PID vị trí động cơ có ở bài điều khiển vị trí động cơ DC cùng cụm chủ đề.

5. ESP32 tự dò thông số cho lò nhiệt

Chương trình đọc nhiệt độ bằng LM35 (10 mV/°C, GPIO34) và điều khiển điện trở nhiệt qua SSR (GPIO26). Lệnh t chạy thử rơ le quanh điểm đặt: đo 6 chu kỳ dao động, lấy 4 chu kỳ cuối để tính biên độ a, chu kỳ Tu, Ku. Từ đó chương trình đặt hệ số Ziegler–Nichols PI và chuyển sang chạy PID, bắt đầu từ công suất nền để không giật.

/*
  PYWORLD.VN - PID điều khiển nhiệt cho ESP32 có chế độ tự dò thông số bằng thử rơ le (Åström–Hägglund) rồi tính Ziegler–Nichols
  PYWORLD.VN - Cảm biến LM35 ở GPIO34 (10 mV/°C), SSR ở GPIO26 điều khiển theo tỉ lệ thời gian trong cửa sổ 2 giây
  PYWORLD.VN - Lệnh (115200, New Line): s 60 = điểm đặt 60 °C, t = tự dò, b 40 = công suất nền khi tự dò, p 2, i 0.05, d 10
*/

#include <math.h>

const int CHAN_CAM_BIEN = 34;
const int CHAN_SSR = 26;
const uint32_t CUA_SO_MS = 2000;            // SSR: bật một phần của mỗi 2 giây theo % công suất
const float DT = 1.0f;                      // chu kỳ PID, giây
const float TRE_ROLE = 0.2f;                // độ trễ (°C) của rơ le khi tự dò, chống đóng ngắt vì nhiễu
const float BIEN_DO_ROLE = 20.0f;           // % công suất cộng, trừ quanh mức nền khi tự dò
const int SO_CHU_KY_DO = 6;                 // số chu kỳ dao động khi tự dò (lấy 4 chu kỳ cuối)

float diemDat = 60.0f, nhietLoc = NAN;
float Kp = 2.0f, Ki = 0.05f, Kd = 0.0f;     // % trên °C; % trên °C·s; % trên °C/s
float tichPhan = 0, nhietTruoc = 0, congSuat = 0, congSuatNen = 40.0f;
bool dangTuDo = false;

// Số liệu thử rơ le
bool roleCao = true;
int soLanDoiLen = 0;
float dinhTam = -1e9, dayTam = 1e9, tongDinh = 0, tongDay = 0;
uint32_t lucDoiLen[SO_CHU_KY_DO + 1];

float docNhiet() {
  return analogReadMilliVolts(CHAN_CAM_BIEN) / 10.0f;   // LM35: 10 mV/°C
}

void batDauTuDo() {
  dangTuDo = true;
  roleCao = true;
  soLanDoiLen = 0;
  tongDinh = tongDay = 0;
  dinhTam = -1e9;
  dayTam = 1e9;
  Serial.printf("# Tu do: role %.0f%% +- %.0f%% quanh %.1f C\n", congSuatNen, BIEN_DO_ROLE, diemDat);
}

// Một bước thử rơ le; trả về công suất. Đủ số chu kỳ thì tính Ku, Tu và hệ số Ziegler–Nichols (PI cho lò nhiệt)
float buocTuDo(float nhiet) {
  dinhTam = fmaxf(dinhTam, nhiet);
  dayTam = fminf(dayTam, nhiet);
  if (roleCao && nhiet > diemDat + TRE_ROLE) {
    roleCao = false;
    if (soLanDoiLen > SO_CHU_KY_DO - 4) tongDinh += dinhTam;
    dayTam = 1e9;
  } else if (!roleCao && nhiet < diemDat - TRE_ROLE) {
    roleCao = true;
    if (soLanDoiLen > SO_CHU_KY_DO - 4) tongDay += dayTam;
    dinhTam = -1e9;
    lucDoiLen[soLanDoiLen++] = millis();
    if (soLanDoiLen > SO_CHU_KY_DO) {
      float a = (tongDinh / 4 - tongDay / 4) / 2;                 // biên độ dao động nhiệt độ
      float tu = (lucDoiLen[SO_CHU_KY_DO] - lucDoiLen[SO_CHU_KY_DO - 4]) / 4000.0f;
      float ku = 4 * BIEN_DO_ROLE / ((float)M_PI * sqrtf(fmaxf(a * a - TRE_ROLE * TRE_ROLE, 1e-6f)));
      Kp = 0.45f * ku;                                            // Ziegler–Nichols PI
      Ki = Kp / (tu / 1.2f);
      Kd = 0;
      tichPhan = (congSuatNen - 0) / Ki;                          // bắt đầu PID từ công suất nền, không giật
      dangTuDo = false;
      Serial.printf("# Ket qua: a %.2f C, Tu %.1f s, Ku %.3f -> Kp %.3f, Ki %.4f, Kd %.3f\n", a, tu, ku, Kp, Ki, Kd);
    }
  }
  return roleCao ? fminf(100, congSuatNen + BIEN_DO_ROLE) : fmaxf(0, congSuatNen - BIEN_DO_ROLE);
}

float buocPID(float nhiet) {
  float e = diemDat - nhiet;
  float u = Kp * e + Ki * tichPhan - Kd * (nhiet - nhietTruoc) / DT;   // D theo giá trị đo: không giật khi đổi điểm đặt
  nhietTruoc = nhiet;
  bool baoHoa = (u >= 100 && e > 0) || (u <= 0 && e < 0);
  if (!baoHoa) tichPhan += e * DT;
  return constrain(u, 0.0f, 100.0f);
}

void xuLyLenh() {
  String s = Serial.readStringUntil('\n');
  s.trim();
  if (s.length() < 1) return;
  char c = tolower(s.charAt(0));
  float so = s.substring(1).toFloat();
  if (c == 's') diemDat = so;
  else if (c == 't') batDauTuDo();
  else if (c == 'b') congSuatNen = constrain(so, 0.0f, 100.0f);
  else if (c == 'p') Kp = so;
  else if (c == 'i') { Ki = so; tichPhan = 0; }
  else if (c == 'd') Kd = so;
  Serial.printf("# dat %.1f C, Kp %.3f, Ki %.4f, Kd %.3f\n", diemDat, Kp, Ki, Kd);
}

void setup() {
  Serial.begin(115200);
  pinMode(CHAN_SSR, INPUT);                   // gói esp32 3.x: ghi sẵn mức tắt rồi mới chuyển sang ngõ ra
  digitalWrite(CHAN_SSR, LOW);
  pinMode(CHAN_SSR, OUTPUT);
  analogSetPinAttenuation(CHAN_CAM_BIEN, ADC_11db);
}

void loop() {
  static uint32_t lucLoc = 0, lucDK = 0, dauCuaSo = 0;
  uint32_t bayGio = millis();
  if (Serial.available()) xuLyLenh();

  if (bayGio - lucLoc >= 100) {                         // đọc 10 lần/giây, lọc EMA
    lucLoc = bayGio;
    float n = docNhiet();
    nhietLoc = isnan(nhietLoc) ? n : nhietLoc + 0.2f * (n - nhietLoc);
    if (dangTuDo) congSuat = buocTuDo(nhietLoc);        // rơ le phản ứng ngay theo số đo
  }
  if (bayGio - lucDK >= (uint32_t)(DT * 1000)) {
    lucDK = bayGio;
    if (!dangTuDo) congSuat = buocPID(nhietLoc);
    Serial.printf("dat:%.1f,nhiet:%.2f,cong_suat:%.1f\n", diemDat, nhietLoc, congSuat);
  }
  if (bayGio - dauCuaSo >= CUA_SO_MS) dauCuaSo += CUA_SO_MS;   // tỉ lệ thời gian: bật phần đầu của cửa sổ
  digitalWrite(CHAN_SSR, (bayGio - dauCuaSo) < congSuat / 100.0f * CUA_SO_MS ? HIGH : LOW);
}

• SSR đóng ngắt chậm, nên công suất xuất theo tỉ lệ thời gian: trong mỗi cửa sổ 2 giây, SSR bật một phần tương ứng với % công suất. Dùng SSR loại đóng ở điểm 0 (zero-cross) cho tải điện trở.
• Rơ le khi tự dò có độ trễ 0,2 °C để nhiễu cảm biến không làm SSR đóng ngắt liên tục; công thức Ku đã tính đến độ trễ này: Ku = 4d / (π √(a² − ε²)).
• Chọn công suất nền (b) gần mức giữ được điểm đặt: rơ le dao động cân đối, kết quả chính xác hơn. Lò thật lần đầu nên đưa nhiệt lên gần điểm đặt bằng công suất vừa phải rồi mới chạy t, vì trễ của lò làm nhiệt vọt khá xa (trong mô phỏng, từ điểm đặt +30 °C lên +56 °C).
• Phần điện 220 V (SSR, điện trở nhiệt, cầu chì, aptomat) do thợ điện lắp; nên có thêm bộ ngắt nhiệt độc lập phòng khi chương trình hoặc cảm biến hỏng.

6. Lỗi hay gặp

7. Xem thêm

• Giữ tốc độ động cơ DC ổn định bằng PID trên ESP32
• Mạch điều khiển bay của flycam: IMU, vòng PID, bộ lọc rung
• ESP32 DevKit V1 30 chân, cổng Type-C

Bài cùng cụm chủ đề: điều khiển vị trí động cơ DC bằng PID; lọc EMA cho cảm biến.