PID tuning... It's broken, but how is best to fix it?

Developers
1

PID autotuning in Klipper currently uses the Ziegler-Nichols method to generate PID parameters, with the critical point determined by the Astrom-Hagglund method but using a relay with hysteresis. It’s all in this file.

When implemented, PID uses a Conditional Integration PID loop.

It appears such a setup has numerous shortcomings, and fixing each one just makes the others more problematic… They are:

  • The hysteresis level in the Astrom-Hagglund method (TUNE_PID_DELTA, a 5C constant) would ideally be zero. By changing that constant, we can see quite how problematic that is:
Send: PID_CALIBRATE HEATER=heater_bed TARGET=50 TUNE_PID_DELTA=0.5
Recv: // PID parameters: pid_Kp=217.643 pid_Ki=24.183 pid_Kd=489.697

Send: PID_CALIBRATE HEATER=heater_bed TARGET=50 TUNE_PID_DELTA=1.0
Recv: // PID parameters: pid_Kp=152.234 pid_Ki=10.356 pid_Kd=559.459

Send: PID_CALIBRATE HEATER=heater_bed TARGET=50 TUNE_PID_DELTA=2.0
Recv: // PID parameters: pid_Kp=80.175 pid_Ki=3.448 pid_Kd=466.015

Send: PID_CALIBRATE HEATER=heater_bed TARGET=50 TUNE_PID_DELTA=5.0
Recv: // PID parameters: pid_Kp=41.185 pid_Ki=0.708 pid_Kd=599.243
  • Yes, the Kp constant really does vary by a factor of 5x depending how its measured!! I think this is the main reason people complain of a slow frequency response and long heatup times. They are using grossly undertuned PID parameters.

The TUNE_PID_DELTA value cannot be set smaller because it is vulnerable to temperature sensor noise. It would be necessary to simulate an ideal relay and use frequency analysis to find the critical frequency, but since the temperature sensors don’t sample at consistent intervals you need some pretty specialized methods to do this (ie. you’d end up needing most of scipy as a dependency). I did this anyway to get the perfect parameters… Which brings me to the next shortcoming:

  • The PID loop uses Conditional Integration - ie. we don’t update the integral term if the output is max power or zero. With undertuned PID constants, this works fine… But when correctly tuned with the Astrom-Hagglund method, the output value frequently hits maximum or zero because the constants are much higher, while the sampling rate and noise are unchanged. That in turn means that we miss enough integral updates to start to affect performance and stability. I tried removing Conditional Integration and stability is now great… but obviously there are now integral windup issues instead.

Here are some ideas for solutions:

  • Run the tuning twice, with TUNE_PID_DELTA=5.0 and then =2.5. Use that to linearly interpolate the values that would have been measured with an ideal realy (ie. delta=0). That should get close to the correct PID values, at the expense of double the tuning time.

  • Modify the PID loop to back-calculate the integrator value as an antiwindup mechanism. Ie. figure 6.7 here.

Ideas welcome.

2

Well, for my setup the PID tuning works good (good enough for me) but you are right, there are edge cases where it does not yield the wanted results. You might want to check:

3

I’m sure it is possible to improve the PID and PID_CALIBRATE process. However, my own experience is that one can spend dozens of hours attempting to tune/code a control system and end up with only minor improvements. In many cases it even subtly reduces robustness. FWIW, I therefore place high value on keeping things simple.

Some random comments:

TUNE_PID_DELTA: As you indicate, it is there to account for signal measurement noise. I’ve read papers that said it is okay to use this type of delta during the relay test, but I don’t have the citations handy. I always felt it was a bit odd to have it, but it is a common mechanism in relay tests.

Conditional integration: As indicated it is a common mechanism to avoid “integral windup”. It’s very common in PID controllers and this is the first I’ve heard that it is a significant issue. I agree that “back calculating” the impact of a full-on/full-off heat setting would be theoretically better - however it would require an accurate model of the heater - which would add significant complexity and risk a loss of robustness.

I don’t know why you say this. The time between temperature measurements should be very precise (handful of microseconds).

-Kevin

4

I tried out this method for a ceramic heater that was jumping around wildly and higher delta was just making the problem worse. I actually used ChatGPT to do the interpolation since I suck at math. It completely fixed the tuning on my hotend, and its hitting my set temperature now without a single overshot. Seems to work much better than a higher tune_pid_delta alone.

Thanks @Hello1024!

5

Just my opinion, scary.

3 Likes
6

@hcet It scares me too honestly. If it helps its not exactly easy to get good results, requires some intuition. Also when I say I used ChatGPT I mean I had ChatGPT write python code that would parse klipper log messages and interpolate results. I can read math I just don’t like to write it.

2 Likes
7

3 posts were split to a new topic: GAN - Generative Artificial Nonsense

8

Sorry for necroing this thread. But I was inspired by what @kins did with ChatGPT, so I used AI to help me write a script which automates the calibrations and the interpolation. So far, it’s working well for me. Please let me know what you think. Here it is:

import time
import shutil
from pathlib import Path
import csv
import subprocess
import pandas as pd
import numpy as np
import math
from scipy.signal import find_peaks
from datetime import datetime

home = Path.home()

# Function to read and update the PID calibration file
def update_pid_calibrate_file(old_delta, new_delta):
    with open(home / 'klipper/klippy/extras/pid_calibrate.py', 'r') as file:
        filedata = file.read()
    
    filedata = filedata.replace(f'TUNE_PID_DELTA = {old_delta}', f'TUNE_PID_DELTA = {new_delta}')
    
    with open(home / 'klipper/klippy/extras/pid_calibrate.py', 'w') as file:
        file.write(filedata)

# Function to restart Klipper using systemctl
def restart_klipper():
    subprocess.run(["sudo", "systemctl", "restart", "klipper"], check=True)

# Function to wait for the log entry to start with "Stats "
def wait_for_stats_entry():
    log_file_path = home / 'printer_data/logs/klippy.log'
    target_string = "Stats "

    print("Waiting for Klipper to finish restarting...")
    time.sleep(10)  # Add delay to wait for klippy to start logging
    while True:
        with open(log_file_path, 'r') as log_file:
            lines = log_file.readlines()
            if lines:
                # print(f"Current last line: {lines[-1].strip()}")  # Debug print statement
                if lines[-1].startswith(target_string):
                    break
        time.sleep(1)  # Wait a bit before checking again

# Function to append command to klippy.serial
def send_calibrate_command(delta, heater, target):
    with open(home / 'printer_data/comms/klippy.serial', 'a') as myFile:
        print(f"PID_CALIBRATE HEATER={heater} TARGET={target} WRITE_FILE=1", file=myFile)
    print(f"Sent calibration command for {heater} with target {target} and delta {delta}")

# Function to wait for the log entry
def wait_for_log_entry(heater):
    log_file_path = home / 'printer_data/logs/klippy.log'
    target_string = f"save_config: set [{heater}] pid_Kd ="

    print("Waiting for calibration to complete...")
    while True:
        with open(log_file_path, 'r') as log_file:
            lines = log_file.readlines()[-10:]  # Get the last 10 lines
            for line in lines:
                if line.startswith(target_string):
                    return
        time.sleep(1)  # Wait a short time before checking again

# Function to copy and parse the output file
def process_output_file(delta, heater, target):
    timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
    source_file = Path('/tmp/heattest.txt')
    destination_dir = home / 'PID' / 'Raw_Data'
    destination_file = destination_dir / f'heattest_{heater}_{target}_{delta}_{timestamp}.txt'

    print(f"Copying output file for delta {delta}...")
    destination_dir.mkdir(parents=True, exist_ok=True)
    shutil.copy(source_file, destination_file)

    # Parse the output file
    with open(destination_file, 'r') as file:
        lines = file.readlines()

    # Filter out lines that start with 'pwm:'
    filtered_lines = [line for line in lines if not line.startswith('pwm:')]

    # Write the filtered data to a new CSV file
    csv_file = destination_dir / f'parsed_heattest_{heater}_{target}_{delta}_{timestamp}.csv'
    print(f"Parsing output file for delta {delta}...")
    with open(csv_file, 'w', newline='') as file:
        writer = csv.writer(file)
        writer.writerow(['Time', 'Temperature'])  # Write the header
        for line in filtered_lines:
            time, temp = line.strip().split()
            writer.writerow([time, temp])

    return csv_file  # Return CSV file path for further processing

# Function to find ultimate period Tu
def calculate_Tu(time, temp):
    peaks, _ = find_peaks(temp,prominence=1)
    peak_times = time[peaks[1:]]  # Ignore the first peak
    Tu = np.diff(peak_times).mean()  # Average period between peaks
    return Tu, temp[peaks[1:]]

# Function to calculate amplitude
def calculate_amplitude(temp):
    peaks, _ = find_peaks(temp,prominence=1)
    valleys, _ = find_peaks(-temp,prominence=1)

    # Ensure equal number of peaks and valleys
    num_cycles = min(len(peaks), len(valleys)) - 1
    peaks = peaks[1:]
    valleys = valleys[1:]

    # Create new lists with sequential indexing
    sequential_peaks = []
    sequential_valleys = []
    for i in range(num_cycles):
        sequential_peaks.append(temp[peaks[i]])
        sequential_valleys.append(temp[valleys[i]])

    amplitudes = []
    for i in range(num_cycles):
        amplitude = (sequential_peaks[i] - sequential_valleys[i]) / 2
        amplitudes.append(amplitude)

    average_amplitude = sum(amplitudes) / len(amplitudes)
    return average_amplitude

# Function to calculate ultimate gain Ku
def calculate_Ku(amplitude, max_power):
    return 4.0 * max_power / (math.pi * amplitude)

# Function to interpolate values
def interpolate(x1, y1, x2, y2, x):
    return y1 + (y2 - y1) / (x2 - x1) * (x - x1)

# Function to calculate PID parameters based on various tuning rules
def calculate_pid_parameters(Ku, Tu, rule):
    if rule == "Classic Ziegler-Nichols":
        Kp = 0.6 * Ku * 255
        Ti = 0.5 * Tu
        Td = 0.125 * Tu
    elif rule == "Pessen Integral Rule":
        Kp = 0.7 * Ku * 255
        Ti = 0.4 * Tu
        Td = 0.15 * Tu
    elif rule == "Some Overshoot":
        Kp = 0.33 * Ku * 255
        Ti = 0.5 * Tu
        Td = 0.33 * Tu
    elif rule == "No Overshoot":
        Kp = 0.2 * Ku * 255
        Ti = 0.5 * Tu
        Td = 0.33 * Tu
    else:
        raise ValueError("Unknown rule")
    
    Ki = Kp / Ti
    Kd = Kp * Td
    return Kp, Ki, Kd

print("\nPID Calibration Menu:")
print("1. Calibrate Extruder")
print("2. Calibrate Bed")
print("3. Exit")

while True:
    choice = input("Enter your choice (1/2/3, press Enter for Extruder): ")
    if choice == '1' or choice == '':
        heater = "extruder"
        while True:
            try:
                target_temp = int(input("Enter target temperature (default 250): ") or 250)
                break
            except ValueError:
                print("Invalid input. Please enter an integer value for the target temperature.")
        break
    elif choice == '2':
        heater = "heater_bed"
        while True:
            try:
                target_temp = int(input("Enter target temperature (default 100): ") or 100)
                break
            except ValueError:
                print("Invalid input. Please enter an integer value for the target temperature.")
        break
    elif choice == '3':
        print("Exiting...")
        exit()
    else:
        print("Invalid choice. Please try again.")

while True:
    try:
        max_power = float(input("Enter maximum power (default 1.0): ") or 1.0)
        break
    except ValueError:
        print("Invalid input. Please enter a float value for maximum power.")

while True:
	try:
		first_delta = float(input("Enter first delta value (default 2.5): ") or 2.5)
		break
	except ValueError:
		print("Invalid input. Please enter a float value for the first delta.")

while True:
	try:
		second_delta = float(input("Enter second delta value (default 5.0, 0 to skip): ") or 5.0)
		break
	except ValueError:
		print("Invalid input. Please enter a float value for the second delta.")

# First run for first_delta
print(f"\nStarting first calibration run with TUNE_PID_DELTA = {first_delta}")
update_pid_calibrate_file(5.0, first_delta)
restart_klipper()
wait_for_stats_entry()
send_calibrate_command(first_delta, heater, target_temp)
wait_for_log_entry(heater)
csv_file_first = process_output_file(first_delta, heater, target_temp)

# Only run second calibration if second_delta is not 0
if second_delta != 0:
    print(f"\nStarting second calibration run with TUNE_PID_DELTA = {second_delta}")
    update_pid_calibrate_file(first_delta, second_delta)
    restart_klipper()
    wait_for_stats_entry()
    print("Waiting 1 minute for heater to cool down")
    time.sleep(60)
    send_calibrate_command(second_delta, heater, target_temp)
    wait_for_log_entry(heater)
    csv_file_second = process_output_file(second_delta, heater, target_temp)
else:
    csv_file_second = None

# Restore the PID calibration file to original value
update_pid_calibrate_file(second_delta if second_delta != 0 else first_delta, 5.0)
restart_klipper()

# Load the parsed CSV data
data_first = pd.read_csv(csv_file_first)

# Calculate Tu and amplitude for the first dataset
Tu_first, peaks_first = calculate_Tu(data_first['Time'], data_first['Temperature'])
amplitude_first = calculate_amplitude(data_first['Temperature'])
valleys_f, _ = find_peaks(-data_first['Temperature'], prominence=1)
valleys_first = data_first['Temperature'][valleys_f[1:]]

# Calculate ultimate gain Ku for the first dataset
Ku_first = calculate_Ku(amplitude_first, max_power)

# Define the tuning rules
rules = ["Classic Ziegler-Nichols", "Pessen Integral Rule", "Some Overshoot", "No Overshoot"]

# Prepare to write results to a text file
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
results_file = home / f'PID/results_{heater}_{target_temp}_{timestamp}.txt'

with open(results_file, 'w') as file:
    file.write(f'Target Temperature: {target_temp}\n')  # Record the target temperature

# Calculate and print PID parameters for each tuning rule and the first dataset
print(f'\n\033[91mTune PID Delta={first_delta}:\033[0m')
print(f'Peak values: {peaks_first.tolist()}')
print(f'Valley values: {valleys_first.tolist()}')
print(f'Max peak value: {max(peaks_first)}')
print(f'Min valley value: {min(valleys_first)}')
print(f'Amplitude: {amplitude_first}')
print(f'Ultimate Period Tu={Tu_first}')
print(f'Ultimate Gain Ku={Ku_first}')

with open(results_file, 'a') as file:
    file.write(f'\nTune PID Delta={first_delta}:\n')
    file.write(f'Peak values: {peaks_first.tolist()}\n')
    file.write(f'Valley values: {valleys_first.tolist()}\n')
    file.write(f'Max peak value: {max(peaks_first)}\n')
    file.write(f'Min valley value: {min(valleys_first)}\n')
    file.write(f'Amplitude: {amplitude_first}\n')
    file.write(f'Ultimate Period Tu={Tu_first}\n')
    file.write(f'Ultimate Gain Ku={Ku_first}\n')

    for rule in rules:
        Kp, Ki, Kd = calculate_pid_parameters(Ku_first, Tu_first, rule)
        pid_output = f"#*# [{heater}]\n#*# control = pid\n#*# pid_kp = {Kp:.3f}\n#*# pid_ki = {Ki:.3f}\n#*# pid_kd = {Kd:.3f}"
        print(f'\033[93m{rule}\033[0m PID parameters:')
        print(pid_output)
        file.write(f'{rule} PID parameters:\n')
        file.write(f'{pid_output}\n')

# Only calculate and print interpolated PID parameters if second_delta is not 0
if second_delta != 0:
    # Load the second parsed CSV data
    data_second = pd.read_csv(csv_file_second)
    
    # Calculate Tu and amplitude for the second dataset
    Tu_second, peaks_second = calculate_Tu(data_second['Time'], data_second['Temperature'])
    amplitude_second = calculate_amplitude(data_second['Temperature'])
    
    # Calculate ultimate gain Ku for the second dataset
    Ku_second = calculate_Ku(amplitude_second, max_power)

    # Define valleys_second
    valleys_s, _ = find_peaks(-data_second['Temperature'], prominence=1)
    valleys_second = data_second['Temperature'][valleys_s[1:]]

    print(f'\n\033[91mTune PID Delta={second_delta}:\033[0m')
    print(f'Peak values: {peaks_second.tolist()}')
    print(f'Valley values: {valleys_second.tolist()}')
    print(f'Max peak value: {max(peaks_second)}')
    print(f'Min valley value: {min(valleys_second)}')
    print(f'Amplitude: {amplitude_second}')
    print(f'Ultimate Period Tu={Tu_second}')
    print(f'Ultimate Gain Ku={Ku_second}')
    
    with open(results_file, 'a') as file:
        file.write(f'\nTune PID Delta={second_delta}:\n')
        file.write(f'Peak values: {peaks_second.tolist()}\n')
        file.write(f'Valley values: {valleys_second.tolist()}\n')
        file.write(f'Max peak value: {max(peaks_second)}\n')
        file.write(f'Min valley value: {min(valleys_second)}\n')
        file.write(f'Amplitude: {amplitude_second}\n')
        file.write(f'Ultimate Period Tu={Tu_second}\n')
        file.write(f'Ultimate Gain Ku={Ku_second}\n')
        
        for rule in rules:
            Kp, Ki, Kd = calculate_pid_parameters(Ku_second, Tu_second, rule)
            pid_output = f"#*# [{heater}]\n#*# control = pid\n#*# pid_kp = {Kp:.3f}\n#*# pid_ki = {Ki:.3f}\n#*# pid_kd = {Kd:.3f}"
            print(f'\033[93m{rule}\033[0m PID parameters:')
            print(pid_output)
            file.write(f'{rule} PID parameters:\n')
            file.write(f'{pid_output}\n')

    print(f'\n\033[91mInterpolated\033[0m PID parameters for \033[91mTUNE_PID_DELTA=0.0:\033[0m')
    with open(results_file, 'a') as file:
        file.write(f'\nInterpolated PID parameters for TUNE_PID_DELTA=0.0:\n')
        for rule in rules:
          Kp_first, Ki_first, Kd_first = calculate_pid_parameters(Ku_first, Tu_first, rule)
          Kp_second, Ki_second, Kd_second = calculate_pid_parameters(Ku_second, Tu_second, rule)

          # Interpolate Ku and Tu
          Ku_0 = interpolate(second_delta, Ku_second, first_delta, Ku_first, 0.0)
          Tu_0 = interpolate(second_delta, Tu_second, first_delta, Tu_first, 0.0)

          # Calculate final PID parameters using interpolated Ku and Tu
          Kp_0, Ki_0, Kd_0 = calculate_pid_parameters(Ku_0, Tu_0, rule)

          print(f'\033[93m{rule}\033[0m PID parameters:')
          pid_output = f"#*# [{heater}]\n#*# control = pid\n#*# pid_kp = {Kp_0:.3f}\n#*# pid_ki = {Ki_0:.3f}\n#*# pid_kd = {Kd_0:.3f}"
          print(pid_output)
          file.write(f'{rule} PID parameters:\n')
          file.write(f'{pid_output}\n')

print(f"Results written to {results_file}")

You might have to install pandas to run the script

pip3 install pandas

I also modified my system to allow the user “pi” to run the command “sudo systemctl restart klipper” without requiring a password in order to automate the whole process. But that’s optional.