Electronic – Why is a system description in control theory important

The question assumes there is one damping factor i.e. the transfer function is dominated by a 2nd order system. In practice there could be several interacting 2nd order systems so care has to be taken here.

Assuming it is dominated by a single 2nd order system, this diagram may be useful: -

This picture allows you to calculate \$\zeta\$ by looking at the peak amplitude value in the bode plot. From this you can calculate \$\omega_n\$ using the 2nd formula. Once \$\omega_n\$ is known you can double check the value of Q (1/2\$\zeta\$) because, for a 2nd order filter like this Q IS the magnitude of the amplitude response at \$\omega_n\$.

Here's the bigger picture showing how the bode plot and pole zero diagram are related: -

Your control is quite good despite if you feel in other way. I do think that a P-control with feed forward LUT (Look Up Table) would solve 98% of your need.

Try this way:

For each temperature setpoint you create a LUT table with output % that maintain the setpoint temperature at the steady state. If this setpoint is just 50 deg. C, then you don't even need a LUT, just feed approx. 20% as from your open loop tests. Then close the loop with a P-controller, the Kp=4 seems very good from your tests.

u=Kp*(T_sp-T_act)+20; // in %
if u>=100
  u=100;
else if u<=0
  u=0;
end if

As the error is subject to the whole loop gain, there will be always a little bit off, but maybe not so big error if you want to avoid complexity with PID. You don't even need sample time interrupt, you will be surprised by this simple control.

Once you get a stable control, then you can experiment with PID as well, then it's up to you to choose.

EDIT 1: You can improve the initial overshoot by using the ramp on a setpoint value. This will allow you to set a higher Kp.

EDIT 2:

Aström controller with tracking anti-windup (Source):

Aström controller with tracking anti-windup and with filtered D-part (Source):

input parameters
Ti ...integral time (s)
Td ...derivative time (s)
Kp ...proportional gain
b ...setpoint weighing factor for proportional part (0..1), default 1.0
Ts ...sampling time (s)
N ...D-part filter 1.5 to 30, default 15
umax ... max output limit
umin ... min output limit

--

float bi, ad, bd, a0, Tt, Kpold, Tiold, Tdold, ysp, y, yold, P, I, D;
float Kp = 7;
float b = 1.0;
float Ti = 0;
float Td = 0;
float N = 15.0;
float umin = 0.0;
float umax = 100.0;
float Ts = 0.1; // 100ms sample time, same for ISR interrupt
int UpdateParam = 0; // flag for parameter update

--

func init(void)
{
 if Td>0.0 
  Tt = sqrt(Ti*Td);
 else
  Tt = Ti/2.0;  
 
 if (Ti>0.0) // switch for disabling integral part 
  {
  bi = Kp*Ts/Ti; 
  a0 = Ts/Tt;
  }
 else
  {
  bi = 0.0;
  a0 = 0.0;
  }
 if Td>0.0 && (N>=0.0)
 {
  ad = (2*Td-N*Ts)/(2*Td+N*Ts);
  bd = (2*Kp*N*Td)/(2*Td+N*Ts);
 }
 else
 {
  ad = 0.0;
  bd = 0.0;
 }
}

--

func ISR(void)
{
if UpdateParam
{
  init();
  I=I+Kpold*(bold*ysp-y)-Kp*(b*ysp-y); //bumpless transfer of new parameters
  bold = bi;
  Kpold = Kp;
  Tiold = Ti;
  Tdold = Td; 
  UpdateParam = 0;
}
ysp = read(setpoint value);
y = read(process value);
P = Kp*(b*ysp-y);  // compute proportional part
D = ad*D - bd*(y-yold);  // update derivative part
v = P + I + D; // compute temporary output
if v>=umax  // limit output
 u=umax;
else if v<=umin
 u=umin;
else
 u=v;
output = write(u);
I = I + bi*(ysp-y) + a0(u-v);
yold = y;

}

--

main
{
  init();
  I = 0.0;
  bold = bi;
  Kpold = Kp;
  Tiold = Ti;
  Tdold = Td;   
  __ISR_Enable();

  while(1)
  {
    if (Kp<>Kpold) || (Ti<>Tiold) || (Td<>Tdold)
     UpdateParam = 1;
  }
}

You start with Ti=0, Td=0 and then execute the ZN tuning procedure. According to ZN you do calculate Kp, Ti, Td. A 20% oversoot is a normal thing, you can avoid the overshoot by introducing a ramp function on the setpoint value. It would gradually increase the temperature.

--

func Ramp(void)  //execution in ISR at every Ts
{
// delta_max ... ramp slope (degrees C/s)
// yspin ...  input setpoint value
step = delta_max * Ts;
if yspin > (ysp + step)
  ysp = ysp + step;
else if yspin < (ysp - step)
  ysp = ysp - step;
else
 ysp = yspin;
}

--