I started playing around with simulating an attempt to automate PID tuning using the Ziegler Nichols method (using this as a reference: https://en.wikipedia.org/wiki/Ziegler%E2%80%93Nichols_method). I'm attaching the code I have been using. 

The gist of what I'm trying to simulate:

-Increase Kp by some increment

-Get step response data from closed-loop system

-Fourier transform step response data to see if there are stable oscillations (currently done by comparing to some reference magnitude that I am guesstimating/trying to determine heuristically although once I do that, it might be able to be related to the Kp of the current iteration)

-If there are stable oscillations, log the oscillation period and then calculate PID values

Issues/thoughts:

There's an envelope around 0Hz because of the DC component, and I'm trying to figure out a way to remove it. What I've come up with so far is heuristically determining how wide it is and then just chopping that section off, but if the oscillation frequency is within that envelope (especially since temperature control is really low frequency) then that kind of ruins the whole point. Also I'm relying on heuristically determining a few things that would be specific to each loop, which I'm not super happy about, but even the more optimized (and complicated) approach to PID tuning I found relied on heuristically determining a few things about the system to set up the tuning algorithm. After Googling this for a while, I don't know if/don't think it's even feasible to totally get around tweaking heuristically determined parameters for tuning each system.

Anyway, I'm going to keep poking at this. If anyone sees this and has suggestions, that would be super appreciated.

Edit 2/5/2017: Whelp, realized I should probably be LombScargling? To be continued...

#!/usr/bin/env python

import numpy, time
from control.matlab import *
import matplotlib.pyplot as plt

def simulate_heater(total_time,timestep,Kp,Ki,Kd):
	num_steps = numpy.ceil(total_time/timestep)
	timestamps = linspace(0,total_time,num=num_steps)