PID using integer math. – 8 days old

[keytapper]

Hello Forum,
I'm struggling to implement a PID controller for my oven. The sampling rate is dictated by the time to read a MAX6675, which is 0.25 seconds to complete a full conversion.
The output should drive a homemade SSR, which using MOC3061 with ZCD. So I can only modulate the number of half sine wave on and off.
As far as my researches I found the AN964, but the parameters are different so I don't get the idea how to convert the 360 degree from that note into 100 pulses per minute in my attempt.
So the Kp, Ki and Kd can be an integer size and then there should be a divider to reduce them into fractional values. But the output should be from 0 to 100, with 250 mS time frame. So I'm puzzled to make the right corrections.

Currently I just implement a proportional method and it seems to be rather enough, as fas as I simulated it.

    Symbol MINTEMP  40                  ' minimum temperature level
    Symbol MAXTEMP 260                  ' maximum temperature level
    Symbol TEMPSPAN MAXTEMP - MINTEMP   ' temperature span

  Clear Result                      ' forget last state
  Clear AlarmFlag                   ' forget The flag
  If currtemp > Setpoint Then       ' when is over the Setpoint
    If currtemp > TEMPSPAN Then     ' may be over the guard limit
      Btfss no_PID                  ' when the flag is clear
      Set AlarmFlag                 ' it will signalize an error
    End If
    Return                          ' Stop there, return result zeroed
  End If
  If no_PID = TRUE Then             ' disable the proportional controller
    Set Result
    Return
  End If
  If currtemp > 0 Then              ' whether is higher than minimum temperature
    _PP2 = currtemp * 100/TEMPSPAN
    y = currtemp * KPD /Setpoint
    If ticks >= y.LowByte Then      ' the y higher the lesser ticks to stay on
      Set Result                    ' let fire the output
    End If
  Else                              ' in case it hasn't got to minimum level
    Set Result                      ' let the output to be on.
  End If
EndProc

I should say that the temperature is reduced to one byte, as far as the range is from 40° to 260°C. When it's simulated I can see the oscillation is limited within 3 degrees, but when I trying the thermocouple it comes with a bigger hysteresis, that's why I suppose the PID would give a better solution.
I found several pages, but all using floating point for the three constant Kp. Ki & Kd. I think I can do similar with integer values as long as using the method to multiply and divide with integers, such is done in that application note.

[John Drew]

Don't be afraid to use floating point. FP maths is sufficiently accurate for many applications and is very easy to use in Proton. Unless you're dealing with a tiny device at very low speed use FP if it makes your maths easier. I wouldn't hesitate. Unfortunately I didn't understand your abbreviations in the post.
John

[keytapper]

Thanks, John.
I realized that the presented code is rather meaningless. The ZCD stands for Zero Crossing Detection.
My intent is to remain on the tiny MCU, so here I'm to find what that application note is able to accomplish within few hundreds of instruction code.
The above procedure evaluate the percentage to keep the TRIAC on or off, according the counter that is counting downward at every zero cross, which is sensed by the external interrupt. KPD is a correction constant to give a little better solution.
Then in the main loop it will set the pin output according what the procedure result was comparing. The task in the main loop is carried out in less than 0.2 milliseconds and then there's plenty of idle time.

The only part that I'm missing is how to scale down the output result in a number that should be from 0 to 100. Whereas the AN964 has a range from -180 to 180. Then I will be able to handle how many semi-waves to skip. There could be even a better way to distribute the OFF time in one second period, E.G. if the result si 10% it will be possible to skip 1 over 10, rather than 10 over 100. But this is another use case that in a tiny MCU will not fit.

One thing I could decipher is that the three constant are read from the assembled potentiometers so the calculation will get the 10 bit ADC value and divide by 1024. The most relevant variables are SWord size, only the output is 32bit wide. But I still have room in my source for the RAM. I don't understand the scaling factor and perhaps the way I should take samples in a 2^n value, so the calculation will be just a bit shifting matter.

[See_Mos]

If I read your post correctly you are detecting every zero crossing?

When using a zero crossing opto you do not need to detect zero crossing at all because the opto will take care of that. While opto LED is on it will drive both half cycles and in a worst case you will only miss one half cycle.

Decide on a gate time. for a small heater it might be one second or for a very big heater it might be ten minutes. Some of the ones I worked on for nitrogen extractors were 65KW and the gate time was 10 minutes.

Then just calculate the on and off times as a percentage of the gate time. If you use an interrupt timer for the gate that leaves plenty of time to take the sample and do the sums.

Interestingly the end user could not understand why they could not run the nitrogen extractors on 110KVA generators. The company that built them had used phase angle control so at 50% heat requirement the generator was hit with a sudden current demand at the peak of the sine wave. The result was that the generator became very unstable, electrically and mechanically! The automatic voltage controller could not cope, the generators were going well over voltage and blowing 1500 volt voltage dependent resistors in the heater control.

[keytapper]

If I read your post correctly you are detecting every zero crossing?

When using a zero crossing opto you do not need to detect zero crossing at all because the opto will take care of that.

I started by using the timer0 to count 10 mSec, the maximum permissible count, so I can use only one byte to measure one second. Because that is the part of the timer that I will use to set a heating period. I was worry that the TRIAC firing wasn't synchronize with the main. So I turned to use the count over the frequency main supply. Even the MOC3061 is with ZCD, so my best PWM is based on skipping some semi-wave.

To have a regulated output I suppose to work on the skipped wave number. As the rectified waves are 100, therefore I consider to adjust the percentage that the TRIAC will stay on in a second.
The MAX6675 needs at least 0.17 second for its ADC conversion. That implies the time frame won't be less than that, because the PID error calculation is possible only when a new sample is taken.
To compare with the AN964 I need to understand what is the delta time to be used to calculate the integral and perhaps the derivative terms.

I found a very informative page, I still need some more efforts to grasp the entire concept.

The integral of a signal is the sum of all the instantaneous values that the signal has been, from whenever you started counting until you stop counting.

This is the point where I should find the way to determine how long should be a time frame for my case. The two incognitos are when is time to take a sample and what is the relationship with 100 pulse per second (50 Hz x 2, from the main).

I may guess that the degree per second might be about 5°C/sec. Then I prepared a thermal switch at 260°C, just to add further safety in case the system may be faulty. But the oven should go up to 250° and not higher. That is a particular condition that the PID controller won't have to exceed the set point for a long period.

[keytapper]

One further step, I converted the AN964 source to proton basic

On_Interrupt GoTo interrupt_Isr
Declare Xtal 8
Device 16F684

Symbol TRUE 1
Symbol FALSE 0
Symbol ADON = ADCON0.0                 ' A/D ON bit
Symbol GO_DONE = ADCON0.1              ' A/D Conversion Enable
Symbol CHS0 = ADCON0.2                 ' Analog Channel Select bits
Symbol CHS1 = ADCON0.3                 ' Analog Channel Select bits
Symbol CHS2 = ADCON0.4                 ' Analog Channel Select bits
Symbol VCFG = ADCON0.6                 ' VCFG: Voltage Reference bit
Symbol ADFM = ADCON0.7                 ' A/D Result Formed Select bit
Symbol PS0 = OPTION_REG.0              ' Prescaler Rate Select
Symbol PS1 = OPTION_REG.1              ' Prescaler Rate Select
Symbol PS2 = OPTION_REG.2              ' Prescaler Rate Select
Symbol PSA = OPTION_REG.3              ' Prescaler Assignment
Symbol T0CS = OPTION_REG.5             ' TMR0 Clock Source Select
Symbol ADCS0 = ADCON1.4             ' A/D Conversion Clock Select bit
Symbol ADCS1 = ADCON1.5             ' A/D Conversion Clock Select bit
Symbol ADCS2 = ADCON1.6             ' A/D Conversion Clock Select bit
Symbol RAIF = INTCON.0               ' RA Port Interrupt Flag
Symbol TMR0IF = INTCON.2             ' TMR0 Overflow Interrupt Flag
Symbol INT0IE = INTCON.4             ' RA0 External Interrupt Enable
Symbol TMR0IE = INTCON.5             ' TMR0 Overflow Interrupt Enable
Symbol GIE = INTCON.7                ' Global Interrupt Enable
Symbol ANS2 = ANSEL.2                 ' Pin Analog / Digital Selection
Symbol ANS4 = ANSEL.4                 ' Pin Analog / Digital Selection
Symbol ANS5 = ANSEL.5                 ' Pin Analog / Digital Selection
Symbol IRCF0 = OSCCON.4             ' Internal Oscillator Frequency Select bits
Symbol IRCF1 = OSCCON.5             ' Internal Oscillator Frequency Select bits
Symbol IRCF2 = OSCCON.6             ' Internal Oscillator Frequency Select bits
Symbol TMR2ON = T2CON.2              ' Timer2 ON
Symbol DC1B0 = CCP1CON.4              ' PWM Duty Cycle Least Significant bits
Symbol DC1B1 = CCP1CON.5              ' PWM Duty Cycle Least Significant bits
Symbol P1M1 = CCP1CON.7               ' PWM Output Configuration bit
Symbol RA4 = PORTA.4                   ' Bi-directional I/O pin (with interrupt on change)
Symbol  SumE_Max = 30000
Symbol  SumE_Min = 1 - SumE_Max

Dim do_PID As Bit
Dim ADchanel As ADRESL.Word
Dim en0 As SByte
Dim en1  As SByte
Dim en2  As SByte
Dim en3  As SByte
Dim off_set  As SByte
Dim temp As Byte
Dim temp_int As SWord
Dim ki As Word
Dim kd As Word
Dim kp As Word
Dim SumE As SWord
Dim integral_term As SWord
Dim derivative_term As SWord
Dim Cn As SDword

main:
  Init()
  Set_Constants()
  Do
    If do_PID = TRUE Then
      PID()
    End If
  Loop
End

Sub Init()
  PORTA = 0
  TRISA = 101101
  PORTC = 0
  TRISC = 000011
  CMCON0 = $07
  IRCF0 = 1
  IRCF1 = 1
  IRCF2 = 1
  CCP1CON = 001100
  ECCPAS = 0
  PR2 = $3F
  T2CON = 0
  CCPR1L = 0
  TMR2ON = 1
  ANSEL = 110101
  VCFG = 0
  ADFM = 1
  ADCS0 = 1
  ADCS1 = 0
  ADCS2 = 1
  CHS0 = 0
  CHS1 = 0
  CHS2 = 0
  ADON = 1

  Clear en0
  Clear en1
  Clear en2
  Clear en3
  Clear ki
  Clear kd
  Clear kp
  Clear off_set
  Clear temp_int
  Clear integral_term
  Clear derivative_term
  do_PID = 1
  T0CS = 0
  TMR0 = 10
  PSA = 0
  PS0 = 0
  PS1 = 0
  PS2 = 1
  INTCON = 0
  PIE1 = 0
  TMR0IE = 1
  GIE = 1
EndSub


Sub PID()
  Clear integral_term
  Clear derivative_term
  SumE = SumE + en0
  If  SumE > SumE_Max Then
    SumE = SumE_Max
  End If
  If  SumE < SumE_Min Then
    SumE = SumE_Min
  End If

  integral_term = SumE / 256
  integral_term = integral_term * ki
  derivative_term = en0 - en3
  If  derivative_term > 120 Then
    derivative_term = 120
  End If
  If  derivative_term < -120 Then
    derivative_term = -120
  End If

  derivative_term = derivative_term * kd
  derivative_term = derivative_term >> 5

  If  derivative_term > 120 Then
    derivative_term = 120
  End If
  If  derivative_term < -120 Then
    derivative_term = -120
  End If

  Cn = en0 + integral_term + derivative_term
  Cn = Cn * kp / 1024

  If  Cn >= 1000 Then
    Cn = 1000
  End If
  If  Cn <= -1000 Then
    Cn = -1000
  End If
  If  Cn == 0 Then
    Clear DC1B1
    Clear DC1B1
    Clear CCPR1L
  End If
  If  Cn > 0 Then
    P1M1 = 0
    temp = Cn
    If  temp ^ 000001 = TRUE Then
      DC1B0 = 1
    Else
      DC1B0 = 0
    End If
    If  temp ^ 000010 = TRUE Then
      DC1B1 = 1
    Else
      DC1B1 = 0
    End If
    CCPR1L = Cn >> 2
    off_set = off_set +1
    If  off_set > 55 Then
      off_set = 55
    End If
  Else
    P1M1 = 1
    temp_int = Abs(Cn)
    temp = temp_int
    If temp ^ 000001 = TRUE Then
      DC1B0 = 1
    Else
      DC1B0 = 0
    End If
    If temp ^ 000010 = TRUE Then
      DC1B1 = 1
    Else
      DC1B1 = 0
    End If
    CCPR1L = temp_int >> 2
    off_set = off_set -1
    If  off_set < -55 Then
      off_set = -55
    End If
  End If
  en3 = en2
  en2 = en1
  en1 = en0
  en0 = 0
  do_PID = 0
  RA4 = 0
EndSub

Sub Set_Constants()
  ANS2 = 1
  ANS4 = 1
  ANS5 = 1
  ADFM = 1
  CHS0 = 0
  CHS1 = 0
  CHS2 = 1
  GO_DONE = 1
  While GO_DONE = TRUE
    DelayUS 20
  Wend
  ki = ADchanel
  CHS0 = 1
  CHS1 = 0
  CHS2 = 1
  GO_DONE = 1
  While GO_DONE = TRUE
    DelayUS 20
  Wend
  kd = ADchanel
  CHS0 = 0
  CHS1 = 1
  CHS2 = 0
  While GO_DONE = TRUE
    DelayUS 20
  Wend
  GO_DONE = 1
  kp = ADchanel
  CHS0 = 0
  CHS1 = 0
  CHS2 = 0
EndSub

interrupt_Isr:
  Context Save
  If TMR0IF = TRUE Then
    If TMR0IE = TRUE Then
      TMR0 = 10
      TMR0IF = 0
      RA4 = 1
      temp_int = ADchanel
      temp_int = temp_int - 512
      en0 = temp_int + off_set/8
      do_PID = 1
      GO_DONE = 1
    Else
      PIR1 = 0
      RAIF = 0
      INT0IE = 0
    End If
  End If
  If temp_int > 180 Then
    GoTo Clrvars
  End If
  If  temp_int < -180 Then
Clrvars:
    Clear DC1B0
    Clear DC1B1
    Clear CCPR1L
    Clear en0
    Clear en1
    Clear en2
    Clear en3
    Clear off_set
    Clear Cn
    Clear integral_term
    Clear derivative_term
    Clear SumE
    Clear RA4
    Clear do_PID
  End If
  Context Restore

So I could figure out the resulting volume of the assembly. Definitely there will be less operations and no Analog input to add. Perhap the three constant will go on the EEPROM.
The original source seems to me that had no much optimization, but it fits in less than 1 Kwords.

[top204]

With most controls, the full PID is not required, and can be slow to respond.

I have, generally, used a PI method, where the error is stored and then the output is altered until the error is 0, or close to it.

I created a PID control for a constant current source when plating copper, but even on a dsPIC33, the initial response from the PID was slow, until it had matched close to what it should be, which casues rough surfaces on the copper plating, so I changed it to a simple error offset mechanism and it worked much better, and the volume of heat and fluid made the current changes slow to alter, so full PID was not required.

[See_Mos]

For a small AC oven heater control using BURST FIRE CONTROL I would start with a gate period of ten seconds. Divide the period by 100 to give tenth of a second control which is equal to a minimum of four and a half or five cycles ON time at 50Hz depending on when the MOC3061 LED is turned on

Any finer resolution would be a waste of processing time.

[keytapper]

Forgive my lack of knowledge.
I'm trying con conceive how to write the formula.
I'm seeing that one degree increase may take a couple of seconds. So I agreed that the time frame may be 10 seconds long. During this period I can keep the ON_time of the TRIAC as the result of the calculation will give.

My concern is about keeping the tolerance within +/- 2 degree, otherwise it would not make sense to set the temperature and see 10% variations.

Maybe Fanie's wife got some experience about carbonized foods .

@ Les
I tried the simple ON/OFF control but I'm getting an horrible 8~10% variations around the set point.
I tried a proportional control which reduces the power factor as soon as the temperature is close to the set point. The simulation promising a 1% variation, but the hardware won't reach the set point at all.

[See_Mos]

Yes, I was thinking again. Ten seconds might be too long so maybe one second but then the MAX6675 is a problem.

[keytapper]

one second but then the MAX6675 is a problem.

The max6675 needs 0.17 seconds (typical) to do a full conversion. I set for 0.25 seconds, so that's the minimum time frame and it rounds up with 4 samples per second.

[See_Mos]

You are correct, after posting I read the data sheet and it only takes about 3.2mS to read the information and typically 170mS while the device does the conversion so if you only take one or two readings a second there is plenty of time to do the maths and do the rest.

Unfortunately maths is my weakness so cannot help there.

[keytapper]

I am studying all day how to whip up the formula. Roughly it cost some 400 words. But I can't get the correct scaling to set the PWM for the TRIAC. In fact the TRIAC should keep on a modulated percentage until a new reading arrive.
Compare to those formulas I saw on the web, the time frame is not necessary whenever it's working in a constant rate. So as I may guess the procedure may be a couple of second long, in that time the temperature doesn't vary that much. Even doing with 0.01 seconds I still have lot of idle time at 20 MHz clock.

A very rough design is as follow:

Device 16F628A
Declare Xtal 20

Symbol Kp 1298
Symbol Ki 987
Symbol Kd 14
Symbol DeltaTime 10
' rough figure in mills. Better in binary value to reduce calculation in a shift right.
$define fractionOf(x) x / 1024 

Dim err As SByte
Dim previous_error As SByte
Dim derivative As SByte
Dim integral As SWord
Dim PIDout As SByte
Dim tmp As SDword
Dim currtemp As Byte
Dim Setpoint As Byte

Do
  err = Setpoint - currtemp
  tmp = err * fractionOf(Kp)
  integral = integral + err*DeltaTime
  If integral > 32000 Then integral = 32000
  If integral < -32000 Then integral = -32000
  tmp = tmp + integral * fractionOf(Ki)
  derivative = (err - previous_error)/DeltaTime
  tmp = tmp + derivative * fractionOf(Kd)
  If tmp > 100 Then tmp = 100
  If tmp < 0 Then tmp = 0
  PIDout = tmp
  previous_error = err
  DelayMS DeltaTime
Loop

If sampling rate is constant, then the DeltaTime is 1, therefore may be eliminated.
Thus it's like the bloggers suggest, but I don't realize how to adjust the output in the correct scaling which will give the T_ON and T_OFF for the TRIAC. I'm facing a huge jump on the integral that it won't recover. Also the output starts to reduce very early.

[keytapper]

I think the algorithm works. It just take in account that the procedure should follow the response time of what is measured. like this page mention, the sampling would be about 10 time smaller than the measure variation. Thus when the variation may be 15 second per degree, the sample shouldn't be smaller than 1.5 second, otherwise the integral variable will reach its limits without any chance to recover.
There's just a nasty situation when starting that cause an overshooting. That seems normal, as many pages reporting.

I just try the algorithm with the simulator, the result is pretty like what I expect. Tomorrow I'll try with my hardware simulation.

[Ken_K]

A most interesting topic, I did not see the C code for the app note 964 so I re posted it here.

http://read.pudn.com/downloads37/sou...ndulum.c__.htm

and the pendulum app note again

http://ww1.microchip.com/downloads/en/Appnotes/00964A.pdf

These are worth saving and are now in my library.