PID and Ziegler-Nichols Tuning

// A PID temperature loop, with the two things every textbook leaves out:
// integral anti-windup, and a derivative that does not scream at the sensor noise.
//
// Hardware: a 10 W power resistor bolted to a small aluminium block, driven by a logic-level
// MOSFET from a PWM pin; a 10k NTC thermistor in the block with a 10k series resistor to A0.
// The block is the plant, and it has exactly the two features that make control interesting:
// a long time constant (the block's heat capacity) and a DEAD TIME (heat takes seconds to
// reach the thermistor from the resistor).
//
// Set MODE to run the same rig as an on-off controller, a P, a PI or a PID, and read the
// four rows of the notebook's table off your own hardware.

const int PIN_HEATER = 9;      // PWM to the MOSFET gate
const int PIN_TEMP   = A0;     // thermistor divider

// 0 = on/off   1 = P only   2 = PI   3 = PID
const int MODE = 3;

// Ziegler-Nichols starts by finding these two on YOUR block. Raise KP with KI and KD at
// zero until the temperature oscillates steadily; that KP is Ku and the period is Tu.
const float Ku = 12.0;         // <-- measure this
const float Tu = 40.0;         // <-- and this, in seconds

const float SETPOINT = 45.0;   // degrees C - warm, safe, and far above the room
const float DT       = 0.5;    // seconds per loop
const float OUT_MAX  = 255.0;
const float N_FILTER = 8.0;    // derivative filter: limits the D gain to N times the P gain

float integral = 0.0, dFilt = 0.0, prevMeas = 0.0;
unsigned long tPrev = 0;

float readTempC() {
  // 10k NTC, beta 3950, in a divider with a 10k fixed resistor to 5 V.
  int raw = analogRead(PIN_TEMP);
  float r = 10000.0 * raw / (1023.0 - raw);
  float invT = 1.0/298.15 + log(r/10000.0)/3950.0;
  return 1.0/invT - 273.15;
}

void setup() {
  Serial.begin(115200);
  pinMode(PIN_HEATER, OUTPUT);
  Serial.println(F("t_s,temp_C,output,integral"));
  prevMeas = readTempC();
  tPrev = millis();
}

void loop() {
  if (millis() - tPrev < (unsigned long)(DT * 1000)) return;
  tPrev = millis();

  float meas  = readTempC();
  float error = SETPOINT - meas;

  // Ziegler-Nichols, ultimate-sensitivity method (1942).
  float Kp = 0.0, Ti = 1e9, Td = 0.0;
  if      (MODE == 1) { Kp = 0.50*Ku; }
  else if (MODE == 2) { Kp = 0.45*Ku; Ti = Tu/1.2; }
  else if (MODE == 3) { Kp = 0.60*Ku; Ti = Tu/2.0; Td = Tu/8.0; }

  float out;
  if (MODE == 0) {
    out = (meas < SETPOINT) ? OUT_MAX : 0.0;         // on/off, no gains at all
  } else {
    // Derivative on the MEASUREMENT, not the error: a setpoint change then does not
    // produce an infinite spike, and the sign works out the same for a disturbance.
    float dMeas = (meas - prevMeas) / DT;
    dFilt += (dMeas - dFilt) * (DT / (Td/N_FILTER + DT));
    float unclamped = Kp*error + Kp/Ti*integral - Kp*Td*dFilt;

    out = constrain(unclamped, 0.0, OUT_MAX);
    // ANTI-WINDUP: only accumulate while the actuator has somewhere to go. Without this
    // the integral fills up during the first slow warm-up and the block sails 15 degrees
    // past the setpoint before the term unwinds.
    if (unclamped == out) integral += error * DT;
  }
  prevMeas = meas;

  analogWrite(PIN_HEATER, (int)out);

  Serial.print(millis()/1000.0, 1); Serial.print(',');
  Serial.print(meas, 2);            Serial.print(',');
  Serial.print(out, 0);             Serial.print(',');
  Serial.println(integral, 1);
}