Merge pull request #26 from adq/ziegler

Ziegler-Nicols kiln-tuner

docs/ziegler_tuning.md

@@ -0,0 +1,89 @@
+# PID Tuning Using Ziegler-Nicols
+
+This uses the Ziegler Nicols method to estimate values for the Kp/Ki/Kd PID control values.
+
+The method implemented here is taken from ["Ziegler–Nichols Tuning Method"](https://www.ias.ac.in/article/fulltext/reso/025/10/1385-1397) by Vishakha Vijay Patel
+
+One issue with Ziegler Nicols is that is a **heuristic**: it generally works quite well, but it might not be the optimal values. Further manual adjustment may be necessary.
+
+## Process Overview
+
+1. First of all, you will record a temperature profile for your kiln.
+2. Next, we use those figures to estimate Kp/Ki/Kd.
+
+## Step 1: Record Temperature Profie
+
+Ensure `kiln-controller` is **stopped** during profile recording: The profile must be recorded without any interference from the actual PID control loop (you also don't want two things changing the same GPIOs at the same time!)
+
+Make sure your kiln is completely cool - we need to record the data starting from room temperature to correctly measure the effect of kiln/heating.
+
+There needs to be no abnormal source of temperature change to the kiln: eg if you normally run with a kiln plug in place - make sure its in place for the test!
+
+To record the profile, run:
+```
+python kiln-tuner.py recordprofile zn.csv
+```
+
+The above will drive your kiln to 400 and record the temperature profile to the file `zn.csv`. The file will look something like this:
+
+```
+time,temperature
+4.025461912,45.5407078
+6.035358906,45.5407078
+8.045399904,45.5407078
+10.05544925,45.59087846
+...
+```
+
+## Step 2: Compute the PID parameters
+
+Once you have your zn.csv profile, run the following:
+
+```
+python kiln-tuner.py zn zn.csv
+```
+
+The values will be output to stdout, for example:
+```
+Kp: 3.853985144980333 1/Ki: 87.78173053095107 Kd: 325.9599328488931
+```
+(Note that the Ki value is already inverted ready for use in config)
+
+------
+
+## Sanity checking the results
+
+If you run
+```
+python kiln-tuner.py zn zn.csv --showplot
+```
+
+It will display a plot of the parameters. It should look simular to this ![kiln-tuner-example.png](kiln-tuner-example.png).
+
+Note: you will need python's `pyplot` installed for this to work.
+
+The smooth linear part of the chart is very important. If it is too short, try increasing the target temperature (see later).
+
+The red diagonal line: this **must** follow the smooth part of your chart closely.
+
+## My diagonal line isn't right
+
+You might need to adjust the line parameters to make it fit your data properly. You can do this as follows:
+
+```
+python kiln-tuner.py zn zn.csv --tangentdivisor 4
+```
+
+`tangentdivisor` modifies which parts of the profile is used to calculate the line.
+
+It is a floating point number >= 2; If necessary, try varying it till you get a better fit.
+
+## Changing the target temperature
+
+By default it is 400. You can change this as follows:
+
+```
+python kiln-tuner.py recordprofile zn.csv --targettemp 500
+```
+
+(where the target temperature has been changed to 500 in the example above)

kiln-tuner.py

@@ -0,0 +1,206 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import csv
+import time
+import argparse
+
+
+def recordprofile(csvfile, targettemp):
+
+    try:
+        sys.dont_write_bytecode = True
+        import config
+        sys.dont_write_bytecode = False
+
+    except ImportError:
+        print("Could not import config file.")
+        print("Copy config.py.EXAMPLE to config.py and adapt it for your setup.")
+        exit(1)
+
+    script_dir = os.path.dirname(os.path.realpath(__file__))
+    sys.path.insert(0, script_dir + '/lib/')
+
+    from oven import RealOven, SimulatedOven
+
+    # open the file to log data to
+    f = open(csvfile, 'w')
+    csvout = csv.writer(f)
+    csvout.writerow(['time', 'temperature'])
+
+    # construct the oven
+    if config.simulate:
+        oven = SimulatedOven()
+    else:
+        oven = RealOven()
+
+    # Main loop:
+    #
+    # * heat the oven to the target temperature at maximum burn.
+    # * when we reach it turn the heating off completely.
+    # * wait for it to decay back to the target again.
+    # * quit
+    #
+    # We record the temperature every second
+    try:
+        stage = 'heating'
+        if not config.simulate:
+            oven.output.heat(1, tuning=True)
+
+        while True:
+            temp = oven.board.temp_sensor.temperature + \
+                config.thermocouple_offset
+
+            csvout.writerow([time.time(), temp])
+            f.flush()
+
+            if stage == 'heating':
+                if temp >= targettemp:
+                    if not config.simulate:
+                        oven.output.heat(0)
+                    stage = 'cooling'
+
+            elif stage == 'cooling':
+                if temp < targettemp:
+                    break
+
+            sys.stdout.write(f"\r{stage} {temp:.2f}/{targettemp}           ")
+            sys.stdout.flush()
+            time.sleep(1)
+
+        f.close()
+
+    finally:
+        # ensure we always shut the oven down!
+        if not config.simulate:
+            oven.output.heat(0)
+
+
+def line(a, b, x):
+    return a * x + b
+
+
+def invline(a, b, y):
+    return (y - b) / a
+
+
+def plot(xdata, ydata,
+         tangent_min, tangent_max, tangent_slope, tangent_offset,
+         lower_crossing_x, upper_crossing_x):
+    from matplotlib import pyplot
+
+    minx = min(xdata)
+    maxx = max(xdata)
+    miny = min(ydata)
+    maxy = max(ydata)
+
+    pyplot.scatter(xdata, ydata)
+
+    pyplot.plot([minx, maxx], [miny, miny], '--', color='purple')
+    pyplot.plot([minx, maxx], [maxy, maxy], '--', color='purple')
+
+    pyplot.plot(tangent_min[0], tangent_min[1], 'v', color='red')
+    pyplot.plot(tangent_max[0], tangent_max[1], 'v', color='red')
+    pyplot.plot([minx, maxx], [line(tangent_slope, tangent_offset, minx), line(tangent_slope, tangent_offset, maxx)], '--', color='red')
+
+    pyplot.plot([lower_crossing_x, lower_crossing_x], [miny, maxy], '--', color='black')
+    pyplot.plot([upper_crossing_x, upper_crossing_x], [miny, maxy], '--', color='black')
+
+    pyplot.show()
+
+
+def calculate(filename, tangentdivisor, showplot):
+    # parse the csv file
+    xdata = []
+    ydata = []
+    filemintime = None
+    with open(filename) as f:
+        for row in csv.DictReader(f):
+            try:
+                time = float(row['time'])
+                temp = float(row['temperature'])
+                if filemintime is None:
+                    filemintime = time
+
+                xdata.append(time - filemintime)
+                ydata.append(temp)
+            except ValueError:
+                continue  # just ignore bad values!
+
+    # gather points for tangent line
+    miny = min(ydata)
+    maxy = max(ydata)
+    midy = (maxy + miny) / 2
+    yoffset = int((maxy - miny) / tangentdivisor)
+    tangent_min = tangent_max = None
+    for i in range(0, len(xdata)):
+        rowx = xdata[i]
+        rowy = ydata[i]
+
+        if rowy >= (midy - yoffset) and tangent_min is None:
+            tangent_min = (rowx, rowy)
+        elif rowy >= (midy + yoffset) and tangent_max is None:
+            tangent_max = (rowx, rowy)
+
+    # calculate tangent line to the main temperature curve
+    tangent_slope = (tangent_max[1] - tangent_min[1]) / (tangent_max[0] - tangent_min[0])
+    tangent_offset = tangent_min[1] - line(tangent_slope, 0, tangent_min[0])
+
+    # determine the point at which the tangent line crosses the min/max temperaturess
+    lower_crossing_x = invline(tangent_slope, tangent_offset, miny)
+    upper_crossing_x = invline(tangent_slope, tangent_offset, maxy)
+
+    # compute parameters
+    L = lower_crossing_x - min(xdata)
+    T = upper_crossing_x - lower_crossing_x
+
+    # Magic Ziegler-Nicols constants ahead!
+    Kp = 1.2 * (T / L)
+    Ti = 2 * L
+    Td = 0.5 * L
+    Ki = Kp / Ti
+    Kd = Kp * Td
+
+    # outut to the user
+    print(f"Kp: {Kp} 1/Ki: {1/ Ki}, Kd: {Kd}")
+
+    if showplot:
+        plot(xdata, ydata,
+             tangent_min, tangent_max, tangent_slope, tangent_offset,
+             lower_crossing_x, upper_crossing_x)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description='Kiln tuner')
+    subparsers = parser.add_subparsers()
+    parser.set_defaults(mode='')
+
+    parser_profile = subparsers.add_parser('recordprofile', help='Record kiln temperature profile')
+    parser_profile.add_argument('csvfile', type=str, help="The CSV file to write to.")
+    parser_profile.add_argument('--targettemp', type=int, default=400, help="The target temperature to drive the kiln to (default 400).")
+    parser_profile.set_defaults(mode='recordprofile')
+
+    parser_zn = subparsers.add_parser('zn', help='Calculate Ziegler-Nicols parameters')
+    parser_zn.add_argument('csvfile', type=str, help="The CSV file to read from. Must contain two columns called time (time in seconds) and temperature (observed temperature)")
+    parser_zn.add_argument('--showplot', action='store_true', help="If set, also plot results (requires pyplot to be pip installed)")
+    parser_zn.add_argument('--tangentdivisor', type=float, default=8, help="Adjust the tangent calculation to fit better. Must be >= 2 (default 8).")
+    parser_zn.set_defaults(mode='zn')
+
+    args = parser.parse_args()
+
+    if args.mode == 'recordprofile':
+        recordprofile(args.csvfile, args.targettemp)
+
+    elif args.mode == 'zn':
+        if args.tangentdivisor < 2:
+            raise ValueError("tangentdivisor must be >= 2")
+
+        calculate(args.csvfile, args.tangentdivisor, args.showplot)
+
+    elif args.mode == '':
+        parser.print_help()
+        exit(1)
+
+    else:
+        raise NotImplementedError(f"Unknown mode {args.mode}")

lib/oven.py

@@ -27,8 +27,10 @@ class Output(object):
             log.warning(msg)
             self.active = False
 
-    def heat(self,sleepfor):
+    def heat(self,sleepfor, tuning=False):
         self.GPIO.output(config.gpio_heat, self.GPIO.HIGH)
+        if tuning:
+            return
         time.sleep(sleepfor)
 
     def cool(self,sleepfor):