diff --git a/docs/kiln-tuner-example.png b/docs/kiln-tuner-example.png new file mode 100644 index 0000000..9f4c28d Binary files /dev/null and b/docs/kiln-tuner-example.png differ diff --git a/docs/ziegler_tuning.md b/docs/ziegler_tuning.md new file mode 100644 index 0000000..89719c8 --- /dev/null +++ b/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) diff --git a/kiln-tuner.py b/kiln-tuner.py new file mode 100755 index 0000000..1fdda65 --- /dev/null +++ b/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}") diff --git a/lib/oven.py b/lib/oven.py index 974c5a0..2d4c05a 100644 --- a/lib/oven.py +++ b/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):