PID with simulated heater or motor

I wrote an example PID simulation for the PID_RT library. In the sketch code, it simulates a heater similar to the extruder block on a 3D printer so you can experiment with a PID with standardized, (simulated) hardware and see how it responds.

I'm not discussing how to tune the parameters, since that depends greatly on your system and is well covered elsewhere in places specific to your system.

I'm trying to demonstrate that:

  • The default kP=2, kI=5, kD=1 PID parameters do respond to a standard test case (a simulated 3d printer heater block).
  • That the Arduino can run code to simulate its own external system so you can do software development before the hardware is complete.

Since the code is in Rob Tillart's PID_RT repository as an example, it will show up in your Arduino IDE as File/Examples/PID_RT/PID_simulated_heater after you've installed the PID_RT library.

The code is in github and is copied here:

//
//    FILE: PID_simulated_heater.ino
//  AUTHOR: drf5n  (based on basic example)
// PURPOSE: demo
//
//  This simulates a 20W heater block driven by the PID
//  Wokwi https://wokwi.com/projects/356437164264235009
//
//  https://github.com/RobTillaart/PID_RT/issues/5


#include "PID_RT.h"

PID_RT PID;

const int PWM_PIN = 3;  // UNO PWM pin

int op = 0;
float input = 0;


void setup()
{
  Serial.begin(115200);
  Serial.println(__FILE__);

  PID.setPoint(125);
  PID.setOutputRange(0, 255);   //  PWM range
  PID.setInterval(50);
  PID.setK(2, 5, 1);
  PID.start();

  op = analogRead(A0);
}


void loop()
{
  float heaterWatts = ((int)op)*20.0/255;   //  20W heater
  float blockTemp = simPlant(heaterWatts);  //  simulate heating
  input = blockTemp;
  //  input = analogRead(A0);
  if (PID.compute(input))
  {
    op = PID.getOutput();
    analogWrite(PWM_PIN, op);

  //  Serial.print(PID.getInput());
  //  Serial.print('\t');
  //  Serial.println(op);
  }
  report();
}


void report(void){
  static uint32_t last = 0;
  const int interval = 1000;
  if (millis() - last > interval){
    last += interval;
//    Serial.print(millis()/1000.0);
    Serial.print(PID.getSetPoint());
    Serial.print(' ');
    Serial.print(input);
    Serial.print(' ');
    Serial.println(op);
  }
}

float simPlant(float Q){ // heat input in W (or J/s)
   //  simulate a 1x1x2cm aluminum block with a heater and passive ambient cooling
   //  next line "C is not used", fails in arduino-build-ci for ESP32
   //  float C = 237;  // W/mK thermal conduction coefficient for Al
   float h    = 5;    //  W/m2K thermal convection coefficient for Al passive
   float Cps  = 0.89; //  J/g°C
   float area = 1e-4; //  m2 area for convection
   float mass = 10;   // g
   float Tamb = 25;   // °C

   static float T       = Tamb;   //   °C
   static uint32_t last = 0;      //  last call
   uint32_t interval    = 100;    //  milliseconds

   if(millis() - last >= interval){
     last += interval;
     T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
   }
   return T;
}


//  -- END OF FILE --

The difference between other PID example is the code that simulates the process that the PID is controlling. The code that simulates the heater block is these lines and this function:

void loop()
{
...
  float heaterWatts = ((int)op)*20.0/255;   //  20W heater
  float blockTemp = simPlant(heaterWatts);  //  simulate heating
...
}
...
float simPlant(float Q){ // heat input in W (or J/s)
   //  simulate a 1x1x2cm aluminum block with a heater and passive ambient cooling
   //  next line "C is not used", fails in arduino-build-ci for ESP32
   //  float C = 237;  // W/mK thermal conduction coefficient for Al
   float h    = 5;    //  W/m2K thermal convection coefficient for Al passive
   float Cps  = 0.89; //  J/g°C
   float area = 1e-4; //  m2 area for convection
   float mass = 10;   // g
   float Tamb = 25;   // °C

   static float T       = Tamb;   //   °C
   static uint32_t last = 0;      //  last call
   uint32_t interval    = 100;    //  milliseconds

   if(millis() - last >= interval){
     last += interval;
     T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
   }
   return T;
}

...which takes the 0-255 output of the PID on PWM pin 3, converts it into the watts of heating that a MOSFET driver and cartridge heater would produce. It then takes those watts as input to a radiant heat transfer model to transform it into the temperature of the block, which the PID can take as its input.

There's a Wokwi simulation of it here:

...and an expansion on it here to give a user-adjustable setpoint:

//    FILE: PID_simulated_heater_pot.ino
//  AUTHOR: drf5n  (based on basic example)
// PURPOSE: demo
//
//  This simulates a 20W heater block driven by the PID
//  Vary the setpoint with the Pot, and watch the heater drive the temperature up
// 
//  Code at https://wokwi.com/projects/357374218559137793
//
//  Based on  
//  Wokwi https://wokwi.com/projects/356437164264235009
//
//  https://github.com/RobTillaart/PID_RT/issues/5
//

#include "PID_RT.h" // https://github.com/RobTillaart/PID_RT

PID_RT PID;  // https://github.com/RobTillaart/PID_RT

const int PWM_PIN = 3;  // UNO PWM pin
const int INPUT_PIN = -1; // Analog pin for Input (set <0 for simulation)
const int SETPOINT_PIN = A1;   // Analog pin for Setpoint Potentiometer
const int SETPOINT_INDICATOR = 6; // PWM pin for indicating setpoint
const int INPUT_INDICATOR = 5; // PWM pin for indicating Input

int op = 0;
float input = 0;


void setup()
{
  Serial.begin(115200);
  Serial.println(__FILE__);

  if (SETPOINT_INDICATOR >= 0) pinMode(SETPOINT_INDICATOR, OUTPUT);
  if (INPUT_INDICATOR >= 0) pinMode(INPUT_INDICATOR, OUTPUT);
  PID.setPoint(125);
  PID.setOutputRange(0, 255);  // PWM range
  PID.setInterval(50);
  PID.setK(2, 5, 1);
  PID.start();
  Serial.println("Setpoint Input Output");
  op = analogRead(A0);
}


void loop()
{
  float heaterWatts = ((int)op) * 20.0 / 255; // 20W heater
  if (INPUT_PIN > 0 ) {
    input = analogRead(INPUT_PIN);
  } else {
    float blockTemp = simPlant(heaterWatts); // simulate heating
    input = blockTemp;   // read input from simulated heater block
  }

  if (PID.compute(input))
  {
    op = PID.getOutput();
    analogWrite(PWM_PIN, op);

    //  Serial.print(PID.getInput());
    //  Serial.print('\t');
    //  Serial.println(op);
    PID.setPoint(analogRead(SETPOINT_PIN) / 4); // Read setpoint from potentiometer
    if (INPUT_INDICATOR >= 0) analogWrite(INPUT_INDICATOR, input);
    if (SETPOINT_INDICATOR >= 0) analogWrite(SETPOINT_INDICATOR, PID.getSetPoint());
  }
  report();
}


void report(void) {
  static uint32_t last = 0;
  const int interval = 1000;
  if (millis() - last > interval) {
    last += interval;
    //    Serial.print(millis()/1000.0);
    Serial.print(PID.getSetPoint());
    Serial.print(' ');
    Serial.print(input);
    Serial.print(' ');
    Serial.println(op);
  }
}

float simPlant(float Q) { // heat input in W (or J/s)
  // simulate a 1x1x2cm aluminum block with a heater and passive ambient cooling
  float C = 237; // W/mK thermal conduction coefficient for Al
  float h = 5 ; // W/m2K thermal convection coefficient for Al passive
  float Cps = 0.89; // J/g°C
  float area = 1e-4; // m2 area for convection
  float mass = 10 ; // g
  float Tamb = 25; // °C
  static float T = Tamb; // °C
  static uint32_t last = 0;
  uint32_t interval = 100; // ms

  if (millis() - last >= interval) {
    last += interval;
    T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
  }

  return T;
}


// -- END OF FILE --

In the Wokwi simulation, you can move the slide potentiometer up and down to adjust the SETPOINT of the PID controller. The graph icon in the lower right toggles between serial monitor or serial plotter mode, and you can see the SETPOINT, INPUT and OUTPUT of the PID as traces or as a table of values.

You could edit this line of code to adjust the PID parameters kP, kI, and kD to change the behavior.

void setup()
{
  ...
  PID.setK(2, 5, 1);
  ...
}

The units of the parameters are:

  • kP: PWM_output_counts/°Cerror
  • kI: PWM_output_counts/(°Cerror *sec)
  • kI: PWM_output_counts/(°Cerror / sec)

Where the PID controller uses these parameters to compute a 0-255 PWM setting to drive the (simulated) 20W cartridge heater.

The (2,5,1) setting is a very common set of PID parameters due to historical reasons: It is given in all sorts of examples. The Kp, Ki, and Kd parameters really should be customized to the system you are trying to control with a PID. The tuning, and even selection of the PID algorithm, is very dependent on the physical process you are trying to control with the PID.

A heater block is very different than a quadcopter, levitating magnet, or cruise control, but PID can help control any of them.

Here's another iteration, modeled after a cheap commercial PID with a display, a manual/auto switch, adjustments to the manual value, and the PID_v1 recommended reset method of PID.SetMode(MANUAL); Output=desired; PID.SetMode(AUTOMATIC).

/********************************************************
   PID Basic simulated heater Example
   Reading simulated analog input 0 to control analog PWM output 3
 ********************************************************/
//  This simulates a 20W heater block driven by the PID
//  Vary the setpoint with the Pot, and watch the heater drive the temperature up
//
//  Simulation at https://wokwi.com/projects/359088752027305985
//
//  Based on
//  Wokwi https://wokwi.com/projects/357374218559137793
//  Wokwi https://wokwi.com/projects/356437164264235009

#include "PID_v1.h" // https://github.com/br3ttb/Arduino-PID-Library
// local copy of .h and .cpp are tweaked to expose the integral per
// https://github.com/br3ttb/Arduino-PID-Library/pull/133
#define USE_HACK   // access the PID.outputSum variable

//Define Variables we'll be connecting to
double Setpoint, Input, Output;

//Specify the links and initial tuning parameters
//double Kp = 20, Ki = .01, Kd = 10; // works reasonably with sim heater block fo 220deg
double Kp = 25.5, Ki = 0.1, Kd = 0; // +/-10°proportional band 
//double Kp = 255, Ki = 0.05, Kd = 0; // works reasonably with sim heater block 
//double Kp = 255, Ki = .0, Kd = 0; // +/-1° proportional band works reasonably with sim heater block 
//double Kp = 10000, Ki = 0.0, Kd = 0.0; // bang-bang
//double Kp = 2, Ki = 0.0, Kd = 0.0; // P-only
//double Kp = 2, Ki = 5, Kd = 1; // commonly used defaults

PID myPID(&Input, &Output, &Setpoint, Kp, Ki, Kd, P_ON_E, DIRECT);

const int PWM_PIN = 3;  // UNO PWM pin for Output
const int INPUT_PIN = -1; // Analog pin for Input (set <0 for simulation)
const int SETPOINT_PIN = A1;   // Analog pin for Setpoint Potentiometer
const int AUTOMATIC_PIN = 8; // Pin for controlling manual/auto mode, NO
const int OVERRIDE_PIN = 12; // Pin for integral override, NO
const int PLUS_PIN = 4; // Pin for integral override, NO
const int MINUS_PIN =7; // Pin for integral override, NO


#include <LiquidCrystal_I2C.h>

#define I2C_ADDR    0x27
#define LCD_COLUMNS 20
#define LCD_LINES   4

LiquidCrystal_I2C lcd(I2C_ADDR, LCD_COLUMNS, LCD_LINES);

void setup()
{
  Serial.begin(115200);
  Serial.println(__FILE__);
  myPID.SetOutputLimits(0, 255); // -4 for 
  pinMode(OVERRIDE_PIN, INPUT_PULLUP);
  pinMode(AUTOMATIC_PIN, INPUT_PULLUP);
  pinMode(MINUS_PIN, INPUT_PULLUP);
  pinMode(PLUS_PIN, INPUT_PULLUP);
  Setpoint = 0;
  //turn the PID on
  myPID.SetMode(AUTOMATIC);
  if(INPUT_PIN>0){
    Input = analogRead(INPUT_PIN);
  }else{
    Input = simPlant(0.0,1.0); // simulate heating
  }
  lcd.init();
  lcd.backlight();

  Serial.println("Setpoint Input Output Integral");
}

void loop()
{
  // gather Input from INPUT_PIN or simulated block
  float heaterWatts = Output * 20.0 / 255; // 20W heater
  if (INPUT_PIN > 0 ) {
    Input = analogRead(INPUT_PIN);
  } else {
    float blockTemp = simPlant(heaterWatts,Output>0?1.0:1-Output); // simulate heating
    Input = blockTemp;   // read input from simulated heater block
  }

  if (myPID.Compute())
  {
    //Output = (int)Output; // Recognize that the output as used is integer
    analogWrite(PWM_PIN, Output);

  }

  Setpoint = analogRead(SETPOINT_PIN) / 4; // Read setpoint from potentiometer
  if(digitalRead(OVERRIDE_PIN)==LOW) mySetIntegral(&myPID,0); // integral override
  if(digitalRead(AUTOMATIC_PIN)==HIGH !=  myPID.GetMode()==AUTOMATIC){
    myPID.SetMode(digitalRead(AUTOMATIC_PIN)==HIGH ? AUTOMATIC :MANUAL);
  }
  static uint32_t lastButton = 0;
  if(myPID.GetMode()==MANUAL && millis() - lastButton > 250){
    if(digitalRead(PLUS_PIN)==LOW){ 
      Output += 1;
      lastButton = millis();
    }
    if(digitalRead(MINUS_PIN)==LOW){
      Output -= 1;
      lastButton = millis();
    }
  }

  report();
  reportLCD();

}

void report(void)
{
  static uint32_t last = 0;
  const int interval = 250;
  if (millis() - last > interval) {
    last += interval;
    //    Serial.print(millis()/1000.0);
    Serial.print("SP:");Serial.print(Setpoint);
    Serial.print(" PV:");
    Serial.print(Input);
    Serial.print(" CV:");
    Serial.print(Output);
    Serial.print(" Int:");
#if defined(USE_HACK)
    Serial.print(myPID.outputSum);
#endif
    Serial.print(' ');
    Serial.println();
  }
}

void reportLCD(void)
{
  static uint32_t last = 0;
  const int interval = 250;
  if (millis() - last > interval) {
    last += interval;
    //    Serial.print(millis()/1000.0);
   // lcd.clear();
    lcd.setCursor(0,0);
    lcd.print("PV:");
    lcd.print(Input,3);
    lcd.print(" CV:");
    lcd.print(Output,3);
    lcd.print("  ");
    lcd.setCursor(0,1);
    lcd.print("SP:");
    lcd.print(Setpoint,3);
    lcd.print(myPID.GetMode()==AUTOMATIC? " Automatic ":" Manual   ");
    lcd.print("  ");
    lcd.setCursor(0,3);
    lcd.print("Int:");
#if defined(USE_HACK)
    lcd.print(myPID.outputSum,4);
#endif
    lcd.print(' ');
    lcd.println();
  }
}

float simPlant(float Q,float hfactor) { // heat input in W (or J/s)
  // simulate a 1x1x2cm aluminum block with a heater and passive ambient cooling
 // float C = 237; // W/mK thermal conduction coefficient for Al
  float h = 5 *hfactor ; // W/m2K thermal convection coefficient for Al passive
  float Cps = 0.89; // J/g°C
  float area = 1e-4; // m2 area for convection
  float mass = 10 ; // g
  float Tamb = 25; // °C
  static float T = Tamb; // °C
  static uint32_t last = 0;
  uint32_t interval = 100; // ms

  if (millis() - last >= interval) {
    last += interval;
    // 0-dimensional heat transfer
    T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
  }
  return T;
}

void  mySetIntegral(PID * ptrPID,double value ){
   ptrPID->SetMode(MANUAL);
   Output = value;
   ptrPID->SetMode(AUTOMATIC);
}

By modifying the code to use different Kp, Ki, Kd variables, you can experiment with their impacts and behavior. For instance:

  • The double Kp = 2, Ki = 5, Kd = 1; that is commonly used as a default is a poor setting for this system. Producing integral windup and overshoot. This is a good demonstration on the need for proper tuning.
  • double Kp = 2, Ki = 0.0, Kd = 0.0; // P-only will approach the setpoint slowly, but will stop before it gets there since it would turn the output down to zero Kp*error to produce any output.
  • double Kp = 10000, Ki = 0.0, Kd = 0.0; With a very large Kp, you can make a PID operate like a bang-bang controller, since a small error will saturate the output. Rapidly approaching the setpoint, and quickly turning off.
  • double Kp = 255, Ki = .0, Kd = 0 produces a "proportional zone" that is +/-1° wide, which works pretty well for this zero-dimensional heater simulation.

There is a slight difference in PID_v1 library the Wokwi code uses. This wokwi sim uses modified local copies of the PID_v1.cpp and PID_v1.h file, which move the PID::outputSum variable from private to public, exposing the integration term to user-space monitoring and adjustment. This enables monitoting the internal integral on the LCD screen. Aside from the exposure of the the PID_v1::privates, it works exactly the same as the stock library.

This simulated 2cm^3 aluminum heater block is fairly easy to tune, because it uses a zero-dimensional heat transfer model, which has no dead time and instantly distributes the heat throughout the block. A more rigorous model would take some time to conduct heat throughout the heater block, but the differences wouldn't be large.

I think this small change to the library makes it easy to see what the PID is doing, and makes it possible to adapt the library to more complicated PID schemes in user-space. I put in this pull request to the PID_v1 library:

1 Like

Yep, that is utterly brilliant. I'm going to spend some proper time reading and understanding that. Thank you!

1 Like

Here's 1-D heat transfer simulation, attached to a couple works-in-process

//1D heat transfer simulation through aluminum rod, heated on one end, 
// https://wokwi.com/projects/359193529297155073
// for https://forum.arduino.cc/t/pid-with-simulated-heater/1093539

#include <BasicLinearAlgebra.h> // https://github.com/tomstewart89/BasicLinearAlgebra
/*
  This models a 1-D heat transfer problem of flow through with n_ele elements

        convection--v                                          v- convection at ends to Tamb
     heater q_in -> [###round aluminum rod(rod_dia, rod_len)###]

  Code at https://wokwi.com/projects/359193529297155073

*/

// All the functions in BasicLinearAlgebra are wrapped up inside the namespace BLA, so specify that we're using it like
// so:
using namespace BLA;

void setup()
{
  Serial.begin(115200);
}

// https://wokwi.com/projects/359235477558390785
// 0-D heat transfer through 1x1x2sm aluminum block
//  T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
float h = 5; // W/m2K thermal convection coefficient for Al passive
float Cps = 0.89; // J/g°C
float C_tcond = 237; // W/mK thermal conduction coefficient for Al
float Tmelt = 660; // °C
float rho = 2.710e6; // kg/m^3 density // Al: 2,710kg/m3
float rod_dia = 0.02; // m
float rod_len = 0.05; //m
float area = pow(rod_dia, 2) * PI / 4 + 1 * rod_len / 3 * PI * rod_dia; // m2 area for convection
float mass = rod_len * pow(rod_dia, 2) * PI / 4 * rho ; // g
float Tamb = 25; // °C


const int n_ele = 4;
const int n_force = 2;
// setup for three-element 1D heat transfer model
BLA::Matrix<n_ele> x; //initial
BLA::Matrix<n_force> u; // Forcings at ends Watts
BLA::Matrix<n_ele, 2> G; // map from forcings to elements, dT=f(W)

//
BLA::Matrix<n_ele, n_ele> Fa ; // transition matrix

float dt = 0.05;
int dtPrint = 500;

void loop() {
  // Or as a more real life example, here's how you might calculate an updated state estimate for a third order state
  // space model:
  static float t = 0;
  static uint32_t last = 0;
  Tamb = analogRead(A1) * 250.0 / 1024;
  float q = analogRead(A0) * 10.0 / 1024;  // 10W heater
  static bool oneShot = true;
  if (oneShot == true) {
    oneShot = false;
    for (int i = 0; i < n_ele; ++i) {
      x(i) = Tamb;
      //
      if (i > 0) { // conduction to the left
        Fa(i, i - 1) = C_tcond * pow(rod_dia, 2) * PI / 4 / (rod_len / n_ele);
        Fa(i, i) -= Fa(i, i - 1);
      }
      if (i < n_ele - 1) { // conduction to the right
        Fa(i, i + 1) = C_tcond * pow(rod_dia, 2) * PI / 4 / (rod_len / n_ele);
        Fa(i, i) -= Fa(i, i + 1);
      }
    }
    Fa = Fa * 1; // *20 extra conductive
    G(0, 0) = Cps;
    G(n_ele - 1, 1) = Cps;
    u(0) = 0;
    u(n_force - 1) = 0;
    Serial <<"F:"<< Fa << "\n";
    Serial <<"G:" << G << "\n";
    Serial.println("Oneshot initialization finished");
  }
  if (millis() - last >= dtPrint) {
    last += dtPrint;
    Serial << "t:" << t << " Tamb=" << Tamb << " q=" << q << " u=" << u << " x: " << x << "\n";
  }
  // convection + heater
  u(0) = ((Tamb - x(0)) * area * h) + q;
  // convection
  u(1) = (Tamb - x(3)) * area * h;
  x += (Fa * x + G * u) * dt ;
  t += dt;

  delay(1000 * dt); // realtime
}

The 1-d model uses a matrix algebra library to model the temperatures and heat transfer between 5 adjacent elements in a simulated aluminum rod. The heat applied at one end through a 10W heater propagates towards the other end, and both ends convect heat into the ambient environment. You can adjust the watts of heat injected into one end, and the ambient temperature. The output lines show the state of the system:

t:139.35 Tamb=0.00 q=0.00 u=[[-0.01],[-0.01]] x: [[2.18],[2.18],[2.18],[2.18]]

Where:

  • t = simulation time in seconds
  • Tamb = ambient temperature
  • q = heat input (watts)
  • u[...] = heat flow into and out of the two rod ends (watts)
  • x[...] = temperature of the 5 elements (°C)

You can see that x[4] element lags the x[0] element as heat is added to the system.

Applying the PID to the 1D model

That 1-D model is connected to the PID setup in the earlier post. You can change which element the temperature sensor is attached to with the ele_sensor variable in the code. In this screenshot, you can see the dynamic spread of temperatures in the different elements by the family of white lines shown in the serial plotter. At a few seconds into the simulation, you can see a spread of 20°C between the x_i0 and x_i4 elements.

/********************************************************
   PID simulated heater Example
   Reading simulated analog input 0 to control analog PWM output 3
 ********************************************************/
//  This simulates a 5W heater block driven by the PID
//  Vary the setpoint with the Pot, and watch the heater drive the temperature up
//
//  Simulation at https://wokwi.com/projects/359088752027305985
//
//  Based on
//  Wokwi https://wokwi.com/projects/357374218559137793
//  Wokwi https://wokwi.com/projects/356437164264235009

#include "PID_v1.h" // https://github.com/br3ttb/Arduino-PID-Library
// local copy of .h and .cpp are tweaked to expose the integral per
// https://github.com/br3ttb/Arduino-PID-Library/pull/133
#define USE_HACK   // access the PID.outputSum variable
// For 1D heat transfer model:
#include <BasicLinearAlgebra.h> // https://github.com/tomstewart89/BasicLinearAlgebra


//Define Variables we'll be connecting to
double Setpoint, Input, Output;

//Specify the links and initial tuning parameters
//double Kp = 20, Ki = .01, Kd = 10; // works reasonably with sim heater block fo 220deg
double Kp = 25.5, Ki = 0.001, Kd = 0; // +/-10°proportional band
//double Kp = 255, Ki = 0.001, Kd = 0; // works reasonably with sim heater block
//double Kp = 255, Ki = .0, Kd = 0; // +/-1° proportional band works reasonably with sim heater block
//double Kp = 10000, Ki = 0.0, Kd = 0.0; // bang-bang
//double Kp = 2, Ki = 0.0, Kd = 0.0; // P-only
//double Kp = 2, Ki = 5, Kd = 1; // commonly used defaults

PID myPID(&Input, &Output, &Setpoint, Kp, Ki, Kd, P_ON_E, DIRECT);

// expose rod temperature vector to other routines
const int n_ele = 5;
const int sensor_element = 2; // 0 for sensor close to heater, easyist
                              // n_ele -1 for furthest, most lag, hardest
static BLA::Matrix<n_ele> x; //initial

const int PWM_PIN = 3;  // UNO PWM pin for Output
const int FAN_PIN = 5;  // UNO PWM pin for Fan Output

const int INPUT_PIN = -1; // Analog pin for Input (set <0 for simulation)
const int SETPOINT_PIN = A1;   // Analog pin for Setpoint Potentiometer
const int AUTOMATIC_PIN = 8; // Pin for controlling manual/auto mode, NO
const int OVERRIDE_PIN = 12; // Pin for integral override, NO
const int PLUS_PIN = 4; // Pin for integral override, NO
const int MINUS_PIN = 7; // Pin for integral override, NO


#include <LiquidCrystal_I2C.h>

#define I2C_ADDR    0x27
#define LCD_COLUMNS 20
#define LCD_LINES   4

LiquidCrystal_I2C lcd(I2C_ADDR, LCD_COLUMNS, LCD_LINES);

void setup()
{
  Serial.begin(115200);
  Serial.println(__FILE__);
  myPID.SetOutputLimits(-4, 255); // -4 for fan to increase cooling
  pinMode(OVERRIDE_PIN, INPUT_PULLUP);
  pinMode(AUTOMATIC_PIN, INPUT_PULLUP);
  pinMode(MINUS_PIN, INPUT_PULLUP);
  pinMode(PLUS_PIN, INPUT_PULLUP);
  pinMode(PWM_PIN,OUTPUT);
  pinMode(FAN_PIN,OUTPUT);
  Setpoint = 0;
  //turn the PID on
  myPID.SetMode(AUTOMATIC);
  if (INPUT_PIN > 0) {
    Input = analogRead(INPUT_PIN);
  } else {
    Input = simPlant1D(0.0, 1.0); // simulate heating
  }
  lcd.init();
  lcd.backlight();

  Serial.println("Setpoint Input Output Integral");
}

void loop()
{
  // gather Input from INPUT_PIN or simulated block
  float heaterWatts = Output * 5.0 / 255; // 5W heater
  if (INPUT_PIN > 0 ) {
    Input = analogRead(INPUT_PIN);
  } else {
    float blockTemp = simPlant1D(heaterWatts, Output > 0 ? 1 : 1 - Output); // simulate heating
    Input = blockTemp;   // read input from simulated heater block
  }

  if (myPID.Compute() || myPID.GetMode()==MANUAL)
  {
    //Output = (int)Output; // Recognize that the output as used is integer
    analogWrite(PWM_PIN, Output >= 0 ? Output : 0); // Heater
    analogWrite(FAN_PIN, Output < 0 ? -Output*255/4.0:0); // Fan

  }

  Setpoint = analogRead(SETPOINT_PIN) / 4; // Read setpoint from potentiometer
  if (digitalRead(OVERRIDE_PIN) == LOW) mySetIntegral(&myPID, 0); // integral override
  if (digitalRead(AUTOMATIC_PIN) == HIGH !=  myPID.GetMode() == AUTOMATIC) {
    myPID.SetMode(digitalRead(AUTOMATIC_PIN) == HIGH ? AUTOMATIC : MANUAL);
  }
  static uint32_t lastButton = 0;
  if (myPID.GetMode() == MANUAL && millis() - lastButton > 250) {
    float incValue = pow(10.0, map(analogRead(A2),0,1023,-3,1));
    if (digitalRead(PLUS_PIN) == LOW) {
      Output += incValue;
      lastButton = millis();
    }
    if (digitalRead(MINUS_PIN) == LOW) {
      Output -= incValue;
      lastButton = millis();
    }
  }

  report();
  reportLCD();

}

void report(void)
{
  static uint32_t last = 0;
  const int interval = 250;
  if (millis() - last > interval) {
    last += interval;
    //    Serial.print(millis()/1000.0);
    Serial.print("SP:"); Serial.print(Setpoint);
    Serial.print(" PV:");
    Serial.print(Input);
    Serial.print(" CV:");
    Serial.print(Output);
    Serial.print(" Int:");
#if defined(USE_HACK)
    Serial.print(myPID.outputSum);
#endif
    Serial.print(' ');
    for (int i = 0; i < n_ele; ++i) {
      Serial.print("x_i");
      Serial.print(i);
      Serial.print(':');
      Serial.print(x(i), 3);
      Serial.print(' ');
    }

    Serial.println();
  }
}

void reportLCD(void)
{
  static uint32_t last = 0;
  const int interval = 250;
  if (millis() - last > interval) {
    last += interval;
    //    Serial.print(millis()/1000.0);
    // lcd.clear();
    lcd.setCursor(0, 0);
    lcd.print("PV:");
    lcd.print(Input, 3);
    lcd.print(" CV:");
    lcd.print(Output, 3);
    lcd.print("  ");
    lcd.setCursor(0, 1);
    lcd.print("SP:");
    lcd.print(Setpoint, 3);
    lcd.print(myPID.GetMode() == AUTOMATIC ? " Automatic " : " Manual   ");
    lcd.print("  ");
    lcd.setCursor(0, 3);
    lcd.print("Int:");
#if defined(USE_HACK)
    lcd.print(myPID.outputSum, 4);
#endif
    lcd.print(' ');
    lcd.println();
  }
}



float simPlant1D(float Q, float hfactor) { // heat input in W (or J/s)
  // simulate a 1x1x2cm aluminum block with a heater and passive ambient cooling
  // float C = 237; // W/mK thermal conduction coefficient for Al
  // Using a 1D heat diffusion rod model

  // Still 0D  See https://wokwi.com/projects/359193529297155073 for
  // what needs to go in here
  using namespace BLA;

  const int n_force = 2;
  // setup for three-element 1D heat transfer model
  // expose temperature vector to other routines
  // move these to global if you want other routines to see them:
 // const int n_ele = 5;
 // static BLA::Matrix<n_ele> x; //initial
  static BLA::Matrix<n_force> u; // Forcings at ends Watts
  static BLA::Matrix<n_ele, 2> G;  // map from forcings to elements, dT=f(W)
  static BLA::Matrix<n_ele, n_ele> Fa ; // transition matrix


  float h = 5  ; // W/m2K thermal convection coefficient for Al passive
  float Cps = 0.89; // J/g°C
  float C_tcond = 237; // W/mK thermal conduction coefficient for Al
  float Tamb = 25; // °C
  float rho = 2.710e6; // kg/m^3 density // Al: 2,710kg/m3
  float rod_dia = 0.02; // m
  float rod_len = 0.10; //m
  float area = pow(rod_dia, 2) * PI / 4 + 0 * rod_len / 3 * PI * rod_dia; // m2 area for convection
  float mass = rod_len * pow(rod_dia, 2) * PI / 4 * rho ; // g

  static float T = Tamb; // °C
  static uint32_t last = 0;
  uint32_t interval = 10; // ms
  float dt = interval / 1000.0;
  static bool oneShot = true;
  if (millis() - last >= interval) {
    last += interval;

    if (oneShot == true) {
      oneShot = false;
      for (int i = 0; i < n_ele; ++i) {
        x(i) = Tamb;
        if (i > 0) { // conduction to the left
          Fa(i, i - 1) = C_tcond * pow(rod_dia, 2) * PI / 4 / (rod_len / n_ele)/(mass/n_ele);
          Fa(i, i) -= Fa(i, i - 1);
        }
        if (i < n_ele - 1) { // conduction to the right
          Fa(i, i + 1) = C_tcond * pow(rod_dia, 2) * PI / 4 / (rod_len / n_ele)/(mass/n_ele);
          Fa(i, i) -= Fa(i, i + 1);
        }
        Fa = Fa * 1; // *20 extra conductive
        G(0, 0) = Cps;
        G(n_ele - 1, 1) = Cps;
        u(0) = 0;
        u(n_force - 1) = 0;
        Serial << "Fa:" << Fa <<'\n';
        Serial << "G:" << G <<'\n';
        Serial.print("Oneshot");
      }
    } // end oneShot
    // 0-dimensional heat transfer
    //T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
    // 1-D explicit:
    u(0) = ((Tamb - x(0)) * area * h * hfactor) + Q;
    u(1) = (Tamb - x(2)) * area * h * hfactor;

    x += (Fa * x + G * u) * dt ;

    //T = x(n_ele - 1); // end element (most lag) Harder to control
    T = x(sensor_element); // element with heater (some lag)
    //T = x(0); // element with heater (least lag)
  }
  return T;

}


void  mySetIntegral(PID * ptrPID, double value ) {
  // per
  ptrPID->SetMode(MANUAL);
  Output = value;
  ptrPID->SetMode(AUTOMATIC);
}


The tuning parameters depend on the closeness of the coupling between the heater and the temperature sensor, with the system being much easier to control when the temperature is measured in ele_sensor=0, as compared to ele_sensor=4.

Here's a PID with a simulated DC motor to make a DIY servo:

  • The top slider sets the setpoint of the process at 0-100mm
  • The simulated motor drives a simulated linear actuator from 0-100mm, with the simulated process measurements and motor variables measurements printed in the serial monitor
  • The LEDs show the PWM signal and the direction signal that is being fed into the simulated motor.

In this case, the simulated process is some of the physics of a DC motor, and the rotation angle of the motor is translated into the motion of a linear actuator.

// Work in progress -- bidirectional DC motor as a servo

// simulate a PWM driven motor controlled by potentiometer
// DaveX 2022-02-21 CC BY-SA
// https://github.com/drf5n/drf5nArduino/tree/main/dc_motor_sim_servo
// This simulates a PID+PWM controlled DC motor with the simulated motor
// control pin 9 OC1A PWM read back from the OCR1A register.
//
// This simulation at https://wokwi.com/projects/362627481524244481
// uses these Uno Connections:
//
//   input: A0: 0-5V potentiometer on A0 sets distance 0-100
//   simulated input: 0-1024 analog distance
//   output: digitalPin9: PWM value to motor
//           digitalPin10 : direction to motor
//   Serial output: of PID state and Motor state
//
// Motor simulated by a state space model in class PmMotor
// This sketch aims at Servo control of a DC motor
// position
//
// https://wokwi.com/arduino/projects/323970337959051859 does speed control
//
// Discussions at:
// https://forum.arduino.cc/t/what-is-the-inverse-of-analogwrite-pin-val/960774/20
// https://forum.arduino.cc/t/why-doesnt-analogwrite-protect-its-writes-to-16bit-registers/961470
// https://forum.arduino.cc/t/motor-rpm-measurement-using-optical-encoder-for-pi-control/960163/13?u=davex
// https://forum.arduino.cc/t/controlling-a-linear-actuator-using-a-pid-controller/1117777/15

const byte motorPwmPin = 9; 
const byte motorDirPin = 10;

#include <PID_v1.h> // https://playground.arduino.cc/Code/PIDLibrary/

class PmMotor { // Permanent Magnet Motor Simulation
  public:
    // DC PM motor state space model per
    // https://ctms.engin.umich.edu/CTMS/index.php?example=MotorSpeed&section=SystemModeling
    // DaveX 2023-04-21
    // input vars
    float V = 5; //  voltage (controlled by pwm fraction)
    int pwm_level; // 0-255 pwm controlled voltage
    // state vars
    float theta = 0;    // rotational position (rads)
    float theta_dot = 0; // rad/sec
    float i_amps = 0 ; // coil current amps
    // outputs
    float theta_dotdot = 0; // rad/sec/sec
    float di_dt;
    // extras;
    float MaxRPM = 1000;
    //constants
    float J = 0.001; // inertia kg*m2
    float B = 0.0001 ; // fluid friction N*m*s
    float L = 0.05; // motor inductance H
    float Ke = 0.001; // speed constant V/rad/sec
    float Kt = 1; // torqe constant N*m/A
    float R = 1.0; // motor resistance ohms
    float i_full = 0.01; // no load current at full speed voltage
    // governing eqns:
    // J*theta_dotdot + B * theta_dot = Kt *i
    // L * di/dt +R*i = V -Ke*theta_dot

    // time variables
    unsigned long updateInterval = 50; // ms
    unsigned long lastUpdate = 0;

    float speed_mrev_sec = 0; //theta_dot
    long millirevs = 0; // theta in millirevolutions
    long encoderTicks = 0;
    const long ticksPerRev = 1000;
    float rpm;

    void setKeMaxRPMVR(float rpm) {
      // set Ke from....
      // assume V voltage, for pwm/255 fraction
      // assume di/dt =0
      // theta_dot_max = rpm/60*2*pi // to rad/sec
      // V*255/255 = R*i  +Ke*theta_dot
      //
      Ke = (V - 0.0 * i_full * R) / (rpm / 60 * 2 * PI);
      MaxRPM = rpm;
    }
    void setInertia(float x) {
      J = x;
    }
    void setResistance(float x) {
      R = x;
    }

    void setKtFromIB() {
      // From J*theta_dotdot + B * theta_dot = Kt *i
      // assume theta_dotdot = 0, full speed, i
      const float theta_dot_max = MaxRPM / 60 * 2 * PI ;
      float myKt = B * theta_dot_max / i_full;

      Serial.print(Kt);
      Serial.print("->");
      Serial.print(myKt);
      Serial.println();

    }
    void update(uint16_t pwm, int dir) {
      // pwm 0-255, as from an analogWrite()
      // dir in HIGH,LOW as from a digitalwrite)
      // One also could change direction or drive voltage with motor.V 
      pwm_level = pwm * (dir == HIGH? 1 : -1);
      update();
    }

    void update() {
      unsigned long now = millis();
      static float i_last = 0;
      if (now - lastUpdate >= updateInterval) {
        lastUpdate += updateInterval;

        i_amps = (V * pwm_level / 255.0 - Ke * theta_dot - L * di_dt) / R;
        di_dt = i_amps - i_last;

        theta += theta_dot * updateInterval / 1000.0;
        theta_dotdot = (Kt * i_amps - B * theta_dot) / J;
        theta_dot += theta_dotdot * updateInterval / 1000.0;
        encoderTicks += theta_dot * updateInterval / 1000.0 / 2 / PI * ticksPerRev;

      }
    }
    void printConfig(void) { // print the configuration vars
      Serial.print("Motor Configuration");
      Serial.print(" MaxRPM=");
      Serial.print(MaxRPM);
      Serial.print(" Ke=");
      Serial.print(Ke,4);
      Serial.print("V/rad/sec, Kt=");
      Serial.print(Kt,4);
      Serial.print("N*m/A B=");
      Serial.print(B,4);
      Serial.print("N*m*s, J=");
      Serial.print(J,4);
      Serial.print("kg*m2, L=");
      Serial.print(L,4);
      Serial.print("H V=");
      Serial.print(V,4);
      Serial.print("V pwm_level=");
      Serial.print(pwm_level);
      Serial.println('\n');
    }
    void printState(void) { // Print the state variables
      Serial.print(" pwm=");
      Serial.print(pwm_level);
      Serial.print(" i_amps=");
      Serial.print(i_amps);
      Serial.print("A theta=");
      Serial.print(theta);
      Serial.print("rad theta_dot=");
      Serial.print(theta_dot);
      Serial.print("rad/s theta_dotdot=");
      Serial.print(theta_dotdot);

      Serial.println("rad/s2");
    }

    void printStateRPM(void) { //print the state in rev & RPM form
      Serial.print(" Motor: pwm=");
      Serial.print(pwm_level);
      Serial.print(" i=");
      Serial.print(i_amps);
      Serial.print("A ");
      Serial.print(theta / 2 / PI);
      Serial.print("rev ");
      Serial.print(theta_dot / 2 / PI * 60, 5);
      Serial.print("RPM ");
      Serial.print(theta_dotdot / 2 / PI * 60, 5);
      Serial.print("RPM/s ");

      Serial.println();
    }
};

PmMotor myMotor;

double myRPM = 0.0;
double mySetpointDisp = 0.0;
double myCV = 0.0;

//PID motorPID(&myRPM, &myCV, &mySetpointRPM,1,0.1,0, DIRECT);

double myDisplacement = 0;
PID motorPID(&myDisplacement, &myCV, &mySetpointDisp, 7, 0.1, 0, DIRECT);

void report (void) {
  const int interval = 500;
  static unsigned long last = - interval;
  unsigned long now = millis();
  if (now - last >= interval) {
    last += interval;
    //motor.printState();
    Serial.print("PID: SP: ");
    Serial.print(mySetpointDisp);
    Serial.print("mm MV: ");
    Serial.print(myDisplacement);
    Serial.print("mm CV: ");
    Serial.print(myCV);
    Serial.print(" counts PWM9: ");
    //    Serial.print(OCR1A,DEC);
    Serial.print(analogGetPwm9(), DEC);
    Serial.println();
    myMotor.printStateRPM();
    //motor.printConfig();
  }
}

void setup() {
  // put your setup code here, to run once:
  Serial.begin(115200);
  myMotor.V = 5;
  myMotor.setResistance(1.2);
  myMotor.setKeMaxRPMVR(30);
  myMotor.i_full = 0.1; // Amps at no-load full speed
  myMotor.B = 0.0001e-0; // friction
  myMotor.setInertia(.1);
  myMotor.setKtFromIB();
  myMotor.printConfig();
  myMotor.printState();

  myRPM = myMotor.theta_dot * 60 / 2 / PI;
  mySetpointDisp = analogRead(A0); // 0-1024 rpm
  //turn the PID on
  motorPID.SetSampleTime(10); //us
  motorPID.SetOutputLimits(-254, 255);
  motorPID.SetMode(AUTOMATIC);
  pinMode(motorPwmPin,OUTPUT);
  pinMode(motorDirPin,OUTPUT);
  
  delay(500);
}

int analogGetPwm9() { // PWM 9 on UNO is OCR1A
  // Retrieve PWM value from pin per https://forum.arduino.cc/t/what-is-the-inverse-of-analogwrite-pin-val/960774/20
  // Uno pin 9 is OC1A per
  // Uno https://docs.arduino.cc/hacking/hardware/PinMapping168
  // Mega https://docs.arduino.cc/hacking/hardware/PinMapping2560

  return ((TCCR1A & bit(COM1A1)) ? OCR1A : digitalRead(9) * 255);
}

void loop() {
  unsigned long loopNow = millis();
  static long bresenhamSlope = 0; // phase error in ticks

  mySetpointDisp = (analogRead(A0) -0) / 10.23; // setpoint of mm based on slider


  myDisplacement = myMotor.theta * 60 / 2 / PI; // measure the motor speed
  //myPhase = bresenhamSlope;  // phase error

  motorPID.Compute();                 // PID to calculate CV from SP and MV
  // control the motor
  analogWrite(motorPwmPin, abs(myCV)); // push PID control value out PWM pin 9
  digitalWrite(motorDirPin,myCV>0? HIGH:LOW);

  // Simulate a motor attached to a PWM pin
  myMotor.update(analogGetPwm9(),digitalRead(motorDirPin));
  report();
}

The tuning is pretty ad-hoc, but it seems to work OK for this simulated servo.

This simulation does DC motor speed, either PID or manual PWM control:

// simulate a PWM driven motor controlled by potentiometer
// DaveX 2022-02-21 CC BY-SA
// https://github.com/drf5n/drf5nArduino/tree/main/dc_motor_sim
// Wokwi: https://wokwi.com/projects/323970337959051859
// This simulates a PID+PWM controlled DC motor with the simulated motor
// control read back from the pin 9 OC1A PWM read back from the OCR1A register
//
// input: 0-5V potentiometer on A0 sets RPM 0-1023
// output PWM value to motor
//
// motor simulated by a state space motor in class PmMotor
//
// This sketch does PID control of speed of rotation
// for PID of a PLL on angular position see https://wokwi.com/arduino/projects/324274590349001300
//
// Discussions at:
// https://forum.arduino.cc/t/what-is-the-inverse-of-analogwrite-pin-val/960774/20
// https://forum.arduino.cc/t/why-doesnt-analogwrite-protect-its-writes-to-16bit-registers/961470
// https://forum.arduino.cc/t/motor-rpm-measurement-using-optical-encoder-for-pi-control/960163/13?u=davex


#include <PID_v1.h> // https://playground.arduino.cc/Code/PIDLibrary/

class PmMotor { // Permanent Magnet Motor
  public:
    // DC PM motor state space model per
    // https://ctms.engin.umich.edu/CTMS/index.php?example=MotorSpeed&section=SystemModeling
    // input vars
    float V = 5; //  voltage (controlled by pwm fraction)
    uint8_t pwm_level; // 0-255 pwm controlled voltage
    // state vars
    float theta = 0;    // rotational position (rads)
    float theta_dot = 0; // rad/sec
    float i_amps = 0 ; // coil current amps
    // outputs
    float theta_dotdot = 0; // rad/sec/sec
    float di_dt;
    // extras;
    float MaxRPM = 1000;
    //constants
    float J = 0.001; // inertia kg*m2
    float B = 0.0001 ; // fluid friction N*m*s
    float L = 0.05; // motor inductance H
    float Ke = 0.001; // speed constant V/rad/sec
    float Kt = 1; // torqe constant N*m/A
    float R = 1.0; // motor resistance ohms
    float i_full = 0.01; // no load current at full speed voltage
    // governing eqns:
    // J*theta_dotdot + B * theta_dot = Kt *i
    // L * di/dt +R*i = V -Ke*theta_dot

    // time variables
    unsigned long updateInterval = 50; // ms
    unsigned long lastUpdate = 0;

    float speed_mrev_sec = 0; //theta_dot
    long millirevs = 0; // theta in millirevolutions
    float rpm;

    void setKeMaxRPMVR(float rpm) {
      // set Ke from....
      // assume V voltage, for pwm/255 fraction
      // assume di/dt =0
      // theta_dot_max = rpm/60*2*pi // to rad/sec
      // V*255/255 = R*i  +Ke*theta_dot
      //
      Ke = (V - 0.0 * i_full * R) / (rpm / 60 * 2 * PI);
      MaxRPM = rpm;
    }
    void setInertia(float x) {
      J = x;
    }
    void setResistance(float x) {
      R = x;
    }

    void setKtFromIB() {
      // From J*theta_dotdot + B * theta_dot = Kt *i
      // assume theta_dotdot = 0, full speed, i
      const float theta_dot_max = MaxRPM / 60 * 2 * PI ;
      float myKt = B * theta_dot_max / i_full;

      Serial.print(Kt);
      Serial.print("->");
      Serial.print(myKt);
      Serial.println();

    }
    void update(uint16_t pwm) {
      pwm_level = pwm;
      update();
    }

    void update() {
      unsigned long now = millis();
      static float i_last = 0;
      if (now - lastUpdate >= updateInterval) {
        lastUpdate += updateInterval;

        i_amps = (V * pwm_level / 255.0 - Ke * theta_dot - L * di_dt) / R;
        di_dt = i_amps - i_last;

        theta += theta_dot * updateInterval / 1000.0;
        theta_dotdot = (Kt * i_amps - B * theta_dot) / J;
        theta_dot += theta_dotdot * updateInterval / 1000.0;

      }
    }
    void printConfig(void) { // print the configuration vars
      Serial.print(" MaxRPM=");
      Serial.print(MaxRPM);
      Serial.print(" Ke=");
      Serial.print(Ke);
      Serial.print(" Kt=");
      Serial.print(Kt);
      Serial.print(" B=");
      Serial.print(B);
      Serial.print(" J=");
      Serial.print(J);
      Serial.print(" L=");
      Serial.print(L);
      Serial.print(" V=");
      Serial.print(V);
      Serial.print(" pwm_level=");
      Serial.print(pwm_level);
      Serial.println();
    }
    void printState(void) { // Print the state variables
      Serial.print(" pwm=");
      Serial.print(pwm_level);
      Serial.print(" i_amps=");
      Serial.print(i_amps);
      Serial.print(" theta=");
      Serial.print(theta);
      Serial.print(" theta_dot=");
      Serial.print(theta_dot);
      Serial.print(" theta_dotdot=");
      Serial.print(theta_dotdot);

      Serial.println();
    }

    void printStateRPM(void) { //print the state in rev & RPM form
      Serial.print("Motor: pwm=");
      Serial.print(pwm_level);
      Serial.print(" i=");
      Serial.print(i_amps);
      Serial.print("A ");
      Serial.print(theta / 2 / PI);
      Serial.print("rev ");
      Serial.print(theta_dot / 2 / PI * 60);
      Serial.print("RPM ");
      Serial.print(theta_dotdot / 2 / PI * 60);
      Serial.print("RPM/s ");

      Serial.println();
    }
};

PmMotor myMotor;

double myRPM = 0.0;
double mySetpointRPM = 0.0;
double myCV = 0.0;
const byte modePin = 2;
PID motorPID(&myRPM, &myCV, &mySetpointRPM, 1, 0.5, 0, DIRECT);

void report (void) {
  const int interval = 500;
  static unsigned long last = - interval;
  unsigned long now = millis();
  if (now - last >= interval) {
    last += interval;
    //motor.printState();
    Serial.print("PID: SP: ");
    Serial.print(mySetpointRPM);
    Serial.print("rpm MV: ");
    Serial.print(myRPM);
    Serial.print("rpm CV: ");
    Serial.print(myCV);
    Serial.print(" counts PWM9: ");
    //    Serial.print(OCR1A,DEC);
    Serial.print(analogGetPwm9(), DEC);
    Serial.println();
    myMotor.printStateRPM();
    //motor.printConfig();
  }
}

void setup() {
  // put your setup code here, to run once:
  Serial.begin(115200);
  myMotor.V = 5;
  myMotor.setResistance(1.2);
  myMotor.setKeMaxRPMVR(1526);
  myMotor.i_full = 0.1; // Amps at no-load full speed
  myMotor.B = 0.0001e-0; //
  myMotor.setInertia(0.1);
  myMotor.setKtFromIB();
  myMotor.printConfig();
  myMotor.printState();

  pinMode(modePin, INPUT_PULLUP);

  myRPM = myMotor.theta_dot * 60 / 2 / PI;
  mySetpointRPM = analogRead(A0); // 0-1024 rpm
  //turn the PID on
  motorPID.SetSampleTime(10); //us
  motorPID.SetMode(AUTOMATIC);

  delay(500);
}

int analogGetPwm9() { // PWM 9 on UNO is OCR1A
  // Retrieve PWM value from pin per https://forum.arduino.cc/t/what-is-the-inverse-of-analogwrite-pin-val/960774/20
  // Uno pin 9 is OC1A per
  // Uno https://docs.arduino.cc/hacking/hardware/PinMapping168
  // Mega https://docs.arduino.cc/hacking/hardware/PinMapping2560

  return ((TCCR1A & bit(COM1A1)) ? OCR1A : digitalRead(9) * 255);
}

void loop() {

  mySetpointRPM = analogRead(A0); // setpoint based on slider
  myRPM = myMotor.theta_dot * 60 / 2 / PI; // measure the motor speed
  if (digitalRead(modePin) == LOW) {
    motorPID.Compute();                 // PID to calculate CV from SP and MV
    analogWrite(9, myCV); // push PID control value out PWM pin 9
  } else {
    analogWrite(9, mySetpointRPM / 4);
  }

  // Simulate a motor attached to a PWM pin
  myMotor.update(analogGetPwm9());
  report();
}

ETA: Wokwi has an interesting way to encapsulate the motor dynamics simulation into a custom chip: See the dc_motor_chip Wokwi Simulation With that setup, you can separate the motor simulation from the Arduino code.

Dave, Just another thank you for this phenomenally brilliant code and extensive detailed explanation. cheers, Giles

1 Like

This topic was automatically closed 180 days after the last reply. New replies are no longer allowed.

This expression seems strange to me.

(T - Tamb) * area * h = W not T

1 Like

You are right. I should rearrange it to take into account time and heat capacity:

float Qconv =  (T - Tamb) * area * h; 
T = T + (Q - Qconv) * interval / 1000 / mass / Cps ;

take a look to this example, it uses the qLibs library and uses a transfer-function to simulate the system. The PID can be auto-tuned by pressing the pushbutton

GitHub - kmilo17pet/qlibs: A collection of useful libraries for embedded systems : signal smoothing, PID control, Fuzzy Logic, fixed-point math and more... (can be installed from the library manager)

1 Like

-- looks pretty cool.

1 Like

documentation here : API Reference: Overview

A simpler PID setup, relying on the user to simulate the process manually by moving a slide potentiometer is here:

This manual simulation might help in learning how the different Kp, Ki, and Kd parameters affect a PID output when process responds to external influences.

/* PID with two PWM outputs 180° out of phase
   for https://forum.arduino.cc/t/closed-loop-with-2-pwm-signals-180-degree-apart/1286865/25
  
  Use: Adjust the Setpoint Slider to indicate the setpoint to the controller
       Simulate the process and process disturbances by adjusting the 
       "ProcessOutput/ControllerInput" pot based on what you think the process will 
       do based on the duty cycle and whatever else.
       Monitor the Scope to see how the PID adjusts the duty cycles

Adjust the:
 Kp variable higher or lower to tmake the duty cycle more or less responsive to errors
 Ki variable higher or lower to increase or reduce pressure on modifying the duty cycle
 Kd variable higher or lower to anticipate more or less strongly the future effects of change 

*/

#include "TimerHelpers.h" // https://www.gammon.com.au/forum/?id=11504
// for https://forum.arduino.cc/t/two-pwm-signals-on-pins-9-and-10-with-input-on-a1-signals-need-to-be-180-degrees-apart/1280072/4?u=davex

#include "TimerHelpers.h"

#include <PID_v1.h>

int pwmPin1 = 9;
int pwmPin2 = 10;
double analogPin0 = A0;
double analogPin1 = A1;
double Out = 0;
double Involtage = 0;

double Setpoint, Input, Output;
double Kp = 0.4, Ki = 0.05, Kd = 0;
PID myPID(&Input, &Output, &Setpoint, Kp, Ki, Kd, DIRECT);

void setup() {
  pinMode(pwmPin1, OUTPUT);
  pinMode(pwmPin2, OUTPUT);

  TCCR1A = 0;        // reset timer 1
  TCCR1B = 0;

  Timer1::setMode (8, Timer1::PRESCALE_8, Timer1::CLEAR_A_ON_COMPARE | Timer1::SET_B_ON_COMPARE);
  ICR1 = 250; // 250 -> 32kHz w PRESCALE_1 Sets frequency with prescaler as f=F_CPU/PRESCALE/ICR/2

  TCNT1 = 0;         // reset counter
  OCR1A =  ICR1 / 4;     // compare A register value
  OCR1B =  ICR1 - OCR1A;   // compliment and invert

  Serial.begin(115200);
  analogWrite(pwmPin1, 120);

  //initialize the variables we're linked to
  Input = analogRead(analogPin0);
  Setpoint = 200; // 175= 10v 220=12V

  //turn the PID on
  myPID.SetMode(AUTOMATIC);
  myPID.SetOutputLimits(0, ICR1);
}

void loop() {

  Out = analogRead(analogPin0);// read the output voltage pin
  Involtage = analogRead(analogPin1);// read the input voltage pin
  Setpoint = Involtage;
  Input = Out;
  if (myPID.Compute()) {

    OCR1A = Output; // duty cycle
    OCR1B = ICR1 - OCR1A; // complement for inversion
    if (0) {
      Serial.print ("input volts    ");
      Serial.print ((Involtage * (5.0 / 1023)) * (12.3)); // 12.3 for resistors
      Serial.print("  ");
      Serial.print (Involtage);

      Serial.print("    output volts   ");
      Serial.print ((Out * (5.0 / 1023)) * (12.9)); // 12.9 for resistors
      Serial.print("  ");
      Serial.println (Out);
    }
  }
  report();
}

void report() {
  static uint32_t last = 0;
  if (millis() - last >= 500) {
    last += 500;
    Serial.print("Input:"); Serial.print(Input);
    Serial.print(" Setpoint:"); Serial.print(Setpoint);
    Serial.print(" Output:"); Serial.print(Output);
    Serial.print('/');
    Serial.print(ICR1);
    Serial.println();
  }
}



void reportTimer1(void) {
  Serial.print("TCCR1A:0b");
  Serial.println(TCCR1A, BIN);
  Serial.print("TCCR1B:0b");
  Serial.println(TCCR1B, BIN);
  Serial.print("ICR1:");
  Serial.println(ICR1);
  Serial.print("OCR1A:");
  Serial.println(OCR1A);
  Serial.print("OCR1B:");
  Serial.println(OCR1B);
  Serial.print("Duty:");
  Serial.println(OCR1A * 100.0 / ICR1, 3);
  Serial.print("Freq:");
  Serial.print(F_CPU / 256.0 / ICR1 / 2, 4);
  Serial.println("Hz\n");
}

Here are the plain Arduino PID example code demos implemented in Wokwi:

I think the simulation where you are the process/plant and use a slide potentiometer to communicate your ProcessVariable to the PID are the most hands-on in understanding those PID examples:

imagehttps://wokwi.com/projects/412997833205485569

and

image https://wokwi.com/projects/410423413759806465